1寫在前面

列線圖，又稱諾莫圖（Nomogram），是一種用于預(yù)測模型的可視化工具，它可以將多個(gè)影響因素和結(jié)局事件的關(guān)系展示在同一平面上。🥳

列線圖最早是由法國工程師Philbert Maurice d'Ocagne于1884年發(fā)明的，他在1880年提出了nomography的概念，用于快速計(jì)算復(fù)雜公式的圖形表示法，曾經(jīng)在工程和科學(xué)領(lǐng)域廣泛使用，直到被計(jì)算機(jī)取代。🥸

列線圖的基本原理是根據(jù)回歸模型的系數(shù)，給每個(gè)影響因素的不同取值賦分，然后將各個(gè)分?jǐn)?shù)相加得到總分，再根據(jù)總分和結(jié)局事件發(fā)生概率的轉(zhuǎn)換關(guān)系，計(jì)算出預(yù)測值。🤓

本期就盤點(diǎn)一下基于R語言的列線圖繪制方法。😜

2用到的包

rm(list = ls())
library(tidyverse)
library(survival)
library(rms)
library(nomogramFormula)
library(DynNom)

3示例數(shù)據(jù)

今天用一下偉大的cancer數(shù)據(jù)集中的colon。😘

data(cancer)

DT::datatable(colon)

4rms包制作列線圖

4.1 構(gòu)建模型

我們先構(gòu)建個(gè)模型吧。🥸

## 以下2步為必須步驟
ddist <- datadist(colon)
options(datadist='ddist')

fit1 <- rms::lrm(status~rx+sex+age+obstruct+perfor+nodes, 
           data = colon)

summary(fit1)

可視化一下吧。😜

## 繪圖
nomogram <- nomogram(fit1, fun = function(x)1/(1+exp(-x)))
plot(nomogram)

4.2 加入時(shí)間變量

我們再試著加入時(shí)間變量，需要用到psm和Surv。🤓

fit2 <- rms::psm(Surv(time,status) ~ rx + sex + age + obstruct + perfor + nodes, 
                 data=colon, dist='lognormal')
summary(fit2)

可視化一下吧。😘

med  <- Quantile(fit2)
surv <- Survival(fit2)

plot(nomogram(fit2, fun=function(x) med(lp=x), funlabel="Median Survival Time"))

再試著把解決變量改成半年和1年的生存率。🤩

nom <- nomogram(fit2, fun=list(function(x) surv(6, x),
                            function(x) surv(12, x)),
            funlabel=c("6-Month Survival Probability", 
                       "12-month Survival Probability"))
plot(nom, xfrac=.7)

5rms包的補(bǔ)充

rms的確很好用，但是很難計(jì)算出所有項(xiàng)目的總分和概率，這里可以用nomogramFormula包作為補(bǔ)充。😂

相關(guān)函數(shù)：👇

• formula_rd()，formula_lp() ?? polynomial regression;
• points_cal() ?? total points；
• prob_cal() ?? lrm(), cph() or psm() regression。

5.1 formula_lp

formula_lp(nomogram = nom)
formula_lp(nomogram = nom,power = 1)
formula_lp(nomogram = nom,power = 3,digits=6)

5.2 formula_rd

formula_rd(nomogram = nom)
formula_rd(nomogram = nom,power = 1)
formula_rd(nomogram = nom,power = 3,digits=6)

5.3 Calculate Total Points

##get the formula by the best power using formula_lp
results <- formula_lp(nom)
points_cal(formula = results$formula,lp=fit3$linear.predictors)

#get the formula by the best power using formula_rd
results <- formula_rd(nomogram = nom)
points_cal(formula = results$formula,rd=df)

5.4 Calculate Probabilities

# lrm() function
f <- lrm(status~rx+sex+age+obstruct+perfor+nodes, 
         data=colon,
         linear.predictors = T)

head(prob_cal(reg = f))

# cph() function
f <- cph(Surv(time,status)~rx+sex+age+obstruct+perfor+nodes, 
         data=colon,
         linear.predictors=T,
         surv=T)

head(prob_cal(reg = f,times = c(365,365*2)))

# psm() function
f <- psm(Surv(time,status)~rx+sex+age+obstruct+perfor+nodes, 
         data=colon)

head(prob_cal(reg = f,times = c(365,365*2)))

5.5 Caculate Total Points for nomogram Picture

TotalPoints.rms(rd = colon,fit = f,nom = nom)

6regplot包制作列線圖

6.1 構(gòu)建模型

glm函數(shù)來構(gòu)建哦。💪

fit3 <- glm(status~rx+sex+age+obstruct+perfor+nodes,
           data=colon,
           family="binomial")
summary(fit3)

可視化一下吧。😏

regplot(fit3, plots = c("density","boxes"), center = T,
        observation = F,points = T,
        dencol = "#EA5455",boxcol = "#002B5B",
        droplines=T
        )

6.2 加入時(shí)間變量

這次我們換成使用coxph函數(shù)。😎

Coxfit<-coxph(Surv(time,status) ~ rx + sex + age + obstruct + perfor + nodes, 
                 data=colon)

summary(Coxfit)

可視化一下吧。🥳

這里我們顯示一下第50個(gè)病例的180天、360天、1080天的生存概率。🧐

regplot(Coxfit, plots=c("violin","bars"), 
        observation = colon[50,],
        points = T,
        droplines=T,
        title="Survival Nomogram", 
        dencol = "#EA5455",boxcol = "#002B5B",
        prfail=T,  # For survival models only
        failtime=c(180,360,1080), # For survival models only
        #clickable=T
        )

7DynNom包制作列線圖

新鮮的shiny app，不過個(gè)人覺得不是特別好用。😔

# 需要glm函數(shù)
DynNom(fit3, colon)

點(diǎn)個(gè)在看吧各位~ ?.???? ??? ?