This analysis applies survival analysis methods to the nafld1 dataset from the survival package. The dataset contains baseline data from a population study of non-alcoholic fatty liver disease (NAFLD), with one observation per subject. The outcome of interest is time-to-death (futime, in days), with status indicating whether the event occurred (1) or was censored (0).
data(nafld, package = "survival")
# Create survival object
surv_object <- Surv(nafld1$futime, nafld1$status)
# Quick look at the data
head(nafld1) |> kable(caption = "First 6 rows of nafld1")
1 The median survival time is the point at which the survival probability crosses 0.50 – i.e., half of subjects have experienced the event by this time.
The log-rank test yields \(\chi^2\) = 11.295 with p = 8^{-4}. This is statistically significant (p < 0.05), indicating a meaningful difference in survival between males and females.
3. Univariate Cox Models for Continuous Covariates
Each continuous covariate – age, weight, height, and BMI – is tested individually using a Cox proportional hazards model to assess its relationship with survival.
covariates <- c("age", "weight", "height", "bmi")
uni_results <- lapply(covariates, function(var) {
formula <- as.formula(paste("surv_object ~", var))
fit <- coxph(formula, data = nafld1)
s <- summary(fit)
coef_row <- s$coefficients[1, ]
conf <- s$conf.int[1, ]
data.frame(
Covariate = var,
HR = round(conf["exp(coef)"], 4),
Lower_95 = round(conf["lower .95"], 4),
Upper_95 = round(conf["upper .95"], 4),
z = round(coef_row["z"], 3),
p_value = round(coef_row["Pr(>|z|)"], 4)
)
})
uni_table <- do.call(rbind, uni_results)
rownames(uni_table) <- NULL
kable(
uni_table,
col.names = c("Covariate", "Hazard Ratio", "95% CI Lower", "95% CI Upper", "z", "p-value"),
caption = "Univariate Cox Model Results"
)
Univariate Cox Model Results
Covariate
Hazard Ratio
95% CI Lower
95% CI Upper
z
p-value
age
1.1033
1.0985
1.1081
44.080
0.0000
weight
0.9946
0.9917
0.9974
-3.702
0.0002
height
0.9776
0.9719
0.9832
-7.718
0.0000
bmi
1.0005
0.9917
1.0094
0.115
0.9083
2 A hazard ratio (HR) > 1 indicates higher values of the covariate are associated with increased risk of death; HR < 1 indicates a protective association. Significance assessed at α = 0.05.
Covariates with p < 0.05 are considered statistically significant predictors of survival at the univariate level and are candidates for inclusion in the multivariate model.
4. Multivariate Cox Model
Rationale
Because BMI is mathematically derived from weight and height (\(\text{BMI} = \text{weight} / \text{height}^2\)), including all three in the same model introduces multicollinearity. To address this, two candidate models are compared:
Model A:male + age + bmi (uses the composite measure)
Model B:male + age + weight + height (uses the raw components)
cox_bmi <- coxph(surv_object ~ male + age + bmi, data = nafld1)
cox_wh <- coxph(surv_object ~ male + age + weight + height, data = nafld1)
cat("=== Model A: male + age + bmi ===\n")
## === Model A: male + age + bmi ===
summary(cox_bmi)
## Call:
## coxph(formula = surv_object ~ male + age + bmi, data = nafld1)
##
## n= 12588, number of events= 1018
## (4961 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## male 0.366140 1.442158 0.062849 5.826 5.69e-09 ***
## age 0.100601 1.105835 0.002650 37.967 < 2e-16 ***
## bmi 0.017022 1.017168 0.004945 3.443 0.000576 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## male 1.442 0.6934 1.275 1.631
## age 1.106 0.9043 1.100 1.112
## bmi 1.017 0.9831 1.007 1.027
##
## Concordance= 0.827 (se = 0.007 )
## Likelihood ratio test= 1750 on 3 df, p=<2e-16
## Wald test = 1472 on 3 df, p=<2e-16
## Score (logrank) test = 1699 on 3 df, p=<2e-16
cat("\n=== Model B: male + age + weight + height ===\n")
##
## === Model B: male + age + weight + height ===
summary(cox_wh)
## Call:
## coxph(formula = surv_object ~ male + age + weight + height, data = nafld1)
##
## n= 12588, number of events= 1018
## (4961 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## male 0.545900 1.726162 0.090392 6.039 1.55e-09 ***
## age 0.098361 1.103361 0.002769 35.521 < 2e-16 ***
## weight 0.005686 1.005702 0.001832 3.104 0.00191 **
## height -0.017951 0.982209 0.004597 -3.905 9.43e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## male 1.7262 0.5793 1.4459 2.0607
## age 1.1034 0.9063 1.0974 1.1094
## weight 1.0057 0.9943 1.0021 1.0093
## height 0.9822 1.0181 0.9734 0.9911
##
## Concordance= 0.828 (se = 0.007 )
## Likelihood ratio test= 1758 on 4 df, p=<2e-16
## Wald test = 1490 on 4 df, p=<2e-16
## Score (logrank) test = 1729 on 4 df, p=<2e-16
AIC Comparison
aic_table <- AIC(cox_bmi, cox_wh)
rownames(aic_table) <- c("Model A (bmi)", "Model B (weight + height)")
kable(aic_table, caption = "AIC Comparison of Candidate Models")
AIC Comparison of Candidate Models
df
AIC
Model A (bmi)
3
15982.66
Model B (weight + height)
4
15977.05
The model with the lower AIC is preferred as the best balance of fit and parsimony.
Proportional Hazards Assumption
# Test on the best model (update to whichever wins AIC)
best_model <- if (AIC(cox_bmi) <= AIC(cox_wh)) cox_bmi else cox_wh
ph_test <- cox.zph(best_model)
print(ph_test)
## chisq df p
## male 0.0274 1 0.8684
## age 2.3751 1 0.1233
## weight 4.3252 1 0.0376
## height 1.4784 1 0.2240
## GLOBAL 13.4690 4 0.0092
3 The proportional hazards (PH) assumption requires that the hazard ratio between groups remains constant over time. A significant p-value (< 0.05) in the Schoenfeld residuals test indicates a violation for that covariate.
