¿Cómo visualizar una interacción significativa entre dos predictores lineales usando el paquete rms?

Question

Dos predictores lineales interactuar de forma significativa (véase a continuación). ¿Cómo puedo visualizar esta interacción en una parcela?

> 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

Una forma en que puedo pensar es en dividir un predictor y graficar los valores altos con los bajos valores como dos líneas de puntos en una sola figura. Sin embargo, no puedo reproducir el ejemplo en la rms paquete en plot.Predict() utilizando los datos del ejemplo anterior.

> 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

ACTUALIZACIÓN #1

Mi tratando de averiguar cómo parcela de ambos grupos de un dividió en dos predictor en una sola figura, me di cuenta de cómo trazar líneas para varios valores de una interacción predictor en una parcela de el otro predictor en contra de la relación de riesgo.

Yo pienso que este tipo de gráfico muestra el efecto de la interacción término en una forma fácil de entender (lo cual es especialmente importante para un médico g).

1. Este tipo de trama tiene un nombre especial? ¿Cómo le llaman a este tipo de argumento?
2. ¿Cómo puedo interpretar esta interacción? Sería correcto decir que el impacto pronóstico de la bilirrubina aumenta con los valores más bajos de fosfatasa alcalina?

#

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")

primera figura: modelo m2 sin interacción

segunda figura: modelo m1 con la interacción

Answer 1

El rms paquete permite modelar las interacciones entre las variables continuas muy flexible. Esta es una demostración de modelado cruzado splines de regresión con los datos:

(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)

Usted puede encontrar un código similar se describe en Frank Harrell "modelos de Regresión de Estrategias" en el capítulo 10; Sección 10.5: "Evaluación del Ajuste del Modelo'. Que el código fue en el Binario de Regresión capítulo, pero no hay ninguna razón por la que no se puede utilizar en un análisis de supervivencia de contexto. Una advertencia que debe ser agregado es que esta trama se extiende a las regiones de la bili-por-alk-phos el espacio donde no hay datos, especialmente en el cuadrante superior derecho. Voy a agregar un nuevo código de ejemplo de la sección y trazado de la sección que se muestra cómo restringir los contornos de las regiones en realidad contiene los datos y se hacen algunas observaciones acerca de la interpretación cuando el tiempo lo permite.

> 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

En mi propio trabajo utilizando cph con decenas de miles de eventos que generalmente aumentan los valores predeterminados de n para el perímetro de la función, pero en este pequeño conjunto de datos que he encontrado que es necesario disminuir el valor para obtener un buen ejemplo. He utilizado una simple llamada a la trama para ver donde están los datos de la realidad extendida.

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)

El gráfico muestra que en la región donde la bilirrubina es de menos de 4 o 5, que el riesgo aumenta principalmente con el aumento de la bilirrubina con riesgos relativos en el rango de 1-3, pero \ cuando la bilirrubina es mayor que 5, que aumenta el riesgo para los riesgos relativos de 4 a 9 bajo el efecto conjunto de la disminución de la fosfatasa alcalina y de la bilirrubina. Interpretación clínica: En los altos de bilirrubina "dominio" un bajo alk.phos tiene lo que podría considerarse un efecto paradójico. Puede haber contribuido a la pérdida de la sustancia que el resto de los tejidos del hígado es tan agotado que es incapaz de producir la mayor cantidad de la enzima. Este es un típico ejemplo de "efecto de la modificación", utilizando la terminología de la epidemiología. El mayor "efecto" asociado con una disminución en el alk.phos es bastante pequeña cuando la bilirrubina es baja, pero el alk.phos-"efecto" se amplía a medida bili-rubin se eleva.