-----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
determine whether the variance associated with B is
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
CI contains zero. Another would be to use the theoretical
distribution for var_b, but that would require the
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
the MCMC simulation (which runs fine), but have not been
figure out how to use the results of MCMC to test the above
hypothesis. Any hints or pointing out materials that
use MCMC for inference on random effects will be very much
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,
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)
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 ...