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

[Dennis Murphy]

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
>