Statistical inference 2st

在上一节中，我们学习了one-way ANOVA，并简单介绍了使用回归模型来做ANOVA。这节重点介绍回归分析。 回归是统计学的核心，它其实是一个广义的概念，通指那些用一个或多个解释变量（也称自变量）来预测响应变量（也称因变量、结果变量或标记变量）的方法。

让我们举一个简单的例子：假设您的导师要求你预测年销售额。影响销售的因素（驱动因素）可能有一百个。在这种情况下，销售额是您的响应变量。影响销售的因素是解释变量。回归分析可以帮助你回答以下问题：

1）哪些驱动因素对销售有重大影响？

2）哪个是最重要的销售推动力？

3）驱动力因素之间如何互动？

4）明年的年销售额是多少？

插播：R软件更新方法

很多同学在安装R包的时候遇到了R版本过低的问题，一下是更新R软件的方法：

install.packages("installr")

library(installr)

updateR()

接下来会弹出一个窗口，提示你是否要下载和安装最新的R软件版本。点击“是”并按照提示完成更新即可。

2 回归

2.1 回归分析的相关术语

极端异常值：假设数据集中有一个观测值与数据中的其他观测值相比具有非常高或非常低的值，即它不属于总体，一个值就能验证影响回归的结果，这样的观测值称为极端异常值。

多重共线性:当解释变量彼此高度相关时，这些变量被称为多重共线性。许多类型的回归技术都假设数据集中不应存在多重共线性。这是因为它会导致根据重要性对变量进行排序出现问题。或者它使得选择最重要的解释变量（因子）的工作变得困难。

欠拟合和过拟合： 当我们使用不必要的解释变量时，可能会导致过度拟合。过度拟合意味着我们的算法在训练集上表现良好，但无法在测试集上表现更好。它也被称为高方差问题。

当我们的算法表现不佳，甚至无法很好地拟合训练集时，就可以说它对数据拟合不足。它也被称为欠拟合或高偏差问题。

在下图中，我们可以看到：

1）左图拟合线性回归的直线会使数据拟合不足，导致较大的误差。

2）中图使用多项式拟合是平衡的，即这种拟合可以很好地适用于训练和测试集。

3）右图中的拟合将导致训练集的误差较低，但在测试集上效果不佳。

误差项（Error term）:用于解释响应变量与解释变量之间未能被模型捕捉到的部分。它代表了所有未被考虑或无法量化的因素对解释变量的影响。如未测量的变量、数据收集误差、测量误差、随机扰动等。

残差（residuals）:实际观测值与回归模型的预测值之间的差异,对于第i个观测值，其残差可以表示为：残差i = 观测值i - 预测值i。

对误差项和残差进行分析可以帮助我们评估回归模型的拟合程度和预测能力。

截距项（Intercept）：也被称为常数项或偏移量。在线性回归中，我们试图通过拟合一条直线来预测因变量的值。截距项表示当解释变量取值为0时，响应变量的预测值。

2.2 回归的类型

回归分析可以分为线性回归和非线性回归。线性回归假设解释变量和响应变量之间是线性关系，而非线性回归则假设二者之间是非线性关系。

同时，又可以分为简单回归和多元回归。简单回归只包含一个解释变量和响应变量，而多元回归则包含多个解释变量和一个响应变量。

在此基础上，有很多的变体，例如：

简单线性回归(Simple Linear Regression Model)： 用一个解释变量预测一个响应变量，适用于线性关系。

多项式回归(Polynomial Regression)：用一个解释变量预测一个响应变量，模型的关系是n阶多项式。

多元线性回归(Multiple Linear Regression)：用2个或以上的解释变量预测一个响应变量。

Logistic回归(Logistic Regression)：有学者翻译成逻辑回归，也有学者说这样翻译不合适。它其实际上是分类模型，常用于二分类。它是用一个或多个的解释变量预测一个类别型响应变量。响应变量本质上是二元的（有两个类别,比如：有/无，0/1，发生/不发生），解释变量可以是连续的或二元的。

