1
$\begingroup$

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
set.seed(800)
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)
summary(fit)
# 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)
summary(fit2)

# 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.

$\endgroup$

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.