0

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!

MarianD
  • 1,493
  • 2
  • 8
  • 17
lj_growl
  • 3
  • 2
  • 1
    The way you specified your model doesn't include an interaction term. Your model implicitly assumes that the effect of smoking doesn't differ with alcohol use, and vice versa. In R you can include an interaction by adding a term `smokef:alcoholf` to your model, or you can just use `smokef*alcoholf` to replace your additive `smokef+alcoholf` and it will be expanded out to include the individual and interaction terms. I suspect that you will still have questions after that, so please edit your question to show the revised model and to state your remaining questions. – EdM Nov 21 '20 at 16:59
  • @EdM I have updated the model! That makes sense and thank you for catching that. I am still not understanding how to extract the estimated cumulative hazards and then survival functions for this model. I guess that I could use the exp(-H) to get the estimated survivals, but I'm not understanding enough to get there. Thanks for any help! – lj_growl Nov 21 '20 at 17:10
  • As this seems to be a homework or similar type of problem, please add the [self-study](https://stats.stackexchange.com/tags/self-study/info) tag to your question and read the information in that link about the type of help we can provide. – EdM Nov 21 '20 at 17:18
  • @EdM added! Thanks. Any nod in the correct direction would be helpful in understanding this model! I understand how to extract the cumhaz functions for coxph models, just not fully understanding the Aalen model. – lj_growl Nov 21 '20 at 17:24
  • There are usually `predict` functions for regression models in R. and if not then `summary` methods. For objects of class:aareg I see neither. The code for plot.aareg is obtained with `getAnywhere(plot.aareg)` – DWin Nov 25 '20 at 00:57

1 Answers1

0

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.

EdM
  • 57,766
  • 7
  • 66
  • 187