[R] Power analysis in hierarchical models

[ONKELINX, Thierry]

Dear Tom,

I think you failed to generate simulated outcome from the correct model. Hence the zero variance of your random effects. Here is a better working example.

library(lme4)

fake2 <- expand.grid(Bleach = c("Control","Med","High"), Temp = c("Cold","Hot"), Rep = factor(seq_len(3)), ID = seq_len(8))
fake2$rep <- fake2$Bleach:fake2$Temp:fake2$Rep

SDnoise <- 0.77
SDrep <- 1
FFBleach <- c(3.27,3.21, 3.64)
RFrep <- rnorm(length(levels(fake2$rep)), sd = SDrep)
fake2$Growth <- with(fake2, FFBleach[Bleach] + RFrep[rep] + rnorm(nrow(fake2), sd = SDnoise))

model2 <- lmer(Growth~Bleach*Temp+(1|rep),data=fake2)
str(summary(model2))
summary(model2)@coefs #to extract the t-values

Best regards,

Thierry

PS R-sig-mixed models is a better mailing list for this kind of questions.

> -----Oorspronkelijk bericht-----
> Van: r-help-bounces op r-project.org [mailto:r-help-bounces op r-project.org]
> Namens Tom Wilding
> Verzonden: maandag 5 september 2011 16:17
> Aan: r-help op r-project.org
> Onderwerp: [R] Power analysis in hierarchical models
> 
> 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);sum
> mary(model1)
> model2=lmer(Growth~Bleach*Temp+(1|rep),data=fake2);summary(model2)#no
> te:
> 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.
> +++++++++++++++++++++++++++++++
> 
> ______________________________________________
> R-help op r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.