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

[Yue Yu]

Thanks a lot, Dennis.

The reason I am not using lme4 is that its estimated variance
component seems not desirable.

If you use my simulated data and run the following line

glm2 <- glmer(y~1+(1|alpha)+(1|beta),family=gaussian(link="log"),data=simu)

the result will be

Generalized linear mixed model fit by the Laplace approximation
Formula: y ~ 1 + (1 | alpha) + (1 | beta)
   Data: simu
 AIC   BIC logLik deviance
 508 529.2   -250      500
Random effects:
 Groups   Name        Variance   Std.Dev.
 alpha    (Intercept) 0.01957449 0.139909
 beta     (Intercept) 0.00080326 0.028342
 Residual             0.06888192 0.262454
Number of obs: 1500, groups: alpha, 100; beta, 5

Fixed effects:
            Estimate Std. Error t value
(Intercept)  1.13506    0.01907   59.52

The Std.Dev of alpha, beta and residual is far way from the true value
in my simulation. While glmmPQL will give a better result.

But I still need to find the variance matrix for variance components
and the dispersion parameter, any suggestions?

Yue Yu


On Fri, Jun 24, 2011 at 09:52, Dennis Murphy <djmuser at gmail.com> wrote:
> 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
>>
>