R语言基础| 广义线性模型

12.1 广义线性模型(Generalized Linear Model)和glm()函数

前面学到的OLS回归和方差分析仅适用于因变量符合正态分布的数据，广义线性模型则可以用于因变量不合符合正态分布的数据：

1）当因变量为类别型，例如二值变量（是/否）和多分类变量（差/良/优）

2）因变量为计数型（如一周内交通事故数、每日酒水消耗的数量）等，这类变量都是非负的有限值，且均值和方差不独立，具有相关性。

一般使用glm()函数，重点为logistic回归（因变量为类别型）和泊松回归（因变量为计数型）。

12.1.1 glm()函数

语法：

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

family为概率分布，分别有以下几种，以及相应默认的连接函数（function):

与lm()函数类似，glm()函数也有一些连用函数，用来查看一些glm()结果的：

12.2 Logistic回归(family=binomial())    

当通过一系列连续性或/和类别型变量来预测二值型因变量时，使用logistic回归。

12.2.1 进行logistic回归时的步骤

1.确保因变量为二值型因变量(yes/no)

2.需要先将因变量转化为二值型因子变量

3.使用glm(因变量～自变量1+自变量2+自变量3+自变量4，data=,family=binomial())

12.2.2举例

用AER包中的Affairs数据集来举例。该数据从601个变量身上收集了9个变量，包括一年来出现感情危机的频率（出轨频率）、参与者性别、年龄、婚龄、是否有子女、宗教信仰程度（1-5，数字越大越信仰）、学历、职业、对婚姻的自我评分（1-5，数值越大，越感觉幸福）。

install.packages("AER")

library(AER)

## Warning: 程辑包'AER'是用R版本4.3.2 来建造的

## 载入需要的程辑包：lmtest

## 载入需要的程辑包：zoo

## 
## 载入程辑包：'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## 载入需要的程辑包：sandwich

## Warning: 程辑包'sandwich'是用R版本4.3.2 来建造的

data(Affairs,package = "AER")

summary(Affairs)

##     affairs          gender         age         yearsmarried    children 
##  Min.   : 0.000   female:315   Min.   :17.50   Min.   : 0.125   no :171  
##  1st Qu.: 0.000   male  :286   1st Qu.:27.00   1st Qu.: 4.000   yes:430  
##  Median : 0.000                Median :32.00   Median : 7.000            
##  Mean   : 1.456                Mean   :32.49   Mean   : 8.178            
##  3rd Qu.: 0.000                3rd Qu.:37.00   3rd Qu.:15.000            
##  Max.   :12.000                Max.   :57.00   Max.   :15.000            
##  religiousness     education       occupation        rating     
##  Min.   :1.000   Min.   : 9.00   Min.   :1.000   Min.   :1.000  
##  1st Qu.:2.000   1st Qu.:14.00   1st Qu.:3.000   1st Qu.:3.000  
##  Median :3.000   Median :16.00   Median :5.000   Median :4.000  
##  Mean   :3.116   Mean   :16.17   Mean   :4.195   Mean   :3.932  
##  3rd Qu.:4.000   3rd Qu.:18.00   3rd Qu.:6.000   3rd Qu.:5.000  
##  Max.   :5.000   Max.   :20.00   Max.   :7.000   Max.   :5.000

现研究以上这些自变量对是否出轨（有过感情危机/没有过感情危机）这一二值型因变量的影响：

1.首先这一因变量满足二值型变量

2.将该二值型变量设置为二值型因子变量。这里的affairs的数值为0（无出轨）和大于0（有出轨），需要手动添加这一二值型变量，再将其变为因子变量：

Affairs$ynaffairs[Affairs$affairs>0]<-1
Affairs$ynaffairs[Affairs$affairs==0]<-0
Affairs$ynaffairs <- factor(Affairs$ynaffairs,
                            levels = c(0,1),
                            labels = c("No","yes"))
table(Affairs$ynaffairs)

## 
##  No yes 
## 451 150

3.代入glm()函数

由于变量较多，怕写错或漏了，可以使用attach()函数调出变量名

attach(Affairs)

## The following object is masked _by_ .GlobalEnv:
## 
##     age

