Statistical inference 1st

1 简介
2 方差分析

1 简介

本章介绍《数据统计与软件应用》课程的第四个核心部分，即统计推论（Statistical inference）。

统计推论是指根据样本数据对总体特征进行推断的过程，它是统计学的核心分支之一，主要包括参数估计和假设检验两部分。

参数估计是基于样本数据来推断总体参数的过程。常见的参数估计方法包括点估计和区间估计。

点估计是指用样本统计量估计总体参数值的方法，换句话说，总体样本中的任何一个统计量就是一个点估计量。

区间估计是指推测出总体参数可能取值范围的方法，例如，我们可以用置信区间来估计总体均值的范围。

假设检验是用来检验总体特征是否符合我们的期望或者假设的过程。在假设检验中，我们先提出一个关于总体特征的假设，然后根据样本数据来判断这个假设是否成立。通常我们会设定一个显著性水平，来决定拒绝还是接受假设。如果拒绝了假设，则意味着我们认为总体特征与假设不符；反之则说明我们没有足够的证据去拒绝这个假设。

通过统计推理，我们可以对总体特征进行推断和检验，从而更加准确地了解事物的本质。

方差分析和回归分析是统计推论中最广泛应用的统计技术之一，它们可以帮助我们进行参数估计和假设检验，从而对数据进行更深入的分析和解释。

方差分析（ANOVA）：是一种用于比较两个或多个样本均值是否存在显著差异的统计方法。它通过将总体方差分解为组间变异和组内变异，从而判断因素对总体均值是否产生显著影响。

回归分析（Regression）：是一种寻找解释变量与因变量之间关系的方法。它通过建立一个数学模型来描述自变量对因变量的影响，并利用样本数据对模型参数进行估计。

这章中我们将重点学习基于R语言的方差分析和回归分析。希望课程结束时，你们已经能够使用R软件、读懂R语言，并且可以结合我们前期学习的实验设计、数据探索性分析和统计推论这三大部分，将其灵活运用到你们即将到来的毕业论文设计和数据分析中，不要等到开始写论文了才发现自己的实验设计和数据收集漏洞百出，这门课的目的也就达到了。

2 方差分析

2.1 术语回顾

方差（Variance）: 衡量数据分布的离散程度。它描述的是一组数据与其平均值之间的偏离程度。

方差分析（Analysis of Variance，ANOVA）：单看术语容易被误导，实际上方差分析并不关心方差的分析，而是研究均值的变异，通过比较组间方差和组内方差，来测试多个组的平均值是否存在差异。

组内方差：同一组内各个样本之间的差异程度，它代表了组内变异，是通过将每个观察值与其组平均值进行比较来计算的。
组间方差：代表了组间变异，它是通过将每个组的平均值与数据的总体平均值进行比较来计算的。
总方差（平方和）：组间平方和+组内平方和。

单因素方差分析（One-way ANOVA）：独立样本t检验的扩展，用于在有两个以上组的情况下比较平均值。这是方差分析最简单的情况，其中数据仅根据一个分组变量（也称为因子变量）组织成多个组。其他同义词有：单向方差分析、单因素方差分析和受试者间方差分析。

双因素方差分析（Two-way ANOVA ）用于同时评估两个不同分组变量对连续结果变量的影响。其他同义词有：两个因子设计、双向方差分析或双向受试者间方差分析。

三因素方差分析（Three-way ANOVA ）用于同时评估三个不同分组变量对连续结果变量的影响。其他同义词有：阶乘方差分析、三向方差分析或三向受试者间方差分析。

残差（residuals）:是指实际观测值与模型预测值（例如回归模型或其他统计模型）之间的差异。在统计学和机器学习中，我们经常使用模型来预测观测值，而残差则代表了这些预测与实际观测之间的误差。如果模型能够完美地拟合数据，那么所有的残差都将接近于零。然而，在现实情况下，模型通常无法完美地拟合数据，因此残差是不可避免的，并且它们可以提供关于模型拟合程度和预测准确性的重要信息。

2.2 ANOVA机制

ANOVA背后的逻辑非常简单，即如果组间方差与组内方差的比值(F值)足够大，那么可以得出结，至少有一个组的平均值不等于其他组的平均值。