时间序列回归(Time Series Regression)：用于建模因变量和时间之间的关系。通常用于预测金融数据、天气模式和其他时间相关现象。

泊松回归（Poisson Regression）：用一个或多个解释变量预测一个代表频数的响应变量，它拟合出的模型可以用于预测未来的计数数据，并进行因素分析和比较。它假设响应变量遵循泊松分布。

## `geom_smooth()` using formula = 'y ~ x'

Cox回归：用一个或多个解释变量预测一个事件随时间发生的概率。例如：从客户开设账户到流失的时间；癌症治疗后直至死亡的时间；从第一次心脏病发作到第二次心脏病发作的时间。

贝叶斯回归(Bayesian regression)：使用贝叶斯推断来建模因变量和自变量之间的关系。它允许将先验知识和不确定性纳入模型中。

非线性回归：当因变量和自变量之间的关系是非线性的时使用的技术。它可以用于建模复杂的关系，如指数或对数关系。

随机森林回归：使用决策树集合来创建强预测器的另一种机器学习算法。通常用于具有复杂变量之间相互作用的高维数据。

广义线性模型：扩展线性回归，允许因变量的分布不是正态分布。其中包括Logistic回归和泊松回归等模型。

混合效应模型：当数据中既有固定效应又有随机效应时使用的技术。通常用于纵向研究或存在重复测量的情况下。

下面我将介绍如何用R函数拟合最小二次回归模型(OLS，Ordinary least squares)、评价拟合优度、检验假设条件以及选择模型。包括简单线性回归、多项式回归和多元线性回归。

2.3 OLS回归

最小二乘回归法（OLS）是一种常用的线性回归分析方法，该方法通过求解残差平方和最小化的方式，来找到解释变量和响应变量之间的线性关系。

为了能够恰当的解释OLS模型系数，数据必须满足一下统计假设：

假设1：对于固定的解释变量值，响应变量值成正太分布

假设2：误差项在不同观测值之间是独立的。这意味着一个观测值的误差项不受其他观测值的误差项影响。

假设3：每组内没有明显的异常值。

假设4：同方差性（Homoscedasticity），响应变量的方差不随解释变量的水平不同而变化。如果存在异方差性（Heteroscedasticity），则回归模型的标准误差和置信区间会变得不准确，从而可能导致对系数的偏误估计。

示例 -随着收入的增加，食品消费的波动性也会增加。较贫穷的人总是吃便宜的食物，因此会花费相当稳定的金额；较富有的人有时可能会购买便宜的食物，有时会吃昂贵的饭菜。收入较高的人的食品消费变化较大。

2.2.1 lm()拟合回归模型

在R中，拟合线性模型最基本的函数就是lm()，格式为：

lmfit <- lm(formula, data)

其中， formula（表达式）指要拟合的模型形式， data是你的数据集，结果对象存储在lmfit，包含了所有拟合模型的大量信息。 formula表达式的形式如下：

Y ~ X1 + X2 + ... + Xk

Y为响应变量，X1、X2、X3、…、Xk代表各个解释变量，预测变量之间用+符号分隔。

如果要查看交互效应，可以在Xk后面 + X1:X2 + X1:Xk + X2:Xk + …. ，也可以使用表达式：Y ~ X1 * X2 * … * Xk。

例如：你有三个解释变量x1, x2, x3, 表达式可以写成：

Y ~ x1 * X2 * X3

等同于

Y ~ x1 + X2 + X3 + X1:X2 + X1:X3 + X2:X3 + X1:X2:X3

2.2.1 简单线性回归

#导入R内置数据
data("women")
#使用lm()函数拟合模型，这里的formula表达式是weight ~height（Y ~X），只有一个解释变量。这个模型是用来预测体重（weight）和身高（height）之间的关系的。

onewayfit <- lm(weight ~height, data = women)
#summary 模型数据
summary(onewayfit)

