[R] Power analysis in hierarchical models

[Tom Wilding]

Dear All
I am attempting some power analyses, based on simulated data.
My experimental set up is thus:
Bleach: main effect, three levels (control, med, high),  Fixed.
Temp: main effect, two levels (cold, hot), Fixed.
Main effect interactions, six levels (fixed)
For each main-effect combination I have three replicates.
Within each replicate I can take varying numbers of measurements
(response variable = Growth (of marine worms)) but, for this example,
assume eight).  (I’m interested in changing this to see if the
experimental power changes much).  
Total size = 3 x 2 x 3 x 8 = 144
The script thus far goes:
=========== start of script =================
library(lme4)
#Data frame structure
Bleach=rep(c("Control","Med","High"),each=48)
Temp=  rep(rep(c("Cold","Hot"),each=24),3)
Rep=  (rep(rep(rep(c("1","2","3"),each=8),2),3))
Ind= (rep(rep(rep(c(1:8),3),2),3))#not required for stats

#Fake data (based on pilot studies), only showing a single main effect
(bleach)
Growth=c( rnorm(48,3.27,0.77),rnorm(48,3.21,0.77),rnorm(48,3.64,1.17))
fake2=data.frame(Bleach,Temp,Rep,Ind,Growth);head(fake2)
#generate factor level for lmer as per Crawley, page 649
fake2$rep=fake2$Bleach:fake2$Temp:fake2$Rep#rep is used in the lmer
model
with(fake2,table(rep))#check that each rep contains 8 measurements

# run alternative (?equivalent) models 
model1=aov(Growth~Bleach*Temp+Error(Bleach*Temp/Rep),data=fake2);summary(model1)
model2=lmer(Growth~Bleach*Temp+(1|rep),data=fake2);summary(model2)#note:
see above, rep<>Rep!
============ end of script ==========
I'd like to get familiar with using lme4 because it is likely that the
final results of the experiment will be unbalanced (which precludes the
use of aov I think).  The df given by model1 seem to make sense.  Any
guidance on any of the following would be much appreciated:
1. Are model1 and model2 equivalent?
2. For model1 - is the random component correctly specified and is
there a (simple) mechanism to get the appropriate F ratios and P
values?
3. For model2 - again, is the random component correct (probably not) 
and why is the random effect (rep) variance and standard deviations so
low (zero in most iterations)?
4. For both models - how do I isolate (so I can tabulate and create
histograms) the appropriate P and/or t values?  (for model2 - the
‘mer’ object doesn’t seem to contain the t values but maybe
I’m missing something).
Direction to any more generic sources of information regarding power
analysis in hierarchical models would be gladly received.  
Thank you
Tom. 


-------------------------------------------------------------------------
Tom Wilding, MSc, PhD, Dip. Stat.
Scottish Association for Marine Science,
Scottish Marine Institute,
OBAN
Argyll.  PA37 1QA
United Kingdom.
Phone (+44) (0) 1631 559214
Fax (+44) (0) 1631 559001
------------------------------------------------------------------------
+++++++++++++++++++++++++++++++
The Scottish Association for Marine Science (SAMS) is registered in
Scotland as a Company Limited by Guarantee (SC009292) and is a
registered charity (9206).  SAMS has an actively trading wholly owned
subsidiary company: SAMS Research Services Ltd a Limited Company
(SC224404). All Companies in the group are registered in Scotland and
share a registered office at Scottish Marine Institute, Oban Argyll PA37
1QA.

The content of this message may contain personal views which are not
the views of SAMS unless specifically stated.

Please note that all email traffic is monitored for purposes of
security and spam filtering. As such individual emails may be examined
in more detail.
+++++++++++++++++++++++++++++++