R语言-建模(广义)线性(加性、混合)模型-阿里云开发者社区

1、R 中的简单线性模型

$$ Y_i = \beta_0 + \beta_1X_i + \epsilon_i , \epsilon_i \sim N(0,\sigma^2)$$

mod <- lm(formula = mpg ~ wt,data = mtcars)
#mod <- lm(formula = mpg ~ wt + 0,data = mtcars) ### 生成不带截距项的简单线性模型
summary(mod)

Multiple R-squared:在一元线性回归中，Pearson相关系数R的平方等于线性回归的R²(Multiple R-squared)。这是因为在一元线性回归中，只有一个自变量和一个因变量，因此两个变量之间的线性相关程度和自变量对因变量的解释程度是一致的。但在多元线性回归中，Pearson相关系数R和线性回归的R²可能不相等，因为多元线性回归中有多个自变量，而Pearson相关系数R只能衡量两个变量之间的线性相关程度。

提取$\beta_1$ 的97.5% 置信区间：

mod$coefficients[[2]] - qnorm(0.975) * summary(mod)$coefficients[2,2]
mod$coefficients[[2]] + qnorm(0.975) * summary(mod)$coefficients[2,2]

预测新的数据： predict(mod,newdata = data.frame(wt = c(1,2,3))) ### 返回一个 vector 数据类型

2、R 中的多元线性模型

2.1 多元线性模型

$$ Y_i = \beta_0 + \beta_1X_i + \beta_2X_j + ... + \epsilon_i , \epsilon_i \sim N(0,\sigma^2)$$

mod <- lm(formula = mpg ~ wt + hp + am,data = mtcars)  #formula = y~.  表示纳入数据框的所有变量进行多元线性回归
summary(mod)

Multiple R-squared: 随着变量的增加，模型的拟合优度始终提升；
Adjusted R-squared: 根据变量数目进行调整R²，在多变量的前提下能更准确地反映模型的拟合优度，同时暗示变量不是越多越好。

由于R²无法回答模型是否统计显著的问题，进而提出了基于模型F统计量检验: $F = \frac {MSE}{MSM}$
其中MSE(mean squared error) = SSE/(n-p) 代表预测误差，MSM(mean squared model) = SSM/(p-1)代表模型误差。
在零假设下(模型无用)，MSE 服从自由度为(n-p)的卡方分布，MSM 服从自由度为(p-1)的卡方分布，从而$F_{modle} \sim F(n-p,p-1)$。当$R^2 \sim 1\ \rightarrow F\sim \infty$。
p为模型参数个数，n为拟合样本数。

修正已有拟合模型的实现:

mod_update <- update(mod, .~. + cyl) # 在已有模型上添加一个新变量，会重新拟合更新模型
mod_update <- update(mod, .~. - wt)  # 在已有模型上删减一个旧变量，会重新拟合更新模型
mod_update <- update(mod, log(.) ~. - am) # 在更新模型的因变量为 log 对数转换

2.2 多项式模型 Polynomial

$$ Y_i = \beta_0 + \beta_1X_i + \beta_2X_i*X_j + \beta_3X_j + ... + \epsilon_i , \epsilon_i \sim N(0,\sigma^2)$$
模型生成方式：mod <- lm(formula = mpg ~ wt + am + wt*am ,data = mtcars) 或 mod_new <- update(mod, .~. + am*wt)

2.3 多元回归的变量共线性问题

如果有部分变量互相相关，极端情况下如两变量本身呈严格线性关系$x_i = k * x_j$，会导致其中一个变量无法估计出其参数。一般情况下虽然并不会出现严格线性相关，但是在模型拟合的时候还是会出现这些变量无法准确估计其系数以及出现极高方差的问题。
如：

mtcars$wt2 <- 2*mtcars$wt
mod <- lm(formula = mpg ~ wt + wt + wt2, data = mtcars)
summary(mod)

mtcars$wt2 <- 2*mtcars$wt + rnorm(nrow(mtcars),0,0.1)
mod <- lm(formula = mpg ~ wt + wt + wt2, data = mtcars)
summary(mod)

2.4 多元线性回归模型的变量修剪- AIC 方法

