stats | 线性回归(一)——模型表达式和输出结果

线性回归是形式最简单的回归模型,这里就不对其概念和含义做过多叙述了。本篇先简单介绍下书写线性回归表达式以及提取模型结果的方法,从中也可以理解R语言进行回归分析的基本语法。

使用的数据集来自R中自带的mtcars,这里从中提取三个变量用于后续分析。首先使用plot.data.frame函数观察三个变量的相互关系:

(不理解plot.data.frame函数的可以点击graphics | 基础绘图系统(五)——plot函数功能再探和低级绘图函数

data <- mtcars[, c("mpg", "wt", "qsec")]
plot(data)
  • 从中可以看出wtqsec两个变量都与mpg大致呈线性关系,其中前者更明显,因此下文就以mpg为因变量,wtqsec为自变量;

  • 这里不考虑变量的实际含义,只做技术分析。

stats包中做线性回归的函数为lm,语法结构如下:

lm(formula, data, subset, weights, na.action,
   method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE,
   singular.ok = TRUE, contrasts = NULL, offset, ...)

1 模型表达式和输出结果

1.1 表达式

formula参数就是模型的表达式,data参数是表达式中变量所在的数据框。表达式的基本格式为y ~ x,具体可见下表[1]

建立以下三个线性模型:

model.1 <- lm(mpg ~ wt, data = data)
model.2 <- lm(mpg ~ wt + qsec, data = data)
model.3 <- lm(mpg ~ wt + qsec + wt:qsec, data = data)
  • model.1:一元线性模型;

  • model.2:多元线性模型;

  • model.3:含交互项的线性模型。

使用as.formula函数构建模型表达式

xnam <- paste0("x", 1:25)
b <- paste("y ~ ", paste(xnam, collapse= "+"))
b

## [1] "y ~  x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11+x12+x13+x14+x15+x16+x17+x18+x19+x20+x21+x22+x23+x24+x25"

as.formula(b)

## y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11 + 
##     x12 + x13 + x14 + x15 + x16 + x17 + x18 + x19 + x20 + x21 + 
##     x22 + x23 + x24 + x25

使用updata函数往已有表达式中添加变量

formula.1 <- mpg ~ wt

formula.2 <- update(formula.1, . ~ . + qsec)
formula.2

## mpg ~ wt + qsec

该函数还能直接对模型使用:

update(model.1, .~. + qsec)

## 
## Call:
## lm(formula = mpg ~ wt + qsec, data = data)
## 
## Coefficients:
## (Intercept)           wt         qsec  
##     19.7462      -5.0480       0.9292

1.2 输出模型结果

使用summary函数即可预览模型结果的主要信息:

summary(model.1)

## 
## Call:
## lm(formula = mpg ~ wt, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5432 -2.3647 -0.1252  1.4096  6.8727 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
## wt           -5.3445     0.5591  -9.559 1.29e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared:  0.7528, Adjusted R-squared:  0.7446 
## F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10

使用fitted函数输出模型对因变量的拟合结果:

fit.1 <- fitted(model.1)

plot(data$wt, data$mpg, ylim = c(8, 35))
lines(data$wt, fit.1, col = "red")
arrows(data$wt, fit.1, data$wt, data$mpg,
       length = 0.1, col = "blue")

1.3 偏移项

偏移项也属于模型表达式的一部分,不同之处在于其系数固定为1,其数学含义是因变量已知的成分(或分量),因此不需要对其系数进行估计。偏移项在lm中使用offset参数表示。

model.11 <- lm(formula.1, data, offset = qsec)
summary(model.11)

## 
## Call:
## lm(formula = formula.1, data = data, offset = qsec)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5322 -2.0295 -0.2541  1.4675  5.7309 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  18.4098     1.5754   11.69 1.08e-12 ***
## wt           -5.0254     0.4691  -10.71 9.00e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.556 on 30 degrees of freedom
## Multiple R-squared:  0.8281, Adjusted R-squared:  0.8223 
## F-statistic: 144.5 on 1 and 30 DF,  p-value: 5.351e-13

在线性模型中,以下表达式等价:

为了验证这一点,比较下面的model.12和前面的model.11的结果:

data$diff <- data$mpg - data$qsec
model.12 <- lm(diff ~ wt, data)
summary(model.12)

## 
## Call:
## lm(formula = diff ~ wt, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5322 -2.0295 -0.2541  1.4675  5.7309 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  18.4098     1.5754   11.69 1.08e-12 ***
## wt           -5.0254     0.4691  -10.71 9.00e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.556 on 30 degrees of freedom
## Multiple R-squared:  0.7928, Adjusted R-squared:  0.7858 
## F-statistic: 114.8 on 1 and 30 DF,  p-value: 8.995e-12

也可以不生成新变量,直接将mpg - qsec当作因变量写在表达式的左半部分:

model.13 <- lm((mpg - qsec) ~ wt, data)
summary(model.13)

## 
## Call:
## lm(formula = (mpg - qsec) ~ wt, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5322 -2.0295 -0.2541  1.4675  5.7309 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  18.4098     1.5754   11.69 1.08e-12 ***
## wt           -5.0254     0.4691  -10.71 9.00e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.556 on 30 degrees of freedom
## Multiple R-squared:  0.7928, Adjusted R-squared:  0.7858 
## F-statistic: 114.8 on 1 and 30 DF,  p-value: 8.995e-12
  • 仔细观察以上三个模型的结果,model.12和model.13的输出结果完全一致,下文就以model.12代替两者,而model.11与它们的模型系数是一样的,但是R方和F检验结果是不一样的,这说明model.11与model.12并不完全等价。

model.11与model.12的残差是一致的:

head(residuals(model.11))

##         Mazda RX4     Mazda RX4 Wag        Datsun 710    Hornet 4 Drive 
##       -0.70328941        0.01818514       -2.56090653       -0.29318213 
## Hornet Sportabout           Valiant 
##        0.55753071       -3.14196148

head(residuals(model.12))

##         Mazda RX4     Mazda RX4 Wag        Datsun 710    Hornet 4 Drive 
##       -0.70328941        0.01818514       -2.56090653       -0.29318213 
## Hornet Sportabout           Valiant 
##        0.55753071       -3.14196148

但模型的拟合值不一样:

head(fitted(model.11))

##         Mazda RX4     Mazda RX4 Wag        Datsun 710    Hornet 4 Drive 
##          21.70329          20.98181          25.36091          21.69318 
## Hornet Sportabout           Valiant 
##          18.14247          21.24196

head(fitted(model.12))

##         Mazda RX4     Mazda RX4 Wag        Datsun 710    Hornet 4 Drive 
##          5.243289          3.961815          6.750907          2.253182 
## Hornet Sportabout           Valiant 
##          1.122469          1.021961
  • 这是因为model.11的因变量是mpg,而model.12的因变量是mpg - qsec的整体。也就是说,二者拟合的对象不同。

因为线性回归的R方等于真实值和拟合值相关系数的平方,可以做以下检验:

(cor(fitted(model.11), data$mpg))^2

## [1] 0.8260324

(cor(fitted(model.12), data$diff))^2

## [1] 0.7927538
  • 可以发现,手动计算的model.12的R方与模型输出结果是一致的,而model.11的计算结果和输出结果存在略微差异。

以下是分割线


为了探究为什么会存在以上差异,考虑到R方与相关系数平方相等只在线性模型中成立,而模型中加入偏移项的线性模型和普通的线性模型可能存在区别,因此可以使用R方的真正定义去计算:

yerror <- residuals(model.11)
ycenter = scale(data$mpg, scale = F)
1 - sum(yerror^2)/sum(ycenter^2)

## [1] 0.8259889

ycenter = scale(data$diff, scale = F)
1 - sum(yerror^2)/sum(ycenter^2)

## [1] 0.7927538

结果显示,计算出model.11的R方还是与输出结果不一致,反而和使用相关系数平方计算的结果更接近,而model.12的计算结果仍然和输出结果一致。

有清楚原因的读者可以后台联系~

参考资料

[1]

表: R语言实战(第2版)


### 多元线性回归残差分析在植被研究中的应用 #### 什么是多元线性回归残差? 多元线性回归种统计建模技术,用于探索多个自变量对单因变量的影响关系[^2]。在这种模型中,残差是指实际观察值与通过回归方程预测的值之间的差异。具体来说,对于每个样本点 \(i\),其残差定义为: \[ e_i = y_i - \hat{y}_i \] 其中,\(y_i\) 是第 \(i\) 个样本的实际观测值,而 \(\hat{y}_i\) 则是由回归模型预测出来的值。 #### 残差的意义及其重要性 残差提供了关于模型拟合质量的重要信息。理想情况下,个好的回归模型应该使得残差满足以下几个条件: - **均值接近于零**:表明模型没有系统性的偏差。 - **恒定方差(同方差性)**:意味着不同水平上的误差波动致。 - **正态分布假设成立**:这有助于后续推断检验的有效性。 - **独立性**:各次测量之间不存在关联[^3]。 这些特性可以通过绘制残差图来直观评估,并进步验证模型是否合理适用当前数据集。 #### 在植被研究中的具体应用场景 当涉及到生态学领域内的植被变化分析时,研究人员经常采用多元线性回归方法来量化环境因子如温度、降水量等因素如何共同作用影响某区域内的植物生长状况或者覆盖率情况等指标[^5]。然而,在构建这样的模型之后,还需要仔细考察所得结果背后隐藏的信息——即那些未能被解释清楚的部分,也就是所谓的“残差”。 例如,在项针对某地区多年间NDVI指数随气候变化模式的研究项目里,科学家们可能会先设定个初步假说认为该地区的年平均气温升高会促进当地草本类作物产量增加的同时减少灌木丛密度;与此同时降雨量增多则有利于提高整体绿色程度。于是他们建立了如下形式的标准线性回归表达式: \[ NDVI_t = β_0 + β_1TAVG_t + β_2PRECIP_t + ε_t \] 这里 NDVI 表示归化差分植被指数作为响应变量 (dependent variable),分别对应时间序列 t 的 TAVG PRECIP 分别代表同期气象记录里的日均温以及累积降水量数值充当解释变量 (independent variables) 。ε 符号标记的是随机扰动项或者说未预期因素造成的偏离现象。 完成上述基本框架搭建以后,下步便是执行正式运算过程并通过 MATLAB 软件包内置功能 regress 函数实现自动化参数估计流程[^4]: ```matlab % 假设 X 包含两列分别为 TAVG PRECIP 数据矩阵 Y 存储着对应的 NDVI 测量系列 [b,bint,r,rint,stats] = regress(Y,[ones(size(X,1),1) X]); disp(stats); % 显示 R-square F-test p-value 等关键性能衡量标准 plot(r,'o'); title('Residuals Plot'); ``` 运行结束后可以获得系列输出成果其中包括但不限于斜率系数 b 向量组描述各个输入维度权重贡献比例大小关系之外还有专门用来刻画剩余部分特性的 r 数组集合代表着每轮迭代过程中产生的个体级别层面的具体差距表现形态特征。 通过对这些残留成分深入剖析可以帮助识别潜在遗漏的关键驱动要素或者是发现异常极端案例的存在可能性从而指导改进原有理论构架使其更加贴近现实世界复杂多变的真实场景需求。 #### 结论 综上所述,多元线性回归不仅能够揭示主要控制变量间的相互联系规律而且还借助残差诊断工具挖掘深层次未知秘密推动科学研究不断向前迈进的步伐特别是在生态环境监测保护工作中发挥不可替代的作用价值所在之处显而易见值得高度重视并广泛推广运用实践当中去造福全人类社会可持续发展目标愿景达成贡献力量。 ---
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值