## 
## Call:
## lm(formula = weight ~ height, data = women)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7333 -1.1333 -0.3833  0.7417  3.1167 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -87.51667    5.93694  -14.74 1.71e-09 ***
## height        3.45000    0.09114   37.85 1.09e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.525 on 13 degrees of freedom
## Multiple R-squared:  0.991,  Adjusted R-squared:  0.9903 
## F-statistic:  1433 on 1 and 13 DF,  p-value: 1.091e-14

#打印weight的观测值
women$weight

##  [1] 115 117 120 123 126 129 132 135 139 142 146 150 154 159 164

#打印拟合值，即weight得预测值
fitted(onewayfit)

##        1        2        3        4        5        6        7        8 
## 112.5833 116.0333 119.4833 122.9333 126.3833 129.8333 133.2833 136.7333 
##        9       10       11       12       13       14       15 
## 140.1833 143.6333 147.0833 150.5333 153.9833 157.4333 160.8833

#打印残差，残差 = 观测值 - 预测值
residuals(onewayfit)

##           1           2           3           4           5           6 
##  2.41666667  0.96666667  0.51666667  0.06666667 -0.38333333 -0.83333333 
##           7           8           9          10          11          12 
## -1.28333333 -1.73333333 -1.18333333 -1.63333333 -1.08333333 -0.53333333 
##          13          14          15 
##  0.01666667  1.56666667  3.11666667

#可视化简单线性回归模型
plot(women$height, women$weight, xlab = "hight (in inches)", ylab ="weight (in inches)")
abline(onewayfit)

summary(onewayfit)

## 
## Call:
## lm(formula = weight ~ height, data = women)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7333 -1.1333 -0.3833  0.7417  3.1167 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -87.51667    5.93694  -14.74 1.71e-09 ***
## height        3.45000    0.09114   37.85 1.09e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.525 on 13 degrees of freedom
## Multiple R-squared:  0.991,  Adjusted R-squared:  0.9903 
## F-statistic:  1433 on 1 and 13 DF,  p-value: 1.091e-14

从summary的结果来看，模型的拟合优度很高。因为Multiple R-squared和Adjusted R-squared都接近于1，说明模型能够很好地解释因变量（体重）的变异。

从Coefficients部分可以看出，截距项（Intercept）的估计值为-87.51667，身高（height）对体重的系数估计值为3.45。这意味着单位身高增加对体重的影响约为3.45。身高和体重线性关系的预测式为：

\(\hat{weight}\) = - \(87.51667\) + \(3.45\) × \(height\)

t统计量和对应的p值表明身高对体重的影响是显著的，因为p值（Pr(>|t|)）远小于0.05。

Residual standard error是评估模型拟合好坏的一个指标，果模型能够完全拟合数据，则残差为0，否则残差不为0。 Residual standard error则表示残差的标准差，即实际数据点和线性回归模型的拟合直线之间的差距的平均值。该拟合模型的Residual standard error值为1.525。

综上所述，这个线性回归模型在解释体重和身高之间的关系方面效果很好，身高对体重有显著影响，模型的拟合优度高。

2.2.2 多项式回归

尽管上面的简单线性回归已经很好的拟合了women数据集中的weight和hight的关系，但从plot图中可以看出，可以添加一个二次项（即\(x^2\)）来提高回归的精确度。

#拟合含二次项的等式
PolynomialFit<- lm(weight ~ height +I(height^2), data=women)

#拟合含二次项的等式
summary(PolynomialFit)

## 
## Call:
## lm(formula = weight ~ height + I(height^2), data = women)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.50941 -0.29611 -0.00941  0.28615  0.59706 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 261.87818   25.19677  10.393 2.36e-07 ***
## height       -7.34832    0.77769  -9.449 6.58e-07 ***
## I(height^2)   0.08306    0.00598  13.891 9.32e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3841 on 12 degrees of freedom
## Multiple R-squared:  0.9995, Adjusted R-squared:  0.9994 
## F-statistic: 1.139e+04 on 2 and 12 DF,  p-value: < 2.2e-16

身高和体重线性关系的新的预测式为：

\(\hat{weight}\)= \(261.87818\) - \(7.34832\) × \(height\) + \(0.08306\) × \(height^2\)

