Interaction plots in Lifereg (parametric regression)

[OTs]

Hello, SAS community,

I am in the final stage of the research revision.

I received needed results where the main effect and interaction are significant.

I want to illustrate the interaction effects on plots similarly to those reported in other regressions.

In SAS 14.1 User Guide, on page 5113 Output 69.8.4 provides examples of plots that I wanted to develop; however, the code was not provided, and in the 14.2 version, it was removed.

I have a moderator which is a binary variable (Gender), and also the independent variables, which are also binary and reflect if the subject got treatment (Pill1 and/or Pill2, or both) at time t or not. My model also includes a dozen continuous control variables.

proc lifereg data=data;

class Gender;

model time*Died(0)=  Gender*Pill1 Gender*Pill2 Pill1 Pill2 Gender controls/dist=weibull;

run;

Gender has a significant negative effect on survival time, and pills have both significant and positive effects on survival time. However, the interaction of Gender with both pills is positive and significant.

I want to illustrate the positive effect on survival of taking pill1 and pill2 and the positive interaction of these predictors with the moderator variable (Gender) on plots.

Please share the best examples available.

[PaigeMiller]

When you run PROC LIFEREG, use the STORE statement to store the fitted model.

Then you can use PROC PLM to create interaction plots via the EFFECTPLOT statement.

[OTs]

Thank you Paige for the suggestion.

The suggested method did not work with lifereg and also phreg, only logistic I was able to get the interaction plot. See the screen below.

Are there any other ways of getting an interaction plot for Lifereg models?

[PaigeMiller]

Looks like my mistake, EFFECTPLOT doesn't work with LIFEREG. Maybe @StatDave or @Rick_SAS can suggest another method.

[StatDave]

Current documentation does not include the EFFECTPLOT statement as valid syntax and does not include that plot. However, you can obtain a plot of predicted quantile values against a predictor by saving the quantile estimates using the OUTPUT statement. Using the data set in that example (see the link at the top of the example to find the full code and data set), the following fits the model and saves the estimated median (quantile=.5) values (PRED) for plotting. To do this,  add a control variable (C) to the data with the value 1 for all observations. Then, CONTROL=C in the OUTPUT statement requests predicted quantiles for all observations. The estimated time medians at each stage and age are plotted using the SERIES statement in PROC SGPLOT. The SCATTER statement adds the observed time values. You could, of course, produce a similar plot using estimates at other quantiles of time.

With additional variables in the model, you would want to add observations to the input data set that sets the other variables at specified values, sets the time response to missing, and allows only the predictor of interest to vary. This is illustrated in the first (Motorette) example. You could then use just those observations in the plot to show the effect of the variable of interest on the time quantile when the other variables are at that setting.

proc lifereg data=Larynx order=data;
   class Stage;
   model Time*Death(0) = age stage;
   output out=out p=pred quantile=.5 control=c;
   run;
proc sort data=out; by stage age; run;
proc sgplot; 
   series  y=pred x=age/group=stage;
   scatter y=time x=age/group=stage;
   run;

[OTs · Accepted solution]

Thank you, @StatDave, for the suggestions; they were helpful. I needed to modify code as I was looking for plotting interaction effects, which had a lot of controls and a huge sample size. The code below allows plotting interaction, accommodates controls, and reduces output data for potting.

proc lifereg data= mydata  order=data;

model time*Died(0)=  Pill|Gender control1avr control2avr control3avr /dist=weibull;
* I had a lot of controls, but they do not effect here ;
output out=out p=pred quantile=.05 control=c;
* as only 20% died in my data, I used such a small quantile;
run; 

* This code can be useful for folks who deal with huge samples based on which SAS files
to build graphs due to limits in my sample, 20 million, so I needed to reduce it;
 PROC SURVEYSELECT DATA=out OUT=out METHOD=SRS
  SAMPSIZE=10000 SEED=1999;
  samplingunit id;
  RUN;
* if your sample smaller then 500k do not use it;

proc sort data=out; by Gender Pill; run;

proc sgplot; 
title "Interaction effect";
label Gender="Gender";
   series  y=pred x=Pill/group=Gender ;
xaxis values=(0 1) label="Pill: Placebo=0 Treatment=1 ";
yaxis label=" Predicted survival time, days";
INSET   ("Quantile"="=.05") / BORDER  ;
run;

* try different quantiles also if you have more than one IV,
replace not engaged in interaction IVs with their average levels;

Result (I use the Journal setting (Tools/Options/Preferences)