[R-sig-ME] Power calculations for random effect models
Christos Hatzis
christos.hatzis at nuverabio.com
Fri Feb 20 21:20:40 CET 2009
Thanks, Dimitri. I'll take a look a these references.
I am still hoping to figure out how to use the results of MCMC since it
seems to be a relevant tool and is built in lme4.
Thanks again.
-Christos
> -----Original Message-----
> From: r-sig-mixed-models-bounces at r-project.org
> [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf
> Of Dimitris Rizopoulos
> Sent: Friday, February 20, 2009 3:01 PM
> To: christos at nuverabio.com
> Cc: r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] Power calculations for random effect models
>
> you could also have a look at the following recent articles:
>
> Greven, S., Crainiceanu. C., Kuchenhoff, H. and Peters, A. (2008).
> Restricted Likelihood Ratio Testing for Zero Variance
> Components in Linear Mixed Models. Journal of Computational
> and Graphical Statistics 17, 870-891.
>
> Scheipl, F., Greven, S. and Kuchenhoff, H. (2008). Size and
> power of tests for a zero random effect variance or
> polynomial regression in additive and linear mixed models.
> Computational Statistics & Data Analysis 52, 3283-3299.
>
> and the associated package:
> http://cran.r-project.org/web/packages/RLRsim/index.html
>
>
> I hope it helps.
>
> Best,
> Dimitris
>
>
> Christos Hatzis wrote:
> > Hello,
> >
> > I have a nested random effects model of the form
> >
> > Yijk = M + Ai + Bj(i) + Ek(ij)
> >
> > where biopsies (B) are nested within persons (A) and arrays (E) are
> > nested within biopsies.
> > I am interested in estimating the power of a study of a
> given size to
> > determine whether the variance associated with B is
> trivial, i.e. H0:
> > var_b = 0 vs Ha: var_b > 0 at a fixed Type-I error rate.
> >
> > I have written a function to simulate data from this model and used
> > lmer to estimate the random effects as shown below. What is the
> > recommended way of going about testing the null hypothesis? One
> > approach would be to use the estimated standard deviation of the
> > random effect and assuming normality to test whether the
> appropriate
> > CI contains zero. Another would be to use the theoretical
> chi-square
> > distribution for var_b, but that would require the
> appropriate degrees
> > of freedom. Or to use mcmc to estimate the distribution of
> var_b and use this distribution for inference.
> >
> > I would think that the mcmc approach is the recommended one, but I
> > would appreciate any advise on this. In this case, I have
> tried to run
> > the MCMC simulation (which runs fine), but have not been
> able yet to
> > figure out how to use the results of MCMC to test the above
> > hypothesis. Any hints or pointing out materials that
> explain how to
> > use MCMC for inference on random effects will be very much
> appreciated.
> >
> > Thank you.
> > -Christos
> >
> > Christos Hatzis, Ph.D.
> > Nuvera Biosciences, Inc.
> > 400 West Cummings Park
> > Suite 5350
> > Woburn, MA 01801
> > Tel: 781-938-3830
> > www.nuverabio.com <http://www.nuverabio.com/>
> >
> >> fake.dt <- tumor.heter.dt(10, 3, 2, k=0, sa=4, sb=.6,
> se=.2) fake.dt
> > person biopsy array resp bioWper
> > 1 1 1 1 4.308434 1:1
> > 2 1 1 2 4.293186 1:1
> > 3 1 2 1 5.503841 1:2
> > 4 1 2 2 5.640362 1:2
> > 5 1 3 1 5.579461 1:3
> > 6 1 3 2 5.416201 1:3
> > 7 2 1 1 12.479513 2:1
> > 8 2 1 2 12.311426 2:1
> > 9 2 2 1 13.283566 2:2
> > 10 2 2 2 13.138276 2:2
> > 11 2 3 1 12.954277 2:3
> > 12 2 3 2 13.059925 2:3
> > 13 3 1 1 11.649726 3:1
> > 14 3 1 2 11.694472 3:1
> > 15 3 2 1 12.342050 3:2
> > 16 3 2 2 12.316214 3:2
> > 17 3 3 1 12.762695 3:3
> > 18 3 3 2 12.451338 3:3
> > 19 4 1 1 17.168248 4:1
> > 20 4 1 2 17.573315 4:1
> > 21 4 2 1 16.885176 4:2
> > 22 4 2 2 16.819536 4:2
> > 23 4 3 1 15.924120 4:3
> > 24 4 3 2 16.038491 4:3
> > 25 5 1 1 7.981279 5:1
> > 26 5 1 2 7.777929 5:1
> > 27 5 2 1 6.701801 5:2
> > 28 5 2 2 6.978284 5:2
> > 29 5 3 1 7.136417 5:3
> > 30 5 3 2 7.378498 5:3
> > 31 6 1 1 21.499796 6:1
> > 32 6 1 2 21.397173 6:1
> > 33 6 2 1 22.551116 6:2
> > 34 6 2 2 22.521006 6:2
> > 35 6 3 1 21.027998 6:3
> > 36 6 3 2 21.416872 6:3
> > 37 7 1 1 13.312857 7:1
> > 38 7 1 2 13.559062 7:1
> > 39 7 2 1 13.753016 7:2
> > 40 7 2 2 13.608642 7:2
> > 41 7 3 1 13.556446 7:3
> > 42 7 3 2 13.400672 7:3
> > 43 8 1 1 12.758548 8:1
> > 44 8 1 2 12.486574 8:1
> > 45 8 2 1 13.388409 8:2
> > 46 8 2 2 13.263029 8:2
> > 47 8 3 1 12.991308 8:3
> > 48 8 3 2 12.962116 8:3
> > 49 9 1 1 11.420214 9:1
> > 50 9 1 2 11.308010 9:1
> > 51 9 2 1 13.186774 9:2
> > 52 9 2 2 12.778966 9:2
> > 53 9 3 1 11.625079 9:3
> > 54 9 3 2 11.637015 9:3
> > 55 10 1 1 9.108635 10:1
> > 56 10 1 2 8.895658 10:1
> > 57 10 2 1 9.718046 10:2
> > 58 10 2 2 9.552510 10:2
> > 59 10 3 1 8.794893 10:3
> > 60 10 3 2 9.047379 10:3
> >> heter.lmer <- lmer(resp ~ (1 | person) + (1 | bioWper), fake.dt)
> >> heter.lmer
> > Linear mixed model fit by REML
> > Formula: resp ~ (1 | person) + (1 | bioWper)
> > Data: fake.dt
> > AIC BIC logLik deviance REMLdev
> > 98.24 106.6 -45.12 92.85 90.24
> > Random effects:
> > Groups Name Variance Std.Dev.
> > bioWper (Intercept) 0.312327 0.55886
> > person (Intercept) 21.780993 4.66701
> > Residual 0.020633 0.14364
> > Number of obs: 60, groups: bioWper, 30; person, 10
> >
> > Fixed effects:
> > Estimate Std. Error t value
> > (Intercept) 12.368 1.479 8.36
> >
> >
> >> VarCorr(heter.lmer)[["bioWper"]]
> > (Intercept)
> > (Intercept) 0.3123273
> > attr(,"stddev")
> > (Intercept)
> > 0.5588625
> > attr(,"correlation")
> > (Intercept)
> > (Intercept) 1
> >
> >> heter.mcmc <- mcmcsamp(heter.lmer, n=1000)
> >> str(heter.mcmc)
> > Formal class 'merMCMC' [package "lme4"] with 9 slots
> > ..@ Gp : int [1:3] 0 30 40
> > ..@ ST : num [1:2, 1:1000] 3.89 32.49 1.44 27.06 1.24 ...
> > ..@ call : language lmer(formula = resp ~ (1 | person)
> + (1 | bioWper),
> > data = fake.dt)
> > ..@ deviance: num [1:1000] 92.9 92.9 114.5 119.1 134 ...
> > ..@ dims : Named int [1:18] 2 60 1 40 1 2 0 1 2 5 ...
> > .. ..- attr(*, "names")= chr [1:18] "nt" "n" "p" "q" ...
> > ..@ fixef : num [1, 1:1000] 12.4 14.1 11.9 12.6 12.6 ...
> > .. ..- attr(*, "dimnames")=List of 2
> > .. .. ..$ : chr "(Intercept)"
> > .. .. ..$ : NULL
> > ..@ nc : int [1:2] 1 1
> > ..@ ranef : num[1:40, 0 ]
> > ..@ sigma : num [1, 1:1000] 0.144 0.128 0.126 0.212 0.228 ...
> >
> > _______________________________________________
> > R-sig-mixed-models at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
>
> --
> Dimitris Rizopoulos
> Assistant Professor
> Department of Biostatistics
> Erasmus University Medical Center
>
> Address: PO Box 2040, 3000 CA Rotterdam, the Netherlands
> Tel: +31/(0)10/7043478
> Fax: +31/(0)10/7043014
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>
More information about the R-sig-mixed-models
mailing list