从summary结果看Residual standard error从1.525降低到了0.3841，Multiple R-squared和Adjusted R-squared: 从0.991和0.9903 升高到了0.9995和0.9994，说明模型拟合的更加精确。从plot的曲线看，该数据集通过多项式回归确实比简单线性模型拟合的好。

2.2.3 多元线性回归

当预测变量不止一个时，简单线性回归就变成了多元线性回归，分析也稍微复杂些。

以基础包中的state.x77为例，我们想探究一个州的犯罪率和其他因素的关系，包括人口、文盲率、收入和结霜天数（温度在冰点一下的平均天数）。

state.x77是个矩阵，需要先转化成数据集格式，用如下转化：

state <- as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])
head(state)

##            Murder Population Illiteracy Income Frost
## Alabama      15.1       3615        2.1   3624    20
## Alaska       11.3        365        1.5   6315   152
## Arizona       7.8       2212        1.8   4530    15
## Arkansas     10.1       2110        1.9   3378    65
## California   10.3      21198        1.1   5114    20
## Colorado      6.8       2541        0.7   4884   166

• Murder：谋杀率
• Population：人口
• Illiteracy：文盲率
• Income：收入
• Frost：霜冻日数

计算各个变量之间的相关系数,它将返回一个与原始数据集变量数目相等的方阵。每个元素表示两个变量之间的相关性，取值范围从-1到1。

cor(state)

##                Murder Population Illiteracy     Income      Frost
## Murder      1.0000000  0.3436428  0.7029752 -0.2300776 -0.5388834
## Population  0.3436428  1.0000000  0.1076224  0.2082276 -0.3321525
## Illiteracy  0.7029752  0.1076224  1.0000000 -0.4370752 -0.6719470
## Income     -0.2300776  0.2082276 -0.4370752  1.0000000  0.2262822
## Frost      -0.5388834 -0.3321525 -0.6719470  0.2262822  1.0000000

创建一个矩阵型散点图

suppressWarnings({
library(car)
scatterplotMatrix(state, spread = FALSE, smoother.args = list(lty=2),main ="Scatter plot matrix")
})

## Loading required package: carData

代码解读：
suppressWarnings({ … })：这部分代码用于抑制在运行过程中产生的警告信息，以确保在执行时不会干扰结果的输出。
spread = FALSE：设置了 spread 参数为 FALSE，表示不展开变量，而是以紧凑方式显示散点图矩阵。
smoother.args = list(lty=2)：指定了平滑线的参数，其中 lty=2 表示使用虚线。

结果解读：
scatterplotMatrix()函数默认在非对角区域绘制变量间的散点图，并添加最佳拟合直线（实线）和最佳拟合曲线（虚线）。对角线区域绘制每个变量的核密度估计图（Kernel Density Estimate，KDE),是一种非参数方法的概率密度估计图，。
从图中可以看到， 谋杀率Murder是双峰的曲线，每个预测变量都一定程度上出现了偏率。谋杀率随着人口和文盲率的增加而增加，随着收入水平和结霜天数增加而下降。同时，越冷的州府文盲率越低，收入水平越高。

使用lm()函数拟合多元线性回归模型

mulFit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = state)
summary(mulFit)

## 
## Call:
## lm(formula = Murder ~ Population + Illiteracy + Income + Frost, 
##     data = state)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.7960 -1.6495 -0.0811  1.4815  7.6210 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.235e+00  3.866e+00   0.319   0.7510    
## Population  2.237e-04  9.052e-05   2.471   0.0173 *  
## Illiteracy  4.143e+00  8.744e-01   4.738 2.19e-05 ***
## Income      6.442e-05  6.837e-04   0.094   0.9253    
## Frost       5.813e-04  1.005e-02   0.058   0.9541    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.535 on 45 degrees of freedom
## Multiple R-squared:  0.567,  Adjusted R-squared:  0.5285 
## F-statistic: 14.73 on 4 and 45 DF,  p-value: 9.133e-08

