I'm a bit puzzled by how to write out additive random effects in glmmPQL. In my situation, I have a factorial design on two (categorical) random factors, A and B. At each combination, I have a binary response, y, and two binary fixed covariates, C and D. If everything were fixed, I would use glm(y ~ A + B + C + D, family = binomial) My first thought was to use glmmPQL(y ~ A + B, random = ~ C + D, family = binomial) but glmmPQL wants to see a grouping variable in the random term. Something like glmmPQL(y ~ A + B, random = ~ C + D | CD, family = binomial) where CD is a a variable combining C and D, eats up all my memory, while glmmPQL(y ~ A + B, random = ~ 1 | CD, family = binomial) doesn't seem like the model I want. Perhaps this model is too hard to fit, but before I quit this approach I want to make sure that I'm not just coding it incorrectly. Thanks, Steve Buyske
glmmPQL and additive random effects?
2 messages · Steve Buyske, Douglas Bates
Steve Buyske <buyske at stat.rutgers.edu> writes:
I'm a bit puzzled by how to write out additive random effects in glmmPQL. In my situation, I have a factorial design on two (categorical) random factors, A and B. At each combination, I have a binary response, y, and two binary fixed covariates, C and D. If everything were fixed, I would use glm(y ~ A + B + C + D, family = binomial) My first thought was to use glmmPQL(y ~ A + B, random = ~ C + D, family = binomial) but glmmPQL wants to see a grouping variable in the random term. Something like glmmPQL(y ~ A + B, random = ~ C + D | CD, family = binomial) where CD is a a variable combining C and D, eats up all my memory, while glmmPQL(y ~ A + B, random = ~ 1 | CD, family = binomial) doesn't seem like the model I want. Perhaps this model is too hard to fit, but before I quit this approach I want to make sure that I'm not just coding it incorrectly.
lme and, by extension, glmmPQL do not handle crossed random effects easily. You must create a factor of the same length as y, A, B, C, and D with a single level const = factor(rep(1, length(y))) then use the non-obvious formulation glmmPQL(y ~ A + B, random = list(const = pdBlocked(pdIdent(~ C - 1), pdIdent(~ D - 1))))