glmmADMB random interaction effects and likelihood estimates
Seth W. Bigelow <seth at ...> writes:
Dear list: I'm (still) working on a model of snag (dead tree) abundance
before, after, & 5 yr after forest treatments using negative binomial
distribution with glmmADMB. The data suggest that variance may differ
sharply before and after treatment (factor variable 'time', values 1/2/3),
so I have tried specifying the random term to model variance within time.
e.g., glmmadmb(snag ~ time + Ecozone + (time|plot)), where Ecozone is a
factor variable which takes values East or West, and 'plot' is a subject
variable. The issue is that log-likelihood does not change whether I use
(1|plot) or (time|plot). Am I missing something obvious, or is there a bug
in the way likelihood is being calculated/reported?
Appreciatively, Seth W. Bigelow
glmmadmb(formula = snag ~ time + Ecozone + (1 | SettingID), data = l2,
family = "nbinom1", zeroInflation = F)
[snip]
Number of observations: total=192, SettingID=64 Random effect variance(s): Group=SettingID
Variance StdDev
(Intercept) 2.37 1.54
Negative binomial dispersion parameter: 1.812 (std. err.: 0.33774)
Log-likelihood: -255.758
glmmadmb(formula = snag ~ time + Ecozone + (time | SettingID),
data = l2, family = "nbinom1", zeroInflation = F)
[snip]
Random effect variance(s):
Group=SettingID
Variance StdDev
(Intercept) 2.370e+00 1.5396428
time2 1.765e-07 0.0004201
time3 8.542e-08 0.0002923
Negative binomial dispersion parameter: 1.812 (std. err.: 0.33774) Log-likelihood: -255.758
My guess here is that the log-likelihood differs only by a tiny amount, because the estimated variances for time2 and time3 are about 7-8 orders of magnitude less than the intercept variance. The basic issue is that in a data set of this size (not tiny, but not huge and quite noisy), you just don't have much power to detect differences in slope among plots. You got a bit lucky getting the models to fit at all; on my platform (Ubuntu 10.04 32-bit) the models with a time:plot interaction (with time either continuous or a factor) failed because the hessian was non-positive-definite. You could compare logLik() for the two models to see if they differ by some tiny amount, but it's also possible (because of the way AD Model Builder returns its results) that the log-likelihoods will have been rounded sufficiently that they appear identical. Ben Bolker