当解释变量不止一个时，回归系数的含义为：一个预测变量增加一个单位，其他预测变量保持不变时，响应变量将要增加的数量。

例如本例中，文盲率的回归系数为4.14， 表示控制人口、收入和温度不变时，文盲率上升1%，谋杀率将会上升4.14%，它的系数在p = 2.19e-05 (P<0.001) 的水平下显著不为0（0代表不相关）,即有文盲率和谋杀率显著相关性。相反，Frost的系数没有显著不为0（P=0.954）,说明放控制其他变量不变时，Frost与Murder不呈线性相关。

总体来看，所有的解释变量解释了各州谋杀率56.7%的方差。换句话说，43.3%的方差还未被解释。这很容易理解，影响犯罪率的因素肯定不止人口、文盲率、收入和结霜天数等变量，若想需求更多的解释，需要收集更多的因素。

2.2.4 有交互项的多元线性回归

许多有趣的研究都会涉及交互项的解释变量，以mtcars数据框中的汽车数据为例，若你对汽车的重量和马力感兴趣，可以把它们作为解释变量，并包含交互项来拟合回归模型。

data("mtcars")
head(mtcars)

##                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
## Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

str(mtcars)

## 'data.frame':    32 obs. of  11 variables:
##  $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
##  $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
##  $ disp: num  160 160 108 258 360 ...
##  $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
##  $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
##  $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
##  $ qsec: num  16.5 17 18.6 19.4 17 ...
##  $ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
##  $ am  : num  1 1 1 0 0 0 0 0 0 0 ...
##  $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
##  $ carb: num  4 4 1 1 2 1 4 2 2 4 ...

• mpg：每加仑的英里数（英里/加仑），表示汽车的燃油经济性能。
• cyl：汽缸数，表示汽车引擎的气缸数量。
• disp：排量（立方英寸），表示发动机的总体积。
• hp：马力，表示发动机的功率。
• rat：后桥齿轮比，表示后轮与引擎转速之间的比例。
• wt：重量（千磅），表示汽车的重量。
• qsec：四分之一英里所需时间（秒），表示汽车从静止加速到四分之一英里所需的时间。
• vs：发动机配置，0表示V字型，1表示直列型。
• am：传动方式，0表示自动传动，1表示手动传动。
• gear：前进档位数，表示汽车的前进档位数量。
• carb：化油器数量，表示汽车的化油器数量。

hpwtfit <- lm(mpg ~ hp + wt +hp:wt, data = mtcars)
summary(hpwtfit)

## 
## Call:
## lm(formula = mpg ~ hp + wt + hp:wt, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0632 -1.6491 -0.7362  1.4211  4.5513 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 49.80842    3.60516  13.816 5.01e-14 ***
## hp          -0.12010    0.02470  -4.863 4.04e-05 ***
## wt          -8.21662    1.26971  -6.471 5.20e-07 ***
## hp:wt        0.02785    0.00742   3.753 0.000811 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.153 on 28 degrees of freedom
## Multiple R-squared:  0.8848, Adjusted R-squared:  0.8724 
## F-statistic: 71.66 on 3 and 28 DF,  p-value: 2.981e-13

结果解读：

该模型用于预测汽车的每加仑英里数（mpg），其中自变量包括马力（hp）、车重（wt）以及马力和车重的交互作用（hp:wt）。

1）拟合的模型方程：

\(\hat{mpg}\) = \(49.80842\) - \(0.12010\) × \(hp\) - \(8.21662\) × \(wt\) + \(0.02785\) × \(hp:wt\)

2）模型拟合度：

• 多重 R-squared（决定系数）为 0.8848，表示模型能够解释因变量（mpg）的总变异的 88.48%。
  
• 调整 R-squared（调整决定系数）为 0.8724，考虑到自变量的数量，该值较为可靠。
  

3）模型系数：

• 截距 (Intercept) 的估计值为 49.80842，表示当 hp 和 wt 的值都为 0 时，mpg 的估计值为 49.80842。
  
• hp 的系数估计值为 -0.12010，表示每增加一个单位的 hp，mpg 会减少 0.12010 个单位。
  
