[R-sig-ME] glmmPQL Help: Random Effect and Dispersion Parameter
Douglas Bates
bates at stat.wisc.edu
Sat Jun 25 18:53:41 CEST 2011
This time with the enclosure :-)
On Sat, Jun 25, 2011 at 11:52 AM, Douglas Bates <bates at stat.wisc.edu> wrote:
> On Fri, Jun 24, 2011 at 2:28 PM, Yue Yu <parn.yy at gmail.com> wrote:
>> 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.
>
> Hmm, in lme4a the results, shown in the enclosed, are much closer to
> the values used for simulation. However, this example does show a
> deficiency in this implementation of glmer in that the estimate of the
> residual standard deviation is not calculated and not shown here.
>
>> But I still need to find the variance matrix for variance components
>> and the dispersion parameter, any suggestions?
>
> The best suggestion is don't do it. Estimates of variance components
> are not symmetrically distributed (think of the simplest case of the
> estimator of the variance from an i.i.d Gaussian sample, which has a
> chi-squared distribution).
>
>> 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
>>>>
>>>
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>
-------------- next part --------------
R version 2.14.0 Under development (unstable) (2011-06-24 r56210)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-unknown-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> 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
>
> set.seed(12345)
> 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)
>
> library(lme4a)
Loading required package: Matrix
Loading required package: lattice
Attaching package: ?Matrix?
The following object(s) are masked from ?package:base?:
det
Loading required package: minqa
Loading required package: Rcpp
Loading required package: MatrixModels
Attaching package: ?MatrixModels?
The following object(s) are masked from ?package:stats?:
getCall
> summary(glm2 <- glmer(y ~ 1 + (1|alpha) + (1|beta), simu, gaussian(link="log"), verbose=2L))
npt = 4 , n = 2
rhobeg = 0.2 , rhoend = 2e-07
0.020: 7: 512.232;0.600458 0.986477
0.0020: 14: 509.499;0.554646 0.887092
0.00020: 51: 490.819;0.536841 0.0779844
2.0e-05: 60: 490.815;0.539243 0.0788170
2.0e-06: 62: 490.815;0.539294 0.0787849
2.0e-07: 66: 490.815;0.539295 0.0787820
At return
72: 490.81485: 0.539295 0.0787822
npt = 5 , n = 3
rhobeg = 0.2 , rhoend = 2e-07
0.020: 6: 495.473;0.539295 0.200000 1.10770
0.0020: 14: 491.648;0.546593 0.109993 1.11128
0.00020: 24: 490.787;0.539576 0.0785051 1.09634
2.0e-05: 28: 490.787;0.539552 0.0788343 1.09663
2.0e-06: 37: 490.787;0.539619 0.0787000 1.09667
2.0e-07: 53: 490.787;0.539613 0.0787062 1.09676
At return
66: 490.78710: 0.539613 0.0787075 1.09676
Generalized linear mixed model fit by maximum likelihood ['summary.mer']
Family: gaussian
Formula: y ~ 1 + (1 | alpha) + (1 | beta)
Data: simu
AIC BIC logLik deviance
496.7871 512.7268 -245.3935 490.7871
Random effects:
Groups Name Variance Std.Dev.
alpha (Intercept) 0.291182 0.53961
beta (Intercept) 0.006195 0.07871
Number of obs: 1500, groups: alpha, 100; beta, 5
Fixed effects:
Estimate Std. Error z value
(Intercept) 1.09676 0.06532 16.79
>
> proc.time()
user system elapsed
11.780 0.270 12.248
More information about the R-sig-mixed-models
mailing list