For a brief description of the Student t-test and its power curve, please refer to this page.
This is the plot of the power vs. normalized $\delta$.
The discrete data points come from a simulation.
The continuous curve is the theoretical curve based on the non-central Student t-distribution.
## $Id: simpower.R,v 1.4 2009/02/17 07:17:46 marek Exp $
## Copyright (C), Marek Rychlik 2009
##
## Note: Contains parts of John Bear's HW #2 solution.
##
## Produce a simulated power curve for a small sample of size 10.
## Load this entire file, and call power.curve.t(). It produces a graph of
## the power curve. It also returns a list of type I and type II
## errors, and power.
power.curve.t <-
function(sample.size = 10,
n = 100, # Simulation length
alpha = .05, # 1 - confidence
max.delta = 3, # Maximum delta (plot)
delta.step = .05 # Delta increment (plot)
)
{
delta.count <- ceiling(max.delta / delta.step)
## Preallocate space for the results
type.I.results <- numeric(delta.count)
type.II.results <- numeric(delta.count)
pwr <- numeric(delta.count)
df <- 2 * sample.size - 2 # Assuming equal sample size (the
# same as number of replications)
t.alpha.half <- qt(alpha/2, df = df, lower.tail = FALSE)
delta.range <- ( 0:(delta.count - 1) ) * delta.step
for (j in 1:delta.count) {
delta <- delta.range[j]
type.I.error.count <- 0;
type.II.error.count <- 0;
## Select a random sample n times
for(i in 1:n) {
## Select two samples with mean zero
x.insignificant <- rnorm(sample.size)
y.insignificant <- rnorm(sample.size)
## Select one sample with a significantly different mean
y.significant <- rnorm(sample.size) + delta
## Difference of the means
d.insignificant <- mean(x.insignificant) - mean(y.insignificant);
## Shifted difference of the means by the resulting mean
d.significant <- mean(x.insignificant) - mean(y.significant);
## Standard deviation of the difference of the means
sd <- sqrt((var(x.insignificant) * (sample.size - 1) + var(y.significant) * (sample.size - 1)) / df) *
sqrt(1/sample.size + 1/sample.size)
## Evaluate the statistic
statistic.insignificant <- d.insignificant / sd;
statistic.significant <- d.significant / sd;
## Check for type I and type II errors
type.I.error <- abs( statistic.insignificant ) >= t.alpha.half;
type.II.error <- abs( statistic.significant ) < t.alpha.half;
if(type.I.error) {
type.I.error.count = type.I.error.count + 1;
}
if(type.II.error) {
type.II.error.count = type.II.error.count + 1;
}
}
## Find the error rates
## type.I.error.rate <- as.double(type.I.error.count) / n;
## type.II.error.rate <- as.double(type.II.error.count) / n;
type.I.results[j] <- type.I.error.count / n
type.II.results[j] <- type.II.error.count / n
pwr[j] <- 1 - type.II.results[j]
}
## Plot the power curve:
plot(delta.range,
pwr,
ylab = "Power (simulated - discrete, exact - continuous)",
xlab = "Noncentrality parameter",
main = c("t-Test power with ", sample.size, " replications"))
## This function computes the theoretical power
theoretical.power <- function(delta) {
ncp <- delta/sqrt( 1/sample.size + 1/sample.size )
probability.of.type.II.error <- pt(t.alpha.half, ncp = ncp, df = df) -
pt(-t.alpha.half, ncp = ncp, df = df);
return(1 - probability.of.type.II.error)
}
## Add the plot of the theoretical power curve
points(delta.range, theoretical.power(delta.range), type = 'l');
## Return a list of type I error rates, type II rates, and power.
list("type I error rate" = type.I.results,
"type II error rate" = type.II.results,
"power" = pwr)
}