R语言时变Cox模型中各区域年龄效应估计及置信区间计算问询

Hey there! Let's walk through how to pull the age effect estimates (and their 95% confidence intervals) for each Region from your time-varying Cox model, so you can build your effect modification table.

Step 1: Break Down Your Model Structure #

First, a quick check: your model includes an interaction between age and strata(Region), plus a time-varying BMI term (BMI.time.var) and clustering by id. The strata(Region) lets the baseline hazard vary by region, while the interaction gives you a distinct age effect for each region—exactly what you need for your effect modification analysis.

Step 2: Extract & Calculate Your Effect Estimates #

You can use base R or the broom package to streamline this process. Here are both approaches:

Option 1: Base R (No Extra Packages Needed)

Start by grabbing the full model summary:

summary_fit <- summary(fit1.1)

The coefficient table (summary_fit$coefficients) will have rows for each region-specific age term (look for labels like age:strata(Region=X)). Now extract these rows and compute HRs + 95% CIs (since Cox model coefficients are log(HRs)):

# Filter to keep only age-region interaction terms
age_coefs <- subset(
  summary_fit$coefficients,
  grepl("age", rownames(summary_fit$coefficients))
)

# Calculate HR, 95% CI, and keep p-values
age_effects <- data.frame(
  Region = gsub("age:strata\\(Region=(.*)\\)", "\\1", rownames(age_coefs)),
  HR = exp(age_coefs[, "coef"]),
  HR_Lower = exp(age_coefs[, "coef"] - 1.96 * age_coefs[, "se(coef)"]),
  HR_Upper = exp(age_coefs[, "coef"] + 1.96 * age_coefs[, "se(coef)"]),
  p_value = age_coefs[, "Pr(>|z|)"]
)

This gives you a clean data frame with all the values you need for each region.

Option 2: Using the broom Package (Cleaner & More Readable)

If you don’t have broom installed yet, set it up first:

install.packages("broom")
library(broom)

Use tidy() to convert the model output into a structured table, then filter and compute your metrics:

# Get a tidy coefficient table (keep log HRs for calculation)
tidy_fit <- tidy(fit1.1, exponentiate = FALSE)

# Filter for age-related interaction terms
age_tidy <- subset(tidy_fit, grepl("age", term))

# Add HR, 95% CI, and clean up region labels
age_effects_tidy <- age_tidy %>%
  mutate(
    Region = gsub("age:strata\\(Region=(.*)\\)", "\\1", term),
    HR = exp(estimate),
    HR_Lower = exp(estimate - 1.96 * std.error),
    HR_Upper = exp(estimate + 1.96 * std.error)
  ) %>%
  select(Region, HR, HR_Lower, HR_Upper, p.value)

This output is ready to directly use for your effect modification table.

Step 3: Verify Your Results #

Quick sanity check: cross-reference the values with summary(fit1.1) to make sure the region labels match your data and the HRs align with your expectations.

Example Table Output #

Your final table might look like this (adjust based on your actual regions):

That’s all you need to build your effect modification table!

内容的提问来源于stack exchange，提问作者Trenton Honda