[R-sig-ME] Why se.fit differ in predict.glm and predict.glmmadmb?

[Ben Bolker]

Xandre <alex at ...> writes:

> 
> I am aware of that, sorry. Hopefully, with the following link, 
> https://www.dropbox.com/s/pptl595chabjtly/datos.csv?dl=0 , the example 
> will be completely reproducible.
> 
> Regards,
> 
> Alex

  It looks like glmmADMB is getting the estimate of the correlation
between the fixed effects wrong (I don't yet know whether this is 
a problem with
the optimization in AD Model Builder or an actual bug in the translation
of results from ADMB output to R):

datos <- read.csv("datos.csv")
M1<-glm(response~explanatory,
        data=datos,
        family="binomial")
library(glmmADMB)
M2 <- glmmadmb(response~explanatory,
             data=datos,
             family="binomial")
## Estimated covariance matrix may not be positive definite
##  4.0211 4.1655
coef(summary(M1))
coef(summary(M2))

newdatos <- with(datos,
          data.frame(explanatory=seq(min(explanatory),
                                     max(explanatory),length.out=10)))
predict(M1,type="link",newdata=newdatos,se.fit=TRUE)
predict(M2,type="link",newdata=newdatos,se.fit=TRUE)
X <- model.matrix(~explanatory, data = newdatos)

## compare var-cov matrices
vcov(M1)
vcov(M2)

## compare SEs
sqrt(diag(vcov(M1)))
sqrt(diag(vcov(M2)))

## compare correlations
cov2cor(vcov(M1))
cov2cor(vcov(M2))

## compare predicted SEs
sqrt(diag(X %*% vcov(M1) %*% t(X)))
sqrt(diag(X %*% vcov(M2) %*% t(X)))

## try with glmmTMB
library(glmmTMB)
M3 <- glmmTMB(response~explanatory,
             data=datos,
             family=binomial)
cov2cor(vcov(M3)$cond)



> >>
> >> Just to check problems of code I tried again with a much more simple
> >> example. I made a subset of my original data base (see attached .csv)
> >> and run a much more simple model as follows:
> >>
> >> *> M1<-glm(response~explanatory, **
> >> **+           data=datos,**
> >> **+           family="binomial")**
> >> **≥ M2<-glmmadmb(response~explanatory, **
> >> **+                 data=datos,**
> >> **+                 family="binomial")**
> >> **≥ **
> >> **≥ summary(M1)*
> >>
> >>

 [snip snip snip snip]