[R-sig-ME] glmmPQL Help: Random Effect and Dispersion Parameter

[Yue Yu]

Dear R users,

I am currently doing a project in generalized mixed model, and I find
the R function glmmPQL in MASS can do this via PQL. But I am a newbie
in R, the input and output of glmmPQL confused me, and I wish I can
get some answers here.

The model I used is a typical two-way generalized mixed model with
random subject (row) and block (column) effects and log link function,
y_{ij} = exp(\mu+\alpha_i+\beta_j)+\epsilon.

I can generate a pseudo data by the following R code

===================================================
k <- 5; # Number of Blocks (Columns)
n <- 100; # Number of Subjects (Rows)
m <- 3; # Number of Replications in Each Cell

sigma.a <- 0.5; # Var of Subjects Effects
sigma.b <- 0.1; # Var of Block Effects
sigma.e <- 0.01; # Var of Errors
mu <- 1; # Overall mean

a <- rep(rnorm(n,0,sigma.a),each=k*m);
b <- rep(rep(rnorm(k,0,sigma.b),each=m),n);

# Simulate responses y_{ij}=exp(\mu+\alpha_i+\beta_j)+e
y <- exp(mu+a+b)+rnorm(1,0,sigma.e);

# Indicator vector of subject effects alpha
alpha <- rep(seq(1,n),each=k*m);

# Indicator vector of block effects beta
beta <- rep(rep(seq(1:k),each=m),n);

simu <- data.frame(y,beta,alpha)
===================================================

And I want to use glmmPQL to estimate the mean and variance
components, but I have several questions.

1. How to write the random effect formula?
I have tried
glm <- glmmPQL(y~1,random=~alpha+beta,family=gaussian(link="log"),data=simu);
but it did not work and got the error message "Invalid formula for groups".

And the command
glm <- glmmPQL(y~1,random=~1|alpha/beta,family=gaussian(link="log"),data=simu)
worked, but the result was the nested "beta %in% alpha" variances,
which was not what I want.

2. How to find the estimated variance-covariance matrix for the
variance components, which should be the inverse of information
matrix.
I notice the output variable glm at apVar will give a similar
variance-covariance matrix, but it has the prefix "reStruct." and
attribute "Pars", what are these stand for? I can't find any
explanation in the help document.

3. I am also wondering if there is a way to calculate the dispersion
parameter or not?

Anyone has some tips? Any suggestions will be greatly appreciated.

Best,

Yue Yu