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]区间外的样本数;

plot(mod_reduc)第二幅图

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 , 找出对回归结果影响较大的点。
plot(mod_reduc)第四幅图红色虚线上标注的样本点对回归结果影响较大

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 参数上选择合适的链接函数。


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

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 203,179评论 5 476
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,229评论 2 380
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 150,032评论 0 336
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,533评论 1 273
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,531评论 5 365
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,539评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,916评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,574评论 0 256
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,813评论 1 296
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,568评论 2 320
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,654评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,354评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,937评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,918评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,152评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,852评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,378评论 2 342

推荐阅读更多精彩内容