Example R programs and commands

		      Student-t test power demo
    

# All lines preceded by the "#" character are my comments.
# All other left-justified lines are my input.
# All other indented lines are the R program output.


# Truth: X is a normal random variable with mean 0.34
#   and standard deviation 2
# Hypothesis:  X is a normal random variable with mean 0
#   and unknown standard deviation
#
# How likely are we to reject the (false) hypothesis if we
# can only afford to take 100 samples of X?

### Preliminary investigation:

# Simulate sampling X 100 times:
x <- rnorm(100, mean=0.34, sd=2 )

# Student-t test of "H0:mu=0" from first principles:
s <- sd(x);             # estimated standard deviation of X from samples
se <- s/sqrt(100);      # estimated standard error for a 100-sample average
t <- abs(mean(x)-0)/se; # t-statistic for null hypothesis "H0:mu=0"
2*pt(t,df=99,lower.tail=FALSE)   # p-value for the test.
# ...so given a significance level alpha (such as alpha=0.05), we
#       Do not reject H0 if p>alpha ;
#       Reject H0        if p<alpha ;

# Same Student-t test using an R function:
t.test(x)
# We may extract just the p-value from the output data structure:
t.test(x)$p.value

# Repeat the t-test with several sets of 100 samples:
t.test(rnorm(100, mean=0.34, sd=2))$p.value
t.test(rnorm(100, mean=0.34, sd=2))$p.value
t.test(rnorm(100, mean=0.34, sd=2))$p.value
t.test(rnorm(100, mean=0.34, sd=2))$p.value
t.test(rnorm(100, mean=0.34, sd=2))$p.value
t.test(rnorm(100, mean=0.34, sd=2))$p.value
t.test(rnorm(100, mean=0.34, sd=2))$p.value
t.test(rnorm(100, mean=0.34, sd=2))$p.value
#  ...and note that we sometimes reject H0 and sometimes do not,
#     using a p-value threshold of alpha=0.05


### Simulate repeating the t-test many times to determine power

# Run 1000 t-tests, keeping the p-values:
ps <- vector(length=1000);
for(i in 1:1000) ps[i]<- t.test(rnorm(100, mean=0.34, sd=2))$p.value

# Plot the histogram of p-values to see how many are significant:
hist(ps)

# The (empirical) cumulative distribution function of the p-values,
# or ecdf(ps), shows what fraction of the p-values is below any alpha:
plot(ecdf(ps))

# The output of ecdf() is the cdf function, which may be evaluated at a
#  particular alpha such as 0.05:
ecdf(ps)(0.05)
# ...which estimates the probability of correctly rejecting the false H0.

# Combine the steps and repeat several times:
for(i in 1:1000){ps[i]<-t.test( rnorm(100,mean=0.34,sd=2) )$p.value;}; ecdf(ps)(0.05)
for(i in 1:1000){ps[i]<-t.test( rnorm(100,mean=0.34,sd=2) )$p.value;}; ecdf(ps)(0.05)
for(i in 1:1000){ps[i]<-t.test( rnorm(100,mean=0.34,sd=2) )$p.value;}; ecdf(ps)(0.05)
for(i in 1:1000){ps[i]<-t.test( rnorm(100,mean=0.34,sd=2) )$p.value;}; ecdf(ps)(0.05)
for(i in 1:1000){ps[i]<-t.test( rnorm(100,mean=0.34,sd=2) )$p.value;}; ecdf(ps)(0.05)
for(i in 1:1000){ps[i]<-t.test( rnorm(100,mean=0.34,sd=2) )$p.value;}; ecdf(ps)(0.05)
for(i in 1:1000){ps[i]<-t.test( rnorm(100,mean=0.34,sd=2) )$p.value;}; ecdf(ps)(0.05)
for(i in 1:1000){ps[i]<-t.test( rnorm(100,mean=0.34,sd=2) )$p.value;}; ecdf(ps)(0.05)
# ... and notice that about 39% of the t-tests correctly reject the false H0.




### R function to compute power

# Setup:
#     One-sample test of the mean of a normal population,
#     n=100 elements in the sample,
#     true mean=0.34 differs from H0:mu=0 by 0.34,
#     true population standard deviation is 2,
#     significance level to reject H0 is alpha=0.05

power.t.test(n=100, delta=0.34, sd=2, sig.level=0.05, type="one.sample")

# ...which gives
#
#         "power = 0.3910995"
#
# indicating that about 39% of such tests should correctly reject H0.


# REMARK:  since the population standard deviation is unknown, but the
#  power.t.test() function needs it, it must be estimated somehow, making this
# latter method less general than the simulation method of estimating power.




### R function to compute sample size needed for a given power:

# To get 39% power:
power.t.test(power=0.39, delta=0.34, sd=2, sig.level=0.05, type="one.sample")

# To get 50% power:
power.t.test(power=0.50, delta=0.34, sd=2, sig.level=0.05, type="one.sample")

# To get 90% power:
power.t.test(power=0.90, delta=0.34, sd=2, sig.level=0.05, type="one.sample")

# To get 95% power:
power.t.test(power=0.95, delta=0.34, sd=2, sig.level=0.05, type="one.sample")

# To get 99% power:
power.t.test(power=0.99, delta=0.34, sd=2, sig.level=0.05, type="one.sample")

# To get 99.9999% power:
power.t.test(power=0.999999, delta=0.34, sd=2, sig.level=0.05, type="one.sample")

# REMARK: see what happens when you specify power=1