R help - Apr 2003 - glmmPQL and additive random effects?

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

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))))