• wt 的系数估计值为 -8.21662，表示每增加一个单位的 wt，mpg 会减少 8.21662 个单位。
  
• hp:wt 的交互项系数估计值为 0.02785，表示 hp 和 wt 的交互作用对 mpg 的影响。
  

4）显著性检验：

• 截距项和所有自变量的系数估计都具有统计显著性，因为它们的 p-value 均小于 0.05。
  
• hp:wt 的交互项系数估计也具有统计显著性，因为其 p-value 小于 0.001。
  

5）残差分析：

• 残差的最小值为 -3.0632，最大值为4.5513。残差是观测值与模型预测值之间的差异，这些值应该接近于正态分布并且没有明显的模式。
  
• 残差标准误差为 2.153，表示预测值与实际观测值之间的平均差异约为 2.153 个单位。
  

综合来看，该回归模型对解释 mpg 的变异性具有较高的拟合度（R-squared = 0.8848），说明自变量 hp 和 wt 可以较好地解释 mpg 的变化。同时，交互项 hp:wt 也对 mpg有显著影响。

2.2.5 选择“最佳”回归模型

尝试获取一个回归方程时，实际上你就面对从众多可能的模型中做出选择的问题。是不是 所有的变量都要包括？还是去掉那个对预测贡献不显著的变量？是否需要添加多项式项和/或交互项来提高拟合度？最终回归模型的选择总是涉及预测精度（模型尽可能的拟合数据）与模型简洁度（一个简单且能复制的模型）的平衡问题。比如：如果有两个几乎相同预测精度的模型，我们肯定喜欢简单的那个。

本小节讨论的问题，就是如何在候选模型中进行筛选。注意，“最佳”是打了引号的，因为没有评价的唯一标准，最终的决定取决于你的评判。

2.2.5.1 模型比较

用anova()函数可以比较两个嵌套模型的拟合优度。所谓嵌套模型，即它的一些项完全包含在另一个模型中。

以state.x77转换的state数据为例。在state的多元回归模型中，我们发现Income和Frost的回归系数不显著，此时你可以分别检验不含和包含这两个变量的模型，然后比较两项模型的效果是否一样好。

使用anova()函数比较

state <- as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])
mulFit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = state)
mulFit2 <- lm(Murder ~ Population + Illiteracy, data = state)
anova(mulFit1, mulFit2)

## Analysis of Variance Table
## 
## Model 1: Murder ~ Population + Illiteracy + Income + Frost
## Model 2: Murder ~ Population + Illiteracy
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     45 289.17                           
## 2     47 289.25 -2 -0.078505 0.0061 0.9939

结果解读：

模型1的残差自由度（Res.Df）为45，残差平方和（RSS）为289.17。而模型2的残差自由度为47，残差平方和为289.25。
接下来的列显示了模型2相对于模型1的差异。差异的平方和（Sum of Sq）为-0.078505，F值为0.0061，p值（Pr(>F)）为0.9939。
根据p值，我们可以得出结论：差异是不显著的，也就是说模型1和模型2在解释谋杀率方面没有显著的差异。

因此，我们可以得出结论：可以将Income 和 Frost变量从模型中删除。

2.2.5.2 变量选择

从大量候选变量中选择最终的预测变量有两种主要的方法：逐步回归法（stepwise method）和全子集回归（all-sub）

1. 逐步回归

逐步回归中，模型会一次添加或删除一个变量，直到达到某个判停准则为止。
MASS()函数可以实现逐步回归模型（向前、向后和向前向后），依据的是精准AIC准则。

library(MASS)
mulFit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = state)
stepAIC(mulFit, direction = "backward") # 可选"both", "backward", "forward"

