Skip to content

power analysis for multi-level models

2 messages · Lee Wurm, Greg Snow

Original
Lee Wurm
# 
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
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,