A simulation-based sample size and Power calculation with potential treatment-sex interaction

Statistical model and test

To estimate treatment-sex interaction effect, we will use a simple linear regression model, adjusted for age, baseline DAS, MSC intervention, sex/gender, and sex/gender interaction with MSC intervention on the simulated data. We will simulate 1000 copies of the study data and fit each simulated data to the above regression model and use F-test at 5% significant level to determine statistical significance of the interaction effect. The simulated power will be calculated as the proportion of significant P-values (P-value < 0.05) out of the 1000 simulated iterations.

We code our power calculation function below including steps to simulate synthetic data and the interaction effect.

pwr_normdat <- function(DAS0_m = 25,
                        DAS_sd = 5, 
                        beta_1 = 1, 
                        beta_trt = 5,
                        beta_age = 0,
                        beta_sex = 0,
                        beta_interaction=5, 
                        N=17, #estimated sample size per treatment group based on primary analysis without interaction
                        reps=1000,
                        alpha=.05){
  
  # initiating empty matrix to save regression coefficient and P-value;
  sim_reg_coef <- matrix(nrow = reps, ncol = 5)
  sim_reg_pvalue <- matrix(nrow = reps, ncol = 5)
  
  for (i in 1:reps){
    
    # step 1. data generating function, 
    n_tot = 2*N #total sample size;
    
    trt <- rep(0:1, N) #randomly assign treated and control at 1:1 ratio;
    sex <- rbinom(n_tot, 1, .5) #randomly assign male and female at 1:1 ratio;
    DAS0 <- rnorm(n_tot, mean = DAS0_m, sd = DAS_sd)
    age <- rtruncnorm(n_tot, a = 12, b=17, mean = 15, sd = 2) #restricted age distribution to be bounded by 12 and 17; 
    
    # overall treatment effect and a differential subgroup effect;
    DAS_mean = beta_1*DAS0 +  beta_trt*trt + beta_age*age + beta_sex*sex + beta_interaction*trt*sex
    DAS = rnorm(n_tot, mean = DAS_mean, sd = DAS_sd)
    
    dat <- data.frame(trt, sex, DAS0, DAS)
    
    #fitting linear regression, following the primary analysis model;
    fit<-lm(DAS ~ DAS0 + age + trt*sex, data = dat)
    
    #F-test
    interaction <-aov(fit)
    
    sim_reg_coef[i,] <- summary(fit)$coef[2:6,1] #all coefficients except intercept;
    sim_reg_pvalue[i,] <- summary(interaction)[[1]][["Pr(>F)"]][-6] #remove the last NA p-value for residual;
    
  }
  
  # Summarize results by coefficient
  sim_results <- data.frame(#regression covariates;
                            covariates = rownames(summary(fit)$coef)[-1],
                            #estimated (effect) coeffient on the mean;
                            mean_coef = colMeans(x=sim_reg_coef), 
                            #number of P-values for each regression coef < 0.05;
                            simpower = apply(sim_reg_pvalue, 2, function(x) (length(x[x< alpha])/reps)*100))
  
  
  # Return table of results
  return(sim_results)
}

Power calculation on interaction effect

For interaction effect power calculation, we will consider a range of interaction effect with effect size proportional to the overall treatment effect at 5 points on the primary outcome. We considered a set of interaction effect from 0.1 up-to 3 times the size of the overall treatment effect with 0.1 increment.

We will conduct our power calculations for a sample size at 17 for each intervention group (following the two-sample t-test sample size calculation).

set.seed(123)
interaction_effect <- c(5*seq(0.1,3,by=0.1))

#initializing vector of simulated powers for interaction;
prop_sig_interaction <- rep(NA, length(interaction_effect))

for (i in 1:length(interaction_effect)){

  beta_interaction <- interaction_effect[i]
  
  pwr_res <- pwr_normdat(DAS0_m = 25,
                         DAS_sd = 5, 
                         beta_1 = 1, 
                         beta_trt = 5, 
                         beta_age = 0,
                         beta_sex = 0,
                         beta_interaction= beta_interaction, 
                         N=17,
                         reps=1000,
                         alpha=.05)

  prop_sig_interaction[i]<- pwr_res$simpower[5]
}

power_data1<-data.frame(size_interaction_prop = seq(0.1,3,by=0.1),
                       size_interaction = 5*seq(0.1,3,by=0.1),
                       sig_interaction = prop_sig_interaction)

knitr::kable(power_data1, 
             digits = 1,
             align='ccc',
             col.names = c("Interaction effect <br/> proportional to treatment effect","Interaction effect","Simulated power"),
             row.names = FALSE,
             escape = FALSE)

ggplot(power_data1, aes(size_interaction_prop, sig_interaction)) +
  geom_line()+
  geom_point()+
  scale_x_continuous(breaks = seq(0,3,by=0.1))+
  scale_y_continuous(breaks = seq(0,100,by=5))+
  geom_hline(yintercept = 80, color="red", linetype="dashed")+
  geom_vline(xintercept = 2.12, color="red", linetype="dashed")+
  labs(x="size of interaction effect relative to overall treatment effect",
       y="% of significant interaction at 5% significance level",
       title = "1000 simulated linear regression model with gender and intervention interaction \n Total sample size = 17 x 2 = 34 \n Overall intervention effect at 5 points on DAS")+
  theme_bw()

Figure 1: Power Curve for interaction effect with 17 study subjects per intervention group

Conclusion from the simulation-based power calculation

For a total sample size of 34 (17 in each treatment group), we can successfully detect a treatment-sex interaction effect that is about 2 times the overall treatment effect at 5*2 = 10 points changes on the 3-month DAS with 80% power (Figure above). Therefore, we do not have enough power at 34 patients to detect any true interaction effect that is less than 10 units change on the outcome.