Estimated survival functions for interaction terms in Aalen Model, Cox Ph Model, and Kaplan Meier

Question

I have this question:

For each of the four interaction groups

    (1)  no smoking, no alcohol,
    (2)  smoking, no alcohol,
    (3)  no smoking, alcohol,
    (4)  smoking, alcohol

– produce plots of the estimated survival function based upon Aalen’s linear model, the Cox model fit with smoking status and alcohol use as covariates, the Cox model fit with smoking status and alcohol use as stratifying variables, and the unadjusted Kaplan-Meier estimates. Overlay the plots from the different approaches in single graphs for each of the 4 groups in different colors, i.e. you should produce 4 graphs, each with 4 estimated survival functions.

Smoke is a binary variable with 1=yes and 0=no. Alcohol is a binary variable with 1=yes and 0=no. I am struggling to get the interaction terms in this way. My output is this when I run the Aalen model with code

feed.amod <- aareg(feed.surv~smokef*alcoholf+poverty+racef):

Call:
aareg(formula = feed.surv ~ smokef * alcoholf + race + poverty)

  n= 927 
    35 out of 48 unique event times used

                                   slope      coef se(coef)      z        p
Intercept                        0.08770  0.001410 1.63e-04  8.610 7.02e-18
smokefsmokeNo                    0.04260  0.000531 1.68e-04  3.160 1.60e-03
alcoholfalcoholNo                0.00756  0.000252 3.30e-04  0.765 4.44e-01
race                             0.02980  0.000318 9.98e-05  3.190 1.42e-03
poverty                         -0.01300 -0.000253 1.74e-04 -1.450 1.46e-01
smokefsmokeNo:alcoholfalcoholNo  0.04170  0.000227 5.38e-04  0.421 6.74e-01

Chisq=21.32 on 5 df, p=0.000704; test weights=aalen

I am just needing help understanding how to extract the cumulative hazards for this Aalen model. From there, I am understanding that I need to use the exp(-H) to extract the estimated survival. Any help or nod in the right direction is appreciated!

Answer 1

First, you are not being asked to plot cumulative hazards. You are being asked to plot the survival function, the survival probability as a function of time. Those are not the same thing. See for example this page on this site, or this explanation from the NIST handbook. The hazard at a particular time is based upon having survived up to that time. Survival functions, in contrast, represent the overall survival probability as a function of time, just based on having entered the study at time = 0.

Now for the types of analysis. Questions about implementation in specific software like R are off-topic on this site, but here are hints about how to proceed. (I expect that you already know much of this, but others who come across your question might not.)

Unadjusted Kaplan-Meier estimates are just raw survival curves broken down into the 4 combinations of alcohol and smoking, without correction for other covariates. In R, that analysis might be done with the survfit() function.*

The Cox models assume proportional hazards around baseline hazard(s). If you include alcohol and smoking as covariates you assume a single baseline hazard for all cases. If you stratify by alcohol and smoking, you get separate baseline hazards associated with alcohol and smoking, with the same hazard ratios for the other covariates around all those baseline hazards. In R, one might use the coxph() function with corresponding specification of strata terms. In this case the plots will be predicted plots for specified values of the covariates.

As the help page for the R aareg() function that you used puts it:

The Aalen model assumes that the cumulative hazard H(t) for a subject can be expressed as a(t) + X B(t), where a(t) is a time-dependent intercept term, X is the vector of covariates for the subject (possibly time-dependent), and B(t) is a time-dependent matrix of coefficients. The estimates are inherently non-parametric; a fit of the model will normally be followed by one or more plots of the estimates.

Software typically provides functions like plot() or predict() to display the survival function for each type of analysis. It's simplest to use the corresponding built-in plotting function to display the results of these models. If you want to see how the prediction and plotting is implemented and you are using open-source software like R, then you can examine the source code.

---

*I vaguely remember that you don't need to specify an interaction with that function to get all 4 combinations, but my memory isn't what it used to be.