[ADMB Users] difference between ADMB-RE and R/mgcv in SEs for smoother coefficients in a GAM fitted by maximum likelihood

[Tim Miller (NOAA Federal)]

Dear ADMB users,

I am fitting some GAMs via maximum likelihood by reparameterizing as a 
mixed effects problem.  The penalty in the penalized likelihood that is 
more commonly the maximized function for gams can be thought of as prior 
so that some of the smoother coefficients are random effects with a 
specific variance structure and the smoothing parameter is then an 
associated variance parameter. Simon Wood's book (Generalized Additive 
Models: An introduction with R) explains this idea better than I ever 
could (see Section 6.6.1). As you can see in that section of Wood's 
book, the mixed effects reparameterization is pretty straightforward. It 
just involves separating the eigenvectors corresponding to the zero and 
non-zero eigenvalues of the cubic spline penalty matrix and then using 
the different sets of eigenvectors to form fixed and random effects 
design matrices from the original smoother basis matrix.

I am using this reparameterization and the ADMB random effects module to 
fit some binomial GAMs and other models with more complicated random 
effects structures and distributional assumptions. I use the mgcv 
package in R to obtain the smoother basis matrix and the penalty matrix 
for this.  For a simple binomial gam with a cubic spline smoother, I can 
replicate results for the maximized likelihood, smoother coefficients 
(after appropriate transformation), and predictions provided by the mgcv 
package in R (with very negligible differences I would attribute to the 
differences in estimation procedures). However, the standard 
errors/covariance matrix of the coefficients (after the appropriate 
transformation) and standard errors of predicted smoother are very 
different from those given by R /mgcv.

What is strange and confounding to me is that if I simply multiply one 
of the eigenvectors corresponding to the fixed effects by -1, I get 
standard errors and covariance matrix for the (now incorrect) 
transformed coefficients very similar to that given by R/mgcv. (I only 
stumbled on this because I realized that eigen() in R does not always 
give the same set of eigenvectors and I was having some issues with that 
initially).  Anyway, I corresponded with Simon Wood initially because I 
didn't think this was an ADMB issue, but he can't seem to figure out 
anything I'm doing wrong technically speaking. We were wondering whether 
this might have something to do with the second term in the ADMB random 
effects variance  (Eq. 3.1 in the ADMB-RE manual). At Hans Skaug's 
instruction I approximated the first term in Eq. 3.1 by fixing the fixed 
effects/hyperparameters at the MLEs and estimating the random effects, 
but these standard errors are still not similar to those given by R/mgcv.

The standard errors for the predicted smoother given by my ADMB program 
are often much larger than what is given by R/mgcv and the trend is very 
odd. I've checked repeatedly that this isn't due to using different 
eigenvectors for fitting vs. predictions and I've even gone so far as to 
do all of the calculations internally in the ADMB program using the 
eigenvalues/eigenvectors functions and I still get the same results.  To 
help illustrate the problem better, attached is a plot (comparison.pdf) 
of the smooth given by plot(gam.object) (black) from R/mgcv with the 
same results produced by ADMB added. Solid red (right on top of R/mgcv 
result) are the ADMB estimates. Dashed red is the 95% CI based on the 
ADMB (correct but "bad"?) covariance matrix. Dashed blue is based on the 
covariance matrix that uses the "wrong" eigenvectors (one vector 
multiplied by -1). Note how much wider the CIs based on the (correct?) 
ADMB covariance matrix can be than those given by R/mgcv.

A question that is also maybe more generally applicable is how exactly 
is the covariance/correlation of the MLEs and random effects provided by 
ADMB defined (in say the .cor or the binary .bgs files)? I cannot seem 
to be able to find an answer to this in the ADMBRE manual.  I understand 
covariance of MLEs is a function of the likelihood surface and 
covariance of random effects as a function of the posterior, but I'm 
unclear about the covariance of the two. I wonder whether this is 
related to the strange standard errors because the transformed 
coefficients are functions of both the fixed effects and random effects 
parts of the smoother.

Also attached are simplified tpl and dat files that will produce the 
results below and which are used to make the plots in comparison.pdf.

The maximized log-likelihoods for R/mgcv and ADMB are -850.1552 and 
-850.154, respectively (ADMB by a hair!).

The intercept (first) and smoother coefficients that R/mgcv produce for 
these data are
1.3216274   0.7395298   0.1553618  -0.6677414  -1.1453418  -0.9579393  
-0.6969166  -0.6534421  -0.9104002  -1.0061178
with corresponding standard errors
0.1066258   0.1351087   0.1673951   0.2136631   0.2421144   0.2080006   
0.1791422   0.2958478   0.5712558   0.9365626

The coefficients produced by the ADMB program are
1.3241e+000 7.4057e-001 1.5283e-001 -6.7285e-001 -1.1514e+000 
-9.6254e-001 -6.9886e-001 -6.5363e-001 -9.1114e-001 -1.0075e+000
The strange ADMB standard errors are
1.0888e-001 3.5051e-001 2.3446e-001 2.5181e-001 4.8167e-001 6.1297e-001 
6.1003e-001 5.3677e-001 3.9467e-001 3.3868e-001
Those provided when one of the eigenvectors is multipled by-1 are
1.0888e-001 1.6316e-001 1.9334e-001 2.1642e-001 2.4441e-001 2.3118e-001 
1.9598e-001 2.9668e-001 6.0026e-001 1.0339e+000

Using AMDB 11, MinGW GCC-4.5.2, and Windows 7.

I appreciate any insights anyone may have on this problem.

Best Regards,
Tim


-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: binom_gam_ex.tpl
URL: <http://lists.admb-project.org/pipermail/users/attachments/20121128/dd76f5c2/attachment.ksh>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: binom_gam_ex.dat
URL: <http://lists.admb-project.org/pipermail/users/attachments/20121128/dd76f5c2/attachment-0001.ksh>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: comparison.pdf
Type: application/x-pdf
Size: 43175 bytes
Desc: not available
URL: <http://lists.admb-project.org/pipermail/users/attachments/20121128/dd76f5c2/attachment.bin>