[R-sig-ME] Fwd: generating contrasts from glmm.admb?

Ben Bolker bbolker at gmail.com
Tue Sep 17 21:19:58 CEST 2013


Maureen Durkin <mmdurkin at ...> writes:

> 
> Hi all,
 
> I am using package glmm.admb to run a negative binomial mixed
> effects model on count data of different disturbance for plots. I
> would like to test for differences in counts of disturbances between
> five site.s Not having any trouble with glmm.admb itself, but I
> can't seem to figure out how to generate a contrast matrix for
> differences between individual sites. It seems that the "contrast"
> package cannot handle objects from glmm.admb? Is this correct, and
> if so does anyone know of other methods to generate contrast methods
> from glmm.admb?

  I started to dive into this, but generating the matrix for all
pairwise contrasts is a bit of a pain.  Let's say you already have the matrix
of contrasts (say, each row is a contrast, so if 'sites' is a fixed
effect and you have a usual treatment-contrasts parameterization,
the contrast matrix for comparing site 2 against site 3 would
be c(0,-1,1,0,...))

set.seed(1)
d <- data.frame(x=rnbinom(1000,mu=3,size=1),
      f=factor(rep(LETTERS[1:5],each=200)))
library(glmmADMB)
m1 <- glmmadmb(x~f,data=d,family="nbinom")
cc <- coef(m1)

## compare site 2 vs 3; 3 vs 4; 4 vs 5
M <- matrix(c(0,-1,1,0,0,
              0,0,-1,1,0,
              0,0,0,-1,1),
            byrow=TRUE,nrow=3)

## The estimated value for each contrast:
cvals <- c(M %*% cc)

## The estimated variance-covariance matrix for the contrasts:
vv <- M %*% vcov(m1) %*% t(M)

## The Z-score for each contrast:
cvals/sqrt(diag(vv))

I'll let someone else worry about how to come up with
the all-pairwise-contrasts matrix.

The other alternative would be to bug the glmmADMB and/or
contrasts package maintainers to extend their packages to
allow them to co-operate.

This all assumes that sites are a *fixed* effect -- otherwise
(if they're random), hypothesis tests on the site-to-site
differences aren't meaningful.



More information about the R-sig-mixed-models mailing list