fit <- glm(ynaffairs~age+children+education+gender+occupation+rating+religiousness+yearsmarried,data=Affairs,family = binomial())
summary(fit)

## 
## Call:
## glm(formula = ynaffairs ~ age + children + education + gender + 
##     occupation + rating + religiousness + yearsmarried, family = binomial(), 
##     data = Affairs)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    1.37726    0.88776   1.551 0.120807    
## age           -0.04426    0.01825  -2.425 0.015301 *  
## childrenyes    0.39767    0.29151   1.364 0.172508    
## education      0.02105    0.05051   0.417 0.676851    
## gendermale     0.28029    0.23909   1.172 0.241083    
## occupation     0.03092    0.07178   0.431 0.666630    
## rating        -0.46845    0.09091  -5.153 2.56e-07 ***
## religiousness -0.32472    0.08975  -3.618 0.000297 ***
## yearsmarried   0.09477    0.03221   2.942 0.003262 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 675.38  on 600  degrees of freedom
## Residual deviance: 609.51  on 592  degrees of freedom
## AIC: 627.51
## 
## Number of Fisher Scoring iterations: 4

结果显示：性别、是否有子女、学历和职业对方程的贡献不显著（P>0.05）。去除这些重新拟合看看是否会拟合的更好：

fit1 <- glm(ynaffairs~age+rating+religiousness+yearsmarried,data=Affairs,family = binomial())
summary(fit1)

## 
## Call:
## glm(formula = ynaffairs ~ age + rating + religiousness + yearsmarried, 
##     family = binomial(), data = Affairs)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    1.93083    0.61032   3.164 0.001558 ** 
## age           -0.03527    0.01736  -2.032 0.042127 *  
## rating        -0.46136    0.08884  -5.193 2.06e-07 ***
## religiousness -0.32902    0.08945  -3.678 0.000235 ***
## yearsmarried   0.10062    0.02921   3.445 0.000571 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 675.38  on 600  degrees of freedom
## Residual deviance: 615.36  on 596  degrees of freedom
## AIC: 625.36
## 
## Number of Fisher Scoring iterations: 4

用anova()函数对两套模型进行比较，对于广义线性回归使用卡方检验：

