70-线性回归与方差分析
1、线性回归
> library(pacman)
> p_load(data.table,dplyr,dbplyr)
> df <- readr::read_csv("./data_set/dailymarket.csv")
>
> df$date <- as.Date(df$date,format = "%d/%m/%y")
> str(df)
## tibble [31 x 4] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ date : Date[1:31], format: "2020-01-01" "2020-01-02" ...
## $ revenue : num [1:31] 5064 6115 5305 3185 4182 ...
## $ customers: num [1:31] 100 67 33 111 34 222 44 56 11 66 ...
## $ profits : num [1:31] 2026 2691 1592 318 3764 ...
## - attr(*, "spec")=
## .. cols(
## .. date = col_character(),
## .. revenue = col_double(),
## .. customers = col_double(),
## .. profits = col_double()
## .. )
1.1 理解回归的汇总结果
> fit <- lm(profits ~ revenue + customers,data = df)
> summary(fit)
##
## Call:
## lm(formula = profits ~ revenue + customers, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1728.75 -470.17 -117.61 47.29 2517.89
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 129.4002 1137.2538 0.114 0.910
## revenue 0.3567 0.2226 1.603 0.120
## customers 5.3795 4.7479 1.133 0.267
##
## Residual standard error: 1077 on 28 degrees of freedom
## Multiple R-squared: 0.1315, Adjusted R-squared: 0.06948
## F-statistic: 2.12 on 2 and 28 DF, p-value: 0.1389
Call:
lm(formula = profits ~ revenue + customers, data = df)
调用：当创建模型时，以上代码表明lm如何被调用。
Residuals:
Min 1Q Median 3Q Max
-1728.75 -470.17 -117.61 47.29 2517.89
残差统计项：理想情况下，回归残差将有一个完美的正态分布。普通最小二乘法在数学上保证产生均值为零的残差，因此，中位数(Median)的符号表示偏斜的方向，而中位数的大小表示偏斜的程度。本例中中位数为-117.61，表示向左偏斜。
如果残差有一个很好的钟形分布，第一四分位和第三四分位应该有大约相同的幅度。本例中有较大的1Q，虽然中位数为负，但数据依然有较大的左偏。
最小残差和最大残差提供了一种快速的方法来检测数据中的极端离群值，因为极端离群值（在响应变量中）产生较大的残差。
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 129.4002 1137.2538 0.114 0.910
revenue 0.3567 0.2226 1.603 0.120
customers 5.3795 4.7479 1.133 0.267
系数：标记为Estimate的列包含由普通最小二乘法计算出的估计回归系数。
从理论上来说，如果一个变量的系数为零，那么该变量是毫无意义的，它对模型毫无贡献。然而这里显示的系数只是估计值，它们不会正好为零，那么从统计的角度而言，真正的系数是0的可能性有多大？这是t统计量和p值的目的，它们在汇总中标记为t value和Pr(>|t|) 。
p值是一个概率，它估计系数不显著的可能性，所以越小越好。本例中所有的系数在传统界限0.05时均不显著。有较大p值的变量是可以从模型中移除的候选变量。
Signif. codes行给出了标记项的含义。
Std.Error列是估计的回归系数的标准误差；标记为t value的列是t统计量，p值是据此计算出来的。
Residual standard error: 1077 on 28 degrees of freedom
显示残差标准误差。
Multiple R-squared: 0.1315, Adjusted R-squared: 0.06948
R2（判定系数）是衡量模型拟合质量的指标，它的值越大越好。在数学上，它是回归模型所能解释的响应变量y的方差比例，剩余的方差是模型不能解释的，必须由其他因素解释，例如未知变量或样本的变化。本例中模型解释y的方差修正后为0.02764，剩余的97%未被解释。
修正R2（Adjusted）相对于Multiple考虑了模型中变量的数目，所以能实际地评估模型的有效性。
F-statistic: 2.12 on 2 and 28 DF, p-value: 0.1389
F统计量告诉你模型是否显著，如果有任何回归系数为非零，则该模型就是显著的，如果所有系数均为零，模式就是不显著的。按照惯例，本例中p值大于0.05，表明模型可能是不显著的。
一般来说，先看F统计量，如果不显著，其他的就都不重要了。
1.2 运行无截距项的线性回归
线性回归通常包括截距项，所以它在R中是默认的。极少数情况下可能需要假设截距项为零来拟合数据，这时在回归公式右侧加上“+0”就行了。这种情况下，我们假设模型在x为0时，y也应该为0。
我们通过上面的回归检验下该假设条件，先求出截距项的置信区间：
> confint(fit)
## 2.5 % 97.5 %
## (Intercept) -2.200159e+03 2458.9590177
## revenue -9.920849e-02 0.8126228
## customers -4.346165e+00 15.1050784
因为截距项的置信区间包含0，所以从统计学来说，截距可能为0。另外，从数据理解上来说，如果收入（revenue）和顾客数量（customers）为0，那么利润（profits）肯定也为0。
> fit2 <- lm(profits ~ revenue + customers + 0,data = df)
> summary(fit2)
##
## Call:
## lm(formula = profits ~ revenue + customers + 0, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1735.75 -460.88 -130.43 49.63 2521.35
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## revenue 0.38081 0.06708 5.677 3.87e-06 ***
## customers 5.46728 4.60432 1.187 0.245
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1059 on 29 degrees of freedom
## Multiple R-squared: 0.8314, Adjusted R-squared: 0.8198
## F-statistic: 71.52 on 2 and 29 DF, p-value: 6.14e-12
从结果中可以看到，当假设截距项为0时，模型显著了，revenue的系数也显著了，修正R2的值也提升到了0.8198。
1.3 运行有交互项的线性回归
在回归中，当两个预测变量的乘积也是一个显著的预测变量时，交互作用就会出现。R语法通过在变量之间插入一个*号来表示。
> fit3 <- lm(profits ~ revenue * customers + 0,data = df)
> summary(fit3)
##
## Call:
## lm(formula = profits ~ revenue * customers + 0, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1748.12 -424.03 -170.53 62.21 2538.09
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## revenue 0.376763 0.067721 5.563 5.96e-06 ***
## customers -5.667860 14.929434 -0.380 0.707
## revenue:customers 0.002201 0.002805 0.785 0.439
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1066 on 28 degrees of freedom
## Multiple R-squared: 0.8351, Adjusted R-squared: 0.8174
## F-statistic: 47.25 on 3 and 28 DF, p-value: 4.397e-11
R会自动在模型中包括预测变量，如本例中的revenue和customers，以及它们的乘积，它基于这样一个推理：如果一个模型包含交互项，那么回归理论告诉我们模型还应该包含交互项组成的变量。
同样，如果有更多变量，它们之间插入*号表示它们之间的所有交互作用。如果不需要每个可能的交互，使用:号可以明确指定一个单一乘积。（两个变量时只有一种交互可能）
1.4 选择最合适的回归变量
使用step()函数可以执行逐步回归，前向或者后向。step()函数使搜索自动化，后向逐步回归是最简单的方法，开始时，它包括所有的预测变量，我们称它为全模型，然后函数逐步移除不显著的变量，直到模型达到最优时，我们称为简化模型。
> fit.reduce <- step(fit3,direction = "backward")
## Start: AIC=435.08
## profits ~ revenue * customers + 0
##
## Df Sum of Sq RSS AIC
## - revenue:customers 1 699428 32510399 433.76
## <none> 31810970 435.08
##
## Step: AIC=433.76
## profits ~ revenue + customers - 1
##
## Df Sum of Sq RSS AIC
## - customers 1 1580652 34091050 433.23
## <none> 32510399 433.76
## - revenue 1 36134185 68644583 454.92
##
## Step: AIC=433.23
## profits ~ revenue - 1
##
## Df Sum of Sq RSS AIC
## <none> 34091050 433.23
## - revenue 1 158773808 192864858 484.95
> summary(fit.reduce)
##
## Call:
## lm(formula = profits ~ revenue - 1, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1920.89 -302.13 -231.13 -39.43 2324.90
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## revenue 0.4468 0.0378 11.82 8.13e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1066 on 30 degrees of freedom
## Multiple R-squared: 0.8232, Adjusted R-squared: 0.8173
## F-statistic: 139.7 on 1 and 30 DF, p-value: 8.127e-13
最后只保留了一个变量revenue。
后向逐步回归很容易运行，但有时开始包含所有的变量是不可行的，因为有太多的候选变量。这种情况下使用前向逐步回归，开始不包含变量，然后逐步增加变量以优化回归，当模型不能再改进时停止。
> fit.min <- lm(profits ~ 0,data = df)
> fit.fwd <- step(fit.min,
+ direction = "forward",
+ scope = (~ revenue + customers + revenue:customers),
+ trace = 1)
## Start: AIC=484.95
## profits ~ 0
##
## Df Sum of Sq RSS AIC
## + revenue 1 158773808 34091050 433.23
## + customers 1 124220275 68644583 454.92
## <none> 192864858 484.95
##
## Step: AIC=433.23
## profits ~ revenue - 1
##
## Df Sum of Sq RSS AIC
## <none> 34091050 433.23
## + customers 1 1580652 32510399 433.76
> summary(fit.fwd)
##
## Call:
## lm(formula = profits ~ revenue - 1, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1920.89 -302.13 -231.13 -39.43 2324.90
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## revenue 0.4468 0.0378 11.82 8.13e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1066 on 30 degrees of freedom
## Multiple R-squared: 0.8232, Adjusted R-squared: 0.8173
## F-statistic: 139.7 on 1 and 30 DF, p-value: 8.127e-13
前向和后向最后达成了相同的模式。在实际应用中，建议尝试前向和后向两种回归方法，然后比较结果。
最后，需要注意的是，逐步回归不能替代你小心明智地选择预测变量的过程，它并不是万能的。
1.5 对数据子集回归
参数subset可以指定一部分数据拟合线性模型，该参数的值可以是任何用于索引数据的表达式。
当使用样本内数据创建模型和样本外数据来检验模型时，可以使用该参数。另外，如果有一个因子标识了数据的来源，可以通过该参数将回归限制在标识了的数据上。
> fit.subset <- lm(profits ~ revenue + 0,
+ data=df,subset = 1:floor(length(revenue)*2/3))
> summary(fit.subset)
##
## Call:
## lm(formula = profits ~ revenue + 0, data = df, subset = 1:floor(length(revenue) *
## 2/3))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1721.90 -251.39 -157.20 -80.35 2146.51
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## revenue 0.43090 0.04343 9.921 5.98e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 983.2 on 19 degrees of freedom
## Multiple R-squared: 0.8382, Adjusted R-squared: 0.8297
## F-statistic: 98.43 on 1 and 19 DF, p-value: 5.975e-09
1.6 在回归公式中使用表达式
回归公式的语法禁止对计算出的值进行回归，因为基本二元运算符（+-*/^）都有特殊的含义，比如profits ~ revenue^2 ，R会将其解释为一个交互项，而不是revenue2。这时候可以使用I()函数将计算嵌入到表达式后再进行回归。
> fit4 <- lm(profits ~ I(revenue^2) + 0,data=df)
> summary(fit4)
##
## Call:
## lm(formula = profits ~ I(revenue^2) + 0, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1961.67 -468.16 -42.05 281.23 2459.01
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## I(revenue^2) 8.200e-05 7.658e-06 10.71 9.09e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1155 on 30 degrees of freedom
## Multiple R-squared: 0.7926, Adjusted R-squared: 0.7857
## F-statistic: 114.7 on 1 and 30 DF, p-value: 9.087e-12
1.7 多项式回归
在公式中调用poly(x,n)函数来对x的n阶多项式进行回归。
设定raw=TRUE是必要的，没有它，函数计算正交多项式而不是简单的多项式。
> fit5 <- lm(profits ~ poly(revenue,2,raw = TRUE) + 0,data=df)
> summary(fit5)
##
## Call:
## lm(formula = profits ~ poly(revenue, 2, raw = TRUE) + 0, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1886.84 -312.14 -260.31 -29.65 2306.10
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## poly(revenue, 2, raw = TRUE)1 5.737e-01 2.483e-01 2.310 0.0282 *
## poly(revenue, 2, raw = TRUE)2 -2.402e-05 4.644e-05 -0.517 0.6089
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1079 on 29 degrees of freedom
## Multiple R-squared: 0.8249, Adjusted R-squared: 0.8128
## F-statistic: 68.29 on 2 and 29 DF, p-value: 1.07e-11
1.8 转换数据的回归
x和y之间没有线性关系，需要转换后再进行回归。
> fit6 <- lm(log(profits) ~ log(revenue) + 0,data=df)
> summary(fit6)
##
## Call:
## lm(formula = log(profits) ~ log(revenue) + 0, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.41881 -0.01443 0.01931 0.11381 0.83064
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## log(revenue) 0.89044 0.01196 74.46 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.566 on 30 degrees of freedom
## Multiple R-squared: 0.9946, Adjusted R-squared: 0.9944
## F-statistic: 5545 on 1 and 30 DF, p-value: < 2.2e-16
转换后的修正R2得到了极大的提升，残差项也有了极大的降低。
我们画个图对比一下：
> p_load(ggplot2,patchwork)
> p1 <- ggplot(df,aes(revenue,profits)) +
+ geom_point() +
+ geom_smooth(method = "lm",se = F,col="#18CC97") +
+ labs(title = "Observed") +
+ theme_bw()
>
> p2 <- ggplot(df,aes(revenue,profits)) +
+ geom_point() +
+ geom_smooth(method = "lm",se = F,col="#18CC97") +
+ scale_y_log10() +
+ scale_x_log10() +
+ labs(y="",title = "Transformed") +
+ theme_bw()
>
> p1 + p2
变换前后对比
1.9 寻找最佳幂变换
通过MASS包中的boxcox函数，对响应变量应用一个幂变换来改进模型。这个函数会确定一个指数λ，然后将y转换为yλ。
> p_load(MASS)
> fit7 <- lm(profits ~ revenue + 0,data=df)
> plot(fit7,which = 1)
残差-拟合值图
残差-拟合值图略微有点抛物线，尝试对profits进行幂变换。
> bc <- boxcox(profits ~ revenue + 0,data=df)
找出lambda
函数没有返回λ的最佳值，但是可以通过（x，y）的值求出。
> which.max(bc$y)
## [1] 65
> lambda <- bc$x[which.max(bc$y)];print(lambda)
## [1] 0.5858586
然后对y应用这个幂变换后再次拟合模型：
> fit8 <- lm(I(profits^lambda) ~ revenue + 0,data=df)
> summary(fit8)
##
## Call:
## lm(formula = I(profits^lambda) ~ revenue + 0, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -56.972 -9.432 -2.471 4.878 50.816
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## revenue 0.0175948 0.0009431 18.66 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 26.6 on 30 degrees of freedom
## Multiple R-squared: 0.9207, Adjusted R-squared: 0.918
## F-statistic: 348.1 on 1 and 30 DF, p-value: < 2.2e-16
> summary(fit7)
##
## Call:
## lm(formula = profits ~ revenue + 0, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1920.89 -302.13 -231.13 -39.43 2324.90
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## revenue 0.4468 0.0378 11.82 8.13e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1066 on 30 degrees of freedom
## Multiple R-squared: 0.8232, Adjusted R-squared: 0.8173
## F-statistic: 139.7 on 1 and 30 DF, p-value: 8.127e-13
变换后的残差有极大的改善，R2也有比较大的改善。
但是Box-Cox的结果只能作为一个参考，并不是一个确切的答案。如果λ的置信区间包含1（上图中虚线区间），那么可能没有幂变换才是实际有用的。
1.10 回归系数的置信区间
综上可知，1.8中对x和y进行log变换时模型的性能最好。
要提取模型回归系数的置信区间，使用confint函数。可以通过level参数改变置信水平。
> confint(fit6,level = 0.95)
## 2.5 % 97.5 %
## log(revenue) 0.8660134 0.9148573
1.11 绘制回归残差
残差平滑线为一条恒为0的直线是最好的状况。
如果平滑线呈抛物线，说明模型缺少一个二次因子；如果呈圆锥形，可能表明y有非齐次性方差。
> plot(fit6,which = 1)
回归残差图
1.12 诊断线性回归
拟合线性回归是很容易的，然而这只是个开始，我们需要确定拟合的模型是否奏效并是否能达到很好的效果。
> p_load(car)
> outlierTest(fit6)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 4 -2.824377 0.0084804 0.26289
找到离群点4，即过度影响的观察值。
> plot(fit6)
残差-拟合值图
正态Q-Q图
大小-位置图
残差-影响图
比较好的回归诊断图有如下要求：
1.在残差-拟合值图中的点没有特别的模式，大多数呈随机分布；
2.在正态Q-Q图中的点基本落在线上，表明残差服从正态分布；
3.在大小-位置图和残差-影响图中，点以小组形式存在并且离中心不远。如果存在离中心很远的点，表明存在离群点，它们对回归有过多的影响。
1.13 识别有影响的观察值
当统计学家说一个观察值是有影响的时候，意味着移除这个观察值很大程度上会改变所拟合的回归模型，需要识别出这些观察值，因为它们可能是扭曲模型的离群值。
函数识别出的强影响值后面会有一个*标记，但我们不应该习惯性的剔除它们，需要做一些判断，因为我们不能确定这些值是改善了模型还是破坏了模型。
本例中有影响的观察值有4、19、28。
> influence.measures(fit6)
## Influence measures of
## lm(formula = log(profits) ~ log(revenue) + 0, data = df) :
##
## dfb.lg.. dffit cov.r cook.d hat inf
## 1 0.005915 0.005915 1.069 3.62e-05 0.0325
## 2 0.044515 0.044515 1.069 2.05e-03 0.0339
## 3 -0.086343 -0.086343 1.062 7.65e-03 0.0328
## 4 -0.488484 -0.488484 0.836 1.94e-01 0.0290 *
## 5 0.264729 0.264729 0.993 6.74e-02 0.0310
## 6 0.246284 0.246284 1.010 5.92e-02 0.0335
## 7 0.011851 0.011851 1.070 1.45e-04 0.0337
## 8 -0.001958 -0.001958 1.067 3.96e-06 0.0308
## 9 -0.002346 -0.002346 1.067 5.69e-06 0.0307
## 10 0.005002 0.005002 1.069 2.59e-05 0.0323
## 11 -0.096332 -0.096332 1.055 9.50e-03 0.0297
## 12 0.011650 0.011650 1.070 1.40e-04 0.0337
## 13 0.006255 0.006255 1.069 4.05e-05 0.0326
## 14 0.009359 0.009359 1.070 9.06e-05 0.0332
## 15 0.010848 0.010848 1.070 1.22e-04 0.0335
## 16 -0.000952 -0.000952 1.068 9.37e-07 0.0310
## 17 0.039828 0.039828 1.068 1.64e-03 0.0331
## 18 -0.083020 -0.083020 1.063 7.08e-03 0.0338
## 19 -0.495636 -0.495636 0.855 2.03e-01 0.0327 *
## 20 0.268513 0.268513 0.992 6.93e-02 0.0314
## 21 0.223798 0.223798 1.013 4.92e-02 0.0310
## 22 -0.006612 -0.006612 1.066 4.52e-05 0.0298
## 23 0.001513 0.001513 1.068 2.37e-06 0.0315
## 24 0.010520 0.010520 1.070 1.14e-04 0.0334
## 25 0.017884 0.017884 1.072 3.31e-04 0.0349
## 26 0.033940 0.033940 1.067 1.19e-03 0.0319
## 27 -0.090770 -0.090770 1.059 8.45e-03 0.0315
## 28 -0.496405 -0.496405 0.858 2.04e-01 0.0332 *
## 29 0.279790 0.279790 0.990 7.50e-02 0.0326
## 30 0.240755 0.240755 1.010 5.66e-02 0.0329
## 31 0.004576 0.004576 1.069 2.17e-05 0.0322
1.14 残差自相关检验
> p_load(lmtest)
> dwtest(fit6)
##
## Durbin-Watson test
##
## data: fit6
## DW = 2.1284, p-value = 0.6411
## alternative hypothesis: true autocorrelation is greater than 0
如果p<0.05，则残差显著相关，如果p>0.05，则表明没有自相关的证据。
1.15 预测新值
> preds <- data.frame(date=2020-02-01,revenue=5555,customers=44)
> predict(fit6,newdata = preds)
## 1
## 7.677737
即log(profits)=7.677737,profits=e^7.677737 ≈ 2.718^7.677737 = 2158.008。
1.16 建立预测区间
> predict(fit6,newdata = preds,interval = "prediction",level = 0.95)
## fit lwr upr
## 1 7.677737 6.502807 8.852668
lwr和upr分别为置信水平0.95时，区间的下限和上限。
1.17 运行单因素方差分析
> clas <- readr::read_csv("./data_set/class.csv")
> fit.clas <- lm(weight ~ height,data = clas)
> oneway.test(clas$height ~ clas$sex)
##
## One-way analysis of means (not assuming equal variances)
##
## data: clas$height and clas$sex
## F = 2.1063, num df = 1.000, denom df = 16.727, p-value = 0.1652
本例中将身高（height）按性别（sex）分成两组，然后检测组间均值是否显著不同。
p<0.05表明两个或多个组可能有显著不同的均值，但它不表明所有组有不同的均值；p>0.05表明没有足够的证据。
默认情况下该检验不假定方差相等（not assuming equal variances），如果知道这些组有相等的方差，可以设定参数var.equal=TRUE得到不保守的检验：
> oneway.test(clas$height ~ clas$sex,var.equal = TRUE)
##
## One-way analysis of means
##
## data: clas$height and clas$sex
## F = 2.1101, num df = 1, denom df = 17, p-value = 0.1645
1.18 创建交互关系图
方差分析是线性回归的一种形式，因此，理想的情况是每一个预测变量和响应变量之间有线性关系，非线性的来源之一是两个预测变量之间的交互关系：当一个预测变量的值发生变化时，其他预测变量和响应变量的关系也随之发生变化。
多因素方差分析：应用两个或多个分类变量作为预测变量的方差分析，需要得到预测变量之间可能的交互关系的可视化检验。
> p_load(faraway)
> data(rats)
> interaction.plot(rats$poison,rats$treat,rats$time)
image.png
上图本例中每条线绘制了time和poison的关系，各条线应该是平行的，但前两条线不完全平行，这预示着我们应该检验可能存在的交互关系。
1.19 找到组间均值的不同
1.17中可以检验组间均值是否显著不同，调用aov函数运行方差分析检验，可以找到均值间的差异。
> m <- aov(clas$height ~ clas$sex)
> TukeyHSD(m)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = clas$height ~ clas$sex)
##
## $`clas$sex`
## diff lwr upr p adj
## M-F 3.321111 -1.502522 8.144744 0.1645363
> plot(TukeyHSD(m))
image.png
p>0.05，另外置信区间跨越0值，表明差异不显著。
1.20 执行稳健方差分析
kruskal.test()函数不同于方差分析检验的地方是它不依靠数据的正态性，但它假设各组具有相同的中位数。
> kruskal.test(clas$height ~ clas$sex)
##
## Kruskal-Wallis rank sum test
##
## data: clas$height by clas$sex
## Kruskal-Wallis chi-squared = 1.9301, df = 1, p-value = 0.1648
p<0.05表示在两个或更多组的中位数之间有显著的差异。
1.21 运用方差分析比较模型
> anova(fit6,fit7)
## Analysis of Variance Table
##
## Response: log(profits)
## Df Sum Sq Mean Sq F value Pr(>F)
## log(revenue) 1 1776.19 1776.19 5544.6 < 2.2e-16 ***
## Residuals 30 9.61 0.32
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p<0.05表示模型是显著不同的。
