5
$\begingroup$

Two linear predictors interact significantly (see below). How can I visualize this interaction in a plot?

> data(pbc)
> d <- pbc
> rm(pbc, pbcseq)
> d$status <- ifelse(d$status != 0, 1, 0)
> 
> dd = datadist(d)
> options(datadist='dd')
> (m <- cph(Surv(time, status) ~ bili * alk.phos, data=d))

Cox Proportional Hazards Model

cph(formula = Surv(time, status) ~ bili * alk.phos, data = d)

Frequencies of Missing Values Due to Each Variable
Surv(time, status)               bili           alk.phos 
                 0                  0                106 


                   Model Tests       Discrimination    
                                        Indexes        
Obs      312    LR chi2     94.76    R2       0.264    
Events   144    d.f.            3    Dxy      0.565    
Center 0.676    Pr(> chi2) 0.0000    g        0.641    
                Score chi2 193.72    gr       1.898    
                Pr(> chi2) 0.0000                      

                Coef   S.E.   Wald Z Pr(>|Z|)
bili            0.2280 0.0300  7.59  <0.0001 
alk.phos        0.0001 0.0000  1.83  0.0667  
bili * alk.phos 0.0000 0.0000 -2.86  0.0043

One way I can think of is to dichotomize one predictor and plot the high values with the low values as two line plots in one figure. However, I cannot reproduce the example in the rms package under plot.Predict() using the example data above.

> d$alk.phos.high <- ifelse(d$alk.phos > 1259, 1, 0)
> (m <- cph(Surv(time, status) ~ bili * alk.phos.high, data=d))

Cox Proportional Hazards Model

cph(formula = Surv(time, status) ~ bili * alk.phos.high, data = d)

Frequencies of Missing Values Due to Each Variable
Surv(time, status)               bili      alk.phos.high 
                 0                  0                106 


                  Model Tests       Discrimination    
                                       Indexes        
Obs     312    LR chi2     97.95    R2       0.272    
Events  144    d.f.            3    Dxy      0.540    
Center 0.81    Pr(> chi2) 0.0000    g        0.727    
               Score chi2 194.00    gr       2.069    
               Pr(> chi2) 0.0000                      

                     Coef    S.E.   Wald Z Pr(>|Z|)
bili                  0.2374 0.0277  8.57  <0.0001 
alk.phos.high         0.5667 0.2214  2.56  0.0105  
bili * alk.phos.high -0.1139 0.0309 -3.69  0.0002

UPDATE #1

My trying to figure out how to plot both groups of a dichotomized predictor in a single figure, I figured out how to plot lines for several values of one interacting predictor in a plot of the other predictor against the hazard ratio.

I thinks this kind of plot shows the effect of the interaction term in an easy to understand way (which is especially important for a physician g).

  1. Does this kind of plot have a special name? How would you call this kind of plot?
  2. How can I interpret this interaction? Would it be correct to say that the prognostic impact of bilirubin increases with lower values for alkaline phosphatase?
#
library(rms)
data(pbc)
d <- pbc
rm(pbc, pbcseq)
d$status <- ifelse(d$status != 0, 1, 0)
dd = datadist(d)
options(datadist='dd')

m1 <- cph(Surv(time, status) ~ bili * alk.phos, data=d)
p1 <- Predict(m1, bili, alk.phos=c(850, 1250, 2000), conf.int=FALSE, fun=exp)
plot(p1, ylab="Hazard Ratio")


m2 <- cph(Surv(time, status) ~ bili + alk.phos, data=d)
p2 <- Predict(m2, bili, alk.phos=c(850, 1250, 2000), conf.int=FALSE, fun=exp)
plot(p2, ylab="Hazard Ratio")

first figure: model m2 without interaction

second figure: model m1 with interaction

mmodel #2 without interaction model #1 with interaction

$\endgroup$
  • $\begingroup$ I think alk_phos by bili interactions are very interesting; I will try to construct a response. I don't think I would use discretized variables, though. It would not be pleasing to rms's author. Better would be to visualize the two-D plot of a cross spline fit. I've done this with my own data research and the results are much more useful if you leave the data as continuous. $\endgroup$ – DWin Aug 4 '15 at 4:26
  • $\begingroup$ Would you please give me a paper reference of your own research where you used this? $\endgroup$ – Gurkenhals Aug 4 '15 at 14:29
  • 2
    $\begingroup$ @Gurkenhals: Your update shows another good way of visualizing the relationships predicted by the model. Again a little explanation would be helpful - & what are you asking about this addition? NB Any mucking around with a good model by dichotomizing predictors to aid visualization is IMO putting the cart before the horse - see What is the benefit of breaking up a continuous predictor variable?. $\endgroup$ – Scortchi - Reinstate Monica Aug 4 '15 at 14:49
  • $\begingroup$ I explained my update further. $\endgroup$ – Gurkenhals Aug 4 '15 at 18:30
  • $\begingroup$ Also related $\endgroup$ – AdamO Aug 5 '15 at 19:22
6
$\begingroup$

The rms package allows you to model interactions between continuous variables very flexibly. This is a demonstration of modeling crossed regression splines with that data:

