Greetings,
I've tried to update the power.t.test function to allow for different sample sizes and
sample variances in the case of a two-sample t test. I'd gladly update
the corresponding help page if the code is to replace the current power.t.test
function. The modified power.t.test code is included below.
The changes are as follows:
- Added three new arguments to the function
* ratio : the ratio between group sizes (defaults to 1)
* sd.ratio : the ratio between group sd's (defaults to 1)
* df.method : the method used to calculate the df (defaults to
Welch's/Satterthwaite as in the t.test function)
The three new arguments are unly used for the unpaired two-sample case (i.e., when
type="two.sample"), and they have default values such that function calls in
"old code" should work. However, I've changed the output for the group sizes and sd's
so they now return a 2-vector with the corresponding group sizes and
sd's. For example
# Old code
> power.t.test(power=.8, delta=1, sd=1)
Two-sample t test power calculation
n = 16.71477
delta = 1
sd = 1
sig.level = 0.05
power = 0.8
alternative = two.sided
NOTE: n is number in *each* group
# New code
> power.t.test(power=.8, delta=1, sd=1)
Two-sample t test power calculation
n = 16.71477, 16.71477
delta = 1
sd = 1, 1
sig.level = 0.05
power = 0.8
alternative = two.sided
df.method = welch
NOTE: n is vector of number in each group
This change of output may break some existing code and is not necessarily a "good
idea". Another (possibly better) solution to this would be to return
both ratios and then let print.power.htest present the result in the correct way.
Also, the df.method in the result list could be removed for the one-sample
and the paired case.
Please comment,
Claus Ekstrøm
power.t.test <-
function (n = NULL, delta = NULL, sd = 1, sig.level = 0.05, power = NULL,
ratio = 1, sd.ratio = 1,
type = c("two.sample", "one.sample", "paired"),
alternative = c("two.sided", "one.sided"),
df.method = c("welch", "classical"),
strict = FALSE)
{
type <- match.arg(type)
if (type == "two.sample") {
if (sum(sapply(list(n, delta, sd, power, sig.level, ratio, sd.ratio), is.null)) !=
1)
stop("exactly one of n, delta, sd, power, sig.level, ratio and sd.ratio must be
NULL")
if (!is.null(ratio) && ratio <= 0)
stop("ratio between group sizes must be positive")
if (!is.null(sd.ratio) && sd.ratio <= 0)
stop("sd.ratio between group sd's must be positive")
}
else {
if (sum(sapply(list(n, delta, sd, power, sig.level), is.null)) != 1)
stop("exactly one of n, delta, sd, power, and sig.level must be NULL")
}
alternative <- match.arg(alternative)
df.method <- match.arg(df.method)
tsample <- switch(type, one.sample = 1, two.sample = 2, paired = 1)
tside <- switch(alternative, one.sided = 1, two.sided = 2)
if (tside == 2 && !is.null(delta))
delta <- abs(delta)
p.body <- quote({
nu <- switch(tsample, n-1, switch(df.method, welch=(sd^2/n +
(sd*sd.ratio)^2/(n*ratio))^2/((sd^2/n)^2/(n-1) +
((sd*sd.ratio)^2/(ratio*n))^2/(n*ratio-1)),
classical=(1+ratio)*n-2))
pt(qt(sig.level/tside, nu, lower = FALSE), nu, ncp = switch(tsample,
sqrt(n/tsample), sqrt(n/(1+sd.ratio/ratio))) * delta/sd, lower = FALSE)
})
if (strict & tside == 2)
p.body <- quote({
nu <- switch(tsample, n-1, switch(df.method, welch=(sd^2/n +
(sd*sd.ratio)^2/(n*ratio))^2/((sd^2/n)^2/(n-1) +
((sd*sd.ratio)^2/(ratio*n))^2/(n*ratio-1)),
classical=(1+ratio)*n-2))
qu <- qt(sig.level/tside, nu, lower = FALSE)
ncp <- switch(tsample, sqrt(n/tsample), sqrt(n/(1+sd.ratio/ratio))) * delta/sd
pt(qu, nu, ncp = ncp, lower = FALSE) +
pt(-qu, nu, ncp = ncp, lower = TRUE)
})
if (is.null(power))
power <- eval(p.body)
else if (is.null(n))
n <- uniroot(function(n) eval(p.body) - power, c(2, 1e+07))$root
else if (is.null(sd))
sd <- uniroot(function(sd) eval(p.body) - power, delta * c(1e-07, 1e+07))$root
else if (is.null(delta))
delta <- uniroot(function(delta) eval(p.body) - power, sd * c(1e-07, 1e+07))$root
else if (is.null(sig.level))
sig.level <- uniroot(function(sig.level) eval(p.body) - power, c(1e-10, 1 -
1e-10))$root
else if (is.null(ratio))
ratio <- uniroot(function(ratio) eval(p.body) - power, c(2/n, 1e+07))$root
else if (is.null(sd.ratio))
sd.ratio <- uniroot(function(sd.ratio) eval(p.body) - power, c(1e-07, 1e+07))$root
else stop("internal error")
NOTE <- switch(type, paired = "n is number of *pairs*, sd is std.dev. of
*differences* within pairs",
two.sample = "n is vector of number in each group", NULL)
n <- switch(type, paired=n, two.sample=c(n, n*ratio), one.sample=n)
sd <- switch(type, paired=sd, two.sample=c(sd, sd*sd.ratio), one.sample=sd)
METHOD <- paste(switch(type, one.sample = "One-sample", two.sample = "Two-sample",
paired = "Paired"), "t test power calculation")
structure(list(n = n, delta = delta, sd = sd, sig.level = sig.level,
power = power, alternative = alternative, note = NOTE,
method = METHOD), class = "power.htest")
}
--
*****************************************
Claus Thorn Ekstrøm <[EMAIL PROTECTED]>
Dept of Mathematics and Physics, KVL
Thorvaldsensvej 40
DK-1871 Frederiksberg C
Denmark
Phone:[+45] 3528 2341
Fax: [+45] 3528 2350
______________________________________________
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-devel