如下，方差分析表:

首先，你需要计算平方和（sum of squares, SS），它展示了组间和组内的方差。正如我们已经知道的那样，如果你将每个平方和除以它们各自的自由度（degree of freedom, df），你将得到组间和组内的均方（mean squares, MS）。F统计量就是组间均方与组内均方的比值。

演示：假设我们有 3 个组进行比较，如下图所示。虚线表示组平均值。下图图A显示了组间的变异，图B展示了组内变异。

ANOVA背后的数学过程：

1)计算组内方差(在ANOVA分析输出的结果中也叫残差（residuals）），这告诉我们每个观测值与他们所在组的平均值有多大不同（图B）。

2)计算组间均值方差（图A）。

3)计算F值，即组间方差/组内方差的比率。

4)然后根据F值与自由度 (df)，可计算得出p值

请注意，较低的 F 值 (F < 1) 表示所比较的样本均值之间没有显着差异。然而，较高的F值意味着组间变化与组内变化相比有大不同。根据F统计量的计算结果，我们可以查找F分布表或使用统计软件来获取对应的p值（aov函数会直接给出结果）。通过与设定的显著性水平（通常为0.05或0.01）比较p值，我们可以判断组间变异是否显著。

ANOVA的假设 经典的One-way ANOVA测试对数据做出以下前提假设：

假设1：因变量为连续变量。

假设2：解释变量包含2个及以上分组。

假设3：每组间的观测值相互独立。

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

假设5：每组内因变量符合正太分布。

假设6：每组的方差相等（方差齐性检测）

因此，在计算方差分析之前，我们需要执行一些初步测试来检查假设是否满足。

请注意，如果不满足上述假设，则可以使用单因素方差分析的非参数替代方案（ Kruskal-Wallis 检验）。

2.3 单因素方差分析

使用R数据集内置的PlantGrowth，它包含一个对照和两个不同处理条件下获得的植物的重量。

data("PlantGrowth")

levels(PlantGrowth$group)

## [1] "ctrl" "trt1" "trt2"

2.3.1 数据概况

查看数据集的基本情况

summary(PlantGrowth)

##      weight       group   
##  Min.   :3.590   ctrl:10  
##  1st Qu.:4.550   trt1:10  
##  Median :5.155   trt2:10  
##  Mean   :5.073            
##  3rd Qu.:5.530            
##  Max.   :6.310

计算每个组内weight变量的数据概况(count, mean and sd)

library(dplyr)
library(rstatix)
PlantGrowth %>%
  group_by(group) %>% #数据按照group变量进行分组，味着后续的操作将针对每个不同的group值进行。
  get_summary_stats(weight, type = "mean_sd") #weight变量的均值和标准差，并以mean_sd的形式返回结果

## # A tibble: 3 × 5
##   group variable     n  mean    sd
##   <fct> <fct>    <dbl> <dbl> <dbl>
## 1 ctrl  weight      10  5.03 0.583
## 2 trt1  weight      10  4.66 0.794
## 3 trt2  weight      10  5.53 0.443

2.3.2 检测异常值

boxplot识别离群点

library(ggpubr)
ggboxplot(PlantGrowth, x = "group", y = "weight")

使用箱线图方法可以轻松识别离群值，而R 函数 recognize_outliers() [rstatix 包] 中把离群点打印出来，并且识别极端异常值（is.extreme）。

PlantGrowth %>% 
  group_by(group) %>%
  identify_outliers(weight)

## # A tibble: 2 × 4
##   group weight is.outlier is.extreme
##   <fct>  <dbl> <lgl>      <lgl>     
## 1 trt1    5.87 TRUE       FALSE     
## 2 trt1    6.03 TRUE       FALSE

从上面的结果来看，图中的两个值虽然术语离群值，但不能算作极端异常值(is.extreme列给的结果是FALSE)。当然在上面的boxplot中，我们也可以得出相同的结论，因为那两个离群点分布在整体样本的变异范围内。

注意：如果出现极端异常值，如何处理？这部分我们在探索性数据分析章节已详细介绍过。异常值的原因可能是由于：数据输入错误、测量错误或异常值。排查异常值的原因后，做出相应的决策，即保留还是删除。如果你不认为结果会受到重大影响，则可以将异常值包含在分析中。这可以通过比较有和没有异常值的方差分析测试结果来评估。

