Note on R vs. Stata tutorial

I have noticed that different answers are found using the R functions compared to the Stata functions. I beleive this is due to the exact methods used in the power calculations, but I have not yet been able to discern the difference from the help documentation yet. I will keep working on that and update this document when I know, but for now do not worry that the answers are slightly different.


This R tutorial is also missing the Part 2 of the Stata tutorial, which includes commands to calculate power or sample size in a study invloving survival analysis. The corresponding R packages are more complex and require different data then the Stata function and so will not be covered here. Power calculations for survival analysis are not needed for the extra credit project. Please talk with me if you are interested in the relevant R functions and I can share some tutorials.
-Andrew




Sample Size Calculations in R

This problem set is meant to introduce sample size calculations in R To calculate sample size, in general, you will need to know:

  • What type of variable is the exposure: continuous or binary?
  • How will your outcome be assessed: continuous, binary, time to event (Cox Proportional Hazards model)?
  • What is your estimate of the effect size? (Note: You are usually doing the study to answer this question to begin with, so you probably won’t know the answer to this one. If there is no prior research in this area, then your best guess is as good as any body’s.)
  • What is your \(\alpha\) and power? Are you doing a 1- or 2-sided test? These must be pre-specified.

Note: If you have estimates of your sample size(s), you can use these R functions to calculate the power of your study. In these exercises, we are in the design phase of the study and are trying to figure our how many people we will need to recruit for a specified level of power.




Overview

Power analysis is an important aspect of experimental design. It allows us to determine the sample size required to detect an effect of a given size with a given degree of confidence. Conversely, it allows us to determine the probability of detecting an effect of a given size with a given level of confidence, under sample size constraints. If the probability is too low, the study may be a waste of time and resources, because you won’t have the confidence that a null result is due to no effect of the treatment, rather than from an underpowered study. You’ll be left with study results that you can’t make conclusions from.

The following four numbers are all mathematically inter-related. If you know 3, you can calculate the 4th (or make R do it for you):

  • Sample size
  • Effect size
  • Significance level = P(Type I error) = probability of finding a significant effect when there is not truly an effect
  • Power = 1 - P(Type II error) = probability of finding a significant effect when there is truly an effec




New R Function: Part 1, the “pwr” Package

You must install the pwr package if you’ve never used it before using the following code:

install.packages("pwr")

And then each time you use R, run the following command;

library("pwr")
## Warning: package 'pwr' was built under R version 3.2.5


pwr package power calculation functions:

Function Use
pwr.2p.test two proportions (equal n)
pwr.2p2n.test two proportions (unequal n)
pwr.anova.test balanced one way ANOVA
pwr.chisq.test chi-square test
pwr.f2.test general linear model
pwr.p.test proportion (one sample)
pwr.r.test correlation
pwr.t.test t-tests (one sample, 2 sample, paired)
pwr.t2n.test t-test (two samples with unequal n)

For each of these functions, you enter three of the four quantities (effect size, sample size, significance level, power) and the fourth is calculated.


The significance level defaults to 0.05. Therefore, to calculate the significance level, given an effect size, sample size, and power, use the option “sig.level=NULL”.




Part 1. Two sample comparison of proportions

You want to conduct a cross-sectional study to investigate the association between kombucha tea consumption and liver function. The exposed are regular kombucha tea drinkers (who report 1+ kombucha tea drinks per week) and the unexposed do not drink kombucha regularly. You will measure levels of the liver enzyme GGT to assess liver function and classify levels as “normal” or “abnormally high”. In a pilot study, you found that the proportion of the exposed with abnormal liver enzymes was 0.40 and the proportion of the unexposed with abnormal liver enzymes was 0.25. In this case, n1=exposed participants, n2= unexposed participants. You want alpha=0.05 and a two sided test. You want 80% power.


(Note: The health benefits of kombucha tea have not been rigorously investigated.) To calculate the sample size needed for 80% in two sample comparison of proportions outlined above, use the following command:

pwr.2p.test(h = -0.15, n = NULL, sig.level = 0.05, power = 0.8, 
    alternative = "two.sided")
## 
##      Difference of proportion power calculation for binomial distribution (arcsine transformation) 
## 
##               h = 0.15
##               n = 697.6765
##       sig.level = 0.05
##           power = 0.8
##     alternative = two.sided
## 
## NOTE: same sample sizes

Note the setup. We know the effect size (h=(0.25-0.4)=-0.15), the desired significance level (sig.level = 0.05), and the desired power (power = 0.8). The sample size is set as n=NULL because that is the value we want to calculate. After running the command, we see a sample size of 549.56 is calculated. So 550 subjects are needed (divided equally between exposed and unexposed) to have 80% power to detect an effect size of 15% increased proportion of abnormal liver enzymes.


We also specify the alternative hypothesis. If we don’t know for certain whether the exposure is harmful or protective, we choose a "two-sided" test. If we know exposure is harmful or protective, we can choose the options "greater" or "less", respectively. What happens when "less" is chosen instead of "two-sided"?

pwr.2p.test(h = -0.15, n = NULL, sig.level = 0.05, power = 0.8, 
    alternative = "less")
## 
##      Difference of proportion power calculation for binomial distribution (arcsine transformation) 
## 
##               h = -0.15
##               n = 549.5606
##       sig.level = 0.05
##           power = 0.8
##     alternative = less
## 
## NOTE: same sample sizes

The sample size needed is smaller, because we no longer need the power to detect an effect in either direction.



Part 1. Questions

  1. Restate the null hypothesis (in words) for the kombucha tea/liver enzymes example.






  1. Complete the following Sample Size Table below. Use alpha= 0.025 and a two-sided test for all calculations.
Proportion outcome in exposed Proportion outcome in the unexposed Power Estimated Num Subjects Needed
0.1 0.5 0.90 ?
0.3 0.5 0.90 ?
0.45 0.5 0.90 ?
0.2 0.8 0.80 ?
0.4 0.8 0.80 ?
0.6 0.8 0.80 ?
0.2 0.8 0.70 ?
0.4 0.8 0.70 ?
0.6 0.8 0.70 ?






  1. Comment on the Table above: What happens to the sample size as the proportions in the two groups get closer and closer together? What happens as you increase the power?