glmmPQL Help: Random Effect and Dispersion Parameter
Hi: glmmPQL has been around a while, and I suspect it was not meant to handle crossed random effects. This was one of the original motivations for the lme4 package, and it seems to work there, although it's using Gauss-Hermite approximations to the likelihood rather than PQL: library(lme4) mod1 <- lmer(y ~ 1 + (1 | beta) + (1 | alpha), data = simu)
summary(mod1)
Linear mixed model fit by REML ['summary.mer']
Formula: y ~ 1 + (1 | beta) + (1 | alpha)
Data: simu
REML criterion at convergence: 584.4204
Random effects:
Groups Name Variance Std.Dev.
alpha (Intercept) 3.11128 1.7639
beta (Intercept) 0.17489 0.4182
Residual 0.05405 0.2325
Number of obs: 1500, groups: alpha, 100; beta, 5
Fixed effects:
Estimate Std. Error t value
(Intercept) 2.9167 0.2572 11.34
Hopefully that's closer to what you had in mind. If not, take a look
at Ben Bolker's GLMM wiki:
http://glmm.wikidot.com/faq
BTW, thank you for the nice reproducible example.
Dennis
On Thu, Jun 23, 2011 at 9:16 PM, Yue Yu <parn.yy at gmail.com> wrote:
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
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models