[R-sig-ME] glmmADMB random interaction effects and likelihood estimates

[Ben Bolker]

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