7.常用统计分析方法——Logistic回归

最新推荐文章于 2025-02-18 20:08:45 发布

Dr_long1996

最新推荐文章于 2025-02-18 20:08:45 发布

阅读量2.6w

点赞数 27

分类专栏： R语言临床统计分析文章标签： r语言

本文链接：https://blog.csdn.net/Dr_long1996/article/details/131090558

版权

R语言同时被 2 个专栏收录

16 篇文章

订阅专栏

临床统计分析

13 篇文章

订阅专栏

前一章介绍了多重线性回归模型，该模型中的因变量Y是定量变量，且给定自变量时，需要服从正态分布。如果需要分析的因变量为分类变量，如复发与未复发、生存和死亡、疗效、肿瘤组织的类型等，logistic回归就是分析该类因变量与自变量关系的方法。

Logistic回归模型

基础概念

logistic回归属于概率型非线性回归

与直线回归的区别：

线性回归的因变量y是连续性数值型变量，不能是分类变量。logistic回归是研究二分类、多分类、有序多分类（等级资料）观察结果（因变量）与一些影响因素（自变量）之间的关系。如因变量：食管癌发生，一个二分类变量；自变量：吸烟、饮酒、不良饮食习惯等危险因素

logistic回归模型的线性形式：ln(P/(1-P))=logitP=β0+β1X1+...+βnXn

基本原理：

如果只研究一个影响因素，可以使用前面介绍的 $\chi ^{2}$ 检验，如果研究多个因素影响，可以使用Mantel-Haenszel分层分析方法，用于校正多个混杂因素，但该方法存在一定的缺陷，无法描述混杂因素作用大小及方向，无法研究混杂因素之间的交互作用，其次该方法要求大样本，可分析因素少，混杂因素的增多，会发哦之分层变细，分析困难。

因变量Y为二分类变量，取值为0或1时，Y取值为1时通常为阳性结果（复发、死亡、有效等），也即是研究者关心的结果。

如果将影响因素与疾病发生概率 $\pi$ 直接建立回归模型，即 $\pi =\beta _{0}+\beta X$ ，则可能出现 $\pi$ 大于1或者小于0的情况。流行病学中将发病的概率 $\pi$ 与未发病的概率（1- $\pi$ ）之比称为优势（Odds）。Odds的取值范围是 $\left (0, +\infty \right )$ ，此时依然无法将取值在（ $-\infty ,+\infty$ ）的自变量的线性组合建立模型。因此对发病概率进行logit变换。

OR(Odd Ratio):优势比，比如x1表示性别，1表示男生，0表示女性，P1/（1-P1）表示男性发生与不发生的比值，P0/（1-P0）表示女性发生与不发生的优势

参数的意义：

1.β0（常数项）表示暴露剂量为0时个体发病与不发病概率之比的自然对数；

2.βi表示某个危险因素Xi增加一个单位时，结果Y优势比OR的对数值（其他危险因素固定）。

β=0，OR=1，X与Y之间无关
β>0，OR>1, X与Y有关，危险因素
β<0，OR<1, X与Y有关，保护因子

该值近似说明暴露与非暴露相比危险增加的倍数

但暴露因素增加一个单位时，增加前后的ln（Odds）相减，正好等于ln（OR），由此可以推导出该暴露因素的比值比OR=exp^βj。

因为在发生率较低的慢性病，比如心血管病、恶心肿瘤等，由于p很小，优势比可以作为相对危险度RR的近似估计值，RR值才是倍数关系，OR只是危险增加倍数
OR=（P1/（1-P1））/（P0/（1-P0））≈P1/P0=RR

logistic回归分析基本步骤：

1. logistic回归模型的建立
方法：最大似然估计
问题：样本量是否足够：一般按自变量个数的5-20倍，一般10倍。

2. 假设检验

2.1 总模型的假设检验

H0：只含有截距项的模型成立

H1：含有截距和m个自变量的模型成立

对数似然函数lnL的数值可以反映模型拟合的效果，对数似然函数值越大越好，或者对数似然函数的负二倍值越小越好。似然比检验（likelihood ratio test）统计量

$G=-2(lnL_{0}-lnL_{1})=-2ln(L_{0}/L_{1})$

(1)似然比检验引入的变量G=-2lnL-（-2lnL），G服从自由度为1的卡方，所以可以对G进行检验，评价新加入的自变量是否有意义
(2)wald检验
(3)计分检验

三种检验方法的比较：在大样本情况下，三种方法得到的结果是一致的
（1）似然比检验结果相对可靠，推荐使用
（2）在小样本情况下，score的分布可能更接近于卡方分布，犯一类错误的概率要小一些
（3）wald检验在计算和使用上更容易一些，到那时结果略偏于保守

2.2 回归系数的假设检验

H0：β1=β2=...=βm=0
H1：各βj不全为0

3. 拟合优度检验

拟合优度检验时用于检验所选模型与实际数据的吻合程度，评价模型预期值与实际观测值的一致性。

检验的零假设H0：实际观察的频数分布与模型预测的频数分布相符合，即是模型拟合观察资料，模型具有较好的拟合优度，可采用Hosmer-Lemeshow检验。拟合优度好，说明模型的预测能力强，适合作预测。

4. Logistic回归自变量的筛选

logistic回归对变量做筛选
筛选算法：前进法、后退法、逐步法
当比较不同自变量的贡献大小时，不能直接使用b进行比较，因为自变量单位不同，需要标准化
标准回归系数bj‘=bj*sj/（Π/根号3），比较自变量对因变量的贡献大小

5. logistic回归模型的评价

拟合模型预测效果：AUC较大，但拟合优度检验有显著性，说明数据中的信息还没有被模型充分解释，当前的模型还不是最好的。

拟合优度检验：拟合优度好，仅仅说明当前数据中的信息能被充分解释，但不能说明模型用于预测效果一定好

拟合优度：当前模型和饱和模型的差值是否具有统计学差异。

饱和模型：若模型中的参数的个数p+1=自变量组合数，则称为饱和模型。实际上，这样的模型必然时纳入自变量的所有主效应、各阶交互效应的模型

对于预测而言，饱和模型的样本符合率最高

6. logistic回归模型的诊断和修正

6.1 残差分析（异常点）
常用Standized残差和Deviance残差，如果残差绝对值>2,提示该条记录在多维空间中可能是异常点。

6.2 多重共线性
自变量之间的相关性，若存在多重共线性，计算自变量的偏回归系数B，矩阵不可逆，导致B的无解
后果：整个模型有统计学意义，但各自自变量偏回归系数没有统计学意义；自变量的偏回归系数取值大小或符号与实际相违背；专业上认为统计学差异的自变量检验结果却无意义；增加或者减少一个自变量或一条记录，自变量偏回归系数发生较大变化

6.3 共线性的诊断指标：
1.容忍度，容忍度越小，说明共线性越严重，若小于0.1表明存在严重的共线性。
2.方差膨胀因子（VIF）：容忍度的倒数，一般VIF<5或10
3.特征根：对常数项和所有自变量计算主成分。若自变量间存在较强的线性相关，则前面几个主成分将拥有较大的特征根，后面几个主成分较小，甚至接近0
4.条件指数：如果几个条件指数较大（如大于30），则提示存在多重共线性
上述四个指标综合评价共线性问题

6.4 多重共线性的处理：
（1）逐步回归：不适合严重共线性
（2）岭回归：有偏估计
（3）主成分回归：丢失一部分信息
（4）路径分析：清楚了解自变量之间的联系

二分类Logistic回归

案例：利用数据集birthwt探索新生儿低体重的影响因素

分析：

新生儿低体重是一个二分类变量，结局之比有两个“是”和“否”。根据数据纳入的影响因素包括age年龄，lwt母亲体重，race种族，smoke吸烟情况，ptl既往早产次数，ht高血压病史，ui子宫激惹状态，ftv怀孕初期探访医生的次数。

1. 数据准备

library(MASS)
data(birthwt)
str(birthwt)

查看分类变量的特征ptl和ftv。

library(epiDisplay)
tab1(birthwt$ptl)
tab1(birthwt$ftv)

根据数据分布特点，超过84.1%的患者没有早产史，可以将先前早产次数的变量转化为二分类变量：有和无。变量ftv可以转换为三个水平的因子。

birthweight <- birthwt %>% 
  mutate(race = factor(race, labels = c("white", "black", "other")),
         smoke = factor(smoke, labels = c("no", "yes")),
         ptl = ifelse(ptl > 0, "1+", ptl),
         ptl = factor(ptl),
         ht = factor(ht, labels = c("no", "yes")),
         ui = factor(ui, labels = c("no", "yes")),
         ftv = ifelse(ftv > 1, "2+", ftv),
         ftv = factor(ftv)
         )

设置因子化，并且利用条件函数巧妙的转化了分类变量。！！！！！好好学习

2. 模型建立

glm1 <- glm(low ~ age + lwt + race + smoke + ptl + ht + ui + ftv,
            family = binomial, 
            data = birthweight)
summary(glm1)

利用广义的回归模型函数glm，使用参数family设置广义线性模型的具体类型，binomial代表logistic回归，类似于线性模型。

第一部分显示了模型所调用的公式，第二部分给出了回归系数的估计、标准误和显著性检验的结果，第三个部分给出了两个离差平方和（类似于线性回归中的残差平方和），其中第一个空离差平方和（Null deviance）是指包含常数项的模型的离差平方和（也即是模型不含任何自变量，只会得到一个概率），第二个是当前模型的离差平方和（Residual deviance），第三个是AIC和模型拟合过程中迭代的次数，AIC可以用来比较模型，迭代次数一般不关注，当超过25，表示对于现有的数据而言，模型太过于复杂。我们关注的是两个离差平方和之差，234.67-195.48=39.19。这个值可以用来检验模型的显著性。离差的平方和近似服从 $\chi ^{2}$ 分布，所以可以计算出p值：

pchisq(39.19, df = 188 - 178, lower.tail = FALSE)

函数pchisq（）用于计算 $\chi ^{2}$ 分布的分布函数值，自由度df为两个模型自由度的差。参数lower.tail默认值为TRUE，这里设置为FALSE是为了计算大于39.19的概率。也即是说么包含自变量的模型比空模型好。还可以通过anova（）得到检验结果。

anova(glm1)

空自变量模型自由度为188，离差平方和等于234.67，全部包含的模型自由度df=178，离差平方和为195.48。与summary的结果一致。

3. 自变量的筛选

drop1(glm1)

根据drop1（）函数的首先剔除的变量为ftv，因为剔除后模型的AIC最小，为214.83。再次使用step（）函数自动完成。

glm2 <- step(glm1, trace = FALSE)
summary(glm2)

4. 模型的比较

似然比检验：

对数似然比LLR=2（lnL1-lnL2）

其中L1为复杂模型（包含较多自变量的模型）的似然值，L2为简单模型的似然值。LLR近似服从卡方分布，自由度为两个模型参数的个数之差。在比较两个模型时要求样本量相同。

函数logLik（）计算模型glm1和glm2的对数似然值

lnL1 <- as.numeric(logLik(glm1))
lnL2 <- as.numeric(logLik(glm2))
pchisq(2 * (lnL1 - lnL2), df = 3, lower.tail = FALSE)

[1] 0.4981082

自由度的计算，age包含一个参数，变量ftv包含两个参数（3个水平需要两个哑变量），所以模型的参数个数差为3.

anova(glm1, glm2, test = "Chisq")

参数test可以设为“LRT”，得到一样的结果。P=0.4981>0.05，说明含有6个自变量的模型与含有8个自变量的模型拟合得一样好。也即是，年龄和孕早期探访医生次数fvt这两个变量不会显著提高模型得预测精度。

函数AIC（）也可以比较两个或多个模型，一般情况下，AIC值越小模型拟合得越好。

AIC(glm1, glm2)

5. 回归系数的解释

coef(glm2)#查看回归系数
exp(coef(glm2))#转化成为优势比
confint(glm2)#回归系数得置信区间
exp(confint(glm2))#优势比的置信区间

结果表明，lwt母亲孕期体重每增加1lb，新生儿低体重得优势将乘以0.98（保持其他变量不变）。黑人孕妇分娩低体重儿的优势约为白种人的3.67倍。其他种族分娩低体重儿的优势约为白人孕妇的2.35倍。类似的，吸烟、有早产史、高血压和子宫应激症的孕妇分娩低体重儿的优势都有增加。因为自变量不可能全为0，所以常数项在这里没有特定的意义。

6. 预测

newdata <- data.frame(lwt = 120, race = "black", 
                      smoke = "yes", ptl = "0", 
                      ht = "yes", ui = "no")
logit <- predict(glm2, newdata = newdata)
logit
exp(logit)/(1 + exp(logit))

对于建立的孕妇样本，其分娩低体重儿的概率为88%，还可以直接通过设置predict（）函数的参数type=“response”直接得到预测值。

predict(glm2, newdata = newdata, type = "response")

7. 模型的检查

利用Hosmer-Lemeshow检验可以用于判断观测值和预测值之间的差异的显著性。使用ResourceSelection包中的函数hoslem.test（）可以执行该检验

library(ResourceSelection)
hoslem.test(birthweight$low, fitted(glm2))

零假设是观测值和预测概率之间的差异无统计学意义。不能拒绝关于模型很好拟合了数据的假设。

8. 结果汇总输出

library(epiDisplay)
logistic.display(glm2)
write.csv(logistic.display(glm2)$table, file = "model.csv")

结果包含了各个自变量在调整前后的优势比及其95%的置信区间，以及针对回归系数的wald检验的p值和针对自变量的似然比检验的p值。

无序多分类Logistic回归

如果因变量有k个类别，可以选定其中一个类别作为参考类别，分别拟合k-1个logistic模型。

案例：来自一个关于人工流产史是否是宫外孕的危险因素的病例对照研究。该数据一共包括723名患者，每名患者有三个观察指标，其中结果变量outc有3个水平：宫外孕患者EP、来做人工流产的孕妇（IA）和来分娩的孕妇（Deli），变量hia有2个水平，分别代表有无人工流产病史，变量gravi有3个水平，表示怀孕的次数。

library(epiDisplay)
data(Ectopic)
des(Ectopic)
str(Ectopic)

下面把分类变量outc的第一个水平“EP（宫外孕患者）”作为对照组，运用多分类Logistic回归分析人工流产史是否是宫外孕的危险因素。

nnet包中的函数multinom（）可以用于建立这种回归模型。

multi1 <- multinom(outc ~ hia, data = Ectopic)
summary(multi1)

从结果可以看到hia变量ever曾经有过人工流产史的妇女，宫内妊娠（并来做人工流产）的风险是未做人工流产的exp（-0.907）=0.404倍，宫外孕的风险增加了1/0.404-1=1.48倍；如果研究对象有人工流产病史，宫外孕的风险比正常分娩的孕妇增加了1/exp（-1.726）-1=4.617

由回归系数及其标准误可获得每个系数的z值：

st <- summary(multi1)$standard.errors
z <- coef(multi1) / st; z
p.values <- pnorm(abs(z), lower.tail = FALSE) * 2
p.values

confint(multi1)#回归系数的95%置信区间
exp(coef(multi1))#优势比
exp(confint(multi1))#优势比的95%置信区间

epiDisplay包中的函数mlogit.display()可以简化上述的命令：

mlogit.display(multi1)

这里的两个组的比较都是IA与宫外与比较，Deli与宫外孕比较，这里的RR值不是宫外孕的风险，而是其倒数。

multi2 <- multinom(outc ~ hia + gravi, data = Ectopic)
mlogit.display(multi2)

虽然gravi两个水平没有统计学差异，但是与之前得模型相比，离差更小，AIC更小。进行似然比检验函数lrtest（）

lrtest(multi1, multi2)
AIC(multi1, multi2)

P=0.001所以模型2是不同于模型1，且模型2的AIC更低。综上所述，调整了怀孕次数的影响后，人工流产史显著增加了宫外孕的风险。人工流产组宫外孕的相对危险度=1/0.33=3.03。

表格数据的Logistic回归

案例：一个研究吸烟、饮酒与食管癌关系的病例对照资料。

dat.array <- array(c(136, 57, 107, 151, 63, 44, 63, 265), 
             dim = c(2, 2, 2), 
             dimnames = list(smoke = c("no", "yes"),
                             drink = c("no", "yes"),
                             outcome = c("control", "case")))

array（）创建一个三维数组，行变量为吸烟smoke，列变量为drink，层变量为outcome。

dat.table <- as.table(dat.array)#as.table把数组转换成为一个表格
dat.table
dat <- as.data.frame(dat.table) #转化为数据框形式
dat

建立回归模型：

#长数据格式
mod1 <- glm(outcome ~ smoke + drink, family = binomial, 
            weights = Freq, data = dat)
summary(mod1)

#宽数据格式
library(tidyr)
dat.wide <- pivot_wider(dat, names_from = outcome, values_from = Freq)
dat.wide
mod2 <- glm(cbind(case, control) ~ smoke + drink, 
            family = binomial, data = dat.wide)
summary(mod2)

#原始数据格式
dat.original <- dat[rep(1:nrow(dat), dat$Freq), -4]
nrow(dat.original)
mod3 <- glm(outcome ~ smoke + drink, family = binomial, 
            data = dat.original)
summary(mod3)

my.data <- data.frame(outcom=c("control","control","case","case"),
                      smoke=c("yes","no","yes","no"),
                      freq=c(3,2,4,1))
my.data
my.data.original <- my.data[rep(1:nrow(my.data),my.data$freq),-3]#-3表示不要频数这一列
my.data.original

实现频数转化为原始数据的方法：

有序多分类Logistic回归

有序多分类的logit模型中最常用的是累积优势logistic回归模型。

思想：

累积logit模型假定各个自变量在各个累积模型中对累积概率的优势比影响相同，即各个logit模型的偏回归系数相同，差别在于常数项。也就是说，根据拟合的累积logit模型，绘制反应变量的累积概率与自变量所对应的曲线，则各条曲线是平行的，只是截距不同。故累积logistic回归模型需要进行平行性检验，也即是检验各个自变量在k-1个不同模型中的回归系数是否相同。

如果p>0.1，说明满足平行性假设，否则不满足平行性假设。不满足平行性假设的资料，不适合用累积logistic回归，可能Y的不同类别有不同的本质，应该采用无序多分类的logistic回归。

案例：研究人员随机选择84名患者进行临床试验，探讨新旧两种治疗方法对该病的疗效的影响。自变量为性别，和治疗方法X，因变量Y为疗效，无效，有效，痊愈三种等级。

分析：

结果变量有3个水平，按照有序将其分为两类：（1）{无效}，{有效，痊愈}；（2）{无效、有效}，{痊愈}。设 $\textit{P}_{k}=\left ( Y=k|X \right )$ ，则需要建立两个方程：

$logit(P_{1})=ln\frac{P_{1}}{1-P_{1}}=ln\frac{P_{1}}{P_{2}+P_{3}}=\alpha _{1}-\beta _{1}X_{1}-\beta _{2}X_{2}$

$logit(P_{1}+P_{2})=ln\frac{P_{1}+P_{2}}{1-(P_{1}+P_{2})}=ln\frac{P_{1}+P_{2}}{P_{3}}=\alpha _{2}-\beta _{1}X_{1}-\beta _{2}X_{2}$

按照累积优势模型的假定，无论对于哪种分法，治疗方法的效应都是相同的。因此每个方程的回归系数都是相同的。

dat <- array(c(10, 7, 19, 6, 0, 2, 7, 5, 1, 5, 6, 16), 
             dim = c(2, 2, 3), 
             dimnames = list(method = c("old", "new"),
                             sex = c("male", "female"),
                             outcome = c("effectless", "effective", "recover")))
dat <- as.table(dat)
data1 <- as.data.frame(dat)
data1
str(data1)

将结局指标定义为有序因子：

data1$outcome <- ordered(data1$outcome)
data1$outcome

MASS包中的函数polr（）可用于建立累积优势logistic回归模型。weight权重输入为freq

library(MASS)
polr1 <- polr(outcome ~ sex + method, weights = Freq, data = data1)
summary(polr1)
#或者原始记录格式
data2 <- data1[rep(1:nrow(data1), data1$Freq), -4]
head(data2)
str(data2)
polr2 <- polr(outcome ~ sex + method, data = data2)
summary(polr2)

由结果可以得到 $\alpha _{1}$ =1.8128， $\alpha _{2}$ =2.6672，β1=1.319，β2=1.797。由此可以得到两个回归方程。

$logit(P_{1})=ln\frac{P_{1}}{1-P_{1}}=ln\frac{P_{1}}{P_{2}+P_{3}}=1.8128-1.319X_{1}-1.797X_{2}$

$logit(P_{1}+P_{2})=ln\frac{P_{1}+P_{2}}{1-(P_{1}+P_{2})}=ln\frac{P_{1}+P_{2}}{P_{3}}=2.6672-1.319X_{1}-1.797X_{2}$

exp(coef(polr1))#计算OR值
exp(confint(polr1))#OR值得置信区间

本例中得OR是疗效中“无效对有效和痊愈”或“无效和有效对痊愈”的优势比。两个变量的OR值都大于1，且他们的95%置信区间都不包含1，说明女性患者疗效显著优于男性患者，新疗法的疗效优于旧疗法。

如果自变量较多，可以使用函数drop1（）和step（）筛选自变量。

平行性假设检验：

library(VGAM)
polr.p <- vglm(outcome ~ sex + method,
               family = cumulative(parallel = TRUE),
               data = data2)
polr.np <- vglm(outcome ~ sex + method,
                family = cumulative(parallel = FALSE),
                data = data2)
VGAM::lrtest(polr.p, polr.np)

函数vglm（）里面的参数parallel为TRUE和FALSE建立了两个模型，并采用似然比检验比较两个模型。P=0.4797说明没有统计学差异，可以认为两个模型满足平行性假设检验条件。为了避免与epiDisplay包中的lrtest混淆，使用VGAM::irtest()

输出保存结果：

library(epiDisplay)
ordinal.or.display(polr1)
write.csv(ordinal.or.display(polr1), file = "polr.csv")

条件logistic回归

对于匹配的病例对照研究，在招募一个病例时，可以根据患者指标如年龄、性别、或其他情况，选择一个或一组对照进行匹配，即1：1匹配或者1：m匹配。对于这一类资料选择条件logistic回归。

匹配资料进行比较时，比较是在每个匹配组内进行，而不是一个组与其他组相比较。

案例：探索吸烟、饮酒和在橡胶行业工作是否是食管癌的危险因素。每个病例匹配一例同性别同年龄的邻居，共26个匹配，52例观察对象。

library(epiDisplay)
data("VC1to1")
str(VC1to1)
head(VC1to1)

建立模型：

条件logistic回归可以使用survival包中的函数clogit（）建立模型

clogit1 <- clogit(case ~ smoking + alcohol + rubber + strata(matset), 
                  data = VC1to1)
clogit1

使用了三个变量，似然比检验显示该模型没有统计学差异，用drop1选择自变量。

drop1(clogit1, test = "Chisq")

smoking的p值最大，先从模型中剔除，重新建立模型。

clogit2 <- clogit(case ~ alcohol + rubber + strata(matset), data = VC1to1)
drop1(clogit2, test = "Chisq")

rubber也没有统计学差异，去掉后再次建立模型3.

clogit3 <- clogit(case ~ alcohol + strata(matset), data = VC1to1)
drop1(clogit3, test = "Chisq")
AIC(clogit1, clogit2, clogit3)

从AIC可以看出，模型3的拟合程度最好。

summary(clogit3)

得到最终模型的优势比，置信区间和p值

类似于前面的线性回归和非条件logistic回归，epiDisplay包中的函数clogistic.display（），用于汇总条件logistic回归模型的结果。

clogistic.display(clogit3)
write.csv(clogistic.display(clogit3)$table, file = "clogit.csv")