Visualizing Power

Primer on Statistical Significance

Null-hypothesis significant testing (NHST) is a controversial approach to social science research. Although I will not visit the concerns in full, it’s important to have an understanding of concepts related to NHST so you can fully understand why this approach is criticized and its flaws are..well… flaws. One ubiquitous, yet, misunderstood concept is p-values. For you own knowledge:

…beginning with the assumption that the true effect is zero (i.e., the null hypothesis is true), a p-value indicates the proportion of test statistics, computed from hypothetical random samples, that are as extreme, or more extreme, then the test statistic observed in the current study.

(Spence & Stanley, 2018)

or, stated another way:

The smaller the p-value, the greater the statistical incompatibility of the data with the null hypothesis, if the underlying assumptions used to calculate the p-value hold. This incompatibility can be interpreted as casting doubt on or providing evidence against the null hypothesis or the underlying assumptions.

(Wesserstein & Lazar, 2016)

Simply, a p-values indicates the probability of your, or more extreme, test statistic, assuming the null: \(p(data | null)\). For ease, consider the following table:

When we conduct NHST, we are assuming that the null hypothesis is indeed true, despite never truly knowing this. As noted, a p-value indicates the likelihood of your or more extreme data given the null.

Simulating Data

We will simulate some data for this visualization. Consider a researcher who is interested in the association substance use (SU) and suicidal behaviors (SB) in Canadian high school students. Let’s assume that the true association between SU and SB is \(r = .3\); we do not know that this is the true value. Let’s simulate the data with the assumption that there are 100,000 high school students

Note: I will include my r script for anyone interested in following along or adjusting it themselves.

library(tidyverse)
library(simstudy)

#Setting the seed will allow you to get the same results, if you are following along.
set.seed(2434)
mat <- genCorMat(2, c(.3))
mean = c(100, 100)
sd = c(15, 15)
data <- genCorData(100000, mu=mean, sigma = sd, corMatrix = mat, cnames=c("SU", "SB"))

Using the above, we have generated a data set of the SU and SB scores of 100,000 high schoolers. Furthermore, we have set the population correlation to .3. First, let’s plot the data to get an idea of how this relationship looks:

ggplot(data, aes(x=SB, y=SU))+
  geom_point(alpha=.05, color="grey31")+
  geom_smooth(method="lm", color="red")+
  labs(x="Suicidal Behaviors", y="Substance Use")+
  coord_cartesian(xlim=c(50, 150), ylim=c(50,150))+
  theme_classic(15)

## `geom_smooth()` using formula 'y ~ x'

So, it appears that as substance use increases, so do suicidal behaviors We can check to make sure that the population correlation actually came out to be = .3, as we intended when we simulated the data:

cor.test(data$SU, data$SB)

## 
##  Pearson's product-moment correlation
## 
## data:  data$SU and data$SB
## t = 99.473, df = 99998, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2944180 0.3056978
## sample estimates:
##       cor 
## 0.3000684

So, the correlation in the population is \(\rho = .3000684\) (close enough). Although we aren’t so omniscient in the real world, we will use this known population correlation to determine the probability that our proposed study will correctly reject a null hypothesis, which, as you remember, states there is no relationship between SU and SB.

Before conducting our study, we will conduct a power analysis to determine an appropriate sample size required to adequately power our study. We want to have a good probability of rejecting the null, if it were false (here, we know this is the case). To calculate our required sample size, we require: i) our \(\alpha\) criterion level, ii) our desired power (\(1-\beta\)), and iii) our predicted effect size. We will use the standard for our criterion and power, i) \(\alpha = .05\), ii) 1-= .8, and we know our true effect, iii) = .3000684. Although we know the population correlation, in real-world research we would need ot find a good estimate of the effect size, which is typically through reading the literature for effect sizes in similar populations. We can use the pwr package to conduct our power analysis. This package is very useful; you can insert any three of the required four pieces of information to calculate the missing piece. For our analysis, we write:

library(pwr)
pwr.r.test(r=.3000684, sig.level=.05, power=.8)

