[R-sig-ME] power analysis for multi-level models

[Greg Snow]

Here is some code to get you started (based on some assumptions that may be way off):

library(lme4)

sim1 <- function(bSex=0, bFreq=0, bSF=0, b0=1000, Vsubj=1, Vword=1, Verror=1) {
	Subject <- rep( 1:60, each=50 )
	Word <- rep( 1:50, 60 )
	Sex <- rep(c('M','F'), each=50*30)
	
	#  assume frequency is constant accross word, random from 1-100
	tmp <- sample( 1:100, 50, replace=TRUE )
	Frequency <- tmp[Word]

	# random effects per subject
	S.re <- rnorm(60, 0, sqrt(Vsubj))
	
	# random effects per word
	W.re <- rnorm(50, 0, sqrt(Vword))

	# epsilons
	eps <- rnorm(50*60, 0, sqrt(Verror))

	# put it all together
	ReactionTime <- b0 + bSex*(Sex=='M') + bFreq*Frequency + bSF*(Sex=='M')*Frequency +
		S.re[Subject] + W.re[Word] + eps

	# put into a data frame
	mydata <- data.frame( Subject = paste('s',Subject, sep=''), 
					Word = paste('w', Word, sep=''), Sex=Sex, Frequency=Frequency,
					ReactionTime = ReactionTime)

	# analyze looking at interaction term with LR test
	fit1 <- lmer( ReactionTime ~ (Sex*Frequency) + (1|Subject) + (1|Word), data=mydata)
	fit2 <- lmer( ReactionTime ~ Sex + Frequency + (1|Subject) + (1|Word), data=mydata)
	anova(fit2,fit1)[2,7]
}
	


pb <- winProgressBar(max=100) # or tkProgressBar or txtProgressbar

setWinProgressBar(pb, 0)
out1 <- replicate( 100, {setWinProgressBar(pb, getWinProgressBar(pb)+1);
				sim1( bSex=10, bFreq=2, bSF=0.25, Vsub=4000, Vword=2500, Verror=10000)})
hist(out1)
mean( out1 < 0.05 )

####

Now edit the sim1 function to match your real situation (in any cases that I guessed wrong) and analysis.  Run the simulation for reasonable values (guesses) and see what the power is.  I usually start with about 100 runs just to get a feel in a reasonable amount of time, change the values and rerun the last 4 lines several times.  Once you have the values that you want to use, up the number of simulations (change the progress bar as well) to 1,000 or maybe even 10,000 (start it running at the end of the day, then go home and let it run over night) to get your final values.

You may want to include a table/graph that shows the power for different effects of the interaction term, etc.

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 Lee Wurm
> Sent: Tuesday, January 20, 2009 9:18 AM
> To: r-sig-mixed-models at r-project.org
> Subject: [R-sig-ME] power analysis for multi-level models
> 
> I'm really happy to see the explosion of mixed models in more and more
> areas of research, but am now faced with the problem that grant
> reviewers (and some journals) insist on seeing effect sizes and power
> calculations for these models. I sent my query to the r-help forum and
> got a valuable idea, along with the suggestion that I try you folks.
> 
> Suppose I want to see if men or women show a stronger word frequency
> effect. I have 50 words of varying frequency that I show to 30 men and
> 30 women, who are supposed to decide as quickly as possible whether
> it's a real word. My data object would end up being 3000 lines long,
> and look like this:
> 
> Subject  Word  Sex  Frequency  ReactionTime
> s1 w1 M 23 2543
> s1 w2 M 67 1438
> s1 w3 M 1 8033
> ...
> s60 w50 F 4 1099
> 
> I analyze with:
> 
> lmer(ReactionTime ~ (Sex*Frequency) + (1|Subject) + (1|Word)
> 
> Can anyone help me get started with calculating effect sizes (or even
> better, with attempting power analyses) for such a model? Or with
> giving an explanation about why it's an ill-formed question, in terms
> that non-experts could understand? Old ways die hard, sometimes, and
> grant reviewers control the purse-strings.
> 
> Thanks.
> 
> --Lee
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models