Power Analysis for Data Scientists
Power analysis can be an intimidating topic for data scientists, but it’s an essential skill for those venturing into the world of research and statistics. In this article, we’ll walk through a power analysis using the National Health and Nutrition Examination Survey (NHANES) dataset. Along the way, we’ll provide a refresher on power analysis and demonstrate two methods of answering power-related questions. So, grab your Python toolkit, and let’s dive into the fascinating world of power analysis!
The NHANES dataset, a treasure trove of health and nutritional data for the United States, is maintained by the Centers for Disease Control and Prevention (CDC). With 197 attributes and 8,366 samples, it offers countless analysis opportunities. We’ll focus on comparing the weights of male and female survey respondents between the ages of 20 and 30.
To streamline our data import process, we can rely on the “nhanes” package maintained by Russell Poldrack of Stanford. Once we have the data, we can separate the groups by gender and age, ultimately resulting in 307 male respondents and 342 female respondents between the ages of 20 and 30. Armed with this information, we can create a quick visual comparing the distribution of weights between the two groups.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import scipy.stats
import statsmodels.stats.power as smp
import matplotlib.pyplot as plt
from nhanes.load import load_NHANES_data
# Download H1B salary data
df = load_NHANES_data()
# Select the attributes we want and remove missing values
df = df[['Gender',
'AgeInYearsAtScreening',
'AnnualHouseholdIncome',
'AnnualFamilyIncome',
'WeightKg']]
# Drop rows with missing values
dfoi = dfoi.dropna()
print(f'The number of rows in the data after dropping rows with missing values is {len(dfoi)}')
"""
The number of rows in the data after dropping
rows with missing values is 7130
"""
# Filter the data by age and separate the data frames into two groups
df20_30 = dfoi[(dfoi['AgeInYearsAtScreening'] >= 20) & (dfoi['AgeInYearsAtScreening'] < 30)]
# Separate the data based on Gender
M_20_30 = df20_30[(df20_30['Gender'] == 'Male')]
F_20_30 = df20_30[(df20_30['Gender'] == 'Female')]
# Print the number of samples in each group
print(f'Number of samples in M_20_30: {len(M_20_30)}')
print(f'Number of samples in F_20_30: {len(F_20_30)}')
"""
Number of samples in M_20_30: 307
Number of samples in F_20_30: 342
"""
# Visualize the two groups (M_20_30 and F_20_30) using an overlapping kde plot
sns.kdeplot(M_20_30['WeightKg'], fill=False, label='M_20_30')
sns.kdeplot(F_20_30['WeightKg'], fill=False, label='F_20_30')
# Add gridlines
plt.grid(alpha = 0.5)
# Add more x-axis ticks
plt.xticks(np.arange(0, 200, 10))
# Draw vertical lines at the mean of each group
plt.axvline(M_20_30['WeightKg'].mean(), color='b', linestyle='--', alpha=.5)
plt.axvline(F_20_30['WeightKg'].mean(), color='r', linestyle='--', alpha=.5)
# Put ticks at 45 degree angle
plt.xticks(rotation=45)
plt.title('Distribution of WeightKg for M_20_30 and F_20_30')
plt.show()
The red and blue vertical lines represent the mean weight of women and men, respectively. So, based on the plot, we can see our data is not perfectly normal, but it will still work well for power analysis. Now that we have our data and can observe the initial distribution, let’s discuss a little bit more about what Power Analysis is.
A Quick Refresher on Power Analysis
Power analysis is a statistical technique that helps us determine the sample size needed to detect an effect of a given size in a study. By performing a power analysis, we, as data scientists and researchers, can ensure our statistical studies avoid Type II errors. To carry out a power analysis, you’ll need a few crucial pieces of information:
- Effect size: The magnitude of the difference between groups or the strength of the relationship between variables, which is estimated from historical data or can be calculated directly for two groups.
- Sample size: The number of participants or observations in the study.
- Significance level (alpha): The threshold used to determine whether an effect is statistically significant, typically set at 0.05.
- Power (1-beta): The probability of correctly rejecting the null hypothesis when it is false. A common target power is 0.8, which means there is an 80% chance of detecting a true effect.
By trading off these four components against one another, we can calculate the missing element. In our example, we’ll assume a significance level of 0.05, a power of 0.8, and we will calculate the effect size. From there, we can determine the number of samples we need to reach our desired power and significance level between the two groups. This allows us to answer the intriguing question:
“How many samples of men and women do we need to confidently predict whether the participants are male or female based on their body weight?”
Calculating and Simulating Power Analysis
We can calculate power in two ways — by calculation and by simulation. A great resource that we draw upon in the following analysis is Python Companion to Statistical Thinking in the 21st Century, by Russell Poldrack of Stanford (also manages the NHANES python package). Before start calculated power, we must first determine our effect size for the example.
Get Fuher Analytics’s stories in your inbox
Join Medium for free to get updates from this writer.
To determine the effect size for our two groups, we’ll use Cohen’s d statistical measure. Cohen’s d represents the standardized difference between two means and is often used as a measure of effect size in experimental and observational studies. Specifically, it represents the difference between the means of two groups in terms of the standard deviation of the data.
# Create a function for calulating the Cohen's d effect size
def Cohen_d(group1, group2):
"""
group1: a pandas series
group2: a pandas series
"""
if isinstance(group1, pd.Series) and isinstance(group2, pd.Series):
# Convert to Pndas Series
group1 = pd.Series(group1)
group2 = pd.Series(group2)
# Calculate the difference in means between the two groups
diff = np.mean(group1) - np.mean(group2)
# Calculate the pooled standard deviation
n1, n2 = len(group1), len(group2)
var1, var2 = np.var(group1, ddof=1), np.var(group2, ddof=1)
pooled_std = np.sqrt((n1 * var1 + n2 * var2) / (n1 + n2))
# Calculate Cohen's d
d = diff / pooled_std
return d
effect_size = Cohen_d(M_20_30['WeightKg'], F_20_30['WeightKg'])
print(f'The effect size between the two groups is {effect_size}')
"""
The effect size between the two groups is 0.43466982759083467
"""
Upon analyzing the data, we found that the effect size between the two groups is roughly 0.435, placing it within the medium-small range. For reference, an effect size of ≤0.2 is considered small, 0.5 is medium, and ≥0.8 is large. A larger effect size means fewer samples are needed to attain a high power level in your study. With the effect size calculated, we can proceed to determine our power. To support our calculations, let’s construct a power curve by maintaining constant significance level and effect size while varying the sample count (independent variable) and calculating the associated power (dependent variable).
# Calculate the power of the test
import statsmodels.stats.power as power_analysis
power_analysis = smp.TTestIndPower()
power_curve_values = []
for i in range(1, 200):
power = power_analysis.solve_power(effect_size=effect_size,
nobs1=i,
alpha=.05)
# print(f'The power of the test is {power} when nobs1={i}')
power_curve_values.append(power)
# Plot the power curve
plt.plot(power_curve_values)
plt.xlabel('Number of Samples')
plt.ylabel('Power')
# put a crosshair at the point where power is 0.8
plt.axhline(0.8, color='r', linestyle='--', alpha=.5)
# find the index of the first value in the list that is greater than 0.8
index = next(i for i, v in enumerate(power_curve_values) if v >= 0.8)
plt.axvline(index, color='r', linestyle='--', alpha=.5)
plt.title('Power Curve for Two-Sample T-Test for Men and Women Weight Ages 50-60')
# Add text to plot "Power = 0.8 at nobs1 =" in a wheat color box in bottom right corner
plt.text(x=index+20, y=0.7, s=f'Power = 0.8 at nobs1 = {index}', bbox=dict(facecolor='wheat', alpha=0.5))
# Add gridlines
plt.grid()
plt.show()
The plot reveals that a power of ≥0.8 is attained once we have 84 samples. This demonstrates how to calculate a study’s sample size using power analysis. The process is relatively straightforward and even generates an informative visual for sharing with colleagues and supervisors. Now, let’s see if simulation yields comparable results.
To perform the power analysis simulation, we can randomly sample 84 (our predicted sample size) from each group (male and female respondents) and conduct a t-test to determine if the two groups are equal. The simulation aspect involves repeating the above 10,000 times consecutively. The proportion of tests we reject (p < 0.05) corresponds to the test’s power. So, let’s run the simulation and observe the outcome.
"""Simulate the results from the power analysis calculation
round up from estimated sample size"""
sampleSize = index
# create a function that will generate samples and test for
# a difference between groups using a two-sample t-test
def get_t_result(group1_data, group2_data, sampleSize, effectSize):
"""
perform a ttest on random data of n=sampSize
"""
# Sample data from group1_data and group2_data
group1 = np.random.choice(group1_data, size=sampleSize)
group2 = np.random.choice(group2_data, size=sampleSize)
ttresult = scipy.stats.ttest_ind(group1, group2)
return(ttresult.pvalue)
# Iterations:
iters = 10_000
# create input data frame for output
power_sim_results = pd.DataFrame({'p_value': np.zeros(iters)})
for run in range(iters):
power_sim_results.loc[run, 'p_value'] = get_t_result(M_20_30['WeightKg'], F_20_30['WeightKg'], sampleSize, effect_size)
p_reject = np.mean(power_sim_results['p_value'] < 0.05)
print(f'The power of the test is {p_reject} when nobs1={sampleSize}')
"""
The power of the test is 0.8026 when nobs1=84
"""
Our simulation results closely align with those obtained from the calculation and power curve. After 10,000 iterations, the test’s power with 84 samples in each group is 0.8026.
I hope this explanation proves helpful for those tasked with conducting power analyses in a professional setting. Keep in mind the importance of estimating a suitable effect size for your study and collaborating with subject matter experts to determine the appropriate significance level and power for your specific test. I’d like to acknowledge and thank Russell Poldrack for providing valuable resources to address these questions. I highly recommend exploring Statistical Thinking in the 21st Century by Russell Poldrack of Stanford, as it is freely accessible online.
Links:
https://statsthinking21.github.io/statsthinking21-core-site/index.html
https://statsthinking21.github.io/statsthinking21-python/09-StatisticalPower.html
