Logistic回归是最基本的统计分析和建模方法。以下为R代码实现和一些注意事项。代码来源见参考资料部分。
1. R代码实现
# install.packages("rms") # 安装包
library(rms) # 加载rms包
# 手动构建数据集
y <- c(0,1,0,0,0,1,1,1,0,0,0,1,1,0,0,1,0,0,0,1,1,0,1,
1,0,1,1,1,0,1,0,1,0,0,1,1,0,0,1,1,0,1,0,0,0,1,
1,1,1,0,1,1,0,0,0,1,1,1,0,1,1,1,1,0,0,1,1,1,0,
0,0,1,0,1,0,1,0,1)
age <- c(28,42,46,45,34,44,48,45,38,45,49,45,41,46,49,46,44,48,
52,48,45,50,53,57,46,52,54,57,47,52,55,59,50,54,57,60,
51,55,46,63,51,59,48,35,53,59,57,37,55,32,60,43,59,37,
30,47,60,38,34,48,32,38,36,49,33,42,38,58,35,43,39,59,
39,43,42,60,40,44)
sex <- c(0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,
0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,
0,1,0,1,0,1,1,1,0,1,1,1,0,1,1,1,0,1,1,1,0,1,1,1,
0,1,1,1,0,1)
ECG <- c(0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,
0,0,1,1,0,0,1,1,0,0,1,1,0,0,2,1,0,0,2,2,0,0,2,2,
0,1,2,2,0,1,0,2,0,1,0,2,1,1,0,2,1,1,0,2,1,1,0,2,
1,1,0,2,1,1)
dt <- data.frame(y,age,sex,ECG) # 把数据集设置成数据框结构
str(dt) ##查看每个变量结构
## 'data.frame': 78 obs. of 4 variables:
## $ y : num 0 1 0 0 0 1 1 1 0 0 ...
## $ age: num 28 42 46 45 34 44 48 45 38 45 ...
## $ sex: num 0 1 0 1 0 1 0 1 0 1 ...
## $ ECG: num 0 0 1 1 0 0 1 1 0 0 ...
Note: 这里应该把ECG处理为三分类变量。
head(dt) ##查看数据框前几行
## y age sex ECG
## 1 0 28 0 0
## 2 1 42 1 0
## 3 0 46 0 1
## 4 0 45 1 1
## 5 0 34 0 0
## 6 1 44 1 0
#设定环境参数#
ddist <- datadist(dt)
options(datadist='ddist')
###logistic
f <- lrm(dt$y ~ age + sex + ECG, data=dt, family = binomial()) ##注意此处使用lrm()函数构建二元LR
summary(f) ##也能用此函数看具体模型情况，模型的系数，置信区间等
## Effects Response : dt$y
##
## Factor Low High Diff. Effect S.E. Lower 0.95 Upper 0.95
## age 41.25 53.75 12.5 1.1606 0.43868 0.30086 2.0204
## Odds Ratio 41.25 53.75 12.5 3.1920 NA 1.35100 7.5416
## sex 0.00 1.00 1.0 1.3564 0.54645 0.28540 2.4274
## Odds Ratio 0.00 1.00 1.0 3.8823 NA 1.33030 11.3300
## ECG 0.00 2.00 2.0 1.7464 0.76867 0.23983 3.2530
## Odds Ratio 0.00 2.00 2.0 5.7339 NA 1.27100 25.8670
### nomogram
par(mgp=c(1.6,0.6,0),mar=c(2,2,2,2)) ## 设置画布
nomogram <- nomogram(f,fun=function(x)1/(1+exp(-x)), ##逻辑回归计算公式
fun.at = c(0.001,0.01,0.05,seq(0.1,0.9,by=0.1),0.95,0.99,0.999),#风险轴刻度
funlabel = "Risk of Death", #风险轴便签
lp=F, ##是否显示系数轴
conf.int = F, ##每个得分的置信度区间，用横线表示，横线越长置信度越
abbrev = F#是否用简称代表因子变量)
plot(nomogram)
image.png
##模型验证
##以原数据集为验证集
f.glm <- glm(y~.,data=dt,family = binomial(link = "logit"))
P1 <- predict(f.glm,type = 'response') ##获得预测概率值
##关键的一步来了
val.prob(P1,y) ##这个函数前面放概率值，后面放结局变量
image.png
## Dxy C (ROC) R2 D D:Chi-sq
## 5.675676e-01 7.837838e-01 3.164825e-01 2.578779e-01 2.111448e+01
## D:p U U:Chi-sq U:p Q
## NA -2.564103e-02 -3.694822e-13 1.000000e+00 2.835189e-01
## Brier Intercept Slope Emax E90
## 1.885480e-01 -4.335687e-09 9.999998e-01 1.157412e-01 6.085456e-02
## Eavg S:z S:p
## 3.492462e-02 -2.762507e-03 9.977958e-01
2. Notes and tips:
- 分类变量应该处理为factor，其中，二分类处理或者不处理成factor不影响结果，但是多分类变量必须处理成factor。
dtECG, levels = c(0,1,2),labels = c("L","M","H"))
- 一般将单因素分析中p < 0.05的变量纳入多因素回归分析中。也可以将纳入标准放宽，比如单因素p < 0.1的变量就纳入多因素模型中。
- 有时我们在做完t检验、卡方检验、秩和检验之后并不会立即构建多因素回归模型，而是先单后多，这样做的好处是：可以在单因素logistic回归里观察到在未校正其它变量的情况下，各个自变量的Crude OR, 也可以观察到多因素logistic回归里的adjustment OR
model <- glm(y ~age, data = dt, family = binomial())
summary(model)$coefficients
exp(cbind("OR"))=coef(model),confint(clin_model)))
- 多因素回归分析会识别并剔除可能的混杂偏倚，这就是为什么在单因素分析时显著的变量到了多因素里就不显著的原因
- 非常规先单后多：如果我们已经确定是研究自变量X与因变量y之间的关系，此时协变量如何筛选进入多因素模型呢？如果某个协变量Z与因变量y之间单因素分析的p值小于0.1，并且同时分析协变量Z与自变量X与因变量y的关系，但是由于Z的存在，使X的系数较单因素分析时变化超过10%，此时协变量Z应该纳入多因素分析中来。
uni_methods <- function(xvar){
model <- glm(y ~ age,data = dt,family = binomial())
coef <- coef(model)[2]
form <- as.formula(paste0("y ~ age +", xvar))
model2 <- glm(form,data = dt,family = binomial())
coef2 <- coef(model2)[2]
ratio <- abs(coef2-coef)/coef > 0.1
if(ratio){
return(xvar)
}
}
xvar <- c("age","sex","ECG")
xvar <- xvar[-which(xvar == "age")] #除去目标变量
xvar
lapply(xvar, uni_methods) #遍历其它变量
- 只要因变量为二分类变量，样本量有比较大，都可以考虑采用Logistic回归分析
- 如果是无序分类变量可将其化成为若干个哑变量；如果是有序变量则既可以按连续变量的方式处理，也可以按无序分类变量的方式进行处理。
- 在分类变量分组时，组数一般要<=5
参考资料：
- 微信公众号@易学统计R语言Logistic回归模型深度验证以及Nomogram绘制
- B站@大鹏统计工作室R语言logistic回归临床预测模型
- 军事医学统计学第1版，军事医学科学出版社
网友评论