## 
##      approximate correlation power calculation (arctangh transformation) 
## 
##               n = 84.03388
##               r = 0.3000684
##       sig.level = 0.05
##           power = 0.8
##     alternative = two.sided

The results of this suggest that we need a sample of about 84 people (n = 84.03388) to achieve our desired power (power = 0.8), using our known popultion correlation (r = 0.3000684) and (sig.level = 0.05) value. What does this mean? It means that if the true population effect/association/relationship was = .3, then 80% of all hypothetical studies, using a sample size of n = 84 drawn from this population, will yield \(p < \alpha\).

Visualizing Power

We can create a histogram that plots the results of many random studies drawn from our population to determine which meet \(p < \alpha\). First, I will create a function that can easily draw many random samples from our population and provide the resultant correlation value.

sample_cor <- function(x){
  data1 <- sample_n(data, x)
cor <- cor(data1$SU, data1$SB)
print(cor)
}

This function will retrieve a sample of x from our population data, which we created above. We can use this function to draw many samples. In the following, we will draw 2000 samples of n = 84 (these are evident in the script: rep(84, each=2000))

set.seed(5664131)
correlations <- data.frame(correlations=lapply(rep(84, each=2000), sample_cor))
correlations <- data.frame(t(correlations))
names(correlations) <- c("correlation")
rownames(correlations) <- 1:2000

We now have a dataframe, correlations, that consist of the 2000 correlation coefficients from our 2000 random samples.

head(correlations)

##   correlation
## 1   0.3182940
## 2   0.5499953
## 3   0.2369162
## 4   0.2818166
## 5   0.4269413
## 6   0.3988922

So, sample 1 resulted in a correlation of .356, sample 2 a correlation of .273, sample 3 a correlation of .225, and so on…

We can plot these as a histogram. However, this doesn’t tell us anything about p-values. We can calculate the required correlation coefficient that would results in \(p < \alpha = .05\). You can use a critical correlation table to find associated correlation values, or a website such as here. Once we fill in our values, we see that a correlation of at least \(r =|.215|\) is required to obtain \(p < \alpha = .05\) when \(n=84\).

First, we will will plot a histogram of all of our calculated p-values.

ggplot(correlations, aes(correlation))+
  geom_histogram(fill="grey", color="black", binwidth = .025)+
  coord_cartesian(xlim=c(-1, 1))+
  scale_x_continuous(breaks=seq(-1, 1, by=.2))+
  scale_y_continuous(breaks=seq(0, 250, by=25))+
  geom_vline(xintercept=0.300684, color="blue", size=1)+
  geom_vline(xintercept=0.215, color="red", size=1, linetype=2)+
  geom_vline(xintercept=-0.215, color="red", size=1, linetype=2)+
  labs(x="Correlation Coefficient", y="Frequency")+
  annotate("rect", xmin = -Inf, xmax = -.215, ymin = -Inf, ymax = Inf, fill = "red", alpha = 0.1) +
  annotate("rect", xmin = .215, xmax = Inf, ymin = -Inf, ymax = Inf, fill = "red", alpha = 0.1) +
  theme_classic(15)

This graph represents the distribution of correlation coefficients for each of our random samples, which were drawn from the 100,000 high schoolers. Notice how they form a seemingly normal distribution around our true population correlation coefficient, \(\rho=.300684\). Although it may seem normal, it is actually skewed because of the bounds of the correlation (i.e., -1 to 1). The red lines and shaded regions represent correlation coefficients that are beyond the r = .215 cut-off, which would result in \(p < .05\). Thus, these are statistically significant results in our test of the null hypothesis (i.e., \(\rho = 0\)). The blue line represents the true population correlation coefficient (\(\rho=.300684\)). Just how many samples scored at or beyond the critical correlation value of \(r = .215\)?

pwr <- correlations %>% 
  mutate(sig=ifelse(correlation>=.215, 1, 
                    ifelse(correlation<=-.215, 1, 0))) %>% 
  mutate(sig=as.factor(sig)) %>% 
  select(sig) 

table(pwr$sig)

