stats | 线性回归（二）——模型假设和模型估计

2 模型假设

线性模型有以下几个基本假设：

线性假设

线性假设是线性模型的一个很自然的要求，但并非是其独有的假设。它是指因变量所服从分布的均值可以使用自变量及其相关项的线性组合表示。

相关项是指自变量的非线性组合，如 、 等。以下表达式都是线性表达式：

• ；

• ；

• 。

正态性假设和同方差假设

正态性假设是线性模型区别于其他模型最根本的假设。它是指各样本对应的因变量都服从正态分布。同方差假设是指各样本对应的正态分布的方差相同。

结合线性假设，线性模型可以写成如下形式：

因为模型的残差 ，所以有 。正态性假设和同方差假设实际是要求各样本对应的模型残差服从同一个正态分布。

• 不一定是正态分布序列；

• 只有在零模型（即模型中没有自变量）的情况下，模型假设等价于所有 的正态分布均值都等于模型的截距 ，此时 是正态分布序列。

对于表达式：

若 不服从正态分布，则不是线性模型。

使用lm函数拟合以下表达式：

则模型假设 服从同方差的正态分布。

因变量的独立性

独立性是指各样本对应的因变量取值是相互独立的， 的取值不会影响 的取值，即要求 、 、...、 独立同分布。

自变量之间不存在多重共线性

3 模型估计

3.1 最小二乘估计和最大似然估计

模型估计中最常用的方法就是最大似然估计，目标是使似然函数取得最大值，而似然函数是由模型因变量的概率分布决定的。用于估计线性模型的最小二乘法实际上是最大似然估计在概率分布为同方差正态分布下的特例。

最小二乘估计的目标函数：

• 优化方向是使 取得最小值。

最大似然估计的目标函数：

• 表示因变量服从的概率密度函数；

• 优化方向是使 取得最大值。

由于线性回归假设 服从正态分布 ，那么，

因为概率密度函数恒为正数，求解最大似然函数的最大值可以对其取对数。所以，

3.2 使用函数输出模型估计结果

summary函数输出的模型结果中，估计结果位于Coefficients部分：

DATA <- mtcars[, c("mpg", "wt", "qsec")]
model <- lm(mpg ~ wt, DATA)
summary(model)

## 
## 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

summary函数输出的模型估计结果包括以下内容：

• Estimate：截距和自变量的系数；

• Std. Error：系数和截距的标准误；

coef或coefficients函数可以只输出系数和截距结果：

coef(model)
coefficients(model)

## (Intercept)          wt 
##   37.285126   -5.344472

confint函数可以输出系数和截距的置信区间：

confint(object, parm, level = 0.95, ...)

• object：模型对象；

• parm：自变量名称，缺省时默认对所有自变量的系数求置信区间；

• level：置信水平，默认为95%置信区间。

confint(model)

##                 2.5 %    97.5 %
## (Intercept) 33.450500 41.119753
## wt          -6.486308 -4.202635

confint(model, "(Intercept)", 0.99)

##                0.5 %   99.5 %
## (Intercept) 32.12166 42.44859

3.3 模型估计的解析解

记

可以求得模型估计的解析解：

据此，可以手动求解model的截距和系数估计值：

X = cbind(c(rep(1, 32)), DATA$wt)
Y = DATA$mpg

solve(t(X) %*% X) %*% t(X) %*% Y 

##           [,1]
## [1,] 37.285126
## [2,] -5.344472

• solve函数可以求矩阵的逆矩阵；

• t函数可以求矩阵的转置矩阵；

• %*%：矩阵乘法符号。

3.3 使用迭代法进行模型估计

线性模型作为一类经典模型，其系数和截距可以通过严格的公式求得解析解，但大多数模型都只能通过迭代的方法求取数值解的。这里使用迭代法来求线性模型的数值解。

线性模型的估计其实是个数学规划问题。决策变量：

目标函数：

约束条件：

stats中的optim函数可以解决以上优化问题：

optim(par, fn, gr = NULL, ...,
      method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN",
                 "Brent"),
      lower = -Inf, upper = Inf,
      control = list(), hessian = FALSE)

• par：决策变量的初始值；

• fn：目标函数。

# 构建目标函数
f <- function(b) sum((DATA$mpg - b[2]*DATA$wt - b[1])^2)

optim(c(0, 0), f)

## $par
## [1] 37.275657 -5.342921
## 
## $value
## [1] 278.3227
## 
## $counts
## function gradient 
##       95       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

• b1和b2的优化结果分别为37.275657和-5.342921，与模型截距和系数比较接近。

截距和系数的概率分布以及标准误可以通过以下蒙特卡洛试验得到：

Y <- fitted(model) # 模型拟合值
E <- residuals(model) # 模型残差
f <- function(b) sum((Y0 - b[2]*DATA$wt - b[1])^2)
R1  <- R2 <- c(1 : 10000) # 存储每次试验结果
set.seed(123)
for(i in c(1 : length(R1))) {
  Y0 <- Y + sample(E, 32) # 对残差进行随机排序后与拟合值相加得到新的因变量取值
  R1[i] <- optim(c(0, 0), f)$par[1] # 拟合新因变量取值与自变量的线性关系
  R2[i] <- optim(c(0, 0), f)$par[2]
}

library(car)
densityPlot(R1, main = "截距的模拟概率密度分布")
sd(R1) # 标准误

## [1] 1.789091

quantile(R1, probs = c(0.025, 0.975)) # 95%置信区间

##     2.5%    97.5% 
## 33.72962 40.74298

densityPlot(R2, main = "wt系数的模拟概率密度分布")
sd(R2)

## [1] 0.5560882

quantile(R1, probs = c(0.025, 0.975))

##     2.5%    97.5% 
## 33.72962 40.74298

• 注意：以上迭代算法并不严谨，仅供参考。

car工具包中的Boot函数使用自助法对模型进行迭代估计：

modelboot <- Boot(model, R = 10000)
summary(modelboot)

## 
## Number of bootstrap replications R = 10000 
##             original  bootBias  bootSE bootMed
## (Intercept)  37.2851  0.174304 2.34725 37.4438
## wt           -5.3445 -0.076318 0.71313 -5.3798

confint(modelboot)

## Bootstrap bca confidence intervals
## 
##                 2.5 %    97.5 %
## (Intercept) 32.709537 41.943489
## wt          -6.773288 -3.987207

---