I have a large dataset of countries, each categorized into one of two regions, and I'm interested in how the effect of a predictor variable 𝑋 differs between the two regions. I plan to fit separate regression models for each region. While there are alternatives, such as incorporating interaction terms between region and 𝑋 in a single model, I have chosen to start with two distinct models.
Here's the issue: The regions have different numbers of countries — Region 1 has 70 countries, while Region 2 has 20. In my analysis, the main effect of 𝑋 is significant in the larger region (70 countries) but not significant in the smaller region (20 countries). I am concerned that this difference might simply be due to the disparity in sample size rather than reflecting a real difference in the effect of 𝑋.
Approach: To assess whether the observed difference in significance is truly meaningful and not a result of the sample size disparity, I’m considering downsampling the larger region (N = 70) to match the smaller region (N = 20) and refitting the model. Since a single downsampling run may be subject to random variation and introduce biases, I would repeat the downsampling process multiple times (e.g., 𝑘 = 100) to get a more reliable picture.
Questions:
- Is this a methodologically sound approach?
- Does downsampling in this way provide a robust way to compare the results across the two regions?
- Are there any alternative approaches that I should consider to ensure the comparability of the results?
- How should one report the results of repeated downsampling?
- Given that I plan to fit 100 regression models for the downsampled data, what are the best practices for summarizing the results across these models? (Tables, Figures)
- Should I report the average coefficients, confidence intervals, and/or p-values across the 100 models? What other metrics could provide meaningful insights?
Minimal Example in R (using mtcars):
For illustration of my attempt, I’ve used the mtcars dataset in R. This dataset includes a binary variable am (denoting 0 = automatic vs. 1 = manual transmission), which I’m treating as equivalent to 'region' in my actual dataframe. My goal is to downsample the cars with am = 1 (representing the larger region) to match the number of cars with am = 0 (representing the smaller region) and fit separate regression models in each case. Then, I repeat the downsampling procedure, in this case 𝑘 = 10 times.
Here's a simplified version of my approach in R, which should theoretically resemble the ideas behind Monte Carlo simulations.
You can copy the code here.
# Load necessary libraries
library(tidyverse)
library(broom)
library(binom)
# Information on dataframe
?mtcars
# Increase the imbalance in mtcars to N = 19 automatic and N = 8 manual cars
df <- mtcars %>%
mutate(car_model = rownames(mtcars)) %>%
filter(!car_model %in% c('Lotus Europa', 'Ford Pantera L', 'Ferrari Dino', 'Maserati Bora', 'Volvo 142E'))
# Fit two models separately for automatic and manual cars
model_big <- lm(mpg ~ disp, data = df %>% filter(am == 0))
model_small <- lm(mpg ~ disp, data = df %>% filter(am == 1))
# Display summary of the original models
summary(model_big)
summary(model_small)
# Get the sample sizes
n_small <- nrow(df %>% filter(am == 1))
n_big <- nrow(df %>% filter(am == 0))
# Set the number of iterations for downsampling
iterations <- 10
# Create an empty list to store the results
results <- list()
# Downsample the big group (am == 0) to the size of the small group (am == 1) and fit OLS models
set.seed(123)
for (i in 1:iterations) {
# Downsample automatic cars (am == 0) to the size of the manual cars group
downsampled_data <- df %>%
filter(am == 0) %>%
sample_n(n_small)
# Fit the model to the downsampled data
refitted_model <- lm(mpg ~ disp, data = downsampled_data)
# Extract model statistics for the disp coefficient
model_summary <- tidy(refitted_model)
disp_info <- model_summary %>% filter(term == 'disp')
# Store results: estimate and p-value for disp
results[[i]] <- data.frame(
iteration = i,
estimate = disp_info$estimate,
p_value = disp_info$p.value
)
}
# Combine results into a single data frame
results_df <- bind_rows(results)
# Summarize the results across iterations
summary_stats <- results_df %>%
summarise(
mean_estimate = mean(estimate),
sd_estimate = sd(estimate),
mean_p_value = mean(p_value),
prop_significant = mean(p_value < 0.05)
)
# Calculate binomial confidence intervals for the proportion of significant p-values
significant_count <- sum(results_df$p_value < 0.05)
total_count <- nrow(results_df)
# Exact binomial confidence interval
binom_ci_exact <- binom.test(significant_count, total_count, conf.level = 0.95)$conf.int
# Wilson confidence interval
binom_ci_wilson <- binom.confint(significant_count, total_count, method = 'wilson')[c('lower', 'upper')]
# Display summary statistics
summary_table <- tibble(
Statistic = c(
'Mean Estimate',
'SD of Estimate',
'Mean p-value',
'Proportion of Significant p-values',
'Confidence Interval (Exact, 95%)',
'Confidence Interval (Wilson, 95%)'
),
Value = c(
round(summary_stats$mean_estimate, 2),
round(summary_stats$sd_estimate, 2),
round(summary_stats$mean_p_value, 2),
round(summary_stats$prop_significant, 2),
paste("[", round(binom_ci_exact[1], 3), ",", round(binom_ci_exact[2], 3), "]"),
paste("[", round(binom_ci_wilson$lower, 3), ",", round(binom_ci_wilson$upper, 3), "]")
)
)
# Display the summary statistics table
print(summary_table)
# Plot a histogram of p-values
ggplot(data = results_df, aes(x = p_value)) +
geom_histogram(aes(y = after_stat(count) / sum(after_stat(count)) * 100),
binwidth = 0.05, fill = 'lightgrey', color = 'darkgray') +
geom_vline(xintercept = 0.05, linetype = 'dashed', color = 'deeppink') +
labs(title = 'Histogram of p-values with α = 0.05',
x = 'P-value',
y = 'Relative frequency') +
scale_y_continuous(labels = scales::percent_format(scale = 1)) +
theme_minimal()
