在上一篇廣義線性模型一（Generalized Linear Models洽糟，GLM） - 簡書 (jianshu.com)臭家，我們大致了解了glm的應(yīng)用范圍呵恢，并從三個方面探索模型構(gòu)建：
1种蘸、如何使用glm構(gòu)建logistics回歸
2墓赴、如何提取模型中的參數(shù)
3竞膳、不同模型之間如何比較
接下來，我們繼續(xù)從四個方面談一談logistics回歸：
1诫硕、模型可視化(列線圖)
2坦辟、使用構(gòu)建的模型預(yù)測結(jié)局
3、評價模型的效能
區(qū)分度(discrimination ability)：C指數(shù)章办、ROC曲線锉走、NRI、IDI
一致性（calibration ability）：校準曲線
4藕届、臨床效能分析（Decision Curve Analysis）：DCA, 生存曲線
以上四個方面并不能完全分開來講挪蹭，因此我將可以使用同一R代碼的部分一起講解，如下：列線圖的構(gòu)建和校準曲線休偶、DCA放在一起梁厉，預(yù)測結(jié)局變量及ROC放在一起。
一踏兜、列線圖及校準曲線
在上一篇文章中词顾，構(gòu)建logistic回歸我們使用了glm()函數(shù)，構(gòu)建方式如下:
fit.reduced <- glm(ynaffair ~ age + yearsmarried + religiousness + rating, data=Affairs, family=binomial())
當(dāng)我們要對其進行可視化時碱妆，就像我在cox回歸模型的展示中說的肉盹，想要使用rms繪制列線圖，就只能使用rms來構(gòu)建模型
survival包和rms包都生成原材料比如都養(yǎng)豬疹尾，同時rms還做香腸上忍，rms還只用自己家的豬。害，所以沒辦法绿渣，你家豬再優(yōu)質(zhì)因痛，想吃香腸就只能去找rms包
survival包學(xué)習(xí)筆記——cox回歸（一） - 簡書 (jianshu.com)
列線圖
繼續(xù)上一篇構(gòu)建的模型，我們來繪制列線圖
library(rms)
ddist <- datadist(Affairs)
options(datadist='ddist')
fit.reduced1 <- lrm(ynaffair ~ age + yearsmarried + religiousness +
rating, data=Affairs,x = T,y = T)
summary(fit.reduced1 )
nomo <- nomogram(fit.reduced1,
lp = F, #是否為線性
fun = plogis, #系數(shù)轉(zhuǎn)換的具體過程
fun.at = c(0.1,0.2,0.3,0.5,0.7,0.9), #最后一行的分段
funlabel = "離婚率" ) #最后一行的名字
plot(nomo)
這里我們每次使用只需要更改數(shù)據(jù)集的名字即可啦它抱！
校準曲線
一致性分析：校準曲線
首先繪制校準曲線，校準曲線所依賴的模型是上步rms::lrm()構(gòu)建的模型朴艰，因此繪制完列線圖可以直接繪制校準曲線观蓄。
library(rms)
ddist <- datadist(Affairs)
options(datadist='ddist')
fit.reduced1 <- lrm(ynaffair ~ age + yearsmarried + religiousness +
rating, data=Affairs,x = T,y = T)
summary(fit.reduced1 )
calb <- calibrate(fit.reduced1,
cmethod='hare', #logistic回歸使用的方法
method='boot', #使用抽樣的辦法，抽樣1000次
B=1000)
plot(calb,xlim=c(0,1.0),ylim=c(0,1.0))
二祠墅、模型預(yù)測及ROC曲線繪制
模型已經(jīng)擬合好了侮穿，那么他是否可以根據(jù)我們的變量來預(yù)測結(jié)局呢？
結(jié)局預(yù)測
那是當(dāng)然毁嗦，使用函數(shù)predict
S3 method for class 'glm'
predict(object, newdata = NULL, type = c("link", "response", "terms"), se.fit = FALSE, dispersion = NULL, terms = NULL, na.action = na.pass, ...)
S3 method for class 'lm'
predict(object, newdata, se.fit = FALSE, scale = NULL, df = Inf, interval = c("none", "confidence", "prediction"), level = 0.95, type = c("response", "terms"), terms = NULL, na.action = na.pass, pred.var = res.var/weights, weights = 1, ...)
參數(shù)介紹
fit: 這里fit就是剛剛擬合的函數(shù)亲茅，
newdata: 可以是多種多樣的
1、可以是原來的數(shù)據(jù)，預(yù)測后可以用來繪制ROC曲線
2克锣、可以是新的數(shù)據(jù)茵肃，來計算敏感性、特異性
3袭祟、可以是自己生成的數(shù)據(jù)验残，來展示某一變量的影響情況
type: 預(yù)測數(shù)據(jù)的類型
1、response巾乳，fit是線性模型時您没，返回值是規(guī)范化后的連續(xù)值; fit是二分類模型時，返回值是發(fā)生結(jié)局變量的可能性（probabilities）
2胆绊、term: 返回一個矩陣氨鹏，里面包含這個模型中的每個參數(shù)項的預(yù)測值，每一項之和或者其他運算得到返回值（formular中的y值）
實戰(zhàn)
這三個參數(shù)辑舷，已經(jīng)夠我們使用的了喻犁，接下來我們實戰(zhàn)看一下
library(rms)
library(xlsx)
library(tidyverse)
fit <- glm(HPD與否~.,data = data,family = binomial())
summary(fit)
a <- predict(fit,newdata = data,type = "terms",se.fit = T)
b <- predict(fit,newdata = data,type = "response")
計算C指數(shù)及ROC曲線繪制
區(qū)分度——C指數(shù)
rms包計算C指數(shù)
library(rms)
ddist <- datadist(Affairs) #將數(shù)據(jù)組裝
options(datadist='ddist')
fit.reduced1 <- lrm(ynaffair ~ age + yearsmarried + religiousness +
rating, data=Affairs,x = T,y = T)
#相對于基礎(chǔ)包中的glm，這里多了 x 和 y 參數(shù)
fit.reduced1
ROCR包繪制ROC曲線
這里使用ROCR必須依賴rms得到的fit何缓，同列線圖
library(rms)
ddist <- datadist(Affairs) #將數(shù)據(jù)組裝
options(datadist='ddist')
fit.reduced1 <- lrm(ynaffair ~ age + yearsmarried + religiousness +
rating, data=Affairs,x = T,y = T)
Affairs$predvalue<-predict(fit.reduced1)
#這里已經(jīng)在前文中包裝了相關(guān)的數(shù)據(jù)肢础，所以無需賦值newdata，type
library(ROCR)
pred <- prediction(Affairs$predvalue, Affairs$ynaffair)
#結(jié)局變量必須是二分類變量（0碌廓，1）传轰，否則結(jié)果會報錯
perf<- performance(pred,"tpr","fpr")
plot(perf)
abline(0,1)
auc <- performance(pred,"auc")
auc@y.values #auc即是C-statistics
ROCR得到perf后，使用plot(perf)繪制的ROC曲線較簡單谷婆，接下來我們使用pROC優(yōu)化ROC曲線
pROC包
輸入變量只需要predict得到的預(yù)測值和準確值
library(pROC)
pdf(file = 'roc.pdf',width = 5,height = 5)
plot.roc(Affairs$ynaffair~Affairs$predvalue, #輸入變量慨蛙，實際變量~預(yù)測變量
data = Affairs,
axes=T, ## 是否顯示xy軸
legacy.axes=T, ## TRUE為1-specificity,FLASE為specificity
main="ROC for HPD與否", ## Title
col= "black", ## 曲線顏色
lty=1, ## 曲線形狀
lwd=3, ## 曲線粗細
identity=T, ## 是否顯示對角線
identity.col="black", ## 對角線顏色
identity.lty=2, ## 對角線形狀
identity.lwd=2, ## 對角線粗細
print.thres=TRUE, ## 是否輸出cut-off值
print.thres.best.method="youden",
print.thres.pch=20, ## cut-off點的形狀
print.thres.col="black", ## cut-off點和文本的顏色
print.thres.cex=1.2, ## 文本放大的倍數(shù)（比起原始值）
#print.thres.pattern="%.3f", ## cut-off文本的格式
#print.thres.pattern.cex=1.2,
print.auc=T, ## 是否顯示AUC
ci=T, ## 是否顯示CI
print.auc.pattern="AUC = 0.704", ## 展示AUC的格式
print.auc.x=0.4, ## AUC值的X位置
print.auc.y=0.15, ## AUC值的Y位置
print.auc.cex=1.2, ## AUC值的放大倍數(shù)
print.auc.col='black', ## ACU值的顏色
auc.polygon=TRUE, ## 是否將ROC下面積上色
auc.polygon.col='white',
auc.polygon.density=NULL,
auc.polygon.lty=2,
auc.polygon.angle=45,
auc.polygon.border='black',
max.auc.polygon=TRUE,
max.auc.polygon.col='white',
max.auc.polygon.lty=2
)
dev.off()
最后
由于篇幅的問題，剩余的內(nèi)容繼續(xù)在下一篇中分享
