0
$\begingroup$

I am using simulation for a power analysis of an experiment tested with a multiple logistic regression model. I find substantial improvements in power for my model with covariate vs without covariate.

# packages
if (!require("pacman")) install.packages("pacman") 
pacman::p_load(broom, tidyverse)

Simulating Data with Binary outcome

set.seed(42)

N <- 1000
B_INT <- -.5  
B_COV_MC <- .1  
B_COND <- 0.43

df <- 
  tibble(
    x_int = 1,
    x_cond = rbinom(n = N, size = 1, pr = .5),
    x_cov = rnorm(n = N, mean = 40, sd = 15), # extract from your baseline data
    x_cov_mc = x_cov - mean(x_cov), # mean centering covariate
    y_lor = B_INT*x_int + B_COND*x_cond + B_COV_MC*x_cov_mc,
    y_pr = 1/(1+exp(-y_lor)),
    y = rbinom(N, 1, y_pr) 
  ) 

cor(df$x_cov_mc, df$y) # strong correlation of covariate with outcome

Fitting 1000 models with covariate

set.seed(42)
N_SAMPLING_DIST <- 1e3
p_vector <- vector(length = N_SAMPLING_DIST, mode = 'numeric')

for (i in 1:N_SAMPLING_DIST) {
df_cond_cov <- 
  tibble(
    x_int = 1,
    x_cond = rbinom(n = N, size = 1, pr = .5),
    x_cov = rnorm(n = N, mean = 40, sd = 15), # extract from your baseline data
    x_cov_mc = x_cov - mean(x_cov),
    y_lor = B_INT*x_int + B_COND*x_cond + B_COV_MC*x_cov_mc,
    y_pr = 1/(1+exp(-y_lor)),
    y = rbinom(N, 1, y_pr) 
  ) 

m <- glm(y ~ x_cond + x_cov, family = 'binomial', data = df_cond_cov)

p <- tidy(m) %>% 
  filter(term == 'x_cond') %>% 
  pull(p.value)

p_vector[i] <- p
}

cond_cov_power <- mean(p_vector < .05)
cond_cov_power

81% of simulated models with N = 1000 observed a significant effect for x_cond when controlling for covariates.

Fitting 1000 models without covariate

set.seed(42)
p_vector <- vector(length = N_SAMPLING_DIST, mode = 'numeric')

for (i in 1:N_SAMPLING_DIST) {
df_cond_cov <- 
  tibble(
    x_int = 1,
    x_cond = rbinom(n = N, size = 1, pr = .5),
    x_cov = rnorm(n = N, mean = 40, sd = 15), # extract from your baseline data
    x_cov_mc = x_cov - mean(x_cov),
    y_lor = B_INT*x_int + B_COND*x_cond + B_COV_MC*x_cov_mc,
    y_pr = 1/(1+exp(-y_lor)),
    y = rbinom(N, 1, y_pr) 
  ) 

m <- glm(y ~ x_cond, family = 'binomial', data = df_cond_cov)

p <- tidy(m) %>% 
  filter(term == 'x_cond') %>% 
  pull(p.value)

p_vector[i] <- p
}

cond_power <- mean(p_vector < .05)
cond_power

67% of simulated models with N = 1000 observed a significant effect for x_cond when NOT controlling for covariates.

This improvement in power is in the expected direction, but I am surprised that I have higher reported power from power.prop.test() in base R (see below).

power.prop.test() results

group_means <- 
    df_cond_cov %>% 
    group_by(x_cond) %>% 
    summarise(mean = mean(y))

group_means

Group means are 41% & 53%.

power.prop.test(p1 = group_means$mean[1], p2 = group_means$mean[2], n = N/2)

How can I have 97% power using just my treatment condition in power.prop.test(), when my logistic regression with covariates only achieved 81% power?

$\endgroup$

1 Answer 1

1
$\begingroup$

The group means you are using (41%, 53%) depend strongly on the seed.

These appear to be not actually the true marginal average probabilities of your model - note that these were averages of only 1000 samples.

I expect if you instead computed group means over a df_cond_cov that had e.g. 1 million rows, this would fix the issue and probably give you a more reasonable-seeming answer

-

$\endgroup$
1
  • $\begingroup$ Thank you @mike, Indeed creating a base sample of 1 million rows estimated more accurate marginal means for the treatment and control groups (e.g., .413 & .488) regardless of seed, and the power.prop.test() output for these marginal means is now in line with my bootstrapped power calculations. $\endgroup$
    – Joe
    Sep 24, 2021 at 13:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.