2.3 广义线性模型

前面我们介绍了单因素方差分析(One-way ANOVA)和基于最小二乘回归算法（OLS）的线性回归，这类分析方法可以通过一些列连续型或类别型解释变量来预测正太分布的响应变量（包括原始数据、转化后的数据近正太分布、或者残差近正太分布）。

然而在许多情下响应变量并非正太分布，例如下面的情况：

像上面这些情况使用最小二乘回归法就太理想，而广义线性模型（Generalized Linear Model，GLM）是线性回归模型的扩展，它可以处理更多类型的响应变量，包括二元、分类、计数、离散等非正态分布的数据。与传统的线性回归模型相比，广义线性模型引入了连接函数（Link Function）和分布族（Distribution Family）的概念，使得模型的应用范围更广泛。

要理解本章内容，我们首先需要简单回忆一下概率分布。 概率分布(probability distribution)是用来描述随机变量在各个取值上的概率的函数。随机变量是指其取值不确定或是不固定的变量。概率分布则是给出了这个随机变量在每个可能取值上出现的概率大小。常见的概率分布有离散型概率分布和连续型概率分布两种形式。离散型概率分布适用于随机变量只能取某些固定值的情况（计数型），连续型概率分布则适用于随机变量可以取到任意实数值的情况。

常见的一些概率分布包括：正态分布、二项分布、泊松分布、指数分布、伯努利分布、卡方分布等等。

回归模型的选择很大程度上都是取决于数据的概率分布类型。比如：

在本节中，我们将首先简要概述广义线性模型，并介绍如何使用glm()函数来进行估计。然后重点关注两种流行的模型：Logisic回归（响应变量为二元或类别型）和泊松回归（响应变量为计数型）。

2.3.1 GLM的glm()函数

R中可以通过glm()函数拟合广义线性模型。它与线性模型的lm()类似，只是多了一些参数。函数的基本形式为：

glm(formula, family=family(link=function), data)

  • formula：指定模型的公式，用于描述响应变量和预测变量之间的关系。
  • family：指定模型使用的概率分布和连接函数。
  • data：指定用于拟合模型的数据集。

2.3.2 Logistic回归

当通过解释变量来预测二元型（binary）结果变量时,Logistic回归是一个非常有用的工具。它用于对二元结果（即响应变量）进行建模，该响应变量只能有两个可能的值：0 或 1、是或否、患病或未患病。

该回归模型不直接返回观测值的类别。它使我们能够估计类别成员的概率 (p)。概率范围在 0 到 1 之间。你需要确定类别从一个类别翻转到另一个类别的阈值概率。默认情况下设置为p = 0.5，但实际情况应根据分析目的来确定。

下面用实例来介绍。

2.3.2.1 安装和加载R包

安装包

install.packages("mlbench")
install.packages("caret")

加载包

library(mlbench)
library(tidyverse)
library(caret)
2.3.2.2 准备数据

加载数据

data("PimaIndiansDiabetes2")

删除缺失值NA

PimaIndiansDiabetes2 <- na.omit(PimaIndiansDiabetes2)
  • na.omit(): 这是一个基本的R函数，用于删除数据框中包含缺失值（NA值）的观测行。它会返回一个新的不包含缺失值的数据框。

检查数据

sample_n(PimaIndiansDiabetes2, 3)
##     pregnant glucose pressure triceps insulin mass pedigree age diabetes
## 403        5     136       84      41      88 35.0    0.286  35      pos
## 761        2      88       58      26      16 28.4    0.766  22      neg
## 477        2     105       80      45     191 33.7    0.711  29      pos
  • sample_n(): 这是dplyr包中的函数，用于从数据框或数据集中进行随机抽样。它可以按照指定的数量来随机抽取观测。

将数据分为训练集和测试集

set.seed(123)

training.index <- PimaIndiansDiabetes2$diabetes %>% createDataPartition(p = 0.8, list = FALSE)

train.data  <- PimaIndiansDiabetes2[training.index, ]

test.data <- PimaIndiansDiabetes2[-training.index, ]

  • set.seed(123): 这是设置随机数种子的函数，通过给定一个种子值（这里是123），它确保每次运行代码时生成的随机数都是相同的。这样做可以使结果可重复，方便调试和比较。

  • training.index 这行代码使用了 createDataPartition() 函数来创建一个名为training.index的训练集的索引，目标变量是diabetes列。p = 0.8表示将80%的观测分配给训练集，list = FALSE表示返回的索引结果是一个向量而不是列表。

  • train.data 这行代码根据训练集的索引从原始数据框PimaIndiansDiabetes2中选择对应的观测行，将其作为训练集数据。

  • test.data 这行代码使用负号 - 将training.index的相反值应用于原始数据框PimaIndiansDiabetes2，从而选择除了训练集以外的观测行作为测试集数据。

