To answer my own question, the notes that I was drawing from made it seem as
through the mean I was using for my sum of squares was unweighted. In
reality, I needed to weight the contribution of each contributing mean by
the sample size. Once I did that, I got a value that agreed with SAS.
So, the solution for computing power of an unbalanced ANOVA is:
within.var = 9 ^ 2
means = c(17.5, 19, 25, 20.5)
J = length(means)
# Quantile of the cutoff point in the central F
central.quart = qf(.05, J - 1, N - J, lower.tail = FALSE)
weighted.means = data.frame(Mean = means, n = c(10, 10, 50, 10))
N = sum(weighted.means$n)
weighted.mean = weighted.mean(weighted.means$Mean, weighted.means$n)
# Noncentrality parameter for unbalanced ANOVA
noncentral.param = sum(weighted.means$n * (weighted.means$Mean -
weighted.mean) ^ 2) / within.var
# Probability of central quantile in noncentral distribution
noncentral.p = pf(central.quart, J - 1, N - J, noncentral.param,
lower.tail= FALSE)
Will
2008/2/21 Will Holcomb <[EMAIL PROTECTED]>:
> 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@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.
