On Apr 4, 2011, at 01:10 , Timothy Spier wrote:
>
> I've been searching for an answer to this for a while but no joy. I have a
> simple 2-way ANOVA with an interaction. I'd like to determine the power of
> this test for each factor (factor A, factor B, and the A*B interaction). How
> can I do this in R? I used to do this with "proc Glmpower" in SAS, but I can
> find no analogue in R.
They're not massively hard to do by hand, if you know what you're doing (which,
admittedly is a bit hard to be sure of in this case). The basic structure can
be lifted from power.anova.test and the name of the game is to work out the
noncentrality parameter of the relevant F tests. E.g., lifting an example from
the SAS manual:
> twoway <-
> cbind(expand.grid(ex=factor(1:2),var=factor(1:3)),x=c(14,10,16,15,21,16))
> with(twoway,tapply(x,list(ex,var),mean))
1 2 3
1 14 16 21
2 10 15 16
Now, you have 10 replicates of this with a specified SD of 5. If we do a
"skeleton analysis" of the above table, we get
> anova(lm(x~ex*var,twoway))
Analysis of Variance Table
Response: x
Df Sum Sq Mean Sq F value Pr(>F)
ex 1 16.667 16.6667
var 2 42.333 21.1667
ex:var 2 4.333 2.1667
Residuals 0 0.000
Warning message:
In anova.lm(lm(x ~ ex * var, twoway)) :
ANOVA F-tests on an essentially perfect fit are unreliable
In a 10-fold replication, the SS would be 10 times bigger, and the residual Df
would be 54; also, we need to take the error variance of 5^2 = 25 into account.
The noncentrality for the interaction term is thus 43.333/25 and you can work
out the power as
> pf(qf(.95,2,54),2,54,ncp=43.333/25,lower=F)
[1] 0.1914457
Similarly, the main effect powers are
> pf(qf(.95,2,54),2,54,423.333/25,lower=F)
[1] 0.956741
> pf(qf(.95,1,54),1,54,166.66667/25,lower=F)
[1] 0.7176535
(whatever that means in the presence of interaction, but that is a different
discussion)
--
Peter Dalgaard
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd....@cbs.dk Priv: pda...@gmail.com
______________________________________________
R-help@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.