后面我们会使用train.data和test.data这两份数据进行模型拟合的训练和测试。


2.3.2.3 拟合logistic回归

用于广义线性模型的 R 函数 glm() 可用于计算逻辑回归。需要指定分布族选项 family = binomial，它告诉 R 我们想要拟合logistic回归。

使用训练集拟合logistic回归模型

model <- glm( diabetes ~., data = train.data, family = binomial)
  • diabetes：本数据集的响应变量，表示糖尿病的发生与否。
  • ~ .表示我们将使用所有其他变量作为自变量（也称为预测变量）来预测diabetes。点号.表示“所有其他变量”。

模型结果总结

summary(model)
## 
## Call:
## glm(formula = diabetes ~ ., family = binomial, data = train.data)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.053e+01  1.440e+00  -7.317 2.54e-13 ***
## pregnant     1.005e-01  6.127e-02   1.640  0.10092    
## glucose      3.710e-02  6.486e-03   5.719 1.07e-08 ***
## pressure    -3.876e-04  1.383e-02  -0.028  0.97764    
## triceps      1.418e-02  1.998e-02   0.710  0.47800    
## insulin      5.940e-04  1.508e-03   0.394  0.69371    
## mass         7.997e-02  3.180e-02   2.515  0.01190 *  
## pedigree     1.329e+00  4.823e-01   2.756  0.00585 ** 
## age          2.718e-02  2.020e-02   1.346  0.17840    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 398.80  on 313  degrees of freedom
## Residual deviance: 267.18  on 305  degrees of freedom
## AIC: 285.18
## 
## Number of Fisher Scoring iterations: 5

使用测试集进行对模型进行预测

probabilities <- model %>% predict(test.data, type = "response")
predicted.classes <- ifelse(probabilities > 0.5, "pos", "neg")
  • probabilities 这行代码使用训练好的logistic回归模型对测试集 test.data 数据进行预测。
  • predicted.classes 这行代码将预测的概率值转换成类别标签。如果概率值大于0.5，则将其标记为”pos”，否则标记为”neg”。

结果解读：

  • Call: 这一行显示了glm()函数的调用参数，指定了模型的公式、使用的家族（binomial）和数据集等信息。

  • Coefficients: 这个表格显示了每个解释变量的估计系数（Estimate）、标准误差（Std. Error）、z值（z value）、以及对应的P值（Pr(>|z|)）。

  • Signif. codes: 这一列展示了统计学上显著性水平的符号代码。

  • Null deviance和Residual deviance: 这两个值是用来评估模型拟合优度的统计量。Null deviance: 这是一个基准模型的偏差，通常是只包含截距项的模型。Residual deviance表示在当前模型中的偏差。通过比较模型的Residual deviance和Null deviance，我们可以评估模型相对于简单基准模型的改进程度。在这个例子中，模型的Residual deviance相对于Null deviance有显著的下降，表明拟合的模型相对于基准模型有了显著的改进。

  • AIC: AIC（赤池信息准则）是一种模型选择的准则，它综合考虑了拟合优度和参数数量。较小的AIC值表示模型更好。AIC值的大小是相对的，一般是在比较多个拟合模型时，越小代表拟合越好。

  • Number of Fisher Scoring iterations: 这个数字表示在拟合模型时使用Fisher Scoring方法进行了多少次迭代。

通过这些结果，我们可以了解到每个解释变量对糖尿病的影响程度以及它们的显著性。例如，glucose变量的估计系数为0.0371，P值非常显著（小于0.001），说明血糖水平对糖尿病的预测有重要作用。而pressure变量的估计系数为-0.0003876，P值大于0.1，说明血压对糖尿病的预测作用不显著。

模型精度

mean(predicted.classes == test.data$diabetes)
## [1] 0.7564103

  • 结果表示测试集预测准确率为75.64%。

  • 这个结果告诉我们，在使用上面拟合的logistic回归模型对测试集进行预测时，有大约75.64%的样本被正确地分类。换句话说，该模型在预测糖尿病患者和非糖尿病患者时，预测正确的比例为75.64%。


选择“最佳模型”

