[R-sig-ME] Power calculations for random effect models
Greg Snow
Greg.Snow at imail.org
Sat Feb 21 02:12:52 CET 2009
Since you already have a function to simulate the data, it should be fairly simple to use simulations to compare the methods you suggest below, or to do a simple test based on the likelihood ratio, but simulating the cutoff rather than depending on the chi-squared approximation.
See:
http://finzi.psych.upenn.edu/R/Rhelp08/archive/156522.html
for some examples of assessing these types of questions using simulation (I have an updated version of the last example if you want it as well).
Hope this helps,
--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.snow at imail.org
801.408.8111
> -----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 Christos Hatzis
> Sent: Friday, February 20, 2009 11:01 AM
> To: r-sig-mixed-models at r-project.org
> Subject: [R-sig-ME] Power calculations for random effect models
>
> 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
More information about the R-sig-mixed-models
mailing list