Grokbase Groups R r-help April 2011
FAQ

[R] power of 2 way ANOVA with interaction

Timothy Spier
Apr 3, 2011 at 11:10 pm
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.
reply

Search Discussions

2 responses

  • Peter dalgaard at Apr 4, 2011 at 9:03 am

    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úctor(1:2),varúctor(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,ncpC.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.mes at cbs.dk Priv: PDalgd at gmail.com
  • Greg Snow at Apr 4, 2011 at 4:02 pm
    You can use simulation:

    1. decide what you think your data will look like
    2. decide how you plan to analyze your data
    3. write a function that simulates a dataset (common arguments include sample size(s) and effect sizes) then analyzes the data in your planned manner and returns the p-value(s) or other statistic(s) of interest
    4. run the function from 3 a bunch (1000 or more) times, the replicate function is useful for this (progress bars can also be useful)
    5. the proportion of times that the results are significant is your estimate of power

    --
    Gregory (Greg) L. Snow Ph.D.
    Statistical Data Center
    Intermountain Healthcare
    greg.snow at imail.org
    801.408.8111

    -----Original Message-----
    From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-
    project.org] On Behalf Of Timothy Spier
    Sent: Sunday, April 03, 2011 5:11 PM
    To: R Help
    Subject: [R] power of 2 way ANOVA with interaction


    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.

    [[alternative HTML version deleted]]

    ______________________________________________
    R-help at 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.

Related Discussions

Discussion Navigation
viewthread | post