2.3.3 正太分布假设检验

可以使用以下两种方法之一来检查正态性假设：

分析ANOVA模型残差以检查所有组的正态性。这种方法更简单，在你有很多组或每组数据点较少时非常方便。
分别检查每个组的正态性。当你只有少数几组但每组有很多数据点时，可以采用这种方法。

下面展示如何继续执行这两种方法

方法1：通过分析ANOVA模型残差来检查正态性假设

建立线性模型

model  <- lm(weight ~ group, data = PlantGrowth)

创建残差的 QQ-plot

ggqqplot(residuals(model))

Shapiro-Wilk 正态性检验

shapiro_test(residuals(model))

## # A tibble: 1 × 3
##   variable         statistic p.value
##   <chr>                <dbl>   <dbl>
## 1 residuals(model)     0.966   0.438

结果解读：在 QQ-plot中，由于所有点都大致沿着参考线分布，因此我们可以假设呈正态性。这一结论得到了 Shapiro-Wilk 检验的支持，p 值不显着 (p = 0.13)，因此我们可以假设数据集中的连续变量呈正态性。

方法2：按组检查正态性假设

计算每个组级别的 Shapiro-Wilk 检验。如果数据呈正态分布，则 p 值应大于 0.05。

代码1：

PlantGrowth %>%
  group_by(group) %>%
  shapiro_test(weight)

## # A tibble: 3 × 4
##   group variable statistic     p
##   <fct> <chr>        <dbl> <dbl>
## 1 ctrl  weight       0.957 0.747
## 2 trt1  weight       0.930 0.452
## 3 trt2  weight       0.941 0.564

代码2：

by(PlantGrowth$weight, PlantGrowth$group, shapiro.test)

## PlantGrowth$group: ctrl
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.95668, p-value = 0.7475
## 
## ------------------------------------------------------------ 
## PlantGrowth$group: trt1
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.93041, p-value = 0.4519
## 
## ------------------------------------------------------------ 
## PlantGrowth$group: trt2
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.94101, p-value = 0.5643

结果解读：根据 Shapiro-Wilk 正态性检验评估，每组的得分呈正态分布 (p > 0.05)。

为每个组级别创建 QQ 图：

ggqqplot(PlantGrowth, "weight", facet.by = "group")

结果解读：对于每个单元格，所有点都大致沿着参考线分布。所以我们可以假设数据呈正态性。

注意：如果你对数据的正态性有疑问，可以使用 Kruskal-Wallis 检验，这是单向方差分析检验的非参数替代方法。

2.3.4 方差齐性假设检验

残差与拟合图可用于检查方差的齐性（homogeneity）。

plot(model, 1)

在上图中，残差和拟合值（fitted values, 每组的平均值）之间没有明显的关系，这很好。因此，我们可以假设方差的齐性。

使用 Levene 检验

library(car)
leveneTest(weight ~ group, data=PlantGrowth)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  1.1192 0.3412
##       27

结果解读：从上面的输出中，我们可以看到 p 值 > 0.05，并不显著。这意味着，各组之间的方差不存在显著差异。因此，我们可以假设不同处组中的方差是齐性的。

注意：在不满足方差同质性假设的情况下，你可以使用函数 welch_anova_test()[rstatix 包] 计算 Welch one-way ANOVA test。该检验不需要假设方差相等。

2.3.5 ANOVA分析

方法1：使用aov（）

aov.model <-aov(weight ~ group, data=PlantGrowth)
summary(aov.model)

##             Df Sum Sq Mean Sq F value Pr(>F)  
## group        2  3.766  1.8832   4.846 0.0159 *
## Residuals   27 10.492  0.3886                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Df表示自由度（degrees of freedom），即数据集中独立变量的个数；
Sum Sq表示平方和（sum of squares），即每个因素的平方差之和；
Mean Sq表示均方（mean square），即平方和除以自由度；
F value表示F值，即组间均方与组内均方的比值；
Pr(>F)表示p值，即假设检验的显著性水平。例如，0.05的显著性水平表示在5%的情况下，拒绝相应的零假设。

