Before going in to the code, allow me to give you a quick primer on trust region methods. Let $f(x)$ be your objective function and $x_k$ be your current iterate. Iteration $k$ of a generic trust region method looks something like this:
- Pick a maximum step size, $\Delta_k$
- Build an model of $f(x)$ at $x = x_k$; call it $Q(x)$.
- Find the step, $s_k$, that minimizes $Q_k(x_k + s_k)$ subject to the constraint $||s_k|| \leq \Delta_k$
- If $s_k$ is "good enough", let $x_{k+1} = x_k + s_k$
- Otherwise, refine your model and try again
One of the the ways you determine if $s_k$ is "good enough" is by comparing the decrease predicted by the model to the actual decrease in the objective function. This is what the portion of the code after 430 does. However, before that, there's a quick sanity check to see if the model predicts a decrease AT ALL. And that's what's happening at 430.
To understand the value of VQUAD
, we first have to understand a few other variables. Fortunately, there are good comments right below the declaration of SUBROUTINE BOBYQB
. The salient variables are:
GOPT
, the gradient of the model
HQ
, the Hessian of the model
D
, the trial step (I called this $s_k$ above)
Beginning a few lines above 410, you'll see DO 410 J=1,N
. This begins a for-loop (and a nested for-loop) that evaluates the change predicted by the model using trial step D
. It accumulates the predicted change in VQUAD
. The first part of the for-loop evaluates the first-order terms and the nested for-loop evaluates the second-order terms. It would probably be easier to read if the loops were indented, like this:
DO 410 J=1,N
VQUAD=VQUAD+D(J)*GOPT(J)
DO 410 I=1,J
IH=IH+1
TEMP=D(I)*D(J)
IF (I .EQ. J) TEMP=HALF*TEMP
410 VQUAD=VQUAD+HQ(IH)*TEMP
There's another for-loop after this to incorporate other parameters in to the model. I have to admit, I don't fully understand this - my best guess is that it's particular to how they build the model.
At the end of all this, VQUAD
holds the change in objective function predicted by the model. So if VQUAD
is non-negative, that's bad. Now this particular solver can use an alternative step computation (probably a line search), which is where NTRITS
comes in to play. So the logic at 430 is saying, "If the last iteration used the alternative step computation AND the model does not predict a decrease AND IPRINT
> 0, print the warning message." Note that the solver is going to terminate regardless of the value of IPRINT
.
Speaking of IPRINT
, that value is passed to BOBYQA
by the calling function. In this case, your R routine is the calling function. There's a verbose
parameter to glmer
- I would be dimes to dollars that same value is passed to BOBYQA
. Try setting verbose
to 0 and you probably won't see the warning. But it won't change what's going on under the hood, of course.
Best Answer
Just for the sake of resolving the question, I'm going to state that this is probably due to a confusion over the modeling framework. It would rarely (if ever) be sensible to model an average egg mass with a Poisson distribution (which only applies to a unitless count variable). If you have average counts, and have a measurement of the total exposure (i.e. you have total counts and the area or time over which they were collected), you can do a Poisson model with an offset.
In this case it would make more sense to use a log-Normal (i.e., by transforming the response and then fitting a linear mixed model,
lmer(log(Avg_egg_mass) ~ ...
) or a Gamma model (the former, log-Normal approach is easier; I would generally only recommend a Gamma model in cases where there is a strong mechanistic or cultural reason to use one).