AIC(Akaike Information Criterion) 是描述模型拟合优度指标,它可以通过删减不同变量的方式来比较不同模型的拟合优度,从而选择最优模型。虽然 Adjusted R-squared也能反映模型优度，不过一般用 AIC 的比较多。AIC越低意味着参数越少/误差越小，模型越好。

AIC(mod) ：直接计算一个模型的拟合优度；
step(mod) ：自动通过删减不同变量的方式(backward)来比较不同模型的AIC拟合优度，直到AIC分数不减反增的情况下保留的变量即最优模型对应的变量。该方式选出来的模型不一定是最好的，但不会至少不会太差。

被修剪变量合理性检验:
通过方差分析anova(mod_reduc,mod) 可以检验被删除掉的这些变量的系数是否等于零，也就是检验这些变量是否是对模型没有意义的。
$$F = \frac {(SSE_{mod\_reduc} - SSE_{mod})/(df_{mod\_reduc} - df_{mod})}{MSE_{mod}}$$

模型比较的 $F$统计量用于反映，当模型参数减少的时候，模型误差的变化情况。$H_0:$变量削减前后的模型的预测误差不变

如果减少了很多变量后，但是两个模型的拟合误差变化不大，则说明这些被修剪的变量确实是没有意义的；
如果减少了很多变量后，但是两个模型的拟合误差变化很大，则说明这些被修剪的变量是有意义，不应该被修剪；

mod <- lm(formula = mpg ~ .,data = mtcars)
mod_reduc <- lm(formula = mpg ~ wt + qsec + am, data = mtcars)
anova(mod,mod_reduc)
---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 
Analysis of Variance Table

Model 1: mpg ~ cyl + disp + hp + drat + wt + qsec + vs + am + gear + carb + 
    wt2
Model 2: mpg ~ wt + qsec + am
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     20 147.49                           
2     28 169.29 -8   -21.793 0.3694 0.9246 ### P值0.9246 接受原假设，可以认为mod中被删减的变量不影响模型效果

2.5 多元线性回归模型规范化-L1/L2项_ PLS(Penalized least square)

关联链接：ML-模型正则化 - 简书 (jianshu.com)

### LASSO mode--- L1
mod <- glmnet::cv.glmnet(x = as.matrix(mtcars[,c('disp','hp','drat')]), y = mtcars[["mpg"]])
### Ridge mode--- L2
mod <- glmnet::cv.glmnet(x = as.matrix(mtcars[,c('disp','hp','drat')]), y = mtcars[["mpg"]],alpha = 0)

2.6 模型的合理性检验

目的：检查模型是否符合假设;
一般线性模型都是基于误差独立同分布 的前提假设下$\epsilon \sim N(0,\sigma^2)$建立的。所以可以通过对模型的误差进行估计来判断模型是否符合假设前提。由于模型的真实误差不可知，所以实际上使用残差 $y_i-\hat y_i$来估计模型误差。

均值为零(残差应该均匀分布在零的两侧);
正态性检验(qq图/shapiro检验);
同分布(方差相同)

检查样本预测值标准化残差的分布 :plot(x = mod_reduc$fitted.values,y = rstandard(mod_reduc)) ，主要关注标准化残差落在[-2,2]区间外的样本数；

shapiro 检验

shapiro.test(mod_reduc$residuals)
----------------------------------------------------------------------------
    Shapiro-Wilk normality test

data:  mod_reduc$residuals
W = 0.9411, p-value = 0.08043 #P值大于0.05，残差较符合正态性

检查样本回归中的异常样本点:

Leverage score(hat values) , 该指标反映样本估计值的波动波动情况。
主要是检查 $hatvalue > 2*(p+1)/n$的样本点，$n$是回归的样本数，$p$是回归的变量个数;
```
which(hatvalues(mod_reduc) > 2*(3+1) / nrow(mtcars))
```
Cook's distance , 找出对回归结果影响较大的点。

3.1 R中的广义线性模型

广义线性模型是对一般线性模型的拓展来适应不同的建模需求，因变量与自变量的线性组合由连接函数进行联系，广义模型 较多应用于建模概率。

