# Combining individual terms from the predict() function

I want to merge several variables' estimate and also calculate the confidence intervals from the survival:::predict.coxph function's output. Just like the the predict.lm the function allows for estimating specific term estimates and sample errors for specific variables. The rms::contrast function has a function that allows contrasting one outcome with another and handles this very nicely, unfortunately under some circumstances I've found that I can't use the rm::contrast and it would be nice to use the predict function for estimating the same (or similar combined estimate). I'm particularly interested in this for survival functions as this becomes a must when a variable has a time-interaction estimate that needs to merge with the base estimate.

Getting the estimate is a simple sum exercise of the coefficient of interest:

$$\sum_{i=1}^{n}{est_i}$$

I originally thought that the variance could be calculated by

$$\sqrt{\sum_{i=1}^{n}{\hat{\sigma}_i^2}}$$

where $n$ is the number of terms. But as this doesn't really handle the complexity of predictor covariance, thus the formula should be:

$$Var(a*X + b*Y) = a^2*Var(X) + b^2*Var(X) + 2*a*b*Cov(X,Y)$$

This does make everything much more complicated, especially when one of the variances is a spline. I was hoping to use the predict function for getting the raw estimates and then merging them in some fashion:

library(survival)
library(magrittr)
library(dplyr)
data("lung", package = "survival")
lung$sex <- factor(lung$sex,
labels = c("Male", "Female"))
model <- coxph(Surv(time, status) ~ sex + pspline(age) + wt.loss, data=lung)

new_data <- data.frame(sex = rep(levels(lung$sex), each = 100), age = rep(seq(40, 80, length.out = 100), times = 2), wt.loss = median(lung$wt.loss, na.rm=TRUE))

m_pred <- predict(model,
newdata = new_data,
type = "terms",
terms = names(model$assign)[grep("sex|age", names(model$assign))],
se.fit = TRUE)

lapply(m_pred, function(x) knitr::kable(head(x), align = "r"))


This gives the result:

$fit | sex| pspline(age)| |---------:|------------:| | 0.2135748| -0.7408014| | 0.2135748| -0.7064531| | 0.2135748| -0.6721696| | 0.2135748| -0.6379620| | 0.2135748| -0.6038413| | 0.2135748| -0.5698185|$se.fit

|       sex| pspline(age)|
|---------:|------------:|
| 0.0705246|    0.6129619|
| 0.0705246|    0.5727195|
| 0.0705246|    0.5346473|
| 0.0705246|    0.4988786|
| 0.0705246|    0.4655234|
| 0.0705246|    0.4346488|


Here's what we get if we simply add these (Note: this is not the correct way):

new_data$fit <- rowSums(m_pred$fit)

new_data$se.fit <- apply(m_pred$se.fit, 1, function(x) x^2) %>%
colSums() %>% # apply flips the matrix, therefore the colSums
sqrt()

library(ggplot2)
new_data %>%
mutate(risk = exp(fit),
upper = exp(fit + 1.96*se.fit),
lower = exp(fit - 1.96*se.fit)) %>%
ggplot(aes(y = risk, x = age, col = sex, fill = sex)) +
geom_ribbon(aes(ymax = upper, ymin = lower), alpha = .5) +
scale_fill_hue(c=45, l=80) +
scale_color_hue(c=45, l=80) +
scale_y_log10(breaks=2^(-2:3)) +
scale_x_continuous(expand = c(0,0)) +
geom_line(lwd = 2, alpha = .5, color = "black") +
theme_bw() +
ylab("Hazard ratio") +
xlab("Age (years)")


This plot is actually not that dissimilar from what I can get from the rms::contrast. It doesn't though have a waist and I'm uncertain of how to interpret the confidence interval, I guess they're the confidence interval for each individual estimated hazard ratio. Here's the rms version of above:

library(rms)
# The pspline gives: Error in X[, mmcolnames, drop = FALSE] : subscript out of bounds
cph_model <- cph(Surv(time, status) ~ sex + rcs(age, 4) + wt.loss, data=lung)

