功效可视化模拟与内置power函数结果不符技术问询

Let's break down why your simulation is giving lower power (~40-55%) than the 80% predicted by power.t.test, and how to fix it.

Key Issues in Your Code #

1. Incorrect Standard Deviations Used for Data Generation #

Your power.t.test call uses sd = sd1 (which is 6, calculated from mean1=40 and cv=15%), but in the simulation, you're generating both groups using pooled_sd (~6.627) instead of their actual individual SDs:

• Group 1 should have sd = sd1 = 6
• Group 2 should have sd = sd2 = 7.2 (since mean2=48, cv=15% → (15*48)/100=7.2)

Using the pooled SD for both groups changes the variability structure that power.t.test assumed, which directly lowers the simulated power.

2. Wrong Method for Determining Statistical Significance #

You're using confidence interval non-overlap to tag groups as "different":

different = ifelse(ci_lower_s2 > ci_upper_s1 ,'yes', 'no')

This is not equivalent to a two-sample t-test with sig.level=0.05. Confidence interval non-overlap is a stricter condition than p<0.05—many cases where the t-test would reject the null hypothesis (and count towards power) will have overlapping confidence intervals. This means you're undercounting true positives, leading to artificially low power.

3. Inefficient Simulation Structure #

You're generating data for all rep_sequence values in every loop iteration, but only care about replicate_n=10. This doesn't affect power accuracy, but it slows down your code unnecessarily.

Fixed Simulation Code #

Here's the corrected version of your simulation that aligns with power.t.test's assumptions:

cv <- 15 # Coefficient of variance
percent_increase <- 20 # Percent improvement to detect
mean1 <- 40
mean2 <- mean1 + (mean1*(percent_increase/100))
sd1 <- (cv*mean1)/100 # SD for group 1: 6
sd2 <- (cv*mean2)/100 # SD for group 2: 7.2
difference <- (percent_increase/100)*mean1

# Get required sample size from power.t.test
pwrt <- power.t.test(delta=difference, sd=sd1, power=0.8, sig.level = .05, type="two.sample", alternative = "two.sided")
print(paste("Number of replicates needed is", round(pwrt$n)))

# Simulate power correctly
set.seed(123) # For reproducibility
record_test <- c()
replicate_n <- round(pwrt$n) # Use the exact n from power.t.test
n_simulations <- 1000

for(i in 1:n_simulations){
  # Generate data with correct group-specific SDs
  d <- rnorm(replicate_n, mean = mean1, sd = sd1)
  d2 <- rnorm(replicate_n, mean = mean2, sd = sd2)
  
  # Perform two-sample t-test (matches power.t.test's assumptions)
  t_result <- t.test(d, d2, var.equal = FALSE) # Use Welch's t-test (default in power.t.test when type="two.sample")
  # Tag as "yes" if p-value < 0.05
  test_result <- ifelse(t_result$p.value < 0.05, "yes", "no")
  record_test <- c(record_test, test_result)
}

# Calculate simulated power
power_sim <- sum(record_test == "yes") / length(record_test)
print(paste("Simulated power:", round(power_sim*100, 1), "%"))
print(table(record_test)/length(record_test))

What Changed? #

• Correct SDs: We now generate each group's data using their own SD (sd1 for group 1, sd2 for group 2) instead of the pooled SD.
• T-test for Significance: We use a two-sample t-test (Welch's, which is the default for power.t.test) and check if the p-value is below 0.05—this directly matches the statistical test that power.t.test uses to calculate power.
• Simplified Simulation: We only generate data for the required sample size, making the code faster and cleaner.

When you run this, you should see simulated power close to the 80% predicted by power.t.test (usually within ±2-3% due to random variation).

Additional Notes #

• If you want to visualize the confidence intervals alongside the t-test results, you can still generate the plots, but remember that CI overlap doesn't equal non-significance.
• Using set.seed() ensures your simulation results are reproducible, which is helpful for debugging.

内容的提问来源于stack exchange，提问作者Silvio Ortiz