[R-sig-ME] Power calculations for random effect models
Dimitris Rizopoulos
d.rizopoulos at erasmusmc.nl
Fri Feb 20 21:00:52 CET 2009
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
More information about the R-sig-mixed-models
mailing list