[R-sig-ME] Power calculations for random effect models

[Christos Hatzis]

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 ...