文献解析：生存数据和分类结局列线图的做法，史上最全

今天给大家分析一篇文献，作者是浙大邵逸夫医院的章仲恒副主任，论文的主要内容就是教大家怎么画列线图。写的比较详细，拜读之后顺手写个粗浅的解析分享给大家，依然还是建议大家去看原文哦，论文题目见文章最后。

全文包括2分类逻辑斯蒂回归，多分类逻辑斯蒂回归，生存分析的列线图的画法和解析，基本上涵盖了临床上常用的预测方法，所以很值得大家收藏。

接下来我一个一个给大家解析

作者先模拟了一个数据集，大概长这样：

 

这个数据集有1000个个案，6个变量，其中age和lac都是正态分布的连续自变量，sex和shock为因子，y是二分类结局，Y是多分类结局（3分类），所以我们如果用y做结局，就是拟合一个二分类逻辑斯蒂回归，用Y做结局就是多分类逻辑斯蒂回归。接下来我们一个一个看：

二分类结局的列线图画法

画列线图用到的包是rms。

第一步我们需要用datadist()函数来得到预测变量的分布，以此来规定我们的图形的尺度，紧接着拟合我们的逻辑回归模型（用到的函数为lrm），代码如下：

<span style="color:#222222"><code>ddist <- datadist(<span style="color:#114ba6">data</span>)
options(datadist=<span style="color:#00753b">'ddist'</span>)
mod.bi<-lrm(y~shock+lac*sex+age,<span style="color:#114ba6">data</span>)</code></span>

模型拟合好了之后，我们把这个模型直接喂给nomogram函数就行，代码如下：

<span style="color:#222222"><code>nom.bi <- nomogram(mod.bi,          
											lp.at=seq(-<span style="color:#a82e2e">3</span>,<span style="color:#a82e2e">4</span>,<span style="color:#114ba6">by</span>=<span style="color:#a82e2e">0.5</span>),
                     <span style="color:#8a7304"><span style="color:#114ba6">fun</span>=<span style="color:#a82e2e">function</span>(x)</span><span style="color:#a82e2e">1</span>/(<span style="color:#a82e2e">1</span>+exp(-x)),
                     <span style="color:#8a7304"><span style="color:#114ba6">fun</span>.at=<span style="color:#a82e2e">c</span>(<span style="color:#a82e2e">.001</span>,<span style="color:#a82e2e">.01</span>,<span style="color:#a82e2e">.05</span>,seq(<span style="color:#a82e2e">.1</span>,<span style="color:#a82e2e">.9</span>,<span style="color:#114ba6">by</span>=<span style="color:#a82e2e">.1</span>)</span>,.<span style="color:#a82e2e">95</span>,.<span style="color:#a82e2e">99</span>,.<span style="color:#a82e2e">999</span>),
                     funlabel=<span style="color:#00753b">"Risk of Death"</span>,
                     conf.int=c(<span style="color:#a82e2e">0.1</span>,<span style="color:#a82e2e">0.7</span>),
                     abbrev=TRUE,
                     minlength=<span style="color:#a82e2e">1</span>,lp=F)

