library(tidyverse)
## ── Attaching packages ────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.2.1     ✔ purrr   0.3.2
## ✔ tibble  2.1.3     ✔ dplyr   0.8.3
## ✔ tidyr   1.1.3     ✔ stringr 1.4.0
## ✔ readr   1.3.1     ✔ forcats 0.4.0
## Warning: package 'tidyr' was built under R version 3.6.2
## ── Conflicts ───────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

Generate a P-Value Here’s how to use R for a two-sample t test and retrieve the p-value. Perform a two-sample t test with a sample size of n=25 when H0 is true:

population <- rnorm(1000, mean=10, sd=2)
A <- sample(population, size=25)
B <- sample(population, size=25)
H <- t.test(A, B, alternative="two.sided")
H$p.value
## [1] 0.6868481

EXAMPLE Generate sampling distributions for p-values under different conditions.

Define a Test Iterator …and a helper function to be used later. The iterator allows you to generate an arbitrarily large sampling distribution of p-values for the two-sample t test.

pvalueSampling <- function(popA, popB, n=25, REPS=1000) {
  dist <- array(data=NA, dim=REPS)
  for (i in 1:REPS){
    A <- sample(popA, size=n)
    B <- sample(popB, size=n)
    H <- t.test(A, B, alternative="two.sided")    
    dist[i] <- H$p.value
  }
  return(dist)
}

p.less <- function(A, limit=0.05) {
  return( length(A[A < limit])/length(A) )
}

Test H0:μA=μB

# and test
pvalue.data <- data.frame(pval=pvalueSampling(population, population, REPS=4096))
ggplot(data=pvalue.data, aes(pval)) + 
  geom_histogram(aes(y=..density..), bins=64, 
                 col="orange", fill="green") +
  geom_density(col="red") +
  labs(title="Sampling Distribution of p-values", x = "p-value", y = "frequency") +
  theme_minimal(12)

p.reject <- p.less(pvalue.data$pval)
names(p.reject) <- "observed Type I error"
p.reject
## observed Type I error 
##            0.04760742

Test Ha:μA≠μB

population.A <- rnorm(1000, mean=10, sd=2)
population.B <- rnorm(1000, mean=11, sd=2) # so Cohen's D = 0.5

pvalue.data <- data.frame(pval=pvalueSampling(population.A, population.B, REPS=4096))
ALPHA <- 0.05

ggplot(data=pvalue.data, aes(pval)) + 
  geom_histogram(aes(y=..density..), bins=64, 
                 col="orange", fill="green") +
  geom_density(col="red") +
  geom_vline(xintercept=ALPHA, color="red") +
  labs(title="Sampling Distribution of p-values", x = "p-value", y = "frequency") +
  theme_minimal(12)

p.reject <- p.less(pvalue.data$pval)
names(p.reject) <- "observed power"
p.reject
## observed power 
##      0.2770996

Exercise 1. What was the estimated power of the test in the second example above? (Remember: it was a two-sided test, so use α/2 rather than α.)

population.A <- rnorm(1000, mean=10, sd=2)
population.B <- rnorm(1000, mean=11, sd=2) # so Cohen's D = 0.5

pvalue.data <- data.frame(pval=pvalueSampling(population.A, population.B, REPS=4096))
ALPHA <- 0.025

ggplot(data=pvalue.data, aes(pval)) + 
  geom_histogram(aes(y=..density..), bins=64, 
                 col="orange", fill="green") +
  geom_density(col="red") +
  geom_vline(xintercept=ALPHA, color="red") +
  labs(title="Sampling Distribution of p-values", x = "p-value", y = "frequency") +
  theme_minimal(12)

p.reject <- p.less(pvalue.data$pval)
names(p.reject) <- "observed power"
p.reject
## observed power 
##      0.4384766

Exercise 2. Modify the example to compare sample means when Cohens’s D is 1.0 and 1.5. How does the observed power change?

