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 <ekstrom@dina.kvl.dk>
Dept of Mathematics and Physics, KVL
Thorvaldsensvej 40
DK-1871 Frederiksberg C
Denmark
Phone:[+45] 3528 2341
Fax: [+45] 3528 2350