请注意，线性回归的许多概念也适用于逻辑回归建模。例如，您需要执行一些回归诊断以确保你的数据满足模型所做的假设。

在有许多解释变量的情况下，你可以在不影响预测准确性的情况下，使用逐步回归和全子集回归来选择“最佳”拟合模型（同线性模型的章节）。

例如执行逐步回归：

library(MASS)
model <- glm(diabetes ~., data = train.data, family = binomial)
stepAIC(model, direction = "backward")
## Start:  AIC=285.18
## diabetes ~ pregnant + glucose + pressure + triceps + insulin + 
##     mass + pedigree + age
## 
##            Df Deviance    AIC
## - pressure  1   267.18 283.18
## - insulin   1   267.34 283.34
## - triceps   1   267.69 283.69
## - age       1   269.03 285.03
## <none>          267.18 285.18
## - pregnant  1   269.92 285.92
## - mass      1   273.82 289.82
## - pedigree  1   275.30 291.30
## - glucose   1   305.60 321.60
## 
## Step:  AIC=283.18
## diabetes ~ pregnant + glucose + triceps + insulin + mass + pedigree + 
##     age
## 
##            Df Deviance    AIC
## - insulin   1   267.34 281.34
## - triceps   1   267.69 281.69
## - age       1   269.13 283.13
## <none>          267.18 283.18
## - pregnant  1   269.92 283.92
## - mass      1   274.09 288.09
## - pedigree  1   275.34 289.34
## - glucose   1   305.95 319.95
## 
## Step:  AIC=281.34
## diabetes ~ pregnant + glucose + triceps + mass + pedigree + age
## 
##            Df Deviance    AIC
## - triceps   1   267.79 279.79
## <none>          267.34 281.34
## - age       1   269.40 281.40
## - pregnant  1   270.01 282.01
## - mass      1   274.91 286.91
## - pedigree  1   275.57 287.57
## - glucose   1   324.26 336.26
## 
## Step:  AIC=279.79
## diabetes ~ pregnant + glucose + mass + pedigree + age
## 
##            Df Deviance    AIC
## <none>          267.79 279.79
## - age       1   270.21 280.21
## - pregnant  1   270.47 280.47
## - pedigree  1   276.48 286.48
## - mass      1   285.17 295.17
## - glucose   1   324.60 334.60
## 
## Call:  glm(formula = diabetes ~ pregnant + glucose + mass + pedigree + 
##     age, family = binomial, data = train.data)
## 
## Coefficients:
## (Intercept)     pregnant      glucose         mass     pedigree          age  
##   -10.77109      0.09903      0.03814      0.09512      1.36220      0.02956  
## 
## Degrees of Freedom: 313 Total (i.e. Null);  308 Residual
## Null Deviance:       398.8 
## Residual Deviance: 267.8     AIC: 279.8

结果显示：formula = diabetes ~ pregnant + glucose + mass + pedigree + age 表达式拟合的模型”最佳”。

此外，您可以在模型中添加交互项，或包含样条项（参照上节中的有交互项的多元线性回归）。

2.3.2.4 广义加性模型

当无法假设预测变量的线性时，你还可以拟合广义加性模型（Generalized Additive Model, GAM）（参照上节中的多项式回归）。

这可以使用以下包来完成mgcv：

install.packages("mgcv")
library(mgcv)
# Fit the model
gam.model <- gam(diabetes ~ s(glucose) + mass + pedigree,
                 data = train.data, family = "binomial")
  • gam()：这是一个函数，用于拟合 GAM 模型。它接受多个参数来指定模型的结构和数据的来源。
  • s(glucose)：这是一个平滑项（smooth term），使用 s() 函数来对变量 glucose 进行光滑处理。这意味着 glucose 的效果不是线性的，而是通过平滑函数来建模。
  • mass 和 pregnant：这是线性项（linear term），指定了两个自变量 mass 和 pregnant。
  • data：这是指定数据集的参数，这里是 train.data，表示训练数据集。
  • family：这是指定模型的误差分布和链接函数的参数，这里是 “binomial”，表示二项分布并使用逻辑斯蒂链接函数。

Summarize model

summary(gam.model )
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## diabetes ~ s(glucose) + mass + pedigree
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -4.36318    0.84957  -5.136 2.81e-07 ***
## mass         0.08282    0.02299   3.602 0.000315 ***
## pedigree     1.27657    0.46337   2.755 0.005870 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##            edf Ref.df Chi.sq p-value    
## s(glucose)   1      1  58.18  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.338   Deviance explained = 29.1%
## UBRE = -0.074402  Scale est. = 1         n = 314