根据这个结果，我们可以看到：
组间均方为1.8832，组内均方为0.3886；
F值为4.846，p值为0.0159；
星号(*)表示显著性水平小于0.05，即在95%的置信水平下，不同肥料处理组别对植物生长的影响是显著的。

注意：通过aov函数进行方差分析时，返回的结果中包含了残差（residuals），用于计算组内平方和。残差是实际观测值与回归模型预测值之间的差异，表示了因未考虑的因素而导致的变异。因此，在方差分析中，我们使用残差来表示组内变异。

方法2：使用lm（）

lm.model <- lm(weight ~ group, data=PlantGrowth)
summary(lm.model)

## 
## Call:
## lm(formula = weight ~ group, data = PlantGrowth)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0710 -0.4180 -0.0060  0.2627  1.3690 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.0320     0.1971  25.527   <2e-16 ***
## grouptrt1    -0.3710     0.2788  -1.331   0.1944    
## grouptrt2     0.4940     0.2788   1.772   0.0877 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6234 on 27 degrees of freedom
## Multiple R-squared:  0.2641, Adjusted R-squared:  0.2096 
## F-statistic: 4.846 on 2 and 27 DF,  p-value: 0.01591

实际上这是一个线性回归模型（lm）的结果输出，用于描述自变量（group）对因变量（weight）的影响。

关于线性回归，我们将在后面的章节详细介绍，这里我们主要介绍为什么我们可以用线性回归分析来代替方差分析。有的课件中会直接讲“单因素的方差分析”与“单因素的回归分析”是相同的。

从lm()回归分析结果的最后一行，我们可以看到其F值=4.846，p值=0.01591。该结果与aov()函数计算的F值和p值完全相同。因此，同样的，因为p<0.05，我们得出结论，即在95%的置信水平下，不同肥料处理组别对植物生长的影响是显著的，达到了方差分析的目的。

2.3.6 后续检验

如果使用aov()函数，我们只能得出结论三个组中至少有两个组是有显著差异的，但是无法得知两两之间是否有差异。

因此，aov()函数进行ANOVA分析后通常会随后进行 Tukey检验，以在组之间进行多个成对比较。

下面示例Tukey检验的两种函数：

函数1：TukeyHSD()

# 执行单因素方差分析
aov.model <-aov(weight ~ group, data=PlantGrowth)
#执行Tukey检验
TukeyHSD(aov.model)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = weight ~ group, data = PlantGrowth)
## 
## $group
##             diff        lwr       upr     p adj
## trt1-ctrl -0.371 -1.0622161 0.3202161 0.3908711
## trt2-ctrl  0.494 -0.1972161 1.1852161 0.1979960
## trt2-trt1  0.865  0.1737839 1.5562161 0.0120064

函数2：rstatix包的tukey_hsd()

library(rstatix)
tukey_hsd(aov.model)

## # A tibble: 3 × 9
##   term  group1 group2 null.value estimate conf.low conf.high p.adj p.adj.signif
## * <chr> <chr>  <chr>       <dbl>    <dbl>    <dbl>     <dbl> <dbl> <chr>       
## 1 group ctrl   trt1            0   -0.371   -1.06      0.320 0.391 ns          
## 2 group ctrl   trt2            0    0.494   -0.197     1.19  0.198 ns          
## 3 group trt1   trt2            0    0.865    0.174     1.56  0.012 *

#pwc <- PlantGrowth %>% tukey_hsd(weight ~ group, data=PlantGrowth)
#pwc

结果解读，比较两种函数：
diff = estimate：表示两组之间的均值差异估计值。例如，trt1组相对于ctrl组的均值差异估计值为-0.371。
lwr = conf.low：表示置信区间的下限。在95%置信水平下，均值差异估计值的下限。
upr = conf.high：表示置信区间的上限。在95%置信水平下，均值差异估计值的上限。
p adj:表示经过多重比较校正后的P值。如果该值小于显著性水平（通常为0.05），则表示两组之间的均值差异具有统计学显著性。
p.adj.signif: 根据调整后的P值得出的显著性水平标识，ns代表不显著。

介绍不同的函数，目的是让同学们了解：同一种分析使用不同的函数时，输出的结果的列名可能很不相同，但是换汤不换药，不要被迷惑，理解名称代表的含义，去解释结果即可。

