minor numerical change/issue in test-AGHQ test

[paciorek]

In the "AGH Quadrature Comparison to LME4 1 RE for Poisson-Normal" test, line 475 of test-AGHQ.R (formerly test-ADaghq.R), we are now getting this optimization "failure":

  [Warning] In maximizing the Laplace/AGHQ approximation,
            `optim` has a non-zero convergence code: 1,
            with the message 'false convergence (8)'.
            The control parameters of `optim` can be adjusted using the `outerOptimControl`
            list component of the `control` list argument of `buildLaplace` or `buildAGHQ`.

When I compare results to nimble 1.3.0, the numerical values at the optimum are essentially identical and in fact the new values are very slightly better than the old ones in terms of the Laplace-approximated logLik and the optimization stops in 9 iterations instead of 10 previously. Also, when using 1.3.0, if I re-run cmLaplace$findMLE()$par immediately again (i.e., starting a new optimization from the optimized value), it gives a lack of convergence warning. So I think that there are just some minor numerical issues with the optimization generally for this model. That said, it's a bit disconcerting to get an optim "failure" when we didn't before.

As far as I can see, the difference seems to be caused by a very slightly different optimization path because the results of calcLogLik2 and gr_logLik2 in the 1D inner AGHQ/Laplace are slightly different with Perry's recent AD/Laplace changes. The differences are not present for all of the 50 Laplaces involved in this model , but for some of them there are very small floating point differences. I'm presuming that these add up enough to affect the optimization path. I should also note that Paul and I did change how gr_logLik/gr_logDens work to handle inclusion/exclusion of the prior. I don't think this would change numerical values, but it's possible I'm being to cavalier in saying that.

I'm not particularly inclined to pursue this further, and I can monkey with the test to deal with the warning, but wanted to get thoughts from @perrydv @paul-vdb if they have any.

[paul-vdb]

@paciorek I don't love that it is enough to throw a non-convergence error. I'll try and find time this week to look at this. Maybe I'll end up at the same conclusion that it's not worth it for now and we can just adapt the test.

[paciorek]

I did tweak the convergence tolerance for this test from 1e-12 to 1e-10, and the warning disappeared. I'll leave this open for potential further investigation, but wanted to note that change.

mQuad <- buildAGHQ(model = m, nQuad = 21, control = list(outerOptimControl = list(reltol = 1e-10)))
mLaplace <- buildAGHQ(model = m, nQuad = 1, control = list(outerOptimControl = list(reltol = 1e-10)))

[perrydv]

It sounds like just numerical hassles, but one thing we should maybe look out for is if the gradients are in fact accurate. The optimizer could get confused if it can't get higher but is getting non-zero gradients, or such kinds of problems.

[paciorek]

One note here is that in issue #64 I was thinking we might want to use reltol=1e-10 as our general default based on the 1-D case there. This is a non-1-D case that suggests the same conclusion. If I change the test discussed here to use our default reltol, the error in the MLE compared to LME4 is rather larger:

With reltol=1e-10:

> abs(mleLaplace- mleLME4_laplace)
[1] 7.077823e-06 3.728186e-05
> abs(mleQuad- mleLME4_nquad21)
[1] 8.172672e-06 2.777471e-05

With default reltol:

> abs(mleLaplace- mleLME4_laplace)
[1] 1.122801e-04 6.369764e-05
> abs(mleQuad- mleLME4_nquad21)
[1] 1.136439e-04 5.438433e-05