I'm trying to understand/replicate an adjusted logistic regression analysis where a treatment effect is estimated separately in a number (>2) of subgroups and estimating an overall P-value for interaction.

I want to estimate the odds ratios (ORs) with 95% confidence intervals (CIs) in each subgroup (i.e., without keeping any subgroup as the reference group with OR 1.00) while adjusting for two additional variables and present a P-value for the overall interaction effect. The results I want to produce are essentially similar to those of presented in Fig S3 here for the first set of subgroups (age) (page 10 in https://www.nejm.org/doi/suppl/10.1056/NEJMoa1808561/suppl_file/nejmoa1808561_appendix.pdf ; full paper at https://www.nejm.org/doi/full/10.1056/NEJMoa1808561).

I have read several posts here and elsewhere, but have not found one that I am able to understand enough to correctly perform the analysis.

Here's code to generate an example dataset in R and perform parts of the analyses:

# Create example dataset with 2000 patients
dta <- data.frame(trt = c(rep(FALSE, 1000), rep(TRUE, 1000)), # Allocated to treatment group (TRUE) or control group (FALSE)
                  subgrp = factor(rep(LETTERS[1:5], 400), levels = LETTERS[1:5]), # Five subgroups
                  y = NA, # Add outcome variable, which is set in the next lines
                  adj_1 = rbinom(2000, 1, 0.5), # Add two variables to adjust for
                  adj_2 = rnorm(2000, 150, 2))
dta$y[dta$trt == FALSE] <- rbinom(n = 1000, 1, 0.30) # 30% chance of outcome in the control group

for (i in 1:5){
  # Different chances of outcome in each subgroup
  dta$y[dta$trt == TRUE & dta$subgrp == LETTERS[i]] <- rbinom(n = 200, 1, 0.22 + 0.04 * i)

# Test the overall treatment effect
fit <- glm(y ~ trt + adj_1 + adj_2, family = binomial(link = "logit"), data = dta)
# From these results I can calculate the 95% CI for the coefficient for trt and
# exponentiate coefficient and CI to OR and corresponding CI

# Fit model with subgroup interaction (and corresponding main effect of subgrp)
fit2 <- glm(y ~ trt * subgrp + adj_1 + adj_2, family = binomial(link = "logit"), data = dta)

# I am uncertain about how to get the desired results from each subgroup here

# Calculate overall P-value for interaction - which approach is correct (if any)?
drop1(fit2, test = "Chisq") # Or test =  "LRT", same results
anova(fit, fit2, test = "Chisq")

I am uncertain about how to calculate the ORs with CIs in each subgroup, and uncertain which approach is correct for calculating the overall P-value for interaction.

Help in explaining how to modify the above code to accomplish this would be greatly appreciated.

Edit: minor re-phrasing for clarity.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.