## 
##    0    1 
##  413 1587

This tells us that 1613 correlation coefficients were at or beyond the critical value. Do you have any guess what proportion of the total samples that was? Recall that power is the probability that any study will have \(p<\alpha\). Approximately eighty percent (80.65%%) of these studies met that criteria: this was our power! Why is it 80.65% and not 80%? Remember, power refers to a hypothetically infinite number of samples drawn from the population, for any given sample size, effect size, and \(\alpha\). Had we drawn \(\infty\) samples, we would have 80% having \(p<\alpha\). Also, our power analysis suggested a sample size of \(n=84.03388\).

What if we couldn’t recruit 84 participants?

Perhaps we sampled from one high school in a small Canadian city and could only recruit 32 participants. How do you think this would affect our power? Consider that smaller sample sizes are likely to give less precise estimations of population parameters (i.e., the histogram above may be more spread out). However, this influences our critical correlation coefficient, so the red regions shift outward. This will reduce power. Let’s rerun our simulation with 2000 random samples of \(n=32\) to see how it affects out power.

set.seed(372837)
correlations_32 <- data.frame(correlations=lapply(rep(32, each=2000), sample_cor))
correlations_32 <- data.frame(t(correlations_32))
names(correlations_32) <- c("correlation")
rownames(correlations_32) <- 1:2000

Next we can plot the resultant correlation coefficients. Note that the critical correlation coefficient value for \(n = 32\) is larger than for \(n=84\). The new critical value is \(r=349\). The following script will adjust for these differences.

ggplot(correlations_32, aes(correlation))+
  geom_histogram(fill="grey", color="black", binwidth = .025)+
  coord_cartesian(xlim=c(-1, 1))+
  scale_x_continuous(breaks=seq(-1, 1, by=.2))+
  scale_y_continuous(breaks=seq(0, 250, by=25))+
  geom_vline(xintercept=0.300684, color="blue", size=1)+
  geom_vline(xintercept=0.349, color="red", size=1, linetype=2)+
  geom_vline(xintercept=-0.349, color="red", size=1, linetype=2)+
  labs(x="Correlation Coefficient", y="Frequency")+
  annotate("rect", xmin = -Inf, xmax = -.349, ymin = -Inf, ymax = Inf, fill = "red", alpha = 0.1) +
  annotate("rect", xmin = .349, xmax = Inf, ymin = -Inf, ymax = Inf, fill = "red", alpha = 0.1) +
  theme_classic(15)

Hopefully, you can see that despite the distribution still centering around the population correlation of \(.3\), the red region is shifted outward due to a smaller \(n\) (which affects out critical correlation coefficient). This would mean that fewer of the samples result in correlation coefficients that fall in the red regions (i.e., power is lower). Would a formal calculation agree? First, let’s find out how many studies resulted in \(p<\alpha=.05\) and then do a formal power analysis to determine if they are equal.

pwr_32 <- correlations_32 %>% 
  mutate(sig=ifelse(correlation>=.349, 1, 
                    ifelse(correlation<=-.349, 1, 0))) %>% 
  mutate(sig=as.factor(sig)) %>% 
  select(sig) 

table(pwr_32$sig)

## 
##    0    1 
## 1208  792

So, 808 studies yeilded statistically significant results (808/2000 = 40.4%). Would our power align?

pwr.r.test(n=32, sig.level = .05, r=.3)

## 
##      approximate correlation power calculation (arctangh transformation) 
## 
##               n = 32
##               r = 0.3
##       sig.level = 0.05
##           power = 0.3932315
##     alternative = two.sided

Whoa! Close enough for me. Let’s take it one step further and assume we could only get 20 participants.

set.seed(372837)
correlations_20 <- data.frame(correlations=lapply(rep(20, each=2000), sample_cor))
correlations_20 <- data.frame(t(correlations_20))
names(correlations_20) <- c("correlation")
rownames(correlations_20) <- 1:2000

Next we can plot the resultant correlation coefficients. Note that the critical correlation coefficient value for \(n = 20\) is larger than for \(n=32\) (remember that the null distribution or samples would be spread out more for smaller sample sizes). The new critical value is \(r=444\). The following script will adjust for these differences.

