MReyes_EHA_PS3
Matthew Reyes
2026-09-15
Multivariate parametric proportional hazards models
Using the MHAS data (mhas0124.csv) again…
Estimate a continuous-time parametric proportional hazards model assuming the shape of the baseline hazard is exponential, and using FEMALE as the only covariate. Obtain coefficients, not hazard ratios.
expmodel <- phreg(Surv(agecensor, died0124) ~ female , dist="weibull", shape = 1, data = mhas)
summary(expmodel)## Covariate Mean Coef Rel.Risk S.E. LR p
## female 0.549 -0.110 0.896 0.028 0.0001
##
## Events 5269
## Total time at risk 735431
## Max. log. likelihood -31283
## LR test statistic 15.80
## Degrees of freedom 1
## Overall p-value 7.05075e-05exp(expmodel$coefficients) ## female log(scale)
## 0.8960819 131.6190851Interpret the coefficient for FEMALE after manually exponentiating it.
- The exponentiating the female coefficient, we see that the risk of death for females is 89.6% of the hazard/risk for males. So, females face less hazard/risk for death than males do throughout the dataset.
Exponentiate the intercept of this model and interpret this figure.
exp(-1*expmodel$coefficients["log(scale)"])## log(scale)
## 0.007597682- From exponentiating the intercept of the model, we see the base hazard for males, which is about 0.0076 per year of age.
Sum the constant and the coefficient for female, and exponentiate this number. What is this equivalent to/how would you interpret it?
exp( -1*expmodel$coefficients["log(scale)"] + expmodel$coefficients["female"] )## log(scale)
## 0.006808146- This number is the actual hazards for females (as opposed to their hazard ratio). This is 0.0068, which is a lower hazard than that of males per year of age.
Estimate another proportional hazards parametric model controlling for FEMALE, SCHOOLING (or categories of EDUCLEVEL, whatever you prefer, justifying your choice), and LOCSIZE01 (categories or dichotomizing it into whatever threshold you think is meaningful). Assuming the baseline hazard is distributed…
I am splitting my location data into a binary where the larger cities are 1s and the smaller cities are 2s. I think this will help us generalize a thoeretical difference between urbanites and rural denizens.
Exponential
mhas$locbinary <- ifelse(mhas$locsize01 %in% c(1, 2), 1, 0)
# note that 1 and 2 are larger localities, so making them a 1 in the binary puts rural localities as the reference
mhas2 <- mhas[complete.cases(mhas[, c(
"female",
"schooling",
"locbinary",
"agecensor",
"died0124"
)]), ]
mhas2$school4 <- cut(
mhas2$schooling,
breaks = c(-Inf, 0, 5, 8, Inf),
labels = c("0 years", "1-5", "6-8", "9+")
)
mhas2$school4 <- factor(mhas2$school4,
levels = c("0 years", "1-5", "6-8", "9+"))
expmodel2 <- phreg(Surv(agecensor, died0124) ~ female + school4 + locbinary, dist="weibull", shape = 1, data = mhas2)
summary(expmodel2)## Covariate Mean Coef Rel.Risk S.E. LR p
## female 0.549 -0.139 0.871 0.028 0.0000
## school4 0.0000
## 0 years 0.266 0 1 (reference)
## 1-5 0.350 -0.131 0.878 0.034
## 6-8 0.211 -0.263 0.769 0.041
## 9+ 0.173 -0.525 0.592 0.048
## locbinary 0.741 0.009 1.009 0.032 0.7707
##
## Events 5258
## Total time at risk 734343
## Max. log. likelihood -31147
## LR test statistic 162.32
## Degrees of freedom 5
## Overall p-value 0# logLik:
LL <- expmodel2$loglik[2] # second element = log-likelihood at convergence
# number of parameters:
k <- expmodel2$df # number of parameters estimated (df)
# sample size (number of people):
n <- expmodel2$n
AIC_value <- 2*k - 2*LL
BIC_value <- log(n)*k - 2*LL
AIC_value## [1] 62304.58BIC_value## [1] 62340.31- Gompertz
gompmodel <- phreg(Surv(agecensor, died0124) ~ female + school4 + locbinary, dist = "gompertz", data = mhas2)
summary(gompmodel)## Covariate Mean Coef Rel.Risk S.E. LR p
## female 0.549 -0.275 0.759 0.028 0.0000
## school4 0.0000
## 0 years 0.266 0 1 (reference)
## 1-5 0.350 0.145 1.156 0.034
## 6-8 0.211 0.254 1.289 0.042
## 9+ 0.173 0.135 1.145 0.049
## locbinary 0.741 0.028 1.029 0.032 0.3816
##
## Events 5258
## Total time at risk 734343
## Max. log. likelihood -22128
## LR test statistic 145.67
## Degrees of freedom 5
## Overall p-value 0# logLik:
LL <- gompmodel$loglik[2] # second element = log-likelihood at convergence
# number of parameters:
k <- gompmodel$df # number of parameters estimated (df)
# sample size (number of people):
n <- gompmodel$n
AIC_value <- 2*k - 2*LL
BIC_value <- log(n)*k - 2*LL
AIC_value## [1] 44266.8BIC_value## [1] 44302.53- Weibull
weibmodel <- phreg(Surv(agecensor, died0124) ~ female + school4 + locbinary, dist = "weibull", data = mhas2)
summary(weibmodel)## Covariate Mean Coef Rel.Risk S.E. LR p
## female 0.549 -0.280 0.756 0.028 0.0000
## school4 0.0000
## 0 years 0.266 0 1 (reference)
## 1-5 0.350 0.115 1.122 0.034
## 6-8 0.211 0.215 1.240 0.041
## 9+ 0.173 0.099 1.105 0.048
## locbinary 0.741 0.037 1.037 0.032 0.2594
##
## Events 5258
## Total time at risk 734343
## Max. log. likelihood -21965
## LR test statistic 136.13
## Degrees of freedom 5
## Overall p-value 0# logLik:
LL <- weibmodel$loglik[2] # second element = log-likelihood at convergence
# number of parameters:
k <- weibmodel$df # number of parameters estimated (df)
# sample size (number of people):
n <- weibmodel$n
AIC_value <- 2*k - 2*LL
BIC_value <- log(n)*k - 2*LL
AIC_value## [1] 43939.55BIC_value## [1] 43975.29- Log-logistic
loglogmodel <- phreg(Surv(agecensor, died0124) ~ female + school4 + locbinary, dist = "loglogistic", data = mhas2)
summary(loglogmodel)## Covariate Mean Coef Rel.Risk S.E. LR p
## female 0.549 1.682 5.378 0.183 0.0000
## school4 0.0001
## 0 years 0.266 0 1 (reference)
## 1-5 0.350 -0.283 0.753 0.028
## 6-8 0.211 0.095 1.100 0.034
## 9+ 0.173 0.190 1.209 0.042
## locbinary 0.741 0.077 1.080 0.049 0.2088
##
## Events 5258
## Total time at risk 734343
## Max. log. likelihood -21941
## LR test statistic 130.77
## Degrees of freedom 5
## Overall p-value 0exp(-1 * loglogmodel$coefficients)## (Intercept) female school41-5 school46-8 school49+ locbinary
## 0.185926098 1.327328363 0.909438539 0.827018068 0.925793991 0.960131579
## log(scale) log(shape)
## 0.009871142 0.095642596# Log-likelihood
LL_loglog <- as.numeric(logLik(loglogmodel))
# Number of parameters
k_loglog <- loglogmodel$df
# Sample size
n_loglog <- loglogmodel$n
# AIC and BIC
AIC_loglog <- 2*k_loglog - 2*LL_loglog
BIC_loglog <- log(n_loglog)*k_loglog - 2*LL_loglog
AIC_loglog## [1] 43892.35BIC_loglog## [1] 43928.09- Log-normal
lognormmodel <- phreg(Surv(agecensor, died0124) ~ female + school4 + locbinary, dist = "lognormal", data = mhas2)
summary(lognormmodel)## Covariate Mean Coef Rel.Risk S.E. LR p
## female 0.549 5.563 260.579 0.681 0.0000
## school4 0.0002
## 0 years 0.266 0 1 (reference)
## 1-5 0.350 -0.280 0.756 0.028
## 6-8 0.211 0.092 1.096 0.034
## 9+ 0.173 0.183 1.200 0.041
## locbinary 0.741 0.068 1.070 0.049 0.2086
##
## Events 5258
## Total time at risk 734343
## Max. log. likelihood -21920
## LR test statistic 126.58
## Degrees of freedom 5
## Overall p-value 0exp(-1 * lognormmodel$coefficients)## (Intercept) female school41-5 school46-8 school49+ locbinary
## 0.003837613 1.323197833 0.912538948 0.833064199 0.934573535 0.960124578
## log(scale) log(shape)
## 0.004819632 0.325920902# Log-likelihood
LL_lognorm <- as.numeric(logLik(lognormmodel))
# Number of parameters
k_lognorm <- lognormmodel$df
# Sample size
n_lognorm <- lognormmodel$n
# AIC and BIC
AIC_lognorm <- 2*k_lognorm - 2*LL_lognorm
BIC_lognorm <- log(n_lognorm)*k_lognorm - 2*LL_lognorm
AIC_lognorm## [1] 43850.77BIC_lognorm## [1] 43886.51Use the table template provided in the accompanying Excel file to show the hazard ratios (HRs) for all three covariates under these different models and discuss/compare them (for log-logistic and lognormal models, obtain e-β instead of eβ, i.e., flip the coefficients or obtain the inverse of the HRs before you add them to the table. Examine other kinds of output given in all models, but especially in the last two).
See other file.
For females, the interpretation is consistent across all proportional hazard (PH) models that women have lower hazard for death at any point in time than males do. While it looks like the accelerated failure time (AFT) models show a higher hazard, my understanding is that this is because of how we obtained e-β, so it’s basically flipped. So, they too illustrate the lower hazard that women face.
While I’m still a little confused as to why some schooling HRs are below zero and some are above zero, they all still show the same directionality that increased education lowers the hazard.
For locality, it looks like the rural/urban divide has very little effect on mortality under any of these models. Still, it was interesting to me that they all consistently show a slightly higher hazard for urban respondents.
With regard to model fit, log‑normal is the best (lowest AIC and BIC). Log‑logistic is second-best, followed by Weibull, then Gompertz, then exponential.
Cox regression
Estimate a non-parametric model via Cox regression. Add these HRs to the table.
coxmod1 <- coxph(Surv(agecensor, died0124) ~ female + school4 + locbinary, method = "efron", data = mhas2)
summary(coxmod1)## Call:
## coxph(formula = Surv(agecensor, died0124) ~ female + school4 +
## locbinary, data = mhas2, method = "efron")
##
## n= 9392, number of events= 5258
##
## coef exp(coef) se(coef) z Pr(>|z|)
## female -0.28209 0.75421 0.02793 -10.101 < 2e-16 ***
## school41-5 0.09364 1.09817 0.03395 2.758 0.00581 **
## school46-8 0.18433 1.20241 0.04151 4.440 8.98e-06 ***
## school49+ 0.06439 1.06651 0.04853 1.327 0.18462
## locbinary 0.04233 1.04324 0.03246 1.304 0.19216
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## female 0.7542 1.3259 0.7140 0.7966
## school41-5 1.0982 0.9106 1.0275 1.1737
## school46-8 1.2024 0.8317 1.1085 1.3043
## school49+ 1.0665 0.9376 0.9697 1.1729
## locbinary 1.0432 0.9586 0.9789 1.1118
##
## Concordance= 0.554 (se = 0.005 )
## Likelihood ratio test= 128.5 on 5 df, p=<2e-16
## Wald test = 129.3 on 5 df, p=<2e-16
## Score (logrank) test = 130 on 5 df, p=<2e-16extractAIC(coxmod1)## [1] 5.00 84435.54LL <- as.numeric(logLik(coxmod1)) # log-likelihood
k <- length(coef(coxmod1)) # number of covariates
n <- coxmod1$n # number of observations
BIC_cox <- -2*LL + log(n)*k
BIC_cox## [1] 84471.28Prediction, diagnostics, time-varying covariates.
Choosing whatever model you prefer from I or II, assess whether the proportionality of hazards assumption is violated for FEMALE and SCHOOLING/EDUCLEVEL (separately). Does it hold for either/both?
- My understanding of the proportionality of hazards (PH) assumption is that it entails that two groups’ hazards can differ in values but not in shape over time. So, their hazards, when graphed, should be evenly spaced throughout the observation period. Looking at the Weibull model, and using the Cox model as a check for proportionality and using the Schoenfeld test, we see that both gender and schooling violate this assumption given that their p-values are smaller than .05. This is visible when graphing the betas for both variables, as they visibly slope as time (age) progresses. (Interestingly, it looks like the assumption holds for the location variable.)
cox_weibull_check <- coxph(Surv(agecensor, died0124) ~ female + school4 + locbinary,
data = mhas2)
ph_test <- cox.zph(cox_weibull_check)
ph_test## chisq df p
## female 11.90 1 0.00056
## school4 8.07 3 0.04457
## locbinary 1.60 1 0.20598
## GLOBAL 20.93 5 0.00084plot(ph_test)
Estimate a “Cox-stratified” model using a new dummy variable that divides people into whether they have less than 9 vs. 9 or more years of formal schooling as the stratification variable. Would your conclusions about the effect of education (linearly anyway) change doing this?
mhas$edubinary <- factor(ifelse(mhas$schooling >= 9, 1, 0),
levels = c(0, 1),
labels = c("LowEd", "HighEd"))
cox_strat_model <- coxph(
Surv(agecensor, died0124) ~ female + locbinary + strata(edubinary),
data = mhas
)
summary(cox_strat_model)## Call:
## coxph(formula = Surv(agecensor, died0124) ~ female + locbinary +
## strata(edubinary), data = mhas)
##
## n= 9392, number of events= 5258
## (244 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## female -0.28820 0.74961 0.02789 -10.335 <2e-16 ***
## locbinary 0.07731 1.08037 0.03150 2.454 0.0141 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## female 0.7496 1.3340 0.7097 0.7917
## locbinary 1.0804 0.9256 1.0157 1.1492
##
## Concordance= 0.543 (se = 0.005 )
## Likelihood ratio test= 107.6 on 2 df, p=<2e-16
## Wald test = 108.3 on 2 df, p=<2e-16
## Score (logrank) test = 109 on 2 df, p=<2e-16- Looking at the model’s output and the graph of the baseline hazards by strata (see below), it is very easy to see that the effect of education is definitely not linear and not constant over time. Interestingly, it really only breaks down between strata around 90 years of age, which makes me think that it’s not as significant as the previous models made me think it was.
“Recover” and graph the baseline hazard under the Cox model.
bh <- basehaz(cox_strat_model, centered = FALSE)
head(bh)bh_low <- subset(bh, strata == "LowEd")
bh_high <- subset(bh, strata == "HighEd")
plot(bh_low$time, bh_low$hazard, type="l", col="red",
xlab="Time", ylab="Baseline hazard", lwd=2)
lines(bh_high$time, bh_high$hazard, col="blue", lwd=2)
legend("topleft", legend=c("LowEd", "HighEd"),
col=c("red","blue"), lwd=2)
Estimate predicted survival curves or smoothed hazards (whichever you prefer or find more intuitive) for:
Men with 0 years of education
coxmod_clean <- coxph(
Surv(agecensor, died0124) ~ female + school4,
data = mhas2)
men0 <- data.frame(
female = 0,
school4 = factor("0 years", levels = levels(mhas2$school4)))
sf_men0 <- survfit(coxmod_clean, newdata = men0)
# Select only ages 55+
keep1 <- sf_men0$time >= 55
plot(sf_men0$time[keep1],
sf_men0$surv[keep1],
type = "l",
xlab = "Age",
ylab = "Predicted survival probability",
main = "Predicted Survival Curve: Men with 0 years of schooling (55+)",
col = "blue",
lwd = 2)
Women with 0 years of education
female0 <- data.frame(
female = 1,
school4 = factor("0 years", levels = levels(mhas2$school4)))
sf_female0 <- survfit(coxmod_clean, newdata = female0)
# Select only ages 55+
keep2 <- sf_female0$time >= 55
plot(sf_female0$time[keep2],
sf_female0$surv[keep2],
type = "l",
xlab = "Age",
ylab = "Predicted survival probability",
main = "Predicted Survival Curve: Females with 0 years of schooling (55+)",
col = "blue",
lwd = 2)
Men with 6 or 6-8 years of education (depending on the schooling variable you used)
men68 <- data.frame(
female = 0,
school4 = factor("6-8", levels = levels(mhas2$school4)))
sf_men68 <- survfit(coxmod_clean, newdata = men68)
keep68m <- sf_men68$time >= 55
plot(sf_men68$time[keep68m],
sf_men68$surv[keep68m],
type = "l",
xlab = "Age",
ylab = "Predicted survival probability",
main = "Predicted Survival Curve: Men with 6-8 Years of Education (55+)",
col = "red",
lwd = 2)
Women with 6 or 6-8 years of education (depending on the schooling variable you used)
female68 <- data.frame(
female = 1,
school4 = factor("6-8", levels = levels(mhas2$school4)))
sf_female68 <- survfit(coxmod_clean, newdata = female68)
keep68f <- sf_female68$time >= 55
plot(sf_men68$time[keep68f],
sf_men68$surv[keep68f],
type = "l",
xlab = "Age",
ylab = "Predicted survival probability",
main = "Predicted Survival Curve: Men with 6-8 Years of Education (55+)",
col = "red",
lwd = 2)
What this all means.
- Putting these all together (see below), we see that women consistently have a higher predicted survival probability than men do. This is a stronger predictor than education, but the curves for the higher education group is actually somehow lower than for those with no education, meaning that those without education have a higher predicted survival probability.
plot(sf_men0$time[keep1], sf_men0$surv[keep1],
type = "l", lwd = 2, col = "blue",
xlab = "Age", ylab = "Predicted survival probability",
main = "Predicted Survival (55+): Gender × Education")
lines(sf_female0$time[keep2], sf_female0$surv[keep2],
lwd = 2, col = "blue", lty = 2)
lines(sf_men68$time[keep68m], sf_men68$surv[keep68m],
lwd = 2, col = "red")
lines(sf_female68$time[keep68f], sf_female68$surv[keep68f],
lwd = 2, col = "red", lty = 2)
legend("bottomleft",
legend = c("Men, 0 years", "Women, 0 years",
"Men, 6–8 years", "Women, 6–8 years"),
col = c("blue","blue","red","red"),
lwd = 2,
lty = c(1,2,1,2))
Write a short write up with a reflection on the problem set, either describing and interpreting the results or doing something more “meta” about your learning of the techniques, whatever is most useful to you.
- Reflecting on this assignment, it makes me realize how incredibly difficult interpreting the relationships between these covariates can be. In some instances, it looks apparent that education is a somewhat significant predictor of mortality and that those with greater education benefit. In others, it looks like that relationship is inverted and that education has no effect, or even a negative effect, on mortality. My main takeaway is that we have to undertake basically every method we can for which assumptions aren’t violated and do our best to be upfront about conflicting outcomes and shortcomings of our conclusions. In sum, this was hard.