# Now here's the contrast funciton in action
# We set the age to the age span of interest
# and then we specify that we're interested
# in the two levels of the sex variable
# As reference we set the age to 60 and the
# sex to female
cntr <-
contrast(cph_model,
a = list(age = seq(40, 80, length.out = 100),
sex = levels(lung$sex)), b = list(sex = "Female", age = 60)) # The contrast function returns a list that needs to be reformatted # Note that the sex is interchanging while age is continuous # - If you're unfamiliar with the coding I recommend extracting # one sex at the time and then merging them into one data.frame(age = rep(seq(40, 80, length.out = 100), each = 2), sex = rep(levels(lung$sex), times = 100),
risk = exp(cntr$Contrast), upper = exp(cntr$Upper),
lower = exp(cntr$Lower)) %>% ggplot(aes(y = risk, x = age, col = sex, fill = sex)) + geom_ribbon(aes(ymax = upper, ymin = lower), alpha = .5) + scale_fill_hue(c=45, l=80) + scale_color_hue(c=45, l=80) + scale_y_log10(breaks=2^(-2:3)) + scale_x_continuous(expand = c(0,0)) + geom_line(lwd = 2, alpha = .5, color = "black") + theme_bw() + ylab("Hazard ratio") + xlab("Age (years)")  Gives the following plot: ## My intuition so far on the subject In my mind the variance combined from different estimates should only be a form of re-scaling, i.e. the variance is dependent on the unit and the unit can either be age in years or age in 20 years when comparing 40 years with 60 years. The 20 years age difference is simply a constant, i.e. I should multiply the variance by: $$Var(c*\hat{\sigma}) = c^2 * Var(\hat{\sigma}^2) = 20^2 * Var(\hat{\sigma}^2)$$ The properties of the variance extends nicely to summing the variance between age and sex should be (b = 1): $$Var(a*Age + b*Sex) = a^2*Var(Age) + Var(Sex) + 2*a*Cov(Age,Sex)$$ I understand that the key is in the variance-covariance matrix. Looking at the low covariance estimates between age and sex the similarity of the curves is not surprising: >vcov(cph_model) sex=Female age age' age'' wt.loss sex=Female 0.0305 0.0004 -0.0006 0.0029 0e+00 age 0.0004 0.0013 -0.0024 0.0101 0e+00 age' -0.0006 -0.0024 0.0054 -0.0250 0e+00 age'' 0.0029 0.0101 -0.0250 0.1261 -1e-04 wt.loss 0.0000 0.0000 0.0000 -0.0001 0e+00  I guess that if my reasoning is correct it is rather straight forward to generate contrast from pure linear variables. As splines are inherently difficult to calculate by hand it seems that the predict with the terms attribute should provide a combined variance for age. Unfortunately I won't be able to get the covariance for the formula above between the two or can I get around this in some clever way? ## By request Here's the variance-covariance matrix for the first model: > round(vcov(model), 4) sexFemale ps(age)3 ps(age)4 ps(age)5 ps(age)6 ps(age)7 ps(age)8 ps(age)9 ps(age)10 ps(age)11 ps(age)12 ps(age)13 sexFemale 0.0308 -0.0024 -0.0046 -0.0059 -0.0054 -0.0023 -0.0004 -0.0016 -0.0030 -0.0022 0.0002 0.0045 ps(age)3 -0.0024 0.5458 0.8321 0.9147 0.9022 0.8801 0.8734 0.8749 0.8764 0.8764 0.8760 0.8755 ps(age)4 -0.0046 0.8321 1.4123 1.5816 1.5576 1.5130 1.4993 1.5023 1.5053 1.5054 1.5047 1.5037 ps(age)5 -0.0059 0.9147 1.5816 1.8774 1.8551 1.7869 1.7640 1.7679 1.7728 1.7731 1.7723 1.7711 ps(age)6 -0.0054 0.9022 1.5576 1.8551 1.9087 1.8295 1.7910 1.7923 1.7988 1.7999 1.7993 1.7983 ps(age)7 -0.0023 0.8801 1.5130 1.7869 1.8295 1.8157 1.7658 1.7549 1.7594 1.7617 1.7622 1.7620 ps(age)8 -0.0004 0.8734 1.4993 1.7640 1.7910 1.7658 1.7743 1.7474 1.7394 1.7405 1.7423 1.7434 ps(age)9 -0.0016 0.8749 1.5023 1.7679 1.7923 1.7549 1.7474 1.7741 1.7497 1.7393 1.7389 1.7402 ps(age)10 -0.0030 0.8764 1.5053 1.7728 1.7988 1.7594 1.7394 1.7497 1.7783 1.7535 1.7394 1.7314 ps(age)11 -0.0022 0.8764 1.5054 1.7731 1.7999 1.7617 1.7405 1.7393 1.7535 1.7859 1.7620 1.7361 ps(age)12 0.0002 0.8760 1.5047 1.7723 1.7993 1.7622 1.7423 1.7389 1.7394 1.7620 1.8293 1.8429 ps(age)13 0.0045 0.8755 1.5037 1.7711 1.7983 1.7620 1.7434 1.7402 1.7314 1.7361 1.8429 2.0834 ps(age)14 0.0089 0.8750 1.5028 1.7699 1.7973 1.7618 1.7446 1.7415 1.7238 1.7102 1.8539 2.3164 wt.loss 0.0000 -0.0001 -0.0002 -0.0002 -0.0003 -0.0004 -0.0003 -0.0003 -0.0003 -0.0002 -0.0003 -0.0003 ps(age)14 wt.loss sexFemale 0.0089 0e+00 ps(age)3 0.8750 -1e-04 ps(age)4 1.5028 -2e-04 ps(age)5 1.7699 -2e-04 ps(age)6 1.7973 -3e-04 ps(age)7 1.7618 -4e-04 ps(age)8 1.7446 -3e-04 ps(age)9 1.7415 -3e-04 ps(age)10 1.7238 -3e-04 ps(age)11 1.7102 -2e-04 ps(age)12 1.8539 -3e-04 ps(age)13 2.3164 -3e-04 ps(age)14 3.0339 -3e-04 wt.loss -0.0003 0e+00  As a minor note I commonly refer to log-hazard as hazard ratio, I know this is sloppy and I'm sorry for any confusion that it causes. ## Trying to implement the$\delta$method as suggested by AdamO. From the comment I would interpret that I should: $$Var(\delta hz = prediction_1 - prediction_2) = Var(prediction_1) * Var(prediction_2)*(logHazard(prediction_1) - logHazard(prediction_2))^2$$ This would be equivalent in R-code to: m40_data <- data.frame(sex = "Male", age = 40, wt.loss = median(lung$wt.loss, na.rm=TRUE))
f60_data <-
data.frame(sex = "Female",
age = 60,
wt.loss = median(lung$wt.loss, na.rm=TRUE)) pre_m40 <- predict(cph_model, newdata = m40_data, se.fit = TRUE) pre_f60 <- predict(cph_model, newdata = f60_data, se.fit = TRUE) log_hzrd <- pre_m40$linear.predictors -
pre_f60$linear.predictors comb_se <- sqrt(pre_m40$se.fit^2 *
pre_f60$se.fit^2 * log_hzrd^2) sprintf( "Log Hazard: %.3f, Combined s.e.: %.3f", log_hzrd, comb_se) %>% cat  Gives: Log Hazard: 0.155, Combined s.e.: 0.010  If we compare this to the contrast function the variance is not even close: cntr <- contrast(cph_model, a = list(age = 40, sex = "Male"), b = list(sex = "Female", age = 60)) sprintf( "Log Hazard: %.3f, Combined s.e.: %.3f", cntr$Contrast, cntr$SE) %>% cat  Gives: Log Hazard: 0.155, Combined s.e.: 0.545  I guess I've misunderstood something in the instructions, looking at the Wikipedia page I realize that I need to somehow get this formula to work: $$Var(h_r)=\sum_{i}{\left(\frac{\partial h_r}{\partial B_i}\right) Var(B_i) + \sum_{i}{\sum_{j \ne i}{\left(\frac{\partial h_r}{\partial B_i}\right)\left(\frac{\partial h_r}{\partial B_j}Cov(B_i,B_j)\right)}}}$$ While the formula looks complicated to an orthopaedic surgeon, I don't quite grasp the concept of the partial derivatives but it shouldn't be that complicated to convert into useful code. I guess the predict does this for the spline variables and the only covariances that actually are missing are the covariance of$Cov(age, sex)$but this of course expands to$Cov(age_1,sex), Cov(age_2, sex), ...$due to the spline. • Are you sure you are trying to predict hazard ratios and not failure times from this analysis? Commented Jan 4, 2016 at 21:23 • @AdamO: Yes, I'm pretty sure. I've used the flexsurv-package for failure times (if you mean by failure time - time to death and doing an adjusted Kaplan-Meier). Commented Jan 4, 2016 at 21:37 • predicted values would not generally be independent Commented Jan 5, 2016 at 15:30 • @Glen_b: I guess that is the underlying reason for the$2 * a * b* Cov(X, Y)$in the formula or am I missing something? Commented Jan 5, 2016 at 21:19 • My comment related to the variance calculation in your second paragraph (just after "the variance should be calculated by"). Possibly I misunderstand what you were saying there. Commented Jan 5, 2016 at 23:07 ## 1 Answer Like a contrast, a hazard ratio is a ratio of instantaneous hazards, or a difference of log hazards. Unlike a contrast, the referent group, the denominator hazard, is fixed. You can go one step further with this information. One can calculate a difference of hazard ratios (on the hazard ratio scale) using the delta method. I think most contrast implementations will give you this. Usually, you can tell this is the case because the contrast CIs are symmetric on the hazard ratio scale. Differences in hazard ratios are nonsensical to me. Differences in log hazard ratios become, themselves, hazard ratios between two groups of interest when exponentiated. This makes sense. When using contrast functionality, it's important to perform all calculations on the linear predictor scale before exponentiating the output. The "bottleneck" in your plot is given by the referent group: i.e. women who are 50 years old. So the predicted hazard ratio for men who are 70 is a ratio of their predicted hazard over that of women who are 50 years old. This is completely different than the phenomenon attributed to the differences between confidence and prediction intervals. Your predicted hazard ratio for 70 year old men would be interpreted as an interval in which, if the study were conducted an infinite number of times, would be expected to contain$100\times(1-\alpha)$% of other fitted hazard ratios. No matter how many times the study is redone, women who are 50 will always be the referent group, so their HR is always 1. That's why you observe the bottleneck in the plot. Therefore, there's no way to get rid of it. I also note that your predicted hazard ratios are symmetric on the hazard ratio scale. I believe that you are not accounting for the fact that the covariance matrix applies to model effects on the linear predictor scale (log-hazard), and the fitted "hazard ratio" is exponentiated. To calculate asymptotically correct confidence intervals on the hazard ratio scale, you will need to apply the$\delta$-method, i.e. $$\widehat{\log \left( \lambda_{t}(X_i)\right) - \log \left( \lambda_t(X_0) \right)} = \mathbf{X} \hat{\beta}$$ But the hazard ratio is given by: $$\frac{\lambda_{t}(X_i)}{\lambda_{t}(X_0)} = \exp \left( \mathbf{X} \hat{\beta} \right)$$ You can easily apply the delta method to obtain what the approximate variance of the transformation of linear predictors is. EDIT On further inspection, I suspect the issue has to do with how pspline codes the continuous variable age. Can you show the output of pspline? I bet it standardizes the variable to some reasonable average or median value, so that age=50 becomes 0 and therefore a referent value. • Thank you for your input. Unfortunately I've actually never studied the$\delta$method and Wikipedia is slightly beyond my grasp. I would greatly appreciate if you could perhaps show me how to apply the delta method in this example. As I understand I'm using the log-hazard and only converting it at the last step into a hazard ratio. I've added the vcov as requested. Commented Jan 5, 2016 at 15:15 • @MaxGordon In order to apply the delta method for an exponentiated outcome, you will calculate the log hazard ratio using the coefficients in the model as you have shown. You will also calculate the variance of the log hazard ratio using the covariance matrix as you have shown. A$\delta\$ method approximation to the variance of the hazard ratio is given by the product of the variance of the log hazard ratio times the hazard ratio squared. This should give the variance for an approximate normal distribution to the hazard ratio contrast. Commented Jan 5, 2016 at 18:52
• As you can see from my updated question I could not get this to work. I really appreciate your attempt at helping me out here, but I feel that I need to have some working code-example in order to understand before I can check this as the correct answer. Commented Jan 9, 2016 at 10:52