population.A <- rnorm(1000, mean=9.5, sd=2)
population.B <- rnorm(1000, mean=11.5, sd=2) # so Cohen's D = 1.0

pvalue.data <- data.frame(pval=pvalueSampling(population.A, population.B, REPS=4096))
ALPHA <- 0.05

ggplot(data=pvalue.data, aes(pval)) + 
  geom_histogram(aes(y=..density..), bins=64, 
                 col="orange", fill="green") +
  geom_density(col="red") +
  geom_vline(xintercept=ALPHA, color="red") +
  labs(title="Sampling Distribution of p-values", x = "p-value", y = "frequency") +
  theme_minimal(12)

p.reject <- p.less(pvalue.data$pval)
names(p.reject) <- "observed power"
p.reject
## observed power 
##       0.940918
population.A <- rnorm(1000, mean=9.5, sd=2)
population.B <- rnorm(1000, mean=11.5, sd=2) # so Cohen's D = 1.5

pvalue.data <- data.frame(pval=pvalueSampling(population.A, population.B, REPS=4096))
ALPHA <- 0.05

ggplot(data=pvalue.data, aes(pval)) + 
  geom_histogram(aes(y=..density..), bins=64, 
                 col="orange", fill="green") +
  geom_density(col="red") +
  geom_vline(xintercept=ALPHA, color="red") +
  labs(title="Sampling Distribution of p-values", x = "p-value", y = "frequency") +
  theme_minimal(12)

p.reject <- p.less(pvalue.data$pval)
names(p.reject) <- "observed power"
p.reject
## observed power 
##      0.9018555

Exercise 3. Reduce the sample sizes to n=9, and compare to one of your previous cases. How does the observed power change?

pvalueSampling <- function(popA, popB, n=9, REPS=1000) {
  dist <- array(data=NA, dim=REPS)
  for (i in 1:REPS){
    A <- sample(popA, size=n)
    B <- sample(popB, size=n)
    H <- t.test(A, B, alternative="two.sided")    
    dist[i] <- H$p.value
  }
  return(dist)
}

p.less <- function(A, limit=0.05) {
  return( length(A[A < limit])/length(A) )
}
population.A <- rnorm(1000, mean=10, sd=2)
population.B <- rnorm(1000, mean=11, sd=2) # so Cohen's D = 0.5

pvalue.data <- data.frame(pval=pvalueSampling(population.A, population.B, REPS=4096))
ALPHA <- 0.05

ggplot(data=pvalue.data, aes(pval)) + 
  geom_histogram(aes(y=..density..), bins=64, 
                 col="orange", fill="green") +
  geom_density(col="red") +
  geom_vline(xintercept=ALPHA, color="red") +
  labs(title="Sampling Distribution of p-values", x = "p-value", y = "frequency") +
  theme_minimal(12)

p.reject <- p.less(pvalue.data$pval)
names(p.reject) <- "observed power"
p.reject
## observed power 
##      0.2016602

Exercise 4. Perform a similar comparison where the underlying populations are not normally distributed. You might try either a uniform distribution( runif(n, min, max)) or a gamma distribution (rgamma(n, shape, rate)).

population.A <- runif(25, min=0, max=.05)
population.B <- runif(25, min=0, max=.05) 

pvalue.data <- data.frame(pval=pvalueSampling(population.A, population.B, REPS=4096))
ALPHA <- 0.05

ggplot(data=pvalue.data, aes(pval)) + 
  geom_histogram(aes(y=..density..), bins=64, 
                 col="orange", fill="green") +
  geom_density(col="red") +
  geom_vline(xintercept=ALPHA, color="red") +
  labs(title="Sampling Distribution of p-values", x = "p-value", y = "frequency") +
  theme_minimal(12)

p.reject <- p.less(pvalue.data$pval)
names(p.reject) <- "observed power"
p.reject
## observed power 
##     0.02148438