ggplot(correlations_20, aes(correlation))+
  geom_histogram(fill="grey", color="black", binwidth = .025)+
  coord_cartesian(xlim=c(-1, 1))+
  scale_x_continuous(breaks=seq(-1, 1, by=.2))+
  scale_y_continuous(breaks=seq(0, 250, by=25))+
  geom_vline(xintercept=0.300684, color="blue", size=1)+
  geom_vline(xintercept=0.444, color="red", size=1, linetype=2)+
  geom_vline(xintercept=-0.444, color="red", size=1, linetype=2)+
  labs(x="Correlation Coefficient", y="Frequency")+
  annotate("rect", xmin = -Inf, xmax = -.444, ymin = -Inf, ymax = Inf, fill = "red", alpha = 0.1) +
  annotate("rect", xmin = .444, xmax = Inf, ymin = -Inf, ymax = Inf, fill = "red", alpha = 0.1) +
  theme_classic(15)

Hopefully, you can see that despite the distribution still centering around our population correlation of .3, the red regions are further outward than the previous two examples. Again, let’s see how many studies resulted in \(p<\alpha=.05\) and then do a formal power analysis to determine if they are equal.

pwr_20 <- correlations_20 %>% 
  mutate(sig=ifelse(correlation>=.444, 1, 
                    ifelse(correlation<=-.444, 1, 0))) %>% 
  mutate(sig=as.factor(sig)) %>% 
  select(sig) 

table(pwr_20$sig)

## 
##    0    1 
## 1504  496

So, 511 studies yielded statistically significant results (511/2000 = 25.5%). Would our power align?

pwr.r.test(n=20, sig.level = .05, r=.3)

## 
##      approximate correlation power calculation (arctangh transformation) 
## 
##               n = 20
##               r = 0.3
##       sig.level = 0.05
##           power = 0.2559237
##     alternative = two.sided

So, out of 2000 random samples from our population, 25.5% had \(p<\alpha=.05\) and our power analysis suggests that 25.6% would. Again, given infinite samples, there would be 25.59237…% will \(p<\alpha=.05\). In sum, as sample size decreases and all other things are help constant, our power decreases. If your sample is too small, you will be unlikely to reject \(H_{o}\) even if a true effect exists. Thus, it is important to ensure you have an adequately powered study. Plan ahead. Otherwise, even if a population effect exists, you may not conclude that through NHST’. Or, maybe a non-NHST approach is best ;).

What if our population effect was larger than \(\rho = .3\)?

If our population correlation was higher, what do you think would happen to power? Specifically, imagine we reran the three simulations above, with \(n = 84, 32, & 20\), but the population correlation was \(\rho=.6\). It’s intuitive to think that the distribution of sample correlation coefficients would spread with smaller sample sizes and our critical correlation coefficient would change, yet we should have more samples with \(p<\alpha=.05\) because they would be centering further away from \(0\). I will quickly re-run all analysis using a different \(\rho\) and plot the results and power. My script is almost identical to above, so it is excluded for the most part.

## 
##    1 
## 2000

pwr.r.test(n=84, sig.level = .05, r=.6)

## 
##      approximate correlation power calculation (arctangh transformation) 
## 
##               n = 84
##               r = 0.6
##       sig.level = 0.05
##           power = 0.9999918
##     alternative = two.sided

## 
##    0    1 
##   70 1930

pwr.r.test(n=32, sig.level = .05, r=.6)

## 
##      approximate correlation power calculation (arctangh transformation) 
## 
##               n = 32
##               r = 0.6
##       sig.level = 0.05
##           power = 0.965675
##     alternative = two.sided

## 
##    0    1 
##  310 1690

pwr.r.test(n=20, sig.level = .05, r=.6)

## 
##      approximate correlation power calculation (arctangh transformation) 
## 
##               n = 20
##               r = 0.6
##       sig.level = 0.05
##           power = 0.8306364
##     alternative = two.sided