机器算法验证 - 如何使用 rms 包可视化两个线性预测变量之间的显着交互？ - 吾爱随笔录

如何使用 rms 包可视化两个线性预测变量之间的显着交互？

机器算法验证 r 相互作用 cox模型回归策略

2022-04-12 20:16:42

两个线性预测变量显着交互（见下文）。如何在绘图中可视化这种交互？

> 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

我能想到的一种方法是将一个预测变量二分法，并将高值和低值绘制为一张图中的两条线图。但是，我无法plot.Predict()使用上面的示例数据重现 rms 包中的示例。

> 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

更新#1

我试图弄清楚如何在一个图中绘制两组二分预测变量，我想出了如何在另一个预测变量与风险比的关系图中绘制一个交互预测变量的多个值的线。

我认为这种图以易于理解的方式显示了交互项的效果（这对医生g尤其重要）。

这种情节有什么特别的名字吗？这种情节怎么称呼？
我如何解释这种互动？说胆红素的预后影响随着碱性磷酸酶值的降低而增加是否正确？

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

第一张图：m2没有交互的模型

第二个图：m1有交互的模型

1个回答

该rms包允许您非常灵活地对连续变量之间的交互进行建模。这是使用该数据对交叉回归样条进行建模的演示：

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

您可以在 Frank Harrell 的第 10 章的“回归建模策略”中找到类似的代码；第 10.5 节：“模型拟合评估”。在那里，他使用了伪 3D 图，这在交互的可视化中也很有用。该代码在二元回归一章中，但没有理由不能在生存分析 ( cph) 上下文中使用它。需要添加的一个警告是，该图扩展到没有数据的 bili-by-alk-phos“空间”区域，尤其是在右上象限。我将添加一个进一步的示例代码部分和绘图部分，展示如何使用 rmsperim函数将轮廓限制为实际包含数据的区域，并将对解释做一些进一步的评论。

> 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

在我自己使用cph数万个事件的工作中，我通常会增加n周界函数的默认值，但在这个非常小的数据集中，我发现有必要减少该值以获得一个好的示例。我使用了一个简单的调用 plot 来查看数据实际扩展的位置。

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))
# perim constructs a set of ranges where subjects are actually present
bplot( Predict(m, bili=seq(.5,25,by=.5), alk.phos=seq(300,10000, by=500),
      fun=exp),  # contours of relative risk
      lfun=contourplot, # call to lattice function
      at=1:20,          # levels for contours
      perim=perim)      # regions to include

该图显示，在胆红素小于 4 或 5 的区域，风险主要随着胆红素的增加而增加，相对风险在 1-3 范围内，但当胆红素大于 5 时，风险增加至相对风险。在碱性磷酸酶和胆红素降低的共同作用下，风险为4至9 。临床解释：在较高胆红素“域”中，低alk.phos具有可能被认为是自相矛盾的效果。可能存在功能性肝组织的大量损失，以致剩余的肝脏如此减少以至于不能产生尽可能多的碱性磷酸酶。这是使用流行病学术语进行“效果修正”的典型例子。与给定减少相关的更高“效果”alk.phos当胆红素低时，它相当小，但随着-rubin 的升高alk.phos，“效应”会被放大。bili

其它你可能感兴趣的问题

上一篇R 函数 arima() 如何计算其残差？下一篇作为参数的 Spearman 相关性