[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)
timothy.j.miller at noaa.gov
Wed Nov 28 11:54:58 PST 2012
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>
More information about the Users
mailing list