logit回归 ，因变量为二分类变量(0|1) ，对应如下预测函数
$$P( y = 1| x,\theta) = \frac {1}{1+e^{-\theta^Tx}} = h_{\theta}x, P( y = 0| x,\theta) = 1- h_{\theta}x \\ logit(P) = log(\frac {P}{1-P}) = \theta^Tx$$
```
mod_logit <- glm(formula = am ~ wt,family = binomial(link = "logit"),data = mtcars)     #binomial 描述影响变量是二分类
plot(x = mtcars$wt, y = mod_logit$fitted.values)
```
使用正态分布的CDF作为概率函数: mod_logit <- glm(formula = am ~ wt,family = binomial(link = "probit"),data = mtcars)

3.2 广义线性混合模型GLMM（Generalized linear mixed model）

这类模型的核心在于解决 随机效应(Random effect/mixed effect)，例如在实验设计中，可能存在不同的实验组或处理组。特别是在生物医学的应用上，特别需要考虑不同个体基本面(体质)的随机因素对于研究模型例如药物作用评估的影响。GLMM可以对这些随机效应进行建模，并考虑它们对数据的影响。注意 GLMM 模型建模需要认真识别 固定效应和随机效应。

$$ \small \bf{简单线性混合模型}： Y_i = \beta_0 + \beta_1X_i + A + \epsilon_i ； \ \\ A\sim N(0,\sigma_{A}^2);\ \ \epsilon_i \sim N(0,\sigma^2)$$

通常建模的随机效应的变量$A$具有以下特点：

类别型变量(批次)；
效应无法被重复（随机性）；
区别于固定效应$X$

R语言实现上 LMM 使用 lme4::lmer :

### 测试数据，研究睡眠干扰时长对反映能力的关系，18位受试者
### y = Reaction, x = Days, random_effect = Subject
mod <- lme4::lmer(Reaction ~ Days + (1|Subject),data = sleepstudy ) ### 模型仅考虑随机截距
mod <- lme4::lmer(Reaction ~ Days + (Days|Subject),data = sleepstudy )  ### 模型考虑随机 系数和截距

R语言实现上 GLMM 使用 lme4::glmer :
在简单混合线性模型中，引入链接函数，将线性预测器与响应变量进行联系建模 GLMM

### cbpp数据集, 建立牛群发病率与时期的关系，考虑牛群本身的随机因素和牛群大小对模型进行修正
mod <- glmer(incidence / size ~ period + (1 | herd), weights = size,family = binomial, data = cbpp)

3.2 广义加性模型GAM（Generalized additive model）

本质类似于多元线性模型，不考虑变量间的交互作用，用来建模Y~X之间的非线性关系。当变量与因变量的相关系数大于cor(X,Y) < 0.7 或无法准确建模线性回归的时候会考虑建模非线性关系。GAM模型中的每一项均为$X_i$的任意函数(非参数模型，不对变量进行参数学习)；
$$Y = f_1(X_1) + f_2(X_2) + .. + f_n(X_n) + \epsilon$$
R语言下创建加性模型 mgcv::gam :

X <- matrix(rnorm(1000*4),1000,4) %>% data.frame()
Y <- X[[1]] + log1p(abs(X[[2]])) + X[[3]]^2 + exp(X[[4]]) + rnorm(1000)
dat <- cbind(X,Y) %>% data.frame()
library("mgcv")
mod <- gam(Y~s(X1) + s(X2) + s(X3) + s(X4),data = dat)
plot(mod) ### 可视化每个变量Xi 对Y的影响
predict.gam(mod,newdata = data_new)

### mod <- gam(Y~s(X1) + s(X2) + s(X3) + s(X4),data = dat,family = binomial) ### 添加family 转换为广义加性模型

对于以上提到的所有线性模型方法，若需要生成广义模型，只需在建模函数的 family 参数上选择合适的链接函数。

Reference

【R语言-统计分析】线性回归模型的诊断与广义线性模型_哔哩哔哩_bilibili
【统计学博士教你用R语言做统计分析！】（PLS & GAM）_哔哩哔哩_bilibili
一文读懂回归模型准确度评价指标：R-square, AIC, BIC, Cp
非线性模型原理与R语言多项式回归、局部平滑样条、广义相加模型GAM分析_哔哩哔哩_bilibili

R语言-建模(广义)线性(加性、混合)模型