plot(nom.bi,lplabel=<span style="color:#00753b">"Linear Predictor"</span>,
     <span style="color:#8a7304"><span style="color:#114ba6">fun</span>.side=<span style="color:#a82e2e">c</span>(<span style="color:#a82e2e">3</span>,<span style="color:#a82e2e">3</span>,<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">3</span>,<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">3</span>,<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">3</span>)</span>,
     col.conf=c(<span style="color:#00753b">'red'</span>,<span style="color:#00753b">'green'</span>),
     conf.space=c(<span style="color:#a82e2e">0.1</span>,<span style="color:#a82e2e">0.5</span>),
     label.every=<span style="color:#a82e2e">3</span>,
     col.grid = gray(c(<span style="color:#a82e2e">0.8</span>, <span style="color:#a82e2e">0.95</span>)),
     which=<span style="color:#00753b">"shock"</span>)</code></span>

运行上面的代码即可得到下面的列线图：

 

上面的代码中我们自行设定了部分参数，需要给大家写写都是什么意思，首先是lp.at=seq(-3,4,by=0.5)表示我们的自变量的线性组合只在-3到4的范围内以步长为0.5进行显示，不过奇怪的是作者设置了这个参数，又把lp设置为F（就是不显示预测变量的线性组合值），所以作者在上面的代码中设置这个完全是没有必要的，如果我们将lp参数改为T则可看到预测变量的线性组合结果，效果如下：

 

然后是fun这个参数，因为我们的结局是二分类的，这个分类是通过连接函数得到的，我们的预测变量的线性组合直接得到的其实是log odds（看不明白的同学请看之前关于逻辑回归机器学习的文章），我们这儿其实想要的是某一类的概率，所以就用了一个反logit转换，目的是为了得到每一类的概率p，具体，代码中我们将参数fun设定为function(x)1/(1+exp(-x))，其实写成fun=plogis也是可以的，这个plogis就是逻辑斯蒂的累计分布函数（Logistic Cumulative Distribution Function (plogis Function)）

为了更好地给大家说明，这儿给大家做个演示，运行如下代码大家就可以看到逻辑斯蒂的累计分布函数长啥样：

<span style="color:#222222"><code>x_plogis <- se<span style="color:#00753b">q(- 10, 10, by = 0.1)</span> 
y_plogis <- plogis(x_plogis) 
plot(y_plogis) </code></span>

 

看到没，纵轴就是我们想要的p，这个就是通过反logit转换得到p的原理，好好体会哈。

其实逻辑回归在做的就是把自变量的线性组合通过连接函数映射到上面这个图上（就是得到p）从而达到对2分类结局分类的目的，在做列线图的时候我们也得把这个东西以fun参数的形式给指出来，不然它就不会对自变量进行转换，出来的概率就是错误的，转换就是通过fun这个参数实现的。

再看fun.at参数，这个参数是设定函数轴的标签的，我们将其设定为fun.at=c(.001,.01,.05,seq(.1,.9,by=.1),.95,.99,.999)就是希望如果自变量的组合会导致这些值的话，就会在相应的地方打上标签，因为我们自变量的线性组合最小为-3其对应的p为0.0474，所以在途中标签是从0.05开始的；funlabel参数是用来规定轴的名字的。conf.int参数用来设定变量的confidence levels，默认不显示。

接下来看看plot函数中的几个参数：lplabel用来设定预测变量线性组合的标签，fun.side用来设定fun刻度值标签的上下位置，防止重叠；col.conf给置信水平上颜色；conf.space规定置信轴的垂直位置；label.every=3意思是轴的每3个值打一个标签；方便我们对齐的垂直轴的颜色可以通过col.grid调节。

有序结局变量的列线图画法

有序逻辑斯蒂回归也是用lrm()函数进行拟合的，代码如下：

<span style="color:#222222"><code>mod<span style="color:#a82e2e">.ord</span> <- lrm(Y ~ age+rcs(lac,4)*sex)</code></span>

上面的代码中拟合了加了限制性样条的sex和lac的交互项（之后会给大家写什么是样条，本篇文章略过）。

为了更好地给大家解释有序逻辑回归，我们先看看有序逻辑回归的模型表达：

我们在做有序逻辑回归的时候，比如我们的结局变量有4个水平，这个时候我们会同时拟合3个二分逻辑模型（下图中的1，2，3），我们始终是有一个参考水平的，解释的时候依然和二分逻辑回归一样的解释：

 

比如我们的例子中，我们就会同时拟合如下3个模型，第1个模型回答了个案是不是大于第1类，第二个模型回答了是不是大于第二类，第3个模型回答了是不是大于第3类，这样就穷尽了所有的可能性。

 

相应地，我们的概率p也得有三个：

<span style="color:#222222"><code>fun2 <- <span style="color:#8a7304"><span style="color:#114ba6">function</span>(x)</span> plogis(x-<span style="color:#114ba6">mod</span>.ord$coef[<span style="color:#a82e2e">1</span>]+<span style="color:#114ba6">mod</span>.ord$coef[<span style="color:#a82e2e">2</span>])
fun3 <- <span style="color:#8a7304"><span style="color:#114ba6">function</span>(x)</span> plogis(x-<span style="color:#114ba6">mod</span>.ord$coef[<span style="color:#a82e2e">1</span>]+<span style="color:#114ba6">mod</span>.ord$coef[<span style="color:#a82e2e">3</span>])</code></span>

以上就规定好了出图时的自变量线性组合的转换方法，接着就可以出图了，代码为：

<span style="color:#222222"><code>nom.ord <- nomogram(f, <span style="color:#8a7304"><span style="color:#114ba6">fun</span>=<span style="color:#a82e2e">list</span>(<span style="color:#00753b">'Prob Y>=1'</span>=plogis,
                                <span style="color:#00753b">'Prob Y>=2'</span>=fun2,
                                <span style="color:#00753b">'Prob Y=3'</span>=fun3)</span>,
                    lp=F,
                    <span style="color:#8a7304"><span style="color:#114ba6">fun</span>.at=<span style="color:#a82e2e">c</span>(<span style="color:#a82e2e">.01</span>,<span style="color:#a82e2e">.05</span>,seq(<span style="color:#a82e2e">.1</span>,<span style="color:#a82e2e">.9</span>,<span style="color:#114ba6">by</span>=<span style="color:#a82e2e">.1</span>)</span>,.<span style="color:#a82e2e">95</span>,.<span style="color:#a82e2e">99</span>))
plot(nom.ord, lmgp=.<span style="color:#a82e2e">2</span>, cex.axis=.<span style="color:#a82e2e">6</span>)</code></span>

因为我们的模型中是有交互的，是lac和sex交互，其中sex是有两个水平，那么我们出来的列线图就会在sex的两个水平分别画出lac，于是我们做出来的列线图就是下面这个样子的：

 

我们有看到，在列线图中结局Y不同顺序的概率都有出现，通过这个列线图就可以很方便地得到某个个案到底处于结局变量的哪个水平的概率最大了。

生存分析的列线图

对于生存数据我们可以用半参数模型和参数模型进行建模（见我之前的文章），在生存分析中我们一般会感兴趣特定时点的生存概率和中位生存期，当然啦，这些模型都是可以通过列线图进行可视化的。

首先我们需要拟合生存模型，用到的是psm函数，本例中我们用到的数据长这样：

 

需要用到的变量是(time,status，ph.ecog，sex和age，其中ph.ecog是病人日常活动功能的一个指标。我们现在就拟合一个模型来研究ph.ecog对病人生存的影响，将sex和age当作协变量，代码如下：

<span style="color:#222222"><code> <span style="color:#114ba6">mod</span>.sur <- psm(Surv(<span style="color:#114ba6">time</span>,<span style="color:#114ba6">status</span>) ~ ph.ecog+sex+age,
 data, dist=<span style="color:#00753b">'weibull'</span>)</code></span>

然后我们可以用下面的代码得到病人的中位生存期和生存概率：

<span style="color:#222222"><code>med <span style="color:#00753b"><- Quantile(mod.sur)</span>
surv <span style="color:#00753b"><- Survival(mod.sur)</span>
ddist <span style="color:#00753b"><- datadist(lung)</span></code></span>

到这儿，我们画图需要的东西就全部准备好了，接下来进行出图，首先是以生存期为结局的图：

<span style="color:#222222"><code>nom.sur1<-nomogram(mod.sur,
                   fun=<span style="color:#114ba6">list</span>(<span style="color:#8a7304"><span style="color:#114ba6">function</span>(x) <span style="color:#a82e2e">med</span>(lp=x, q=<span style="color:#a82e2e">0.5</span>),
                            <span style="color:#a82e2e">function</span>(x) <span style="color:#a82e2e">med</span>(lp=x,q=<span style="color:#a82e2e">0.25</span>)),
                   <span style="color:#a82e2e">funlabel</span>=<span style="color:#a82e2e">c</span>(<span style="color:#00753b">"Median Survival Time"</span>,
                              <span style="color:#00753b">"1Q Survival Time"</span>),
                   <span style="color:#a82e2e">lp</span>=<span style="color:#a82e2e">F</span>)
<span style="color:#a82e2e">plot</span>(nom.sur1,
     fun.side=list(c(rep(<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">7</span>),<span style="color:#a82e2e">3</span>,<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">3</span>,<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">3</span>),rep(<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">7</span>)),
     col.grid = c(<span style="color:#00753b">"red"</span>,<span style="color:#00753b">"green"</span>))</span></code></span>

运行代码就可以得到下面的图啦：

 

可以看到我们的结局变量有两，分别是中位生存期和第一四分位生存期，因为我们定义fun的时候是定义了两嘛。

接下来是以生存概率为结局进行出图，代码如下：

<span style="color:#222222"><code>nom.sur2 <- nomogram(<span style="color:#114ba6">mod</span>.sur, fun=list(<span style="color:#8a7304"><span style="color:#114ba6">function</span>(x)</span>surv(<span style="color:#a82e2e">200</span>, x),
  <span style="color:#8a7304"><span style="color:#114ba6">function</span>(x)</span> surv(<span style="color:#a82e2e">400</span>, x)),
  funlabel=c(<span style="color:#00753b">"200-Day Survival Probability"</span>,
             <span style="color:#00753b">"400-Day Survival Probability"</span>),
  lp=F)</code></span>

上面的代码分别定义了随访200天和400天时病人的生存概率，也就是我们的结局，接下来就用plot函数出图了，代码如下：

<span style="color:#222222"><code>plot(nom.sur2,
     <span style="color:#8a7304"><span style="color:#114ba6">fun</span>.side=<span style="color:#a82e2e">list</span>(c(rep(c(<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">3</span>)</span>,<span style="color:#a82e2e">5</span>),<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">1</span>),
                   c(<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">1</span>,rep(c(<span style="color:#a82e2e">3</span>,<span style="color:#a82e2e">1</span>),<span style="color:#a82e2e">6</span>))),
     xfrac=.<span style="color:#a82e2e">7</span>,
     col.grid = c(<span style="color:#00753b">"red"</span>,<span style="color:#00753b">"green"</span>))</code></span>

代码中的xfrac参数是用来规定轴标签和图距离的哈，运行代码就可以得到下面的列线图：

 

半参生存分析

半参生存分析其实就是比例风险模型，为什么叫半参呢，这儿再给大家补习一下：

A parametric survival model is one in which survival time (the outcome) is assumed to follow a known distribution. Examples of distributions that are commonly used for survival time are: the Weibull, the exponential (a special case of the Weibull), the log-logistic, the log-normal, etc.

The Cox proportional hazards model, by contrast, is not a fully parametric model. Rather it is a semi-parametric model because even if the regression parameters (the betas) are known, the distribution of the outcome remains unknown. The baseline survival (or hazard) function is not specified in a Cox model (we do not assume any shape or form).

我们依然是先拟合模型，以ph.ecog为i变量并控制sex和age，比例风险模型的拟合用的是cph()函数，本例中代码如下：

<span style="color:#222222"><code><span style="color:#114ba6">mod</span>.cox <- cph(Surv(<span style="color:#114ba6">time</span>,<span style="color:#114ba6">status</span>) ~ ph.ecog+sex+age,
               lung, surv=TRUE)</code></span>

接下来是出图前的准备（注意看fun参数）：

<span style="color:#222222"><code>ddist <- datadist(lung)
options(datadist=<span style="color:#00753b">'ddist'</span>)
surv.cox <- Survival(<span style="color:#114ba6">mod</span>.cox)
nom.cox <- nomogram(<span style="color:#114ba6">mod</span>.cox, fun=list(<span style="color:#8a7304"><span style="color:#114ba6">function</span>(x)</span>
  surv.cox(<span style="color:#a82e2e">200</span>, x),
  <span style="color:#8a7304"><span style="color:#114ba6">function</span>(x)</span> surv.cox(<span style="color:#a82e2e">400</span>, x)),
  funlabel=c(<span style="color:#00753b">"200-Day Sur. Prob."</span>,
             <span style="color:#00753b">"400-Day Sur. Prob."</span>),
  lp=F)</code></span>

然后出图：

<span style="color:#222222"><code>plot(nom.cox,
<span style="color:#8a7304"><span style="color:#114ba6">fun</span>.side=<span style="color:#a82e2e">list</span>(c(rep(c(<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">3</span>)</span>,<span style="color:#a82e2e">5</span>),<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">1</span>),
c(<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">1</span>,<span style="color:#a82e2e">1</span>,rep(c(<span style="color:#a82e2e">3</span>,<span style="color:#a82e2e">1</span>),<span style="color:#a82e2e">6</span>))))</code></span>

运行代码即可做出如下的列线图：

 

 

通过上图就可以很方便地得到某个病人200天和400天时的生存概率了，好了，以上就是论文的所有内容，就给大家解析到这儿，希望对做列线图的同学有用。

论文原文：Zhang, Z., & Kattan, M. W. (2017). Drawing Nomograms with R: applications to categorical outcome and survival data. Annals of translational medicine, 5(10).

小结

今天给大家写了临床常见预测模型的列线图的做法，感谢大家耐心看完，自己的文章都写的很细，代码都在原文中，希望大家都可以自己做一做，请关注后私信回复“数据链接”获取所有数据和本人收集的学习资料。如果对您有用请先收藏，再点赞转发。

也欢迎大家的意见和建议，大家想了解什么统计方法都可以在文章下留言，说不定我看见了就会给你写教程哦，另，咨询代做请私信。