2.3.7 如何在研究论文中报告结果

根据上面的示例，我们的目的是使用方差分析来评估 3 个不同处理组的植物生长是否不同：ctr (n = 10)、trt1 (n = 10) 和 trt2 (n = 10)。数据以平均值+/-标准差表示。

文字报告

ANOVA分析显示不同处理组之间植物生长有统计学显着差异，F(2, 27) = 4.85，p = 0.016。

计算组内均值和标准差显示：与 ctr 组（5.03 +/- 0.58）相比，trt1 组（4.66 +/- 0.79）植物生长下降；与trt1和ctr组相比，trt2组（5.53+/-0.44）的植物生长增加。

Tukey事后分析显示，从trt1到trt2的增加（0.87，95% CI（0.17至1.56））具有统计学显着性（p = 0.012），但其他组间差异没有统计学显着性。

图形报告

library(ggplot2)
library(rstatix)

#anova分析
aov.model <- PlantGrowth %>% anova_test(weight ~ group)
#tukey检验
tuk <- PlantGrowth %>% tukey_hsd(weight ~ group)
#可视化
tuk <- tuk %>% add_xy_position(x = "group")
ggboxplot(PlantGrowth, x = "group", y = "weight") +
  stat_pvalue_manual(tuk, hide.ns = TRUE) +
  labs(
    subtitle = get_test_label(aov.model, detailed = TRUE),
    caption = get_pwc_label(tuk)
    )

comparisons <- list(c("ctrl","trt1"),c("trt1","trt2"),c("ctrl","trt2"))
ggplot(PlantGrowth, aes(x=group, y=weight,color=group)) +
  geom_boxplot(notch=FALSE)+
  scale_fill_hue(l=40, c=35)+
  labs(title = "boxplot")+
  scale_x_discrete(limits=c("ctrl","trt1","trt2"))+
  geom_jitter(shape=16,position=position_jitter(0.2)) +
  stat_compare_means(comparisons = comparisons)

2.3.8 放宽方差齐性假设

标准的one-way ANOVA检验要求假设所有组的方差相等。在我们的示例中，方差齐性假设结果很好：Levene 检验并不显著。

假如数据违反方差齐性假设的情况下，我们如何进行方差分析测试？

在无法假设方差同质性的情况下（即 Levene 检验显着），Welch one-way test是标准 one-way ANOVA的一种替代方法。在这种情况下，Games-Howell 事后检验或成对 t 检验（不假设方差相等）可用于比较两两组差异的所有可能组合。

# Welch One way ANOVA test
res.aov2 <- PlantGrowth %>% welch_anova_test(weight ~ group)
# Pairwise comparisons (Games-Howell)
pwc2 <- PlantGrowth %>% games_howell_test(weight ~ group)
# Visualization: box plots with p-values
pwc2 <- pwc2 %>% add_xy_position(x = "group", step.increase = 1)
ggboxplot(PlantGrowth, x = "group", y = "weight") +
  stat_pvalue_manual(pwc2, hide.ns = TRUE) +
  labs(
    subtitle = get_test_label(res.aov2, detailed = TRUE),
    caption = get_pwc_label(pwc2)
    )

2.3.9 练习

使用iris数据集练习单因素方差分析，分别检验不同基因型（species）对四个变量的影响，最终形成一份可用于论文中的文字和图形报告。要求检验数据是否符合one-way anova的前提假设，然后根据结果选择后续的分析方法。

2.4 更多因素的方差分析

在本节中，我们只学习单因素方差分析，而不在此涉及更多因素的方差分析。因为更多因素的方差分析不如直接做一个多因素的回归分析，因为在做多因素的方差分析时，它只会给你一个overall的F值，这并不是我们想要的，而多因素的回归分析会单独给你每个因素下的t值和p值，更方便查看每个因素对因变量的影响是否显著，并且可以看两两因素间是否有交互作用，一般我们把这种交互作用称为ANCOVA (analysis of covariance)。ANCOVA（协方差分析），是一种结合方差分析（ANOVA）和回归分析的统计技术。本节内容作为回归分析的一个“引子”，在下一节中我将系统的介绍线性回归。