## Start:  AIC=97.75
## Murder ~ Population + Illiteracy + Income + Frost
## 
##              Df Sum of Sq    RSS     AIC
## - Frost       1     0.021 289.19  95.753
## - Income      1     0.057 289.22  95.759
## <none>                    289.17  97.749
## - Population  1    39.238 328.41 102.111
## - Illiteracy  1   144.264 433.43 115.986
## 
## Step:  AIC=95.75
## Murder ~ Population + Illiteracy + Income
## 
##              Df Sum of Sq    RSS     AIC
## - Income      1     0.057 289.25  93.763
## <none>                    289.19  95.753
## - Population  1    43.658 332.85 100.783
## - Illiteracy  1   236.196 525.38 123.605
## 
## Step:  AIC=93.76
## Murder ~ Population + Illiteracy
## 
##              Df Sum of Sq    RSS     AIC
## <none>                    289.25  93.763
## - Population  1    48.517 337.76  99.516
## - Illiteracy  1   299.646 588.89 127.311

## 
## Call:
## lm(formula = Murder ~ Population + Illiteracy, data = state)
## 
## Coefficients:
## (Intercept)   Population   Illiteracy  
##   1.6515497    0.0002242    4.0807366

开始时模型包含4个（全部）变量，然后每一步中删除一个变量，删除的变量选择是根据你输入的表达式中变量的排序。

每个模型中各自计算了：
- AIC值（赤池信息准则）：用于比较不同模型的拟合质量，越小表示模型拟合得越好。
- Df（自由度）表示每个模型中的自由度数量。
- Sum of Sq（平方和）表示每个自变量被移除后的残差平方和。
- RSS（残差平方和）表示模型的残差平方和。

在第一步中，“Frost”被移除，因为它对模型的改进微乎其微。在第二步中，“Income”被移除，因为它对模型的改进也很小。最终的模型只包含”Population”和”Illiteracy”两个自变量。

最终模型的回归方程为：

Murder = 1.6515497 + 0.0002242 × Population + 4.0807366 × Illiteracy

逐步回归法其实存在争议，虽然它可能会找到一个好的模型，但是不能保证模型就是“最佳”模型，因为不是每一个可能的模型都被评价了。为了克服这个限制，便有了全子集回归法。

2. 全子集回归

全自己回归就是所有可能的模型都会被检验。

我们可以选择展示所有可能的结果，也可以展示n个不同子集大小的最佳模型。例如，若nbest=2， 先展示两个最佳的单预测变量模型，然后展示两个最佳的双预测变量模型，以此类推，直到包含所有的预测变量。

library(leaps)
leaps <- regsubsets(Murder ~ Population + Illiteracy + Income + Frost, data = state, nbest = 4)
plot(leaps, scale = "adjr2") 

• “adjr2”（默认值）：调整的P平方（Adjusted R-squared）。

• R平方（R-squared）是一种常用的度量指标，用于衡量回归模型对因变量的解释程度。R平方的取值范围在0到1之间，越接近1表示模型对数据的拟合程度越好。

• 调整的R平方（Adjusted R-squared）是一种修正后的R平方值，它考虑了回归模型中使用的自变量数量对模型拟合度的影响。调整的R平方通过惩罚使用过多自变量的模型，以防止过拟合现象。调整的R平方在模型选择中更常用，因为它能够更准确地衡量模型的拟合程度和复杂性之间的平衡。

• 非调整的R平方（Non-adjusted R-squared）是指在回归模型中使用的解释变量数量不考虑的情况下计算得出的R平方值。它测量了所有自变量对因变量解释的总体效果，但没有考虑到使用的自变量数量可能对模型拟合度的影响。

结果解读：

乍看图比较费解。第一行（图底部开始）代表含intercept(截距项)和income的模型调整R平方和为0.033， 第二行代表含intercept和Popluation的模型调整R平方和为0.1。跳到第十二行（顶部第二行），你会看到含intercept、Popluation、Illiteracy、income的模型调整R平方和为0.54，而含intercept、Popluation、Illiteracy的模R平方为0.55。 此处，你会发现含预测变量越少的模型调R平方越大（对于非调整的R平方，这是不可能的）。

图形表明，双预测变量模型（Popluation、Illiteracy）是最佳模型。

2.3 广义线性回归

2.3.1 广义线性模型和glm()函数

2.3.2 Logistic回归

2.3.3 泊松回归