Make predictions

probabilities <- gam.model %>% predict(test.data, type = "response")
predicted.classes <- ifelse(probabilities> 0.5, "pos", "neg")

Model Accuracy

mean(predicted.classes == test.data$diabetes)
## [1] 0.7692308

2.3.3 泊松回归

当通过解释变量预测计数型结果变量时，泊松回归是一个非常有用的工具。

我们将使用robust包中的Breslow癫痫病数据来介绍泊松回归。

数据描述：患有简单或复杂部分性癫痫的患者被随机分配接受抗癫痫药物普罗加比特或安慰剂治疗。在四次连续的随机后门诊访问中，报告了前两周内发生的癫痫发作次数。

2.3.3.1 准备数据

安装包

install.packages("robust")

加载包

library(robust)

加载数据

data("breslow.dat")
names(breslow.dat)
##  [1] "ID"    "Y1"    "Y2"    "Y3"    "Y4"    "Base"  "Age"   "Trt"   "Ysum" 
## [10] "sumY"  "Age10" "Base4"

  • ID：一个整数值，指定患者的识别号码。

  • Y1：一个整数值，表示第一个两周期间癫痫发作的次数。

  • Y2：一个整数值，表示第二个两周期间癫痫发作的次数。

  • Y3：一个整数值，表示第三个两周期间癫痫发作的次数。

  • Y4：一个整数值，表示第四个两周期间癫痫发作的次数。

  • Base：一个整数值，给出八周基线癫痫发作的次数。

  • Age：一个整数值，给出患者的年龄。

  • Trt：治疗方式，一个因子，包括安慰剂和普罗加比特两个水平。

  • Ysum：一个整数值，表示Y1、Y2、Y3和Y4的总和。

  • sumY：一个整数值，表示Y1、Y2、Y3和Y4的总和。

  • Age10：一个数值，表示年龄除以10得到的结果。

  • Base4：一个数值，表示Base除以4得到的结果。

数据摘要

summary(breslow.dat[c(6,7,8,10)])
##       Base             Age               Trt          sumY       
##  Min.   :  6.00   Min.   :18.00   placebo  :28   Min.   :  0.00  
##  1st Qu.: 12.00   1st Qu.:23.00   progabide:31   1st Qu.: 11.50  
##  Median : 22.00   Median :28.00                  Median : 16.00  
##  Mean   : 31.22   Mean   :28.34                  Mean   : 33.05  
##  3rd Qu.: 41.00   3rd Qu.:32.00                  3rd Qu.: 36.00  
##  Max.   :151.00   Max.   :42.00                  Max.   :302.00

注意：虽然数据集中有12个变量，但我们暂时只关注治3个解释变量（疗条件（Trt）、年龄（Age）和前八周的基础发病数（Base））对响应变量sumY（随机化后八周内的癫痫发病率）。

数据可视化

par(mfrow=c(1,2))
attach(breslow.dat)
hist(sumY, breaks = 50, xlab  ="Seizure Count", main = "Distribution of Seizures") 
#Seizure Count癫痫发作计数
boxplot(sumY ~ Trt, xlab ="Treatment", main ="Group Comparisions")

  • Placebo和Progabide是Trt变量中两个药物的名称。

图中可以清楚的看到响应变量的偏倚分布特性和可能的离群点，注意排查异常值的原因，方法我们不再次复述。

2.3.3.3 拟合泊松回归
model.ps <- glm( sumY ~ Base + Age + Trt, data = breslow.dat, family = poisson)
summary(model.ps)
## 
## Call:
## glm(formula = sumY ~ Base + Age + Trt, family = poisson, data = breslow.dat)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.9488259  0.1356191  14.370  < 2e-16 ***
## Base          0.0226517  0.0005093  44.476  < 2e-16 ***
## Age           0.0227401  0.0040240   5.651 1.59e-08 ***
## Trtprogabide -0.1527009  0.0478051  -3.194   0.0014 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2122.73  on 58  degrees of freedom
## Residual deviance:  559.44  on 55  degrees of freedom
## AIC: 850.71
## 
## Number of Fisher Scoring iterations: 5

结果解读方法同上（logistic回归）

同样的，拟合模型的表达式可以使用部分解释变量或全部解释变量，亦或考虑交互效应。 如何找到适合你数据的“最佳”模型，参考上节的线性回归。