(m <- cph(Surv(time, status) ~ rcs(bili,3) * rcs(alk.phos,3), data=d))
bplot(Predict(m, bili=seq(.5,25,by=.5), alk.phos=seq(300,10000, by=500), fun=exp),
      lfun=contourplot, at=1:20)

You can find similar code described in Frank Harrell's "Regression Modeling Strategies" in chapter 10; Section 10.5:'Assessment of Model Fit'. That code was in the Binary Regression chapter but there is no reason it cannot be used in a survival analysis context. One caveat that needs to be added is that this plot extends to regions of the bili-by-alk-phos space where there is no data, especially in the upper right quadrant. I will add a further example code section and plot section that shows how to restrict the contours to regions actually containing data and will make some further comments on interpretation when time allows.

enter image description here

> anova(m)
                Wald Statistics          Response: Surv(time, status) 

 Factor                                         Chi-Square d.f. P     
 bili  (Factor+Higher Order Factors)            141.86     6    <.0001
  All Interactions                                7.48     4    0.1128
  Nonlinear (Factor+Higher Order Factors)        36.61     3    <.0001
 alk.phos  (Factor+Higher Order Factors)          8.17     6    0.2261
  All Interactions                                7.48     4    0.1128
  Nonlinear (Factor+Higher Order Factors)         3.04     3    0.3854
 bili * alk.phos  (Factor+Higher Order Factors)   7.48     4    0.1128
  Nonlinear                                       6.95     3    0.0735
  Nonlinear Interaction : f(A,B) vs. AB           6.95     3    0.0735
  f(A,B) vs. Af(B) + Bg(A)                        0.13     1    0.7195
  Nonlinear Interaction in bili vs. Af(B)         4.42     2    0.1096
  Nonlinear Interaction in alk.phos vs. Bg(A)     2.75     2    0.2534
 TOTAL NONLINEAR                                 53.97     5    <.0001
 TOTAL NONLINEAR + INTERACTION                   61.75     6    <.0001
 TOTAL                                          147.36     8    <.0001

In my own work using cph with tens of thousands of events I generally increase the default values of n for the perimeter function, but in this very small dataset I found it necessary to decrease the value to get a good example. I used a simple call to plot to see where the data actually extended.

with(d, plot(bili, alk.phos, col=status+1)) 
 # needed to ad one to status to get both deaths and censored.
perim <- with(d, perimeter(bili, alk.phos, n=2))
bplot(Predict(m, bili=seq(.5,25,by=.5), alk.phos=seq(300,10000, by=500), fun=exp),
      lfun=contourplot, at=1:20, perim=perim)

enter image description here

The plot shows that in the region where bilirubin is less than 4 or 5, that the risk increases primarily with increasing bilirubin with relative risks in the range of 1-3, but \ when bilirubin is greater than 5, that the risk increases to relative risks of 4 to 9 under the joint effect of decreasing alkaline phosphatase and bilirubin. Clinical interpretation: In the higher bilirubin "domain" a low alk.phos has what might be considered a paradoxical effect. There may be such substantial loss of hepatic substance that the remaining liver tissue is so depleted that it is unable to produce as much of the enzyme. This is a typical example of "effect modification" using the terminology of epidemiology. The higher "effect" associated with a given decrease in alk.phos is fairly small when bilirubin is low, but the alk.phos-"effect" gets magnified as bili-rubin rises.

I also answered a similar question on StackOverflow albeit not with a focus on the rms facilities: https://stackoverflow.com/questions/14928716/how-to-plot-a-regression-plane-with-an-interaction-in-rgl/14932045#14932045

$\endgroup$
  • $\begingroup$ (+1) But I think it would be very helpful to non-users of rms/R if you could add a little explanation. $\endgroup$ – Scortchi - Reinstate Monica Aug 4 '15 at 11:34
  • $\begingroup$ Thanks for your answer; but yes a little bit more explanation how to read and interpret the plot and the anove would be very helpful. I just updated my question showing two plots with and without interaction (does this kind of plot has a specific name?) $\endgroup$ – Gurkenhals Aug 4 '15 at 14:28
  • $\begingroup$ @Scortchi: I see from your response to the linked question about "binning" that we are both committed disciples of Frank Harrell's wisdom regarding the value of leaving continuous variables continuous. Your answer is getting a well-deserved response from the audience. Hope these additions rise to your challenge for better explanation. Feel free to suggest areas needing clarification. $\endgroup$ – DWin Aug 4 '15 at 19:10
  • 1
    $\begingroup$ @Gurkenhals: I think these plots would be considered "interaction plots". Your individual lines can be thought of as "slices" through a 3D (2D-covariate by 1D risk-level) plot of risk. The bplot function also allows you to specify the wireframe function for plotting, which I find useful when presenting to audiences that have not seen this before. $\endgroup$ – DWin Aug 4 '15 at 19:13
  • 1
    $\begingroup$ anova computes Wald statistics, with optional use of robust sandwich and bootstrap covariance estimates, or estimates adjusted for multiple imputation. The LR $\chi^2$ is the silver standard and should be used for key results. It doesn't suffer from the Hauck-Donner effect, among other things. $\endgroup$ – Frank Harrell Aug 5 '15 at 11:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.