[ADMB Users] Integrating R and ADMB to generalize TPL files

[Ben Bolker]

On 12-11-07 01:17 PM, dave fournier wrote:
> But that is not right.  The model is to assume that the
> random effects in different time periods are independent and that their
> variances are different.  You can do this by figuring out what the
> correct covariance is for the current design matrix or by
> writing the model properly which is trivial to do in ADMB.  As usual it
> is the R interface and these silly design matrices which introduces the
> problem.
> _______________________________________________
> Users mailing list
> Users at admb-project.org
> http://lists.admb-project.org/mailman/listinfo/users


  Just to continue this a little bit:

If I do this:

l2 <- read.csv("bigelow_dat.txt")

## unnecessary picture
library(ggplot2)
ggplot(l2,aes(factor(time),snag))+
    geom_dotplot(binaxis="y",stackdir="center")+
    facet_wrap(~Ecozone+time,scale="free_x")

with(l2,table(table(plot))) ## check exp. design: 3 samples per plot
with(l2,table(table(plot,time))) ## check exp. design: 1 sample per
plot/time
l2$ftime <- factor(l2$time) ## factor-ize time

library(glmmADMB)
## note zero-inflation is FALSE/off throughout (default)
g1 <- glmmadmb(formula = snag ~ ftime + Ecozone + (1 | plot), data = l2,
family = "nbinom1")
VarCorr(g1)
log(sqrt(VarCorr(g1)[[1]]))  ## log std. dev. =0.4315505  ## matches DF
logLik(g1)
##  'log Lik.' -255.758 (df=6)  ## matches DF
coef(g1)
## (Intercept)      ftime2      ftime3    EcozoneW
##   -0.7635323  -0.5884655  -0.6037133   1.2493021

## ?? The fixed effects don't match DF results: did he orthogonalize the
design matrix ... ?
## never mind

## besides being the wrong model, this fails on my system (32-bit Linux)
try(g2 <- glmmadmb(formula = snag ~ ftime + Ecozone + (ftime | plot),
data = l2, family = "nbinom1"))

## this works and is (I think) the right model ...
g3 <- glmmadmb(formula = snag ~ ftime + Ecozone + (ftime-1 | plot), data
= l2, family = "nbinom1",
               save.dir="glmmadmb_tmp")
## the Z matrix looks like I think it should (an indicator variable for
time period)
VarCorr(g3)

## Random effects:
## Structure: Diagonal matrix
## Group=plot
##       Variance StdDev
## time1   0.7957 0.8920
## time2   0.5986 0.7737
## time3   0.9928 0.9964

log(sqrt(diag(VarCorr(g3)[[1]])))
##      ftime1      ftime2      ftime3
## -0.11426024 -0.25658084 -0.00360295

## This doesn't match DF: he gets the log-sd of the variances as

## tmpL:
##  0.352458239856 0.489012745656 0.496221176099

## i.e. increasing from period 1 to periods 2 and 3

For the overdispersion model (which I haven't run yet -- glmmADMB is not
flexible enough) -- DF gets

# tmpL:
 0.436593520224
# tmpL1:
 0.000100000000000
# log_alpha:
 0.826338727816 0.0677913813376 0.657156858701

i.e. the log-sd among plots (tmpL) is about the same as before, and
log_alpha is the log of the NB scale parameter.  This suggests *lower*
variance (at least when scaled by the mean) in the second period,
*higher* variance in the first period -- which seems (??) to match the
(ftime-1|plot) glmmADMB results rather than the DF results (??)

  I'd rather be arguing over what the data say than over whose tools are
better
or worse, or who's an idiot.  I'm happy to concede that I'm an idiot, so we
can get that over right at the beginning ...

  cheers
    Ben