anova(fit,fit1,test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: ynaffairs ~ age + children + education + gender + occupation + 
##     rating + religiousness + yearsmarried
## Model 2: ynaffairs ~ age + rating + religiousness + yearsmarried
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       592     609.51                     
## 2       596     615.36 -4  -5.8474   0.2108

结果显示Pr(chi)为0.21，故无法拒绝原假设，两套模型无显著差异。

补充：在卡方检验中，原假设（null hypothesis）通常是假设两个或多个分类变量之间不存在相关性或差异。具体来说，原假设可以表示为“变量之间独立”或“观察到的频数与期望的频数之间没有显著差异”。

12.2.3 解释模型参数

在logistic回归中，因变量是Y=1的对数优势比(log)。回归系数的含义是当其他自变量不变时，一单位自变量的变化可引起的因变量对数优势比的变化。

回归系数：

coef(fit)

##   (Intercept)           age   childrenyes     education    gendermale 
##    1.37725816   -0.04425502    0.39767213    0.02105086    0.28028665 
##    occupation        rating religiousness  yearsmarried 
##    0.03091971   -0.46845426   -0.32472063    0.09477302

由于对数优势比解释性差，我们可对结果进行指数化：

exp(coef(fit))

##   (Intercept)           age   childrenyes     education    gendermale 
##     3.9640180     0.9567099     1.4883560     1.0212740     1.3235091 
##    occupation        rating religiousness  yearsmarried 
##     1.0314027     0.6259691     0.7227292     1.0994093

可以看到婚龄增加一年，感情危机优势比将乘以1.10（保持年龄、宗教信仰程度和婚姻的自我评分不变）。

12.2.4 以概率的方式思考问题 predict()函数

对于大多数人来说，以概率的方式思考比使用优势比更直观，使用predict()函数可以让我们观察某个自变量在各个水平时对结果概率的影响。换言之，控制其他不变，观察某一个变量对结果的影响的概率是多少。

predict()函数的用法：

#线性回归模型（lm）：
predicted_values <- predict(lm_model, newdata = new_data)
#逻辑回归模型（glm）：
predicted_probabilities <- predict(glm_model, newdata = new_data, type = "response") #type = "response" 参数用于返回概率预测结果。

现在我们使用该方法评价对婚姻的自我评分对感情危机概率的影响。首先，创建一个虚拟数据集，设定年龄、婚龄和宗教信仰程度为它们的均值，对婚姻的自我评分的范围为1-5.

testdata <- data.frame(rating=c(1,2,3,4,5),age=mean(Affairs$age),yearsmarried=mean(Affairs$yearsmarried),religiousness=mean(Affairs$religiousness))
testdata

##   rating      age yearsmarried religiousness
## 1      1 32.48752     8.177696      3.116473
## 2      2 32.48752     8.177696      3.116473
## 3      3 32.48752     8.177696      3.116473
## 4      4 32.48752     8.177696      3.116473
## 5      5 32.48752     8.177696      3.116473

接着，使用虚拟数据集预测相应的概率：

testdata$prob <- predict(fit1,newdata = testdata,type = "response")
testdata

##   rating      age yearsmarried religiousness      prob
## 1      1 32.48752     8.177696      3.116473 0.5302296
## 2      2 32.48752     8.177696      3.116473 0.4157377
## 3      3 32.48752     8.177696      3.116473 0.3096712
## 4      4 32.48752     8.177696      3.116473 0.2204547
## 5      5 32.48752     8.177696      3.116473 0.1513079

12.3 泊松回归

当通过一系列连续和/或类别型自变量来预测计数型因变量时，泊松回归是一个有用的工具。

12.3.1 语法

fit <- glm(因变量～自变量1+自变量2+自变量3，data=,family=poisson())
summary(fit)

12.3.2举例

用robustbase包中的epilepsy癫痫数据为例，研究治疗条件(Trt)、年龄(Age)和前8周内的基础癫痫发病数(Base)对后8周内癫痫发病数(Ysum)的影响。

library(robustbase)

## Warning: 程辑包'robustbase'是用R版本4.3.2 来建造的

## 
## 载入程辑包：'robustbase'

## The following object is masked from 'Affairs':
## 
##     education

## The following object is masked from 'package:survival':
## 
##     heart

data(epilepsy,package = "robustbase")
head(epilepsy)

##    ID Y1 Y2 Y3 Y4 Base Age     Trt Ysum Age10 Base4
## 1 104  5  3  3  3   11  31 placebo   14   3.1  2.75
## 2 106  3  5  3  3   11  30 placebo   14   3.0  2.75
## 3 107  2  4  0  5    6  25 placebo   11   2.5  1.50
## 4 114  4  4  1  4    8  36 placebo   13   3.6  2.00
## 5 116  7 18  9 21   66  22 placebo   55   2.2 16.50
## 6 118  5  2  8  7   27  29 placebo   22   2.9  6.75

fit <- glm(Ysum~Base+Age+Trt,data =epilepsy,family = poisson())
summary(fit)

## 
## Call:
## glm(formula = Ysum ~ Base + Age + Trt, family = poisson(), data = epilepsy)
## 
## 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

12.3.3 解释模型参数

使用coef()函数获得模型参数，或者summary()函数输出结果中的coefficients表格。返回的数字表示的是模型中各个预测变量的系数估计值。这些数字表示了每个预测变量对于响应变量（泊松分布的计数数据）的影响程度。这里是loge(r)的对数形式。

coef(fit)

##  (Intercept)         Base          Age Trtprogabide 
##   1.94882593   0.02265174   0.02274013  -0.15270095

转化为指数形式：

exp(coef(fit))

##  (Intercept)         Base          Age Trtprogabide 
##    7.0204403    1.0229102    1.0230007    0.8583864

可以看到保持其他变量不变，年龄增加一岁，期望的癫痫发病数将乘以1.023.这说明年龄对癫痫病的发病影响较高。而trt的变化（从安慰剂到药物组）则癫痫发病乘以0.858，也就是说保持年龄不变，药物组相对于安慰剂组癫痫发病数大约降低了14%（1-0.858）。

完整教程请查看

R语言基础学习手册