[R] anova power calculations

From: Will Holcomb <wholcomb_at_gmail.com>
Date: Thu, 21 Feb 2008 12:25:42 -0600


I sent a message a couple days ago about doing calculations for power of the ANOVA. Several people got back to me very quickly which I really appreciated.

I'm working now on a similar problem, but instead of a balanced ANOVA, I have an unbalanced one. The first part of the question was:

You assume that the within-population standard deviations all equal 9. You set the Type 1 error rate at = .05. You presume that the population means will have the following values: uA = 17.5, uB = 19, uC = 25, and uD = 20.5. You intend to run 80 subjects in all, with equal n's across all 4 groups. You plan on conducting a one-way ANOVA. Compute your power to reject the null hypothesis under these conditions.

I did:

within.var = 9 ^ 2
means = c(17.5, 19, 25, 20.5)
N = 80
J = length(means)
power.anova.test(groups = J, n = N / J,

                 between.var = var(means),
                 within.var = within.var,
                 sig.level = 0.05)

This gives me 0.6155 which agrees with SAS. The next problem though is:

You have the same Type 1 error rate and make the same assumptions about the population standard deviation and the population means as in part a. You still have 80 subjects in all but now you want to know how power might change by running 10 subjects in groups A, B, and D and 50 subjects in group C. Determine the power under this subject allocation scheme.

For this one I am doing:

# Quantile of the cutoff point in the central F
central.quant = qf(.05, J - 1, N - J, lower.tail = FALSE) weighted.means = data.frame(Mean = means, Weight = c(10, 10, 50, 10))
# Noncentrality parameter for unbalanced ANOVA
noncentral.param = 0
for(i in 1:length(weighted.means$Mean)) {   noncentral.param = (noncentral.param + weighted.means$Weight[i] *

                      (weighted.means$Mean[i] - mean(weighted.means$Mean)) ^
2)
}
noncentral.param = noncentral.param / within.var
# Probability of central quantile in noncentral distribution
noncentral.p = pf(central.quant, J - 1, N - J, noncentral.param, lower.tail= FALSE)
noncentral.p

The logic behind this is in my assignment at:

http://odin.himinbi.org/classes/psy304b/homework_2.xhtml#p2b

This works for a balanced ANOVA and gives the same result as power.anova.test (and SAS). For the unbalanced ANOVA though it is giving me a different result though than SAS, 0.8759455 versus 0.680.

So is there a straightforward way to compute the power of an unbalanced ANOVA? If there isn't, does anyone have any idea what is wrong with my code? The SAS I am comparing it to is:

Data Dep;
Input cue $ mean uneven_weight;
datalines;
A 17.5 1
B 19 1
C 25 5
D 20.5 1
;

proc glmpower;
class cue;
model mean = cue;
weight uneven_weight;
power

        stddev = 9
        alpha = 0.05
        ntotal= 80
        power = .;

run;

Any help would be much appreciated.

Will

        [[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. Received on Thu 21 Feb 2008 - 18:37:01 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Thu 21 Feb 2008 - 20:30:16 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive