Chapter 7 7.1 二分类Logistic回归

7.1.1 Logistic回归模型

rm(list = ls())
# 从0到1，每次间隔0.01，分成100份
p <- seq(from = 0, to = 1, by = .01)
# 优势定义为：p/(1 - p)
odds <- p/(1 - p)
# 绘制曲线图，具体参考https://blog.csdn.net/sinat_41805381/article/details/80164958
# type表示线型，col表示线的颜色
plot(log(odds), p, type = "l", col = "blue", ylab = "Probability", las = 1)
# abline常和plot函数一起使用，具体参考https://www.cnblogs.com/xudongliang/p/6755727.html
# h和v分别表示水平和竖直，lty表示线型
abline(h = 0.5, lty = "dashed")
abline(v = 0, lty = "dashed")

7.1.2 Logistic回归实例

# 1.数据准备
# data函数用法：https://blog.csdn.net/weixin_43369563/article/details/82986495
# str()函数具体用法：https://www.cnblogs.com/geeksongs/p/12426151.html
# str函数来查看这个数据框当中有哪些变量，以及变量的类型。可以简单看数据里包含的内容，不会全部显示出来。
# tab1函数用法：注意此处是1，不是英文字母L,'tab1' 是一种高级的单向制表，提供了一个很好的频率表和条形图
# mutate函数用法：增加新的变量并保持原有的变量。具体参考：https://blog.csdn.net/weixin_57198917/article/details/123149874
# factor函数：将数据转换为因子型。具体参考：https://blog.csdn.net/hsdcc217/article/details/78510087
library(MASS)
data(birthwt)
str(birthwt)
library(epiDisplay)
tab1(birthwt$ptl)
tab1(birthwt$ftv)
library(dplyr)
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)
         )
str(birthweight)
# 2.模型的建立
# glm函数用于拟合广义线性模型，通过给出线性预测器的符号描述和误差分布的描述来指定。官方文档http://127.0.0.1:27513/library/stats/html/glm.html
# 以low为结果变量，其他参数作为自变量，binomial表示Logistic回归
# summary函数获取描述性统计量，可以提供最小值、最大值、四分位数和数值型变量的均值，以及因子向量和逻辑型向量的频数统计等。
# pchisq进行卡方检验，结果是一个p值。默认参数是lower.tail是true，这里修改为false，表示计算大于39.19的概率。
glm1 <- glm(low ~ age + lwt + race + smoke + ptl + ht + ui + ftv,
            family = binomial, 
            data = birthweight)
summary(glm1)
pchisq(39.19, df = 188 - 178, lower.tail = FALSE)

# 3.自变量筛选
# drop1函数用法：剔除一个变量后模型的AIC值。获得到AIC值越小，表示这个自变量影响较小，可以剔除
# step函数表示自动完成drop的剔除过程。不断进行drop
# logLik()函数计算样本的对数似然值
# anova函数可以用于评估不同模型的效应。可参考：https://zhuanlan.zhihu.com/p/128088652
# AIC函数和也可以用于比较两个或者多个模型
drop1(glm1)
glm2 <- step(glm1, trace = FALSE)
summary(glm2)

# 4.模型比较
# logLik()函数计算样本的对数似然值
# anova函数可以用于评估不同模型的效应。可参考：https://zhuanlan.zhihu.com/p/128088652
# AIC函数和也可以用于比较两个或者多个模型
lnL1 <- as.numeric(logLik(glm1))
lnL2 <- as.numeric(logLik(glm2))
pchisq(2 * (lnL1 - lnL2), df = 3, lower.tail = FALSE)

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

# 5.回归系数的解释
# coef(glm2)查看回归系数，回归系数在回归方程中表示自变量对于因变量影响大小的参数。https://baike.baidu.com/item/%E5%9B%9E%E5%BD%92%E7%B3%BB%E6%95%B0/10840879
# exp(coef(glm2))将结果转化为优势比
# conflint得到回归系数的置信空间，置信水平95%。置信区间是指由样本统计量所构造的总体参数的估计区间。https://baike.baidu.com/item/%E7%BD%AE%E4%BF%A1%E5%8C%BA%E9%97%B4/7442583
# exp(confint(glm2))表示优势比的置信空间
coef(glm2)
exp(coef(glm2))
confint(glm2)
exp(confint(glm2))

# 6.预测
# 模型建立好之后，可以对新的样本输入进行预测
# predict函数进行预测，第一个参数是指模型，第二个是需要进行预测的数据
# 得到的logit值是优势的对数值，需要进行exp(P)/(1-exp(P))转化为概率,在predict中加入type = "response"可以直接得到概率值
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))
predict(glm2, newdata = newdata, type = "response")

# 7.模型的检查
# hoslem.test对观测值和预测值差异进行检验，p值小于0.05表示具有统计学意义
install.packages("ResourceSelection")
library(ResourceSelection)
hoslem.test(birthweight$low, fitted(glm2))

# 8.模型结果的汇总输出
# logistic.display汇总回归模型的主要结果
# write.csv将模型输出导入到一个表中
library(epiDisplay)
logistic.display(glm2)
write.csv(logistic.display(glm2)$table, file = "model.csv")

7.1.3 表格数据的Logistic回归

# 先录入表格数据，再将其转换为三维数组
# as.data.frame生成一个含有三个分类变量的全部组合下频数的***数据框***。
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")))
dat.table <- as.table(dat.array)
dat.table
dat <- as.data.frame(dat.table)
dat
# 方法1：使用glm函数建立模型，采用Freq作为频数
mod1 <- glm(outcome ~ smoke + drink, family = binomial, 
            weights = Freq, data = dat)
summary(mod1)

# 方法2：对数据框格式建立logistic回归模型
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)

# 方法3：将表格数据转化为原始数据后进行模型建立
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)

7.2 条件Logistic回归

# head用于查看数据的前几行
# clogit是survival库中的函数，用于建立条件Logistic回归模型。通过最大化条件似然来估计逻辑回归模型。具体描述：http://127.0.0.1:27513/library/survival/html/clogit.html
# 使用drop1对模型中自变量进行剔除
# 使用epiDisplay中的clogistic.display输出汇总条件Logistic回归模型的结果
library(epiDisplay)
data("VC1to1")
str(VC1to1)
head(VC1to1)
clogit1 <- clogit(case ~ smoking + alcohol + rubber + strata(matset), 
                  data = VC1to1)
clogit1
drop1(clogit1, test = "Chisq")
clogit2 <- clogit(case ~ alcohol + rubber + strata(matset), data = VC1to1)
drop1(clogit2, test = "Chisq")
clogit3 <- clogit(case ~ alcohol + strata(matset), data = VC1to1)
drop1(clogit3, test = "Chisq")
AIC(clogit1, clogit2, clogit3)
summary(clogit3)
clogistic.display(clogit3)
write.csv(clogistic.display(clogit3)$table, file = "clogit.csv")

7.3 无序多分类Logistic回归

# des函数用于描述数据框的数据类型和基本描述
# confint计算拟合模型中一个或多个参数的置信区间，默认为 2.5% 和 97.5%
# mlogit.display是epiDisplay中输出回归模型汇总结果的一个函数，与上面的display函数类似
# lrtest是“glm”类对象的似然比检验，
library(epiDisplay)
data(Ectopic)
des(Ectopic)
str(Ectopic)
multi1 <- multinom(outc ~ hia, data = Ectopic)
summary(multi1)
st <- summary(multi1)$standard.errors
z <- coef(multi1) / st; z
p.values <- pnorm(abs(z), lower.tail = FALSE) * 2
p.values
confint(multi1)
exp(coef(multi1))
exp(confint(multi1))
mlogit.display(multi1)
multi2 <- multinom(outc ~ hia + gravi, data = Ectopic)
mlogit.display(multi2)
lrtest(multi1, multi2)
AIC(multi1, multi2)

7.4 有序Logistic回归

# ordered默认按照字母顺序进行排序，后面可以添加条件设定排序类型
# polr1是MASS程序包中的函数，用于累计优势Logistic回归模型
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
# 方法1：使用原始数据建立模型
library(MASS)
polr1 <- polr(outcome ~ sex + method, weights = Freq, data = data1)
summary(polr1)
# 方法2：转化为逐条记录数据建立模型
data2 <- data1[rep(1:nrow(data1), data1$Freq), -4]
head(data2)
str(data2)
polr2 <- polr(outcome ~ sex + method, data = data2)
summary(polr2)
exp(coef(polr1))
exp(confint(polr1))
# VGAM包可以对所有类型的Logistic回归计算，并能够进行累积优势的Logistic回归模型的平行性驾驶检验
# VGAM::lrtest函数表示使用的是VGAM中的lrtest函数，::表示作用域
# ordinal.or.display汇总模型的主要结果
install.packages("VGAM")
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)
library(epiDisplay)
ordinal.or.display(polr1)
write.csv(ordinal.or.display(polr1), file = "polr.csv")