常用统计检验方法总结

基于R语言，整理统计学中常用的假设检验方法，包括基本原理、目的、适用、R语言代码等。

环境设定

# 加载包
library(ggplot2)
library(patchwork)
library(data.table)
library(magrittr)  

1. 数据类型

数据可以分为如下四种测量尺度（Stevens, 1946）：

1. 分类数据：

   1.1 定类数据：类别值、无顺序、数值本身无意义、无法比较、无法运算，例如，1表示男，2表示女

   1.2 定序数据：类别值、有数据、数值本身无意义、可比较、无法运算，例如，第1名、第2名

2. 数值数据：

   2.1 定距数据：数值、可加减、不可乘除、只有相对意义没有绝对意义、没有绝对0点，例如，温度

   2.2 定比数据：数值、可加减乘除、既有相对意义也有绝对意义

2. 统计分布

三大统计分布

2.1 正态分布

2.2 卡方分布

2.3 t分布

2.4 F分布

统计分布总结

在R语言中，统计分布基本都对应四个函数，例如对于正态分布，四个函数包括：

• dnorm: 概率密度函数PDF，正态分布的概率密度值，自变量是正无穷到负无穷的数x，因变量是x对应的正态分布概率密度，即\(P(X=x)\)
• pnorm: 累积分布函数CDF，正态分布的累积分布函数值，自变量是正无穷到负无穷的数x，因变量是x在正态分布上的累计概率密度，即\(P(X\leq x)\)
• qnorm: 分位点函数Quantile，正态分布的分位数，自变量是01区间的概率p，因变量是p的分位数值，即\(P(X\leq x)=p\)对应的x值，含义是，“在这个分布中，多大的值x，使得小于等于x的概率恰好是p”。本质上，pnorm和qnorm互为反函数
• rnorm: 随机数生成函数，生成满足正态分布的一系列随机数

前三个的函数曲线，以及rnorm生成的随机数（point）如下图所示：

# 散点数据
x <- rnorm(100)
dx <- dnorm(x)
px <- pnorm(x)

# 函数图像
p1 <- ggplot()+
  geom_function(fun=dnorm)+
  scale_x_continuous(limits = c(-10,10))+
  labs(x="x", y="概率密度", title="dnorm")+
  theme_light()
p2 <- ggplot()+
  geom_function(fun=dnorm)+
  annotate("point", x=x, y=dx, color="red")+
  scale_x_continuous(limits = c(-10,10))+
  labs(x="x", y="概率密度", title="dnorm + rnorm")+
  theme_light()
p3 <- ggplot()+
  geom_function(fun=pnorm)+
  scale_x_continuous(limits = c(-10,10))+
  labs(x="x", y="累计概率密度（分位点）", title="pnorm")+
  theme_light()
p4 <- ggplot()+
  geom_function(fun=pnorm)+
  annotate("point", x=x, y=px, color="red")+
  scale_x_continuous(limits = c(-10,10))+
  labs(x="x", y="累计概率密度（分位点）", title="pnorm + rnorm")+
  theme_light()
p5 <- ggplot()+
  geom_function(fun=qnorm)+
  scale_x_continuous(limits = c(0,1))+
  labs(x="累计密度（分位点）", y="x", title="qnorm")+
  theme_light()
p6 <- ggplot()+
  geom_function(fun=qnorm)+
  annotate("point", x=px, y=x, color="red")+
  scale_x_continuous(limits = c(0,1))+
  labs(x="累计密度（分位点）", y="x", title="qnorm + rnorm")+
  theme_light()

# 合并
(p1 + p2) / (p3 + p4) / (p5 + p6)

R语言中的分布函数，主要有以下几种（以r*函数为利）

连续分布

分布名	随机数生成函数	描述
正态分布	rnorm()	常见自然现象的分布
均匀分布	runif()	所有值等概率
t分布	rt()	小样本估计均值差异时常用
F分布	rf()	方差分析中常用
卡方分布	rchisq()	检验统计量的分布之一
指数分布	rexp()	表示事件发生的等待时间
伽马分布	rgamma()	等待多个事件的时间
beta分布	rbeta()	定义在 [0,1] 区间的概率
对数正态分布	rlnorm()	对数符合正态的分布
Cauchy分布	rcauchy()	重尾分布
Weibull分布	rweibull()	寿命分析常用
logistic分布	rlogis()	比正态分布胖尾
Pareto分布	VGAM::rpareto()	少数占多数的分布，如财富等

离散分布

分布名	随机数生成函数	描述
二项分布	rbinom()	成功/失败的次数
泊松分布	rpois()	单位时间内事件次数
几何分布	rgeom()	第一次成功前的失败次数
负二项分布	rnbinom()	成功r次前的失败次数
超几何分布	rhyper()	抽卡问题/不放回抽样
多项分布	rmultinom()	多类事件的分布

3. 基础统计量

常用统计量和计算方法包括：

集中趋势

• 均值：mean(x)
• 中位数：median(x)
• 众数：DescTools::Mode(x)

离散程度

• 极差：diff(range(x))，最大值与最小值之差
• 四分位差：IQR(x)，第3四分位数 - 第1四分位数
• 方差：var(x)
• 标准差：sd(x)
• 标准误：sd(x)/sqrt(n)
• 变异系数：sd(x)/mean(x)

分布形态

• 偏度：e1071::skewness()，分布的对称性（偏左或偏右）
• 峰度：e1071::kurtosis()，分布的“尖峭”程度（比正态分布更尖或更平）

位置指标

• 最值：min(x), max(x)
• 百分位数：quantile(x), quantile(x, 0.3)，默认5个分位数，或者自定义哪个分位数

统计量工具

• summary(x)，返回最小值、一三分位数、均值、中位数、最大值
• psych::describe(x)，返回多种描述统计量指标
• pastecs::stat.desc(x)，返回多种描述统计量指标
• skimr::skim(x)，返回多种描述统计量指标

注意

在R语言中，均值都是除以n，而方差、标准差都是样本统计量，即除以(n-1)

4. 统计检验类型

概述统计检验有哪些类别，以及其中包含哪些检验，而检验的具体细节，见第5节。

5. 统计检验方法

介绍常用统计检验方法的目的、适用性、计算公式、结果判定、案例（R语言方法）

注意

1. 统计检验的判定标准，主要有以下几种，因此对于假设检验结果判定，除了统计量，最重要的是计算p值和拒绝域临界值：

   • p值 + 显著性水平α： 若\(p \leq \alpha\)，则拒绝H0

   • 统计量 + 拒绝域： 若统计量落在拒绝域内，则拒绝H0

   • 数据 + 置信区间CI： 若数据落在CI内，则不拒绝H0（不常用）

2. 对于p值和α，可以理解为分别是统计量和拒绝域的累计概率（在密度函数中体现为下围面积），因此如果\(p \leq \alpha\)，说明统计量的面积比α的还小，因此统计量必然落在拒绝域内。或者将p值理解为一种显著性水平，并且是“能拒绝原假设的最小显著性水平”，如果显著性水平更小的话就不能拒绝原假设了，因此比p大的α可以拒绝原假设，比p小的α不能拒绝原假设。

3. 对于p值、拒绝域（或CI）的计算，与原假设和备择假设、统计分布特征有关。而根据定义，拒绝域（拒绝H0）的方向和H0的方向相反、和H1的方向相同：

   • 若H0是=，则H1是≠，拒绝域在两侧。若统计分布中心对称，则p值是大于等于统计量绝对值的概率的2倍，否则是大于等于统计量概率和小于等于统计量概率中，小的那个的2倍。

   • 若H0是≤，则H1是>，拒绝域在右侧，p值是大于等于统计量的概率。

   • 若H0是≥，则H1是<，拒绝域在左侧，p值是小于等于统计量的概率。

4. 对假设检验的结果，建议使用“拒绝原假设”或“接受原假设”进行评判。“通过检验”的说法并不严谨，因为有些检验我们渴望拒绝原假设（如回归系数t检验的原假设为\(\beta=0\)），有些时候我们又渴望接受原假设（如Shapiro-Wilk检验的原假设为数据服从正态分布），“通过”暗示检验结果符合意愿，但对于混淆了原假设是什么，也给p值判定方法造成困扰。

5.1 z检验

也叫u检验，注意并非U检验（非参数的Mann–Whitney U检验）

目的

• 检验单样本：样本均值、比例是否等于总体均值、比例
• 检验两独立样本：均值、比例是否相等

适用

• 检验类型：参数检验（需假设数据服从正态分布，且方差已知）
• 数据类型：数值数据（连续型）
• 数据体量：大样本（通常n≥30，或np≥10且n(1-p)≥10，或小样本但总体方差已知）

注意：

• z检验最为理想化，总体方差（或标准差）已知的条件太过苛刻，实际中，当样本量足够大时，将样本方差近似为总体方差时，才可用，一般情况可使用t检验替代。

• 比例数据实际上是二项分布，当np和n(1-p)较大时，二项分布的抽样趋近于正态分布，因此大样本下可以用z比例检验，也因此，小样本情况下比例数据的t检验就不适用了。

计算公式

1. 单样本均值z检验：\[z=\frac{\bar X - \mu_0}{\sigma / \sqrt n}\]

2. 单样本比例z检验：\[z=\frac{\hat p - p_0}{\sqrt{p_0(1-p_0)/n}}\]

3. 两独立样本均值z检验：\[z=\frac{\bar X_1 - \bar X_2}{\sqrt{\sigma^2_1/n_1 + \sigma^2_2/n_2}}\]

4. 两独立样本比例z检验：\[z=\frac{\hat p_1 - \hat p_2}{\sqrt{\hat p_p(1-\hat p_p)(1/n_1+1/n_2)}}, 其中\hat p_p=\frac{\hat p_1n_1+\hat p_2n_2}{n_1+n_2}\]

5. 单样本配对均值z检验，技术上可以，但理论上，配对差异的总体标准差不可能已知，因此建议用t检验。

上式中的符号含义如下：

• 样本的：\(\bar X\)，样本均值；\(n\)样本数量；\(\hat p\)，样本比例。

• 总体的：\(\mu_0\)，总体均值；\(\sigma\)，总体标准差；\(p_0\)，总体比例。

结果判定

统计量z服从标准正态分布，该分布关于0点左右对称。设定显著性水平为α后，检验结果判定步骤如下：

1. 若H0为=，则双侧检验，\(p=2P(X \geq |z|)\)，拒绝域为\(|z| \geq z_{1-\alpha/2}\) （由于对称，只取右侧）

2. 若H0为≤，则右侧检验，\(p=P(X \geq z)\)，拒绝域为\(z \geq z_{1-\alpha}\)

3. 若H0为≥，则左侧检验，\(p=P(X \leq z)\)，拒绝域为\(z \leq z_{\alpha}\)

注意：

1. z服从标准正态分布，密度函数关于X=0左右对称，因此，\(|z|\)位于密度函数x轴的右侧

2. \(z_{\alpha}\)表示自左至右累计密度为\(\alpha\)时对应的\(z\)值，即等式\(\alpha=P(X\leq z)\)的解。由于设定的\(\alpha\)不超过0.5，因此，\(z_{\alpha}\)位于密度函数x轴的左侧，与\(z_{1-\alpha}\)关于x=0对称，\(z_{1-\alpha/2}\)位于右侧，大于\(z_{1-\alpha/2}\)的面积乘以2等于\(\alpha\)

3. 一个结果判定的例子：经计算\(z=2.12\)，另\(\alpha = 0.05\)，此时：

   • 若H0为=，用双侧检验：p值为\(p=2P(X \geq |2.12|)=0.034\)，p值在R语言中可用2*(1-pnorm(2.12))或2*(pnorm(-2.12))计算得出。拒绝域临界值为\(z_{1-\alpha /2}\)，对应公式为\(P(X \geq z_{1-\alpha/2}) = 0.025\)，可根据R语言qnorm(1-0.05/2)得出\(z_{0.975}=1.960\)。因此，可根据0.034<0.05或2.12>1.960，拒绝H0，认为≠。

   • 若H0为≤，用右侧检验：p值为\(p=P(X \geq 2.12)=0.017\)（即1-pnorm(2.12)），临界值为\(z_{1-\alpha} = z_{0.95}=1.645\)（即qnorm(1-0.05)）。因此拒绝H0，认为＞。

   • 若H0为≥，用左侧检验：p值为\(p=P(X \leq 2.12)=0.983\)（即pnorm(2.12)），临界值为\(z_{\alpha} = z_{0.05}=-1.645\)（即qnorm(0.05)）。因此不拒绝H0，认为≥。

4. 上述判断过程如下图：

案例

z检验只有理论价值，实际不常用，此处案例较为简单，并没有使用检验函数，都是手算的。

1. 单样本均值检验

已知某品牌灯泡寿命的总体均值为1000小时，标准差100小时。为检测质量抽取50个，测的平均寿命为990小时，问是否宣传造假。

计算步骤如下（z计算公式：\(z=\frac{\bar X - \mu_0}{\sigma / \sqrt n}\)）：

# 定义H0：样本均值大于等于总体均值，H1：样本均值小于总体均值，因此用左侧检验

# 设定显著性水平α=0.05
a <- 0.05

# 计算z统计量
z <- (990-1000)/(100/sqrt(50)) #z = -0.7071068

# 定义α=0.05，则对应单（左）侧判定参数：
a <- 0.05
za <- qnorm(0.05)            #zα = -1.644854
p <- pnorm(z)                #p =  0.2397501

# 结论
# 由于z = -0.7071068 > za = -1.644854  （左侧检验）
# 由于p = 0.2397501 > α = 0.05         （p检验不分左右）
# 因此不能决绝H0，接受H0，认为样本均值大于等于总体均值，即没有造假

2. 单样本比例z检验：

某电商宣称订单取消率为5%，现观测500个订单，取消率为6.5%，检验是订单取消率否显著增加。

计算步骤如下（z计算公式：\(z=\frac{\hat p - p_0}{\sqrt{p_0(1-p_0)/n}}\)）：

# 定义H0：取消率减少或没变，H1：取消率增加了，用右侧检验

# 设定显著性水平α=0.05
a <- 0.05

# 计算统计量
z <- (0.065-0.05)/sqrt(0.05*(1-0.05)/500) #z = 1.538968

# 定义α=0.05，则对应单（右）侧判定参数：
a <- 0.05
z1a <- qnorm(1-0.05)        #z1α = 1.644854    #对右侧检验，qnorm是pnorm的逆，1-0.05才是右侧面积
p <- 1-pnorm(z)             #p =  0.06190611  #对右侧检验，pnorm是累计概率，1-pnorm才是右侧面积

# 结论
# 由于z = 1.538968 < z1a = 1.644854  （右侧检验）
# 由于p = 0.06190611 > α = 0.05     （p检验不分左右）
# 因此不能决绝H0，接受H0，认为订单取消量没有显著增加

3. 两独立样本均值z检验

A工厂的产品性能方差为2，B工厂产品为3，现在抽取A产品40个测得性能均值为50，B产品50个测得均值为48，问两厂产品性能均值是否相同。

计算步骤如下（z计算公式：\(z=\frac{\bar X_1 - \bar X_2}{\sqrt{\sigma^2_1/n_1 + \sigma^2_2/n_2}}\)）：

# 定义H0：性能相同，H1：性能不同，用双侧检验

# 设定显著性水平α=0.05
a <- 0.05

# 计算统计量
z <- (50-48)/sqrt(2^2/40+3^2/50) #z = 3.779645

# 定义α=0.05，则双侧判定参数：
a <- 0.05
z1a2 <- qnorm(1-0.05/2)         #z1α2 = 1.959964    #双侧检验即是右侧α除以2
p <- (1-pnorm(abs(z)))*2        #p =  0.0001570523 #这里要注意z的符号，双侧最好取绝对值

# 结论
# 由于z = 3.779645 > z1a2 = 1.959964    （双侧检验，取右侧进行判断）
# 由于p = 0.0001570523 < α = 0.05       （p检验不分左右）
# 因此决绝H0，接受H1，认为两厂产品新能不同

4. 两独立样本比例z检验

对A产品广告，抽取500次点击率为12%，B产品广告抽取500次点击率为8%，判断两广告点击率是否有区别

计算步骤如下（z计算公式：\(z=\frac{\hat p_1 - \hat p_2}{\sqrt{\hat p_p(1-\hat p_p)(1/n_1+1/n_2)}}, 其中\hat p_p=\frac{\hat p_1n_1+\hat p_2n_2}{n_1+n_2}\)）：

# 定义H0：点击率相同，H1：点击率不同，用双侧检验

# 设定显著性水平α=0.05
a <- 0.05

# 计算统计量
pp <- (0.12*500+0.08*500)/(500+500)
z <- (0.12-0.08)/sqrt(pp*(1-pp)*(1/500+1/500))  #z = 2.108185

# 定义α=0.05，则双侧判定参数：
a <- 0.05
z1a2 <- qnorm(1-0.05/2)         #z1α2 = 1.959964    #双侧检验即是右侧除以2
p <- (1-pnorm(abs(z)))*2        #p =  0.03501498    #这里要注意z的符号，双侧最好取绝对值

# 结论
# 由于z = 2.108185 > z1a2 = 1.959964      （双侧检验，取右侧进行判断）
# 由于p = 0.03501498 < a = 0.05           （p检验不分左右）
# 因此决绝H0，接受H1，认为广告点击率不同

5.2 t检验

目的

• 检验单样本：均值是否等于总体均值
• 检验两独立样本：均值是否相等
• 检验配对样本：均值是否相等

适用

• 检验类型：参数检验（需假设数据正态或近似正态，方差齐性）
• 数据类型：数值数据（连续型）
• 数据体量：小样本（n<30 时常用），但大样本也可用（此时接近z检验）

注意：

• 比例检验服从二项分布，大样本下服从正态分布，因此适用于z检验，而不适用于t检验

计算公式

1. 单样本t检验：\[z=\frac{\bar X - \mu_0}{s / \sqrt n}, df=n-1\]

2. 方差齐次的独立样本t检验：\[z=\frac{\bar X_1 - \bar X_2}{\sqrt{s^2_p(1/n_1+1/n_2)}}, df=n_1+n_2-2\]

其中合并方差为\[s^2_p=\frac{(n_1-1)s_1^2+(n_2-1)s_2^2}{n_1+n_2-2}\]

3. 方差不齐次的独立样本t检验（Welch’s t test）：\[z=\frac{\bar X_1 - \bar X_2}{\sqrt{(s_1^2/n_1+s_2^2/n_2)}}\]

\[df = \frac{(s_1^2/n_1+s_2^2/n_2)^2}{\frac{(s_1^2/n_1)^2}{n_1-1}+\frac{(s_2^2/n_2)^2}{n_2-1}}\]

4. 配对样本t检验：\[z=\frac{\bar D}{s_D / \sqrt n}, df=n-1\]

上式中的符号含义如下：

• \(\bar X\)，样本均值；\(\mu_0\)，假设的总体均值；\(s\)，样本标准差；\(n\)，样本量；
• \(D\)，配对数据差，例如\(D_i=X_{i1}-X_{i2}\)；\(\bar D\)，配对差的均值；\(s_D\)，配对差的样本标准差；

注意：

• Welch’s t检验，在方差齐次时仍然是有效的，只是效率较低，尤其是小样本时，因此它更为稳健，是R中的默认方法，也因此，在实操中，进行独立样本t检验前不必进行方差齐次性检验，直接用welch’s t检验就行了

结果判定

统计量t服从t分布（t符号省略自由度），该分布关于0点左右对称。设定显著性水平为α后，检验结果判定步骤如下：

1. 若H0为=，则双侧检验，\(p=2P(X \geq |t|)\)，拒绝域为\(|t| \geq t_{1-\alpha/2}\)

2. 若H0为≤，则右侧检验，\(p=P(X \geq t)\)，拒绝域为\(t \geq t_{1-\alpha}\)

3. 若H0为≥，则左侧检验，\(p=P(X \leq t)\)，拒绝域为\(t \leq t_{\alpha}\)

4. 计显著性水平\(\alpha\)下的置信区间CI，但CI是针对数据的HO的“接受域”。而CI的计算方法与统计量形式有关，其临界值本质与上面的\(t_{\alpha}\)一样，只不过是将判断对象从统计量转换到数据上（通常，单样本对应样本均值\(\bar X\)，或者独立样本对应t统计量的分子部分\(\bar X_1 - \bar X_2\)）。例如对于单样本t检验，其t统计量公式为\(z=\frac{\bar X - \mu_0}{s / \sqrt n}\)，对于双侧检验，t拒绝域为\(|t| > t_{1-\alpha/2}\)，而CI范围是\(\bar X \pm t_{1-\alpha/2} \frac {s} {\sqrt n}\)。

案例

R语言中，使用stats::t.test和rstatix::t_test执行t检验

对于stats::t.test，主要参数和含义如下：

• x: 第一组数据（数值向量或公式）
• y: 第二组数据（第二个样本，如果是独立样本检验的话）
• alternative: 备择假设的方向，包括”two.sided”（双侧），“less”（左侧）和”greater”（右侧）
  
• mu: 零假设下的比较值（例如单样本检验中μ₀）
• paired: 是否为配对样本（TRUE/ FALSE）
• var.equal: 如果是两个独立样本，是否假设方差相等（即方差齐性）
• conf.level: 置信区间CI的显著性水平（默认 0.95）

对于rstatix::t_test，主要参数和含义如下：

• data: 数据框，包含你要比较的变量和分组因子
• formula: 比较公式，例如 value ~ group，如果是单样本t检验，可设置为value ~ 1
• comparisons: 明确指定要比较的组对，比如 list(c(“A”, “B”), c(“A”, “C”))
• ref.group: 指定参考组，和所有其他组两两比较
• p.adjust.method: 多重比较校正方法，如 “holm”、“bonferroni”、“BH” 等
• paired: 是否是配对样本
• var.equal: 是否假设两个组方差相等（用于独立样本）
• alternative: 备择假设方向，“two.sided”、“less”、“greater”
• mu: 零假设的均值差值，默认 0（比如 H0: μ1 - μ2 = 0）
• conf.level: 置信区间水平，默认 0.95
• detailed: 是否返回更详细结果（含CI、方法、备择假设方向等）

注意:

• t.test的结果更为直观，适合单次计算后直接下定论，t_test的结果为tibble格式，给出适用的诊断数据，适合进一步绘图或处理。

• t.test由于使用向量传递参数，所以两个总体时顺序是自定义的，而t_test使用data.frame传递参数，所以默认使用字母顺序排序。在某些时候（如下面的案例3中）他们的原假设是不同的，应根据结果注意区分。

1. 单样本t检验

抽样一个班的学生身高分别是168, 172, 175, 169, 171, 159, 180，问班级总体平均值是否等于180 （或小于180、或大于180）

stats::t.test方法：

# 设定显著性水平α=0.05
a <- 0.05

# 身高样本
height <- c(168, 172, 175, 169, 171, 159, 180)

# 单样本t检验，H0：等于180，双侧检验
res.1 <- t.test(height, mu = 180, alternative = "two.sided")
# 单样本t检验，H0：小于等于180，右侧检验（注意alternative是H1）
res.2 <- t.test(height, mu = 180, alternative = "greater")
# 单样本t检验，H0：大于等于180，左侧检验（注意alternative是H1）
res.3 <- t.test(height, mu = 180, alternative = "less")

# 结果1：说明0.05时，拒绝H0，认为总体均值不等于180
#   One Sample t-test                                    #检验类型，单样本t检验     
# data:  height                                          #数据名称：height
# t = -3.8362, df = 6, p-value = 0.008597                #统计量相关指标计算结果
# alternative hypothesis: true mean is not equal to 180  #备择假设H1
# 95 percent confidence interval:                        #置信区间CI95%，H0的“接受域”
#  164.5574 176.5855
# sample estimates:                                      #样本均值，即mean(height)
# mean of x 
#  170.5714 

# 结果2：说明0.05时，不拒绝H0，认为总体均值小于等于180
#   One Sample t-test
# data:  height
# t = -3.8362, df = 6, p-value = 0.9957
# alternative hypothesis: true mean is greater than 180
# 95 percent confidence interval:                        #数据180在CI中，所以接受H0
#  165.7955      Inf
# sample estimates:
# mean of x 
#  170.5714 

# 结果3：说明0.05时，拒绝H0，认为总体均值小于180
#   One Sample t-test
# data:  height
# t = -3.8362, df = 6, p-value = 0.004299
# alternative hypothesis: true mean is less than 180
# 95 percent confidence interval:                        #数据180没在CI中，所以不接受H0
#      -Inf 175.3474
# sample estimates:
# mean of x 
#  170.5714 

rstatix::t_test方法：

# 加载包 
library(rstatix)
library(magrittr)

# 设定显著性水平α=0.05
a <- 0.05

# 身高样本
dat <- data.frame(height = c(168, 172, 175, 169, 171, 159, 180))

# 单样本t检验，H0：等于180，双侧检验
res_1 <- dat %>% t_test(height~1, mu = 180)
# 单样本t检验，H0：小于180，右侧检验（注意alternative是H1）
res_2 <- dat %>% t_test(height~1, mu = 180, alternative = "greater")
# 单样本t检验，H0：大于180，左侧检验（注意alternative是H1）
res_3 <- dat %>% t_test(height~1, mu = 180, alternative = "less", detailed = T)

# 结果1：
# A tibble: 1 × 7
#   .y.    group1 group2         n statistic    df      p
# * <chr>  <chr>  <chr>      <int>     <dbl> <dbl>  <dbl>
# 1 height 1      null model     7     -3.84     6 0.0086

# 结果2：
# A tibble: 1 × 7
#   .y.    group1 group2         n statistic    df     p
# * <chr>  <chr>  <chr>      <int>     <dbl> <dbl> <dbl>
# 1 height 1      null model     7     -3.84     6 0.996

# 结果3：
# A tibble: 1 × 12
# estimate .y.    group1 group2   n statistic      p    df conf.low conf.high method alternative
# *    <dbl> <chr>  <chr>  <chr>  <int>     <dbl>  <dbl> <dbl>    <dbl>     <dbl> <chr>  <chr>  
# 1     171. height 1    null model     7     -3.84 0.0043     6     -Inf      175. T-test less 

2. 独立样本t检验

比较两种教学方法下学生的成绩(78, 85, 82, 98)和(75, 80, 77, 86, 92)是否无差异

# 加载包 
library(rstatix)

# 设定显著性水平α=0.05
a <- 0.05

# 数据
method_A <- c(78, 85, 82, 98)
method_B <- c(75, 80, 77, 86, 92)
dat <- data.frame(value=c(method_A, method_B),
                  type=c(rep("A", length(method_A)),rep("B", length(method_B))))

# 检验1，按方差齐次算
res.1 <- t.test(method_A, method_B, var.equal = F)

# 检验2，按方差非齐次计算（Welch's t test）
res_2 <- dat %>% t_test(value~type, var.equal = T)

# res.1结果：说明在0.05下，不拒绝H0，认为两总体均值相等
# Two Sample t-test
# data:  method_A and method_B
# t = 0.72276, df = 7, p-value = 0.4933
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:               #CI的意思是barX1-barX2在这个范围内接受H0
#  -8.51865 16.01865 
# sample estimates:
# mean of x mean of y 
#     85.75     82.00 

# res_2结果：说明在0.05下，不拒绝H0，
# A tibble: 1 × 8
#   .y.   group1 group2    n1    n2 statistic    df     p
# * <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl>
# 1 value A      B          4     5     0.723     7 0.493

3. 配对样本t检验

5名患者治疗前的血压为(120, 125, 130, 123, 121)，治疗完变为(115, 108, 111, 107, 112)，问是否起效果

注意：

下面两个检验函数，都设置了左侧检验，即原假设为H0:D≥0，H1:D<0。但问题是t.test和t_test对bef和aft的顺序定义是不同的。t.test由于使用向量传递参数，所以“D=bef-aft”，而t_test使用data.frame传递参数，所以默认使用字母顺序排序，因此“D=aft-bef”。因此当定义H1为“less”时，t.test的H0是“bef-aft≥0”，而t_test的H0是“aft-bef≥0”。因此在本例中，虽然alternative相同、p值是否拒绝H0不同，但是结论是相同的。

# 加载包 
library(rstatix)

# 设定显著性水平α=0.05
a <- 0.05

# 数据
bef <- c(120, 125, 130, 123, 121)
aft <- c(115, 108, 111, 107, 112)
dat <- data.frame(value=c(bef,aft),
                  type=c(rep("bef",length(bef)), rep("aft",length(aft))))

# stats包：H0：bef ≥ aft，用左侧检验
test.1 <- t.test(bef, aft, paired = TRUE, alternative = "less")

# rstatix包：H0：aft ≥ bef，用左侧检验
test_2 <- dat %>% t_test(value~type, paired = T, alternative = "less")

# res.1结果：说明0.05时，不拒绝H0，认为bef ≥ aft
#   Paired t-test
# 
# data:  bef and aft
# t = 4.9749, df = 4, p-value = 0.9962
# alternative hypothesis: true mean difference is less than 0
# 95 percent confidence interval:
#      -Inf 18.85643
# sample estimates:
# mean difference 
#            13.2 

# res_2结果，说明0.05时，拒绝H0，认为aft < bef
# A tibble: 1 × 8
#   .y.   group1 group2    n1    n2 statistic    df       p
# * <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl>   <dbl>
# 1 value aft    bef        5     5     -4.97     4 0.00381 

5.3 F检验

目的

F检验常用以下三种情景；

• 检验两独立样本：总体方差是否相等（方差齐次性检验）
• 检验三组或以上样本：总体均值是否相等（方差分析）
• 检验所有自变量：回归系数是否全部为 0（多元回归方程线性性检验）

适用

• 检验类型：参数检验（用于方差分析或方差齐性检验，需正态分布假设）
• 数据类型：数值数据（连续型）
• 数据体量：小样本或大样本均可，但需满足正态性和方差齐性

计算公式

1. 方差齐次性检验

   • 统计量：\(F=\frac{s^2_1}{s^2_2}\)

   • 自由度：分子自由度\(df_1=n_1-1\)，分母自由度\(df_2=n_2-1\)

   • 注意：

     • 这种用法中，为了手算好查表，会要求分子样本方差要大于分母样本方差，即\(s^2_1>s^2_1\)，因此有些教材中，也定义这种情况下F统计量计算公式为\(F=\frac{max(s^2_1, s^2_2)}{min(s^2_1,s^2_2)}\)
     • F统计量计算时，分子分母都要除以自由度，只不过在这里，样本标准差\(s=\frac {\sum (x_i-\bar x)^2}{n-1}\)，已经包含除以自由度的n-1了。

2. 方差分析

   • 统计量：\(F=\frac{MSB}{MSW}=\frac{SSB/(k-1)}{SSW/(N-k)}\)

   • 组间均方（Mean Square Between, MSB）：\(MSB=\frac{\sum^{k}_{i=1}{n_i(\bar X_i-\bar X)^2}}{k-1}\)

   • 组内均方（Mean Square Within, MSW）：\(MSW=\frac{\sum_{i=1}^{k}{\sum_{j=1}^{n_i}{(X_{ij}-\bar X_i)^2}}}{N-k}\)

   • 自由度：组间（分子）\(df_1=k-1\)，组内（分母）\(df_2=N-k\)

3. 回归分析

   • 统计量：\(F=\frac{SSR/k}{SSE/(n-k-1)}\)

   • 回归平方和SSR（Regression Sum of Squares）或ESS（Explained Sum of Squares），即预测值与样本均值的差异：\(SSR=\sum_{i=1}^{n}{(\hat y_i-\bar y)^2}\)

   • 残差平方和SSE（Error Sum of Squares）或RSS（Residual Sum of Squares），即实际值与预测值的差异：\(SSE=\sum_{i=1}^{n}{(y_i-\hat y_i)^2}\)

   • 总平方和SST（Total Sum of Squares）或SST（Total Sum of Squares），即实际值与样本均值的差异：\(SST=SSR+SSE=\sum_{i=1}^{n}{(y_i-\bar y)^2}\)

   • 自由度：分子（回归）\(df_1=k\)，分母（残差）\(df_2=n-k-1\)

上式中的符号含义如下：

• \(s\)，样本方差

• \(k\)，组数：\(n_i\)，第i组的样本量；\(N=\sum n_i\)，总样本量；\(\bar X_i\)，第i组的样本均值，\(\bar X\)，总体均值

• \(y_i\)，第i个实际观测值；\(\hat y_i\)，第i个预测值；\(\bar y\)，样本均值，\(n\)，样本量

注意：

• 计算MSB时，乘以\(n_i\)是为了反映每组在总体中的代表性或权重，即组的样本量越大，它的均值对整体的贡献越大，偏离总体均值所造成的变异也越大。

结果判定

统计量F服从F(df1,df2)分布，该分布左右不对称。因此，设定显著性水平为α后，检验结果判定步骤如下：

1. 方差齐次性检验

   • 若H0为=，则双侧检验，\(p=2 \times min(P(X \geq F), P(X \leq F))\)，拒绝域为\(F \leq F_{\alpha /2}(n_1-1,n_2-1)\) 或 \(F \geq F_{1-\alpha/2}(n_1-1,n_2-1)\)

   • 若H0为≤，则右侧检验，\(p=P(X \geq F)\)，拒绝域为\(F \geq F_{1-\alpha}(n_1-1,n_2-1)\)

   • 若H0为≥，则左侧检验，\(p=P(X \leq F)\)，拒绝域为\(F \leq F_{\alpha}(n_1-1,n_2-1)\)

2. 方差分析

   • H0: \(\mu_1=\mu_2=...=\mu_k\) vs. H1：\(\mu_i\)不全相等，即双侧检验

   • \(p=2 \times min(P(X \geq F), P(X \leq F))\)，拒绝域为\(F \leq F_{\alpha /2}(k-1,N-k)\) 或 \(F \geq F_{1-\alpha/2}(k-1,N-k)\)

3. 多元回归方程线性性检验

   • H0: \(\beta_1=\beta_2=...=\beta_k=0\) vs. H1：\(\beta_i\)不全为0，即双侧检验

   • \(p=2 \times min(P(X \geq F), P(X \leq F))\)，拒绝域为：\(F \leq F_{\alpha /2}(k,n-k-1)\) 或 \(F \geq F_{1-\alpha/2}(k,n-k-1)\)

注意：

• F分布左右不对称，因此在双侧检验计算统计量p值时，需要判定F位于哪侧，然后将小的那一侧面积乘以2，否则另一侧乘以2结果将大于1。

案例

R语言中，使用stats::var.test执行方差齐次性检验，其他两种情况在方差分析和多元回归内部。对于var.test，参数都包括：

• x：第一个样本的数值向量
• y：第二个样本的数值向量
• ratio：原假设中两个总体方差的比值，默认是 1（即检验是否相等）
• alternative：备择假设的方向，包括two.sided、less、greater
• conf.level：置信区间水平，默认是 95%（即 0.95）

1. 方差齐次性检验中的F检验

# 设定显著性水平α=0.05
a <- 0.05

# 生成示例数据
set.seed(123)
group1 <- rnorm(30, mean = 100, sd = 10)  # 标准差=10
group2 <- rnorm(30, mean = 100, sd = 15)  # 标准差=15

# 函数计算..............................................................................
# F检验要求数据正态分布，这里进行正态性检验，其中shapiro.test原假设正态分布。
# shapiro.test(group1)  #p > 0.05 说明正态
# shapiro.test(group2)  #p > 0.05 说明正态

# F检验（方差齐性检验）
test.1 <- var.test(group1, group2, alternative = "two.side")

# 结果：p<a，说明不拒绝H0，认为两者方差相等
#   F test to compare two variances
# data:  group1 and group2
# F = 0.61331, num df = 29, denom df = 29, p-value = 0.194
# alternative hypothesis: true ratio of variances is not equal to 1
# 95 percent confidence interval:
#  0.2919117 1.2885518
# sample estimates:
# ratio of variances 
#          0.6133053 

# 手工计算..............................................................................
alpha <- 0.05
df1 <- length(group1)-1
df2 <- length(group2)-1

F_value <- var(group1)/var(group2)            #结果为0.6133053，注意这里不要占用关键字F，会影响画图

# 双侧检验的p值，没有小于α，因此不拒绝H0
pvalue <- 2*min(pf(F_value, df1, df2), 1-pf(F_value, df1, df2))  #结果为0.1940367

# 统计量临界值f=0.6133053，没有落在拒绝域内，因此不拒绝H0
Flow <- qf(alpha / 2, df1, df2)         #左临界值0.4759648
Fupp <- qf(1 - alpha / 2, df1, df2)     #右临界值2.100996

2. 方差分析中的F检验

# 加载包
library(rstatix)

# 设定显著性水平α=0.05
a <- 0.05

# 生成四组数据
set.seed(123)
groupA <- rnorm(30, mean = 50, sd = 5)
groupB <- rnorm(30, mean = 55, sd = 5)
groupC <- rnorm(30, mean = 60, sd = 5)
groupD <- rnorm(30, mean = 64, sd = 5)

# 数据框整理
dat <- data.frame(
  value = c(groupA, groupB, groupC, groupD),
  group = factor(rep(c("A", "B", "C", "D"), each = 30))
)

# 函数计算..............................................................................

# 1. 正态性检验（每组），shapiro.test检验原假设正态分布，这里全不拒绝H0
normtest <- tapply(dat$value, dat$group, shapiro.test)

# 2. 方差齐性检验（Levene检验），Levene检验原假设方差齐次，这里全不拒绝H0
homogvar <- levene_test(value ~ group, data = dat)

# 3. 满足条件，进行ANOVA
model <- aov(value ~ group, data = dat)
# summary(model)

# 结果如下（summary的有近似，直接print(model)是精确数值）
# > summary(model)
#              Df Sum Sq Mean Sq F value Pr(>F)    
# group         3   3166  1055.5   52.14 <2e-16 ***
# Residuals   116   2348    20.2                   
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# 对结果的解释：
# 1. 自由度表明，组数为4，样本总数为120（3=4-1,116=120-4）
# 2. Residuals表示残差，在ANOVA中，由于要验证组间差异，因此组内的差异就是残差
# 3. 组间：平方和Sum Sq为3166.455，除以自由度DF，等于Mean Sq1055.5
# 4. 组内：平方和为2348.143，自由度116，均方为20.2
# 5. F统计量为52.14，p值小于2e-16，说明在0.05显著性水平下，拒绝原假设，认为各组均值不全相等
# 6. 至此可以判断至少有一组均值不等，但具体是哪一组，需进一步用TukeyHSD(model)进行检查，结果如下
# 
#   Tukey multiple comparisons of means
#     95% family-wise confidence level
# Fit: aov(formula = value ~ group, data = dat)
# $group
#          diff        lwr       upr     p adj
# B-A  6.127210  3.0990917  9.155329 0.0000037
# C-A 10.357621  7.3295020 13.385739 0.0000000
# D-A 13.766074 10.7379554 16.794193 0.0000000
# C-B  4.230410  1.2022915  7.258529 0.0022695
# D-B  7.638864  4.6107449 10.666982 0.0000000
# D-C  3.408453  0.3803346  6.436572 0.0207466

# 手工计算..............................................................................
# 1.基础指标
a <- 0.05
k <- 4   #组数
na <- length(groupA)            #各组样本数 
nb <- length(groupB)
nc <- length(groupC)
nd <- length(groupD)
N <- na+nb+nc+nd                #样本总数

# 2. 组内均值和总均值
mXA <- mean(groupA)
mXB <- mean(groupB)
mXC <- mean(groupC)
mXD <- mean(groupD)
mX <- mean(c(groupA, groupB, groupC, groupD))

# 3. 组间平方和、组间均方
SSB <- na*(mXA-mX)^2+nb*(mXB-mX)^2+nc*(mXC-mX)^2+nd*(mXD-mX)^2     #结果同上，3166.455
MSB <- SSB/(k-1)                                                   #结果同上，1055.485

# 4. 组内平方和、组内均方
SSW <- sum((groupA-mXA)^2+(groupB-mXB)^2+(groupC-mXC)^2+(groupD-mXD)^2)  #结果同上，2348.143
MSW <- SSW/(N-k)                                                         #结果同上，20.24262

# 5. 计算统计量
F_value <- MSB/MSW     #结果同上，52.14173，注意这里不要占用关键字F，会影响画图

# 6. 双侧检验的p值，小于a，因此拒绝H0
pvalue <- 2*min(pf(F_value, k-1, N-k), 1-pf(F_value, k-1, N-k))  #结果为0，非常小

# 7. 统计量为f=52.14173，落在拒绝域内（落在了右侧），因此拒绝H0
Flow <- qf(alpha / 2, k-1, N-k)         #左临界值0.07169057
Fupp <- qf(1 - alpha / 2, k-1, N-k)     #右临界值3.230795

3. 多元回归中的F检验

# 设定显著性水平α=0.05
a <- 0.05

# 生成数据
set.seed(123)
x1 <- rnorm(100)
x2 <- rnorm(100)
x3 <- rnorm(100)
y <- 5 + 2*x1 + 3*x2 + 3*x3 + rnorm(100)

# 拟合线性模型
model <- lm(y ~ x1 + x2 + x3)

# 函数计算..............................................................................
# 查看F检验结果（模型整体显著性）
# summary(model)

# 结果如下：F统计量为547.5，对应p值< 2.2e-16， p<a，因此拒绝H0
# Call:
# lm(formula = y ~ x1 + x2 + x3)
# 
# Residuals:
#      Min       1Q   Median       3Q      Max 
# -2.49138 -0.65392  0.05664  0.67033  2.53210 
# 
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept)   4.9807     0.1073   46.40   <2e-16 ***
# x1            1.9445     0.1169   16.64   <2e-16 ***
# x2            3.0462     0.1095   27.83   <2e-16 ***
# x3            2.9426     0.1122   26.22   <2e-16 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 1.052 on 96 degrees of freedom
# Multiple R-squared:  0.9448,  Adjusted R-squared:  0.9431 
# F-statistic: 547.5 on 3 and 96 DF,  p-value: < 2.2e-16

# 手工计算..............................................................................
# 1. 基础变量
n <- length(y)    #样本量
k <- 3            #变量个数

# 2. 回归值、样本均值
y_hat <- predict(model)   #数值上用这个是一样的：y_hat <- fitted(model)
y_bar <- mean(y) 

# 3. 回归平方和SSR、残差平方和SSE
SSR <- sum((y_hat-y_bar)^2)
SSE <- sum((y-y_hat)^2)

# 4. F 统计量
F_value <- (SSR/(k))/(SSE/(n-k-1))      #结果同上，547.4956，注意这里不要占用关键字F，会影响画图

# 5. 双侧检验的p值，小于a，因此拒绝H0
pvalue <- 2*min(pf(F_value, k, n-k-1), 1-pf(F_value, k, n-k-1))  #结果为0，非常小

# 6. 统计量为f=547.4956，落在拒绝域内（落在了右侧），因此拒绝H0
Flow <- qf(alpha / 2, k, n-k-1)         #左临界值0.0716408
Fupp <- qf(1 - alpha / 2, k, n-k-1)     #右临界值3.255332

5.4 卡方检验

又叫χ2检验、chi-square检验。

目的

用于分析频数数据（分类变量）之间的关系或分布差异。主要有以下几类经典应用：

1. 检验一个分类变量（k长观测向量）：实际分布与理论分布是否相同（拟合优度检验）

2. 检验两个分类变量（r行c列频数表）：是否独立（独立性检验)

3. 检验不同总体（r行c列频数表）：在某一分类变量上的分布是否相同（同质性检验）

4. 检验分类数据的配对数据检验，处理前后是否一致，分以下几种使用场景：

   • 数据=2分类、测量=2次（2行2列频数表）：McNemar检验，麦克尼马尔检验

   • 数据>2分类、测量=2次（k行k列频数表）：Bowker检验

   • 数据=2分类、测量>2次（n行k列是否表）：Cochran’s Q检验

5. 分层列联表分析

   • 数据=2分类、测量=2次、多层次（2行2列k页频数表）：控制变量后的分类变量是否独立，CMH检验（Cochran–Mantel–Haenszel检验）

   • 数据=2分类、测量=2次、多层次（2行2列k页频数表）：联表层次之间优势比是否同质，Breslow–Day检验

注意：

• 此处所谓的拟合优度，是指样本分布与理论分布的相似程度，和回归分析中说的的R2不是一回事。

• 按统计量计算方法来说，上述1、2、3中，还可分为Pearson卡方（χ2）和似然比卡方（或对数似然卡方，G2），两者计算公式不同，但在大样本下结果相近。而4、5是独有的统计量，没有这种区分。

• Breslow–Day检验是CMH检验的前提检验之一，用以检验CMH检验的一个关键前提条件：分层之间是否存在效应一致性（优势比是否一致，即层变量是否可以作为控制变量）。

适用

• 检验类型：非参数检验（不依赖数据分布）
• 数据类型：分类数据（频数或列联表）
• 数据体量：大样本（期望频数需≥5，需合并类别、连续性校正、蒙特卡洛模拟或使用Fisher精确检验）
• 样本需独立随机抽样，且观测值之间互不干扰（配对数据检验除外）。

计算公式

1. 拟合优度检验

   • 皮尔逊卡方\(\chi^2=\sum_{i=1}^{k}{\frac{(O_i-E_i)^2}{E_i}}\)

   • 似然比卡方\(G^2=2\sum_{i=1}^k{O_i ln(\frac{O_i}{E_i})}\)

   • 自由度\(df=k-1-m\)

   • 其中，\(i\)，类别；\(k\)，类别总数；\(O_i\)，第i类的观测频数；\(E_i\)，第i类的期望频数；\(m\)，使用的参数，通常为0。注意这里的\(E_i\)，是你期望的第i类的频数，例如掷色子，1~6你期望的概率都相等，即均匀分布，则掷60次的话，期望频数都是10次。

2. 独立性检验/同质性检验

   • 皮尔逊卡方\(\chi^2=\sum_{i=1}^{r}{\frac{(O_{ij}-E_{ij})^2}{E_{ij}}}\)

   • 似然比卡方\(G^2=2\sum_{i=1}^{r}{\sum_{j=1}^{c}{O_{ij}ln\frac{O_{ij}}{E_{ij}}}}\)

   • 参数\(E_{ij}=\frac{R_iCj}{N}\)

   • 自由度\(df=(r-1)(c-1)\)

   • 其中；\(i\)代表行，共\(r\)行；\(j\)代表列，共\(c\)列；\(O_{ij}\)，实际观测到的第i行、第j列的频数；\(E_{ij}\)，在原假设（独立）下预期的频数，\(R_i\)是第i行的总和（行边际和），\(C_j\)是第j列的总和（列边际和），\(N\)是样本总数。

   • 注意，独立性检验和同质性检验的计算公式完全相同，但是检验目的（两个分类变量是否独立 vs. 多个总体的某个分类变量分布是否相同）、列联表含义（行是分类变量A列是分类变量B vs. 行是各组列是统一分类变量的取值）和H0（两个变量相互独立 vs. 所有总体分布完全一致）完全不同。

3. McNemar检验

   • 标准大样本（\(b+c\geq25\)）：统计量\(\chi^2=\frac{(b-c)^2}{b+c}\)

   • 小样本连续修正（Edwards修正、Yates连续性修正，\(b+c\leq25\)）：统计量\(\chi^2=\frac{(|b-c|-1)^2}{b+c}\)

   • 自由度\(df=1\)

   • 其中，abcd分别是2乘2配对表第1行第1列、第1行第2列、第2行第1列、第2行第2列中的频数

   • 注意， 减1是连续性修正的核心，目的是降低偏差，有修正时趋向于更保守（更难拒绝原假设）

4. Bowker检验

   • 统计量\(\chi^2=\sum_{i<j}{\frac{(a_{ij}-a_{ji})^2}{a_{ij}+a_{ji}}}\)

   • 自由度\(df=\frac{k(k-1)}{2}-n_0\)

   • 其中，\(a_{ij}\)是k*k频数表中第i行第j列的频数；\(k\)是数据分类数；\(n_0\)是\(a_{ij}+a_{ji}=0\)的个数。

5. Cochran’s Q检验

   • 统计量\(Q=\frac{(k-1)(k\sum_{j=1}^{k}{C_j^2}-T^2)}{k\sum_{i=1}^{n}{R_j}-\sum_{i=1}^n{R_i^2}}\)

   • 参数\(R_i=\sum_{j=1}^k{X_{ij}}\)

   • 参数\(C_j=\sum_{i=1}^n{X_{ij}}\)

   • 参数\(T=\sum_{i=1}^n{\sum_{j=1}^k{X_ij}}\)

   • 自由度\(df=k-1\)

   • 其中，\(k\)，测量次数（或者处理组数、列数）；\(n\)，受试者个数（行数）；\(X_{ij}\)，第i个个体在第j个测量下的二元相应（0或1）；\(R_i\)，第i个对象在所有处理下的成功次数（每行和），\(C_{j}\)第j个测试下成功的对象个数（每列和）；\(T\)，成功的总数。

   • 注意，Q的分子部分，\((k-1)(k\sum_{j=1}^{k}{C_j^2}-T^2)\)，后面是先求和、再乘以k、然后再减去T^2，然后整体再乘以前面的k-1

   • 注意，Q可以看作是方差（处理组之间成功人数）除以个体内部变异（每个人是否一致）

6. CMH检验

   • 无修正统计量\(\chi^2=\frac{[\sum_{i=1}^{k}{(A_i-E_i)}]^2}{\sum_{i=1}^{k}{V_i}}\)

   • 连续修正正统计量\(\chi^2=\frac{[\sum_{i=1}^{k}{(A_i-E_i)}-0.5]^2}{\sum_{i=1}^{k}{V_i}}\)

   • 参数\(E_i=\frac{N1_iM1_i}{T_i}\)

   • 参数\(V_i=\frac{N1_iN2_iM1_iM2_i}{T_i^2(T_i-1)}\)

   • 自由度\(df=1\)

   • 其中符号如下：定义第\(i\)层数据结构如下，所有数据共\(k\)层

     	Col 1	Col 2	Col total
     Row 1	\(A_i\)	\(B_i\)	\(N1_i\)
     Row 2	\(C_i\)	\(D_i\)	\(N2_i\)
     Row Total	\(M1_i\)	\(M2_i\)	\(T_i\)

7. Breslow–Day检验

   • 统计量\(\chi^2=\sum_{i=1}^k{\frac{(A_i-E_i)^2}{V_i}}\)

   • 参数\(\hat\psi_{MH}=\frac{\sum_{i=1}^k{(A_iD_i/Ti)}}{\sum_{i=1}^k{(B_iC_i/T_i)}}\)

   • 参数方程（用以计算\(E_i\)）\(E_i(N1_i-E_i)=\hat\psi_{MH}(M1_i-E_i)(M2_i-N1_i+E_i)\)

   • 参数\(V_i=(\frac{1}{E_i}+\frac{1}{M1_i-E_i}+\frac{1}{N1_i-E_i}+\frac{1}{M2_i-N1_i+E_i})^{-1}\)

   • 自由度\(df=k-1\)

   • 其中，数据第\(i\)层数据结构同上，所有数据共\(k\)层；\(E_i\)，原假设下\(A_i\)的期望值；\(V_i\)，\(A_i\)的渐进方差；\(\hat\psi_{MH}\)，使用CMH估计量计算的公共优势比（OR），通常使用CMH或自定义

注意

• McNemar检验和Bowker检验中，只关注非对角线上的非对称数据，例如bc和i≠j。

• 似然比卡方计算时，如果有某个\(O=0\)，则规定那一项为0，因为\(\lim_{x \to 0^+} xln(\frac x E)=0\)

结果判定

注意，与t检验和t、z检验不同，卡方检验的统计量始终是正数，因此都是右侧检验，即判断统计量是否偏大，而没有所谓的双侧或单侧之分。

上述卡方检验的种类、H0、统计量、自由度

检验类型	原假设 H₀	统计量	自由度
拟合优度检验	实际分布 = 理论分布	X2或G2	k−1−m
独立性检验	两个变量相互独立	X2或G2	(r-1)(c-1)
同质性检验	不同组的分布规律一致	X2或G2	(r-1)(c-1)
McNemar检验	配对前后的判断一致（b≈c）	X2	1
Bowker检验	配对多分类判断是对称的	X2	k(k-1)/2
Cochran’s Q检验	配对二分类多处理条件下比例一致	Q	k-1
CMH检验	控制变量后两分类变量独立	X2	1
Breslow–Day检验	各层同质（优势比相同）	X2	k-1

统计量（这里以\(\chi^2\)代替全部统计量符号）服从卡方分布，该分布左右不对称。设定显著性水平为α后，检验结果判定步骤如下：

1. \(p=P(X\geq \chi^2)\)

2. 拒绝域\(\chi^2 \geq \chi^2_{1-\alpha,df}\)

案例

1. 拟合优度检验（k长观测向量，H0：实际分布 = 理论分布）

R语言中，使用stats::chisq.test和DescTools::GTest执行拟合优度检验，两个函数差不多，其中的参数，含义和用法如下：

• x， 向量（用于拟合优度检验）或表格（用于独立性/同质性检验）

• y， 第二个向量（通常不使用）

• stats::chisq.test中的correct，是否使用连续性修正（Yates correction），仅对 2×2 表有效

• DescTools::GTest中的correct，是否对G2统计量进行修正，选项有：“none”（默认）、“williams”（推荐）、“yates”（连续性修正，仅 2×2）

• p，拟合优度检验中的理论概率分布（默认均匀），默认传入期望频率

• rescale.p， 是否将p归一化成概率，如果定义为T，则p就可以传入频数了

• simulate.p.value，是否使用模拟（Monte Carlo）法估计p值（用于样本太小或分布太偏时）

• B，如果 simulate.p.value = TRUE，模拟的重复次数（默认2000次）

注意：

• 当m≠0时，额外估计的参数，并不能传入stats::chisq.test和DescTools::GTest，此时：

• 如果给定的是“真正外部”的理论分布，可以直接用df=k-1

• 如果基于数据估计 ，必须手动扣除参数数量，另df=k-1-m

• 若期望频数<5，使用模拟simulate.p.value=T或Fisher精确检验fisher.test(observed)

案例1，一个骰子掷了120次，1~6点对应的频数分别为c(20,22,16,18,22,22)，问这个骰子是否是均匀的。

# 设定显著性水平α=0.05
a <- 0.05

# 观察频数
observed <- c(20, 22, 16, 18, 22, 22)

# 函数计算，皮尔逊卡方......................................................................
# 拟合优度检验
test.1 <- chisq.test(observed, p=rep(c(1/6), 6))   #默认p是期望概率，而非频数，可定义为其他频率

# 结果如下：p-value > a，表明不能拒绝H0，认为观测分布与理论分布一致
#   Chi-squared test for given probabilities
# 
# data:  observed
# X-squared = 1.6, df = 5, p-value = 0.9012

# 函数计算，似然比卡方......................................................................
# 加载包
library(DescTools)
test.2 <- GTest(observed, p=rep(c(1/6), 6))     #这里p是期望概率，而非频数，可定义为其他频率

# 结果如下：p-value > a，表明不能拒绝H0，认为观测分布与理论分布一致
#   Log likelihood ratio (G-test) goodness of fit test
# data:  observed
# G = 1.6474, X-squared df = 5, p-value = 0.8955

# 手工计算，皮尔逊卡方及似然比卡方..........................................................
observed <- c(20, 22, 16, 18, 22, 22)      #观测频数
expected <- 120*c(1/6,1/6,1/6,1/6,1/6,1/6) #理论频数
df <- length(observed)-1

# 计算统计量
X2 <- sum((observed-expected)^2/expected)    #X2=1.6
G2 <- 2*sum(observed*log(observed/expected)) #G2=1.647372

# 计算p值和拒绝域临界值
pv_X2 <- 1-pchisq(X2, df)              #pv_X2=0.9012493
pv_G2 <- 1-pchisq(G2, df)              #pv_G2=0.8954622
thres <- qchisq(1-0.05, df)            #thres=11.0705

# 结果表明：p>a，X2<thres，G2<thres，因此不能拒绝H0，认为观测分布与理论分布一致

案例2，有参数估计的拟合优度检验，你测量了一组连续变量，将其分为若干区间（分组频数数据），想检验这组数据是否符合正态分布。你测量的区间为c(-Inf,60,70,80,90,Inf)，对应的频数为c(5,18,42,27,8)。此外，你还估计了分布参数，包括样本均值为75，样本标准差10。

# 设定显著性水平α=0.05
a <- 0.05

# 观察频数
observed <- c(5, 18, 42, 27, 8)
# 区间边界
breaks <- c(-Inf, 60, 70, 80, 90, Inf)
# 样本总数
n <- sum(observed)
# 样本均值与标准差（假设已估计）
mu <- 75
sigma <- 10
# 计算正态分布下的区间概率
probs <- pnorm(breaks[-1], mean=mu, sd=sigma)-pnorm(breaks[-length(breaks)], mean=mu, sd= sigma)

# 函数计算，皮尔逊卡方......................................................................
# 执行拟合优度检验
test.1 <- chisq.test(observed, p = probs, rescale.p = TRUE)

# 结果表明：p-value > a，不能拒绝原假设，认为观测数据确实符合均值为75，标准差为10的正态分布
#   Chi-squared test for given probabilities
# data:  observed
# X-squared = 2.9493, df = 4, p-value = 0.5663           #注意这里df=4，就是k-1

# 手工计算，皮尔逊卡方及似然比卡方..........................................................
# 扣除参数的自由度
df <- length(c(5, 18, 42, 27, 8))-1-2

# 计算频数
expected <- n*probs

# 计算统计量
X2 <- sum((observed-expected)^2/expected)    #X2=2.949326
G2 <- 2*sum(observed*log(observed/expected)) #G2=3.105678

# 计算p值和拒绝域临界值
pv_X2 <- 1-pchisq(X2, df)              #pv_X2=0.2288559
pv_G2 <- 1-pchisq(G2, df)              #pv_G2=0.2116462
thres <- qchisq(1-0.05, df)            #thres=5.991465

# 结果表明：p>a，X2<thres，G2<thres，因此不能拒绝H0，认为观测分布与理论分布一致

2. 独立性检验（r行c列频数表，H0：两个变量相互独立）

某调查机构的调查结果如下，问职业与吸烟行为是否相互独立（即吸烟行为是否与职业相关）？：

职业	吸烟	不吸烟	合计
白领	90	210	300
蓝领	150	150	300
教师	60	240	300
合计	300	600	900

# 设定显著性水平α=0.05
a <- 0.05

# 构造列联表
observed <- matrix(c(90,  210,      #白领
                     150, 150,      #蓝领
                     60,  240),     #教师
                   nrow = 3, byrow = TRUE)

# 函数计算，似然比卡方......................................................................
# 检验期望频数，如果所有格子的期望频数≥5，那么Pearson卡方是稳健的。
# numb_expected <- chisq.test(observed)$expected  #满足大样本
# 如果是小样本，应考虑Fisher精确检验：fisher.test(observed)

# 进行卡方独立性检验
test.1 <- chisq.test(observed)

# 结果如下：p-value < a，表明拒绝H0，认为吸烟行为与职业不独立（有关系）
#   Pearson's Chi-squared test
# data:  observed
# X-squared = 63, df = 2, p-value = 2.088e-14

# 函数计算，似然比卡方......................................................................
# 加载包
library(DescTools)
test.2 <- GTest(observed)

# 结果如下：p-value < a，表明拒绝H0，认为吸烟行为与职业不独立（有关系）
#   Log likelihood ratio (G-test) test of independence without correction
# data:  observed
# G = 63.077, X-squared df = 2, p-value = 2.01e-14

# 手工计算，皮尔逊卡方及似然比卡方..........................................................
# 自由度
df <- (nrow(observed)-1)*(ncol(observed)-1)

# 行、列合计
row_totals <- rowSums(observed)
col_totals <- colSums(observed)
grand_total <- sum(observed)

# 计算期望频数矩阵
expected <- outer(row_totals, col_totals) / grand_total

# 计算统计量
X2 <- sum((observed - expected)^2 / expected)                                #X2=63
G2 <- 2 * sum(ifelse(observed == 0, 0, observed * log(observed / expected))) #G2=63.07716

# 计算p值和拒绝域临界值
pv_X2 <- 1-pchisq(X2, df)              #pv_X2=2.087219e-14
pv_G2 <- 1-pchisq(G2, df)              #pv_G2=2.009504e-14
thres <- qchisq(1-0.05, df)            #thres=5.991465

# 结果表明：p<a，X2>thres，G2>thres，因此拒绝H0，认为吸烟行为与职业不独立（有关系）

3. 同质性检验（r行c列频数表，H0：不同组分布规律一致）

某机构调查了三所学校（A、B、C），并记录了各自满意、一般、不满意的学生人数，问3所学校的学生对自己学校的满意度是否相同。

	满意	一般	不满意	总计
学校A	40	35	25	100
学校B	50	30	20	100
学校C	60	25	15	100

# 设定显著性水平α=0.05
a <- 0.05

# 构造列联表
observed <- matrix(c(40, 35, 25,
                     50, 30, 20,
                     60, 25, 15),
                   nrow = 3, byrow = TRUE)

# 函数计算，似然比卡方......................................................................
# 检验期望频数，如果所有格子的期望频数≥5，那么Pearson卡方是稳健的。
# numb_expected <- chisq.test(observed)$expected  #满足大样本
# 如果是小样本，应考虑Fisher精确检验：fisher.test(observed)

# 进行卡方同质性检验
test.1 <- chisq.test(observed)

# 结果如下：p-value > a，表明不能拒绝H0，认为三所学校的满意度分布一致
#   Pearson's Chi-squared test
# data:  observed
# X-squared = 8.1667, df = 4, p-value = 0.08566
 
# 函数计算，似然比卡方......................................................................
# 加载包
library(DescTools)
test.2 <- GTest(observed)

# 结果如下：p-value > a，表明不能拒绝H0，认为三所学校的满意度分布一致
#   Log likelihood ratio (G-test) test of independence without correction
# data:  observed
# G = 8.2283, X-squared df = 4, p-value = 0.08356

# 手工计算，皮尔逊卡方及似然比卡方..........................................................
# 自由度
df <- (nrow(observed)-1)*(ncol(observed)-1)

# 行、列合计
row_totals <- rowSums(observed)
col_totals <- colSums(observed)
grand_total <- sum(observed)

# 计算期望频数矩阵
expected <- outer(row_totals, col_totals) / grand_total

# 计算统计量
X2 <- sum((observed - expected)^2 / expected)                                #X2=8.166667
G2 <- 2 * sum(ifelse(observed == 0, 0, observed * log(observed / expected))) #G2=8.228288

# 计算p值和拒绝域临界值
pv_X2 <- 1-pchisq(X2, df)              #pv_X2=0.08566027
pv_G2 <- 1-pchisq(G2, df)              #pv_G2=0.08356472
thres <- qchisq(1-0.05, df)            #thres=9.487729

# 结果表明：p>a，X2<thres，G2<thres，表明不能拒绝H0，认为三所学校的满意度分布一致

4. McNemar检验（数据=2分类、测量=2次，2行2列频数表，H0：配对前后的判断一致（b≈c））

R语言中，使用stats::mcnemar.test执行McNemar检验，参数如下：

• x, 一个二维条件表矩阵，或者因子向量

• y, 一个因子向量，如果x是一个矩阵，则为NULL

• correct，是否连续性修正（在小样本下建议修正），默认是做修正

注意：当b+c=0时，McNemar检验无统计量返回。

一项研究调查，访问了100名顾客，问他们在看过某产品的广告之后，购买意向是否改变，统计表格如下：

2乘2条件表	有购买意向（看广告后）	无购买意向（看广告后）	合计
有购买意向（看广告前）	40	10	50
无购买意向（看广告前）	30	20	50
合计	70	30	100

# 设定显著性水平α=0.05
a <- 0.05

# 构造 2x2 表  
observed <- matrix(c(40, 10, 
                     30, 20), 
                   nrow = 2, byrow = TRUE)

# 判断是否需要进行修正
b <- observed[1,2]
c <- observed[2,1]
iscorr <- b+c  #这里大于25，是否修正影响不大，但是否修正都计算了一遍

# 函数计算.................................................................................
# 执行 McNemar 检验
test.1 <- mcnemar.test(observed, correct=F)
test.2 <- mcnemar.test(observed, correct=T)

# 结果如下（有修正时更保守）：

# 对test.1，即大样本无修正时，p-value < a，表明拒绝H0，认为观看广告前后意愿不一致
#   McNemar's Chi-squared test
# data:  observed
# McNemar's chi-squared = 10, df = 1, p-value = 0.001565

# 对test.2，即小样本有修正时，p-value < a，表明拒绝H0，认为观看广告前后意愿不一致
#   McNemar's Chi-squared test with continuity correction
# data:  observed
# McNemar's chi-squared = 9.025, df = 1, p-value = 0.002663

# 手工计算，皮尔逊卡方及似然比卡方..........................................................
df <- 1

# 自由度
X2_big <- (b-c)^2/(b+c)           #大样本无修正时，X2_big=10
X2_lit <- (abs(b-c)-1)^2/(b+c)    #小样本有修正时，X2_lit=9.025

# 计算p值和拒绝域临界值
pv_Fcor <- 1-pchisq(X2_big, df)        #pv_Fcor=0.001565402
pv_Tcor <- 1-pchisq(X2_lit, df)        #pv_Tcor=0.002663119
thres <- qchisq(1-0.05, df)            #thres=3.841459

# 结果表明：pv_Fcor<a，pv_Fcor<a，X2_big>thres，X2_lit>thres，表明能拒绝H0，认为观看广告前后意愿不一致

5. Bowker检验（数据>2分类、测量=2次，k行k列频数表，H0：配对多分类判断是对称的）

注意，在aij+aji全都不为0的时候，mcnemar.test也可以得到相同的结果，貌似是对McNemar检验做了扩展。然而，当存在0时，mcnemar.test的统计量就为NaN了（因为分母有0了，这个问题在McNemar检验中也存在），因此这里自定义了一个函数计算Bowker检验。

研究交通噪音对人的情绪的影响，招募了72人，对某个噪音使人恼怒的程度评价两次（尽可能排除前后两次的关联性）结果如下，问两次评价是否有差异。

第一次＼第二次	很低	低	中等	高	很高	合计
很低	51	28	3	0	0	82
低	15	68	40	5	1	129
中等	0	29	77	21	1	128
高	0	4	19	80	14	117
很高	0	1	5	26	88	120
合计	66	130	144	132	104	576

# 设定显著性水平α=0.05
a <- 0.05

# 构造配对列联表
observed <- matrix(c(51, 28, 3, 0, 0,
                     15, 68, 40, 5, 1,
                     0, 29, 77, 21, 1,
                     0, 4, 19, 80, 14,
                     0, 1, 5, 26, 88),
                   nrow = 5, byrow = TRUE)

# 定义一个bowker检验函数
BowkerTest <- function(tbl) {
  if (!is.matrix(tbl) || nrow(tbl) != ncol(tbl)) {
    stop("输入必须是一个方阵矩阵（k x k 对称列联表）")
  }
  # 初始化
  k <- nrow(tbl)
  chisq <- 0
  df <- 0
  # 计算统计量和自由度
  for (i in 1:(k - 1)) {
    for (j in (i + 1):k) {
      b <- tbl[i, j]
      c <- tbl[j, i]
      if (b + c > 0) {
        chisq <- chisq + (b - c)^2 / (b + c)
        df <- df + 1
      }
    }
  }
  # 计算p值
  p.value <- pchisq(chisq, df = df, lower.tail = FALSE)
  # 整理并输出结果
  result <- list(
    statistic = c("X-squared" = chisq),
    parameter = c("df" = df),
    p.value = p.value,
    method = "Bowker's Test of Symmetry",
    data.name = deparse(substitute(tbl))
  )
  class(result) <- "htest"
  return(result)
}
# 执行检验
test.1 <- BowkerTest(observed)

# 结果如下：p-value > a，表明不能拒绝H0，认为配对多分类判断是对称的，即两次评价无差别
#   Bowker's Test of Symmetry
# data:  observed
# X-squared = 15.162, df = 8, p-value = 0.05608

# 计算拒绝域临界值
thres <- qchisq(1-0.05, df)            #thres=15.50731
# 结果表明：X2<thres，表明不能拒绝H0，认为配对多分类判断是对称的，即两次评价无差别

# mcnemar.test检验，因为有配对和为0，因此无法计算出统计量
test.2 <- mcnemar.test(observed)
#   McNemar's Chi-squared test
# data:  observed
# McNemar's chi-squared = NaN, df = 10, p-value = NA

6. Cochran’s Q检验（数据=2分类、测量>2次，n行k列是否表，H0：配对二分类多处理条件下比例一致）

R语言中，使用DescTools::CochranQTest执行Cochran’s Q检验

比较3种广告样式 A、B、C 对用户点击率的影响，让10名用户每人都看了这3个广告，并记录他们是否点击（1表示点击，0表示未点击），问三种广告的点击率是否一致

用户编号	样式A	样式B	样式C
1	1	1	0
2	1	1	1
3	0	1	0
4	1	0	1
5	1	1	1
6	0	0	1
7	0	1	0
8	1	1	1
9	1	0	0
10	0	1	1

# 加载包
library(DescTools)
library(magrittr)

# 设定显著性水平α=0.05
a <- 0.05

# 构建二元响应数据
observed <- data.frame(
  A = c(1, 1, 0, 1, 1, 0, 0, 1, 1, 0),
  B = c(1, 1, 1, 0, 1, 0, 1, 1, 0, 1),
  C = c(0, 1, 0, 1, 1, 1, 0, 1, 0, 1)
) %>% as.matrix()

# 函数计算..............................................................................
test.1 <- CochranQTest(observed)

# 结果如下：p-value > a，表明不能拒绝H0，认为三组广告样式无差异
#   Cochran's Q test
# data:  y
# Q = 0.28571, df = 2, p-value = 0.8669

# 手工计算..............................................................................
k   <- ncol(observed)            #测试组数
n   <- nrow(observed)            #受试者个数
Ri  <- apply(observed, 1, sum)   #个体i在所有处理中的综合（每行和）
Cj  <- apply(observed, 2, sum)   #第j个处理下成功的个体数（每列和）
Tij <- sum(observed)             #总体成功数
df  <- k-1                       #自由度

# 计算统计量
Q <- (k-1)*(k*sum(Cj^2)-Tij^2)/(k*sum(Ri)-sum(Ri^2))  #Q=0.2857143

# 计算p值和拒绝域临界值
pv <- 1-pchisq(Q, df)        #pv=0.8668779
thres <- qchisq(1-0.05, df)  #thres=5.991465

# 结果表明：p>a，Q<thres，表明不能拒绝H0，认为三组广告样式无差异

7. CMH检验（数据=2分类、测量=2次、多层次，2行2列k页频数表，H0：控制变量后两分类变量独立）

R语言中，使用stats::mantelhaen.test执行CMH检验

现有3个医院，每个医院分别进行一项药物试验，我们想知道，该药物是否比安慰剂更有效。这种情况下，直接比较可能有偏差，因为不同医院的测试可能条件不同，比如医生水平、病人情况、药物保存方式等，因此如果不控制医院这个变量，直接合并数据来比较治愈率，会混入新的不确定性。因此，CMH检验把每家医院当作一个分层，在每个医院内部比较药物与安慰剂的效果，然后再综合所有医院的比较结果。这种“控制医院后”的检验，才能更客观地反映药物的作用。

此时，H0表示控制医院后治疗方式与治愈率独立（排除医院的影响后，药物无治疗效果，即不考虑医院差异，服药无治疗效果），H1表示，控制医院后药物与治疗率不独立（排除医院的影响后，药物有治疗效果，即不考虑医院差异，服药有治疗效果）

按医院分层的药物试验数据如下

医院 1	治愈	未愈	合计
药物组	20	10	30
安慰剂组	12	18	30

医院 2	治愈	未愈	合计
药物组	15	15	30
安慰剂组	10	20	30

医院 3	治愈	未愈	合计
药物组	25	5	30
安慰剂组	15	25	30

# 设定显著性水平α=0.05
a <- 0.05

# 观测数据，三维矩阵
observed <- array(
  c(
    # 医院 1
    20, 10,   # 药物组（治愈，未愈）
    12, 18,   # 安慰剂组

    # 医院 2
    15, 15,
    10, 20,

    # 医院 3
    25, 5,
    15, 15
  ),
  dim = c(2, 2, 3)
)

# 函数计算.................................................................................
# 执行检验
test.1 <- mantelhaen.test(observed, correct = T)    #使用矫正的检验

# 结果如下：p-value < a，表明拒绝H0，认为配控制医院后药物与治疗率不独立，即药物有治疗效果
#   Mantel-Haenszel chi-squared test with continuity correction
# data:  observed
# Mantel-Haenszel X-squared = 11.107, df = 1, p-value = 0.00086
# alternative hypothesis: true common odds ratio is not equal to 1
# 95 percent confidence interval:
#  1.606887 5.600890
# sample estimates:
# common odds ratio 
#                 3 

# 手工计算.................................................................................
k <- dim(observed)[3]       #层数
Ai <- Ei <- Vi <- rep(0,k)  #初始化k层的实际频数、XY期望值、方差
df <- 1

# 按层计算参数
for(i in 1:k){
  obseri <- observed[,,i]
  N1i <- obseri[1,1]+obseri[1,2]
  M1i <- obseri[1,1]+obseri[2,1]
  N2i <- obseri[2,1]+obseri[2,2]
  M2i <- obseri[1,2]+obseri[2,2]
  Ti  <- sum(obseri)
  Ai[i] <- obseri[1,1]
  Ei[i] <- N1i*M1i/Ti
  Vi[i] <- N1i*N2i*M1i*M2i/Ti^2/(Ti-1)
}

# 计算统计量
X2 <- (sum(Ai-Ei)-0.5)^2/sum(Vi)       #X2=11.10696，矫正后的

# 计算p值和拒绝域临界值
pv <- 1-pchisq(X2, df)                 #pv=0.0008600424
thres <- qchisq(1-0.05, df)            #thres=3.841459

# 结果表明：p<a，X2>thres，因此拒绝H0，认为配控制医院后药物与治疗率不独立，即药物有治疗效果

8. Breslow–Day检验（数据=2分类、测量=2次、多层次，2行2列k页频数表，H0：各层同质（优势比相同））

R语言中，使用DescTools::BreslowDayTest执行Breslow–Day检验：

数据同上，这里用于检验不同医院之间是同质，即其优势比是否相同。这里计算繁琐，直接使用BreslowDayTest函数计算，其中，OR优势比未提供，使用CMH法计算，Ei计算使用三次方程求解。具体过程可直接查看该函数源码。

医院 1	治愈	未愈	合计
药物组	20	10	30
安慰剂组	12	18	30

医院 2	治愈	未愈	合计
药物组	15	15	30
安慰剂组	10	20	30

医院 3	治愈	未愈	合计
药物组	25	5	30
安慰剂组	15	25	30

# 设定显著性水平α=0.05
a <- 0.05

# 观测数据，三维矩阵
observed <- array(
  c(
    # 医院 1
    20, 10,   # 药物组（治愈，未愈）
    12, 18,   # 安慰剂组

    # 医院 2
    15, 15,
    10, 20,

    # 医院 3
    25, 5,
    15, 15
  ),
  dim = c(2, 2, 3)
)

# 函数计算.................................................................................
# 执行检验
library(DescTools)
test.1 <- BreslowDayTest(observed)

# 结果如下：p-value > a，表明不能拒绝H0，认为三所医院无差异（优势比同质）
#   Breslow-Day test on Homogeneity of Odds Ratios
# data:  observed
# X-squared = 1.2876, df = 2, p-value = 0.5253

5.5 Fisher检验

Fisher提出了很多检验方法，这里特指Fisher精确检验（Fisher’s Exact Test）。

目的

• 检验两个分类变量：是否独立

适用

• 检验类型：非参数检验（无分布假设）
• 数据类型：分类数据（2×2列联表）
• 数据体量：小样本（期望频数<5时可替代卡方检验）

计算公式

Fisher精确检验本身并没有一个统计量的计算公式，而是直接通过超几何分布概率计算当前观察到的2×2列联表的精确概率。

定义列联表中的符号如下

F＼R	Result1	Result2	Col total
Factor1	\(a\)	\(b\)	\(N1\)
Factor2	\(c\)	\(d\)	\(N2\)
Row Total	\(M1\)	\(M2\)	\(T\)

超极核概率\(P\)的计算公式为：\[P=\frac{\binom{a+b}{a}\binom{c+d}{c}}{\binom{n}{a+c}}=\frac{(a+b)!(c+d)!(a+c)!(b+d)!}{a!b!c!d!n!}\]

基于此，统计检验计算步骤如下：

1. 写出2×2表格并记下a,b,c,d,n

2. 计算当前表格的超几何概率\(P(a)\)

3. 枚举在表格边际和固定的情况下所有可能的a值，其中，\(a \in [max(0, T-M1-N1), min(N1, M1)]\)

4. 计算所有可能a值对应的超几何概率\(P_i\)

结果判定

Fisher精确检验的结果判定，设定显著性水平α后，检验结果的判定步骤如下：

1. 若H0：Factor1和Factor2对于Result1概率相同（即行变量与列变量不相关）

   • 双侧检验

   • 累加法：计算所有小于\(P(a)\)的\(Pi\)的和，即\(p=\sum_{Pi\leq P(a)}{P_i}\)

   • 对称法（R语言默认）：计算\(P_i\)序列中，与\(P(a)\)对称的和，即\(p=P(a)+P(max(i)-a)\)

   • 平滑法（minp法，不常用）：\(p=\sum_{x<a}{P(x)}+\frac 1 2 P(a)\)

2. 若H0：Factor1引发Result1的概率 ≤ Factor2引发Result1的概率 vs. H1：Factor1引发Result1的概率 > Factor2引发Result1的概率

   • 右侧检验

   • \(p=\sum_{x\geq a}{P(x)}\)

3. 若H0：Factor1引发Result1的概率 ≥ Factor2引发Result1的概率 vs. H1：Factor1引发Result2的概率 < Factor1引发Result1的概率

   • 左侧检验

   • \(p=\sum_{x\leq a}{P(x)}\)

4. 若\(p\leq\alpha\)则拒绝H0

案例

R语言中，使用stats::fisher.test执行Fisher精确检验，其中的参数和用法如下（Fisher精确检验计算繁琐，这里只使用函数进行计算，未验证上述公式是否正确）：

• x：矩阵2乘2频数表，或者向量

• y：当x是向量时，y是另一个向量

• alternative：H1方向（等同与检验方向），可选two.sided、greater或less

• or：原假设下的优势比（OR），检验OR是否等于设定值，默认是1

• workspace：分配的内存空间，默认约200kb，当列联表较大时，增加该值以避免计算错误

• hybrid：是否使用混合算法（还是精确概率），默认F，对大于2×2的列联表使用混合近似算法可节省计算资源

• hybridPars：长度为3的向量，用于控制混合算法的阈值

• simulate.p.value：是否使用蒙特卡洛模拟计算p值（适用于大表或复杂情况），默认F

• B：蒙特卡洛模拟次数，默认2000，次数越多结果越稳定。

• conf.int：是否计算OR的置信区间（默认TRUE）

• conf.level：置信水平（默认0.95，即95%置信区间）

案例：研究地区对于环保政策的态度是否有差异，调研数据如下

地区	支持环保政策	不支持	合计
城市居民	30	10	40
乡村居民	20	20	40
合计	50	30	80

# 设定显著性水平α=0.05
a <- 0.05

# 创建表格
observed <- matrix(c(30, 20, 10, 20), nrow = 2)

# H0：地区对于环保政策无差异
test.1 <- fisher.test(observed, alternative = "two.side")

# 结果如下：证明p-value < a，说明拒绝原假设，认为地区对于环保政策有影响
#   Fisher's Exact Test for Count Data
# data:  observed
# p-value = 0.03683
# alternative hypothesis: true odds ratio is not equal to 1
# 95 percent confidence interval:
#  1.060077 8.704448
# sample estimates:
# odds ratio 
#   2.957367 

# 双侧检验.................................................................................
# H0：地区对于环保政策无差异
test.1 <- fisher.test(observed, alternative = "two.side")
# 结果如下：证明p-value < a，说明拒绝原假设，认为地区对于环保政策有影响
#   Fisher's Exact Test for Count Data
# data:  observed
# p-value = 0.03683
# alternative hypothesis: true odds ratio is not equal to 1
# 95 percent confidence interval:
#  1.060077 8.704448
# sample estimates:
# odds ratio 
#   2.957367 

# 右侧检验.................................................................................
# H0：城市居民支持环保政策的概率 ≤ 农村居民支持环保政策的概率
test.2 <- fisher.test(observed, alternative = "greater")
# 结果如下：证明p-value < a，说明拒绝原假设，认为城市居民比农村居民更支持环保政策
#   Fisher's Exact Test for Count Data
# data:  observed
# p-value = 0.01842
# alternative hypothesis: true odds ratio is greater than 1
# 95 percent confidence interval:
#  1.22458     Inf
# sample estimates:
# odds ratio 
#   2.957367 

# 左侧检验.................................................................................
# H0：城市居民支持环保政策的概率 ≥ 农村居民支持环保政策的概率
test.3 <- fisher.test(observed, alternative = "less")
# 结果如下：证明p-value > a，说明拒绝原假设，认为城市居民比农村居民更支持环保政策
#   Fisher's Exact Test for Count Data
# data:  observed
# p-value = 0.9948
# alternative hypothesis: true odds ratio is less than 1
# 95 percent confidence interval:
#  0.000000 7.409392
# sample estimates:
# odds ratio 
#   2.957367

5.6 Wilcoxon检验

包括Wilcoxon符号秩检验和Wilcoxon秩和检验（也称Mann-Whitney U检验）

目的

符号秩检验：

• 检验单样本数据：中位数是否等于某个值，如0

• 检验两组配对数据：中位数差是否为0

秩和检验：

• 检验两组独立样本：分布是否相同

注意

• 秩和检验中，由于是非参数且不依赖分布，因此

  • 当两组数据分布已知且相同时，用于检验两数据中位数的：=、≤、≥

  • 当两组数据分布未知时，只能用于检验前一组数据对后者的：分布相同、分布整体向左、分布整体向右

适用

• 检验类型：非参数检验（不要求正态分布）
• 数据类型：数值数据（连续或序数）
• 数据体量：小样本（配对样本或两独立样本，n<30 时常用）
• 对异常值比参数检验（如t检验）更稳健

计算公式

这里的公式，是大样本下的近似估计值，此时最终的统计量Z服从标准正态分布

• 符号秩检验：

  1. 单样本观测数据为\(X_i\)，共\(n\)个，要检验其中位数是否为\(\mu_0\)，对应差值为\(d_i=X_i-\mu_0\)；或者配对数据观测样本对为(\(X_i\), \(Y_i\))，共\(n\)对，对应差值为\(d_i=X_i-Y_i\)

  2. 排除差值为0的数据对（即\(di=0\)），剩余有效样本为\(n'\)

  3. 求绝对值，分配秩（相同数值赋为平均秩）：\(r_i=rank(|d_i|)\)

  4. 计算正负差值的秩和：\(W_+=\sum_{di>0}r_i\)，\(W_-=\sum_{di<0}r_i\)

  5. 期望：\(\mu_W=\frac{n'(n'+1)}{4}\)

  6. 标准差：\(\sigma_W=\sqrt{\frac{n'(n'+1)(2n'+1)}{24}}\)

  7. 如果剔除\(d_i=0\)后，\(|d_i|\)中有相同值（ties），需要修正标准差。记\(|d_i|,i=1,2,3,...n'\)中，不重复的元素分别为\(d_{i}^{valid},i=1,2,3,...,j\)，对应的频数为\(t_1,t_2,...,t_j\)，则校正标准差为：\(\sigma'_W=\sqrt{\sigma_W-\frac{1}{48}\sum_{j=1}^{n'}{(t_j^3-t_j)}}\)

  8. 标准化统计量（连续性修正，左侧+0.5，右侧或双侧-0.5）：\(Z=\frac{W_+-\mu_W\pm0.5}{\sigma'_W}\)

• 秩和检验

  1. 设两组独立样本分别为\(X_1,X_2,X_3,...,X_m\)，\(Y_1,Y_2,Y_3,...,Y_n\)

  2. 将两组数据合并一起成为新的集合，然后计算秩（相同数值赋为平均秩）

  3. W统计量：取出X组在这个合并集合中的秩，并求和\(W=\sum_{i=1}^m{rank(X')}\)

  4. 对于W：期望\(\mu_E=\frac{m(m+n+1)}{2}\)，标准差\(\sigma_W=\sqrt{\frac{mn(m+n+1)}{12}}\)

  5. U统计量：\(U=\sum_{i=1}^m{\sum_{j=1}^n{I(X_i>Y_i)}}+0.5\sum_{i=1}^m{\sum_{j=1}^n{I(X_i=Y_i)}}\)，其中函数\(I\)是指示函数，满足条件为1，否则为0。或者U统计量还可以由W求出：\(U=W-\frac{m(m+1)}{2}\)

  6. 对于U：期望\(\mu_U=\frac{mn}{2}\)，标准差\(\sigma_U=\sqrt{\frac{mn(m+n+1)}{12}}\)

  7. 关于\(\sigma\)的ties修正：\(\sigma_{cor}=\sqrt{\frac{mn(m+n+1)}{12}-\frac{mn}{12}\sum_{j}{\frac{t_j^3-t_j}{(m+n)(m+n-1)}}]}\)

  8. 标准化统计量（连续性修正，左侧+0.5，右侧或双侧-0.5）：\(Z=\frac{W-\mu_W\pm0.5}{\sigma_{cor}}\)或\(Z=\frac{U-\mu_U\pm0.5}{\sigma_{cor}}\)

注意：

• 上述两个检验的计算方法，都是大样本下的满足渐进正态的估算值，如果样本量较小，则用穷举计算精确解

• 秩和检验中，U统计量是一致的，无论检验是双侧、右侧还是左侧，U的计算公式都是上面的第5点，即x>y。此时分别对应p值的计算方式为标准的(1-pnorm(abs(Z)))*2、1-pnorm(Z)、pnorm(Z)。然而，由于U近似服从正态分布，因此还有个中间对称值，用x<y计算，或者从整体减去另一个\(U'=mn-U\)。ChatGPT说要根据H0选择U的形式，比如H0为=对应min(U,U’)，<对应U，>对应U’，这是乱的，因为后面计算p值的公式，仍然需要根据方向变化，并不是上面标准的方法。

• 秩和检验中，总秩和的计算公式为\(\frac{(m+n)(m+n+1)}{2}\)，当\(m=n\)时，期望秩和\(\mu_W=\frac{m(m+n+1)}{2}\)正好是总秩和的1/2。统计量U的公式为\(U=W-\frac{m(m+1)}{2}\)

• 秩和检验中，统计量W的中心为

• 穷举精确法计算法，直接使用wilcox.test函数算吧，和网上手算的公式不同，对不上结果。

• 对于使用近似估计还是精确估计，通常根据样本量（阈值10~20），当小样本时可用精确估计，大样本用近似估计

• 当有ties时，秩和检验的精确估计失效，此时只能用近似估计

• 当样本量很大时，精确估计需要穷举，会非常耗时

结果判定

统计量Z近似服从标准正态分布。设定显著性水平为α后，符号秩检验或秩和检验结果的判定方式如下：

近似解估计（大样本）时，判定标准同z检验：

1. 若H0为\(d_i=0\)（符号秩检验）或H0为两分布相同（秩和检验），则双侧检验：

   • \(p=2P(X\geq|Z|)\)，拒绝域为\(Z \geq Z_{1-\alpha/2}\)

2. 若H0为\(d_i\leq0\)（符号秩检验）或H0为X整体比Y偏左（秩和检验），则右侧检验：

   • \(p=P(X\geq Z)\)，拒绝域为\(Z \geq Z_{1-\alpha}\)

3. 若H0为\(d_i\geq0\)（符号秩检验）或H0为X整体比Y偏右（秩和检验），则左侧检验：

   • \(p=P(X\leq Z)\)，拒绝域为\(Z \leq Z_{\alpha}\)

案例

1. 符号秩检验

R语言中，使用stats::wilcox.test函数执行Wilcoxon符号秩检验或Wilcoxon秩和检验，其中的参数含义如下：

• x：数值向量，表示第一组样本数据
• y：数值向量，表示第二组样本数据
  • 如果y未提供（NULL），则执行单样本Wilcoxon符号秩检验（检验x的中位数是否等于mu）
  • 如果y提供，且paired=FALSE（默认），执行独立样本Wilcoxon秩和检验（Mann-Whitney U检验）
  • 如果y提供，且paired=TRUE，执行配对样本Wilcoxon符号秩检验
• alternative：指定备择假设的方向，默认为two.sided，可选less、greater
• mu：单样本或配对检验中，假设的中位数差值（默认 0）
• paired：是否执行配对检验，默认为F（单样本符号秩检验或多样本秩和检验），可选T（配对符号秩检验）
• exact：是否计算精确的p值，默认为 NULL（自动选择），可选T强制使用，或F强制不使用
• correct：是否对统计量进行连续性修正（默认T）
• conf.int：是否计算中位数差值的置信区间，默认为FALSE，仅适用于配对检验或单样本检验（独立样本秩和检验不支持）
• conf.level：置信区间的置信水平（默认 0.95，即95%）
• tol.root：计算置信区间时的数值容差（默认 1e-4），一般无需调整
• digits.rank：秩次计算时的数值精度（默认 Inf），一般无需修改

注意：

• 在符号秩检验中，wilcox.test函数的近似解，考虑了ties。

• 使用Wilcoxon秩和检验时，虽然结果显示的是W，但其计算的是U统计量，\(U=W-\frac{m(m+1)}{2}\)

服用某药物前后，15个病人血压有变化如下：服药前c(121,117,107,123,128,107,107,112,100,122,102,115,115,115,129)，服药后c(105,98,97,118,100,111,106,93,100,95,92,103,99,97,107)，问药物对血压是否有影响。

# 设定显著性水平α=0.05
a <- 0.05

# 数据
x <- c(121,117,107,123,128,107,107,112,100,122,102,115,115,115,129)
y <- c(105,98,97,118,100,111,106,93,100,95,92,103,99,97,107)

# 函数计算.................................................................................
# H0：x=y，双侧检验
test.1 <- wilcox.test(x, y, paired = TRUE, alternative = "two.sided", exact = T)
# 结果如下，p-value < a，表明拒绝H0，认为x≠y，服药前后存在差异
#   Wilcoxon signed rank test with continuity correction
# data:  x and y
# V = 103, p-value = 0.001683
# alternative hypothesis: true location shift is not equal to 0

# H0：x≤y，右侧检验
test.2 <- wilcox.test(x, y, paired = TRUE, alternative = "greater", exact = F)
# 结果如下，p-value < a，表明拒绝H0，认为x>y，服药前血压更高
#   Wilcoxon signed rank test with continuity correction
# data:  x and y
# V = 103, p-value = 0.0008414
# alternative hypothesis: true location shift is greater than 0

# H0：x≥y，左侧检验
test.3 <- wilcox.test(x, y, paired = TRUE, alternative = "less", exact = F)
# 结果如下，p-value > a，表明不能拒绝H0，认为x≥y，服药前血压更高
#   Wilcoxon signed rank test with continuity correction
# data:  x and y
# V = 103, p-value = 0.9993
# alternative hypothesis: true location shift is less than 0

# 手工计算.................................................................................
# 数据
x <- c(121,117,107,123,128,107,107,112,100,122,102,115,115,115,129)
y <- c(105,98,97,118,100,111,106,93,100,95,92,103,99,97,107)
# 基础参数
d <- x-y
n  <- length(d)
i_nz <- d!=0                #记录d!=0的位置
d_valid <- d[i_nz]          #有效样本
n_valid <- length(d_valid)     
dp_p <- d_valid>0           #记录大于0、小于0样本位置 
dp_d <- d_valid<0 
W <- rank(abs(d_valid))
W_p <- sum(W[dp_p])         #计算W+

# 大样本近似计算
muW <- n_valid*(n_valid+1)/4                                                  #期望秩和
sigmaW <- sqrt(n_valid*(n_valid+1)*(2*n_valid+1)/24)                          #标准差
ties <- table(abs(d_valid))                                                   #记录ties
sigmaW_cor <- sqrt(n_valid*(n_valid+1)*(2*n_valid+1)/24-sum(ties^3-ties)/48)  #ties校正的标准差
Z_tcor_ccor <- (W_p-muW-0.5)/sigmaW_cor                     #带ties+连续性校正的统计量：3.141146
Z_ccor <- (W_p-muW-0.5)/sigmaW                              #带连续性校正的统计量：3.138824
Z_tcor <- (W_p-muW)/sigmaW_cor                              #带ties校正的统计量：3.172557
Z <- (W_p-muW)/sigmaW                                       #不带校正的统计量：3.170212

# 计算双侧检验p值
p1 <- (1-pnorm(abs(Z_tcor_ccor)))*2   #p1=0.001682881
p2 <- (1-pnorm(abs(Z_ccor)))*2        #p2=0.001696272
p3 <- (1-pnorm(abs(Z_tcor)))*2        #p3=0.001511027
p4 <- (1-pnorm(abs(Z)))*2             #p4=0.001523276

# 计算拒绝域临界值
thes <- qnorm(1-0.05/2)               #thes=1.959964

# 结果表明：p<a，Z>thres，因此拒绝H0，认为x≠y，服药前后存在差异

# 函数验证
# t1 <- wilcox.test(x,y,paired=TRUE,alternative="two.side",exact=F,correct=T) #p=0.001683
# t2 <- wilcox.test(x,y,paired=TRUE,alternative="two.side",exact=F,correct=F) #p=0.001511
# t1的结果同p1，t2的结果同p3，这说明wilcox.test中的标准差考虑了ties问题

2. 秩和检验

将上面例子的描述稍作修改（数据不变）：两个病人的血压分别为：A(121,117,107,123,128,107,107,112,100,122,102,115,115,115,129)，B(105,98,97,118,100,111,106,93,100,95,92,103,99,97,107)，问两人的血压分布是否相同。

# 设定显著性水平α=0.05
a <- 0.05

# 数据
x <- c(121,117,107,123,128,107,107,112,100,122,102,115,115,115,129)
y <- c(105,98,97,118,100,111,106,93,100,95,92,103,99,97,107)

# 函数计算.................................................................................
# H0：分布相同，双侧检验
test.1 <- wilcox.test(x, y, paired = F, alternative = "two.sided", exact = F)
# 结果如下，p-value < a，表明拒绝H0，认为两者分布不同
#   Wilcoxon rank sum test with continuity correction
# data:  x and y
# W = 199.5, p-value = 0.0003238
# alternative hypothesis: true location shift is not equal to 0

# H0：前者分布偏左，右侧检验
test.2 <- wilcox.test(x, y, paired = F, alternative = "greater", exact = F)
# 结果如下，p-value < a，表明拒绝H0，认为前者分布偏右
#   Wilcoxon rank sum test with continuity correction
# data:  x and y
# W = 199.5, p-value = 0.0001619
# alternative hypothesis: true location shift is greater than 0

# H0：前者分布偏右，左侧检验
test.3 <- wilcox.test(x, y, paired = F, alternative = "less", exact = F)
# 结果如下，p-value > a，表明不能拒绝H0，认为前者分布偏右
#   Wilcoxon rank sum test with continuity correction
# data:  x and y
# W = 199.5, p-value = 0.9999
# alternative hypothesis: true location shift is less than 0

# 手工计算.................................................................................
# 数据
x <- c(121,117,107,123,128,107,107,112,100,122,102,115,115,115,129)
y <- c(105,98,97,118,100,111,106,93,100,95,92,103,99,97,107)
m <- length(x)
n <- length(y)
xy <- c(x, y)
W <- sum(rank(xy)[1:m])    #求秩和
muW <- m*n/2
sigmaW <- sqrt(m*n*(m+n+1)/12)
U1 <- sum(outer(x, y, ">")) + 0.5 * sum(outer(x, y, "=="))   #wilcox.test中默认>
U2 <- W-m*(m+1)/2                                            #wilcox.test中的W其实是U2来的 
muU <- m*n/2
sigmaU <- sqrt(m*n*(m+n+1)/12)
ties <- table(xy)                      #记录ties
sigma_cor <- sqrt(m*n*(m+n+1)/12-sum((ties^3-ties)/(m+n)/(m+n-1))*m*n/12)
Z <- (U2-muU)/sigma_cor

# 计算双侧检验p值
p <- (1-pnorm(abs(Z)))*2              #p=0.0002989101

# 计算拒绝域临界值
thes <- qnorm(1-0.05/2)               #thes=1.959964

# 结果表明：p<a，Z>thres，因此拒绝H0，认为x≠y，服药前后存在差异

# 函数验证
# t1 <- wilcox.test(x,y,paired=F,alternative="two.sided",exact=F,correct=F) #p=0.0002989
# t1的结果同p，这说明wilcox.test中的标准差考虑了ties问题
# 注意这组数据不能使用函数计算时，不能correct=T，汇报错，而手算的话没有问题

5.7 Kruskal-Wallis检验

目的

• 检验多个（≥2）独立样本：是否来自相同的总体，或者说这些样本在中位数或总体位置上是否相同

适用

• 检验类型：非参数检验（多组独立样本比较，无需正态分布）
• 数据类型：数值数据或序数数据
• 数据体量：小样本或多组比较（ANOVA的非参数替代）

计算公式

• 统计量\(H=\frac{12}{N(N+1)}\sum_{i=1}^k\frac{R_i^2}{n_i}-3(N+1)\)

• ties修正统计量\(H_{cor}=\frac{H}{1-\frac{\sum_j{(t_j^3-t_j)}}{N^3-N}}\)

• 自由度\(df=k-1\)

其中：

• \(k\)，组数

• \(N\)，所有样本数（将多组数据汇总为一个数据后的）

• \(n_i\)，第i组的样本书来给你

• \(R_i\)，第i组的总秩和（多组数据汇总后求秩，然后将秩按组提取）

• \(t_1,t_2,...,t_j\)为ties，对应合并数据中非重复值对应的频数

结果判定

统计量H近似服从自由度为k-1的卡方分布，设定显著性水平为α，检验结果判定步骤如下（Kruskal-Wallis检验原假设只有≠一种形式）：

H0：所有组的总体分布相同（当已知同分布时，为中位数全相等）vs. H1：至少有一个组的总体分布与其他不同（已知同分布时，中位数不全相等）

1. \(p=P(X\geq H)\)

2. 拒绝域\(H \geq \chi^2_{1-\alpha,df}\)

案例

R语言中，使用stats::kruskal.test函数执行Kruskal-Wallis检验。

比较三种不同教学方法（方法A、B、C）对学生考试成绩的影响，收集的数据为各方法下学生的考试成绩，这些成绩可能不满足正态分布，所以不能用ANOVA，而采用Kruskal-Wallis检验。

A:c(70,75,80,65,78,73,76,81), B:c(88,85,90,87,91,80,79), C:c(60,63,58,66,62,90,88)

# 加载包
library(magrittr)

# 设定显著性水平α=0.05
a <- 0.05

# 数据
A <- c(70,75,80,65,78,73,76,81)
B <- c(88,85,90,87,91,80,79)
C <- c(60,63,58,66,62,90,88)

dat <- data.frame(
  score = c(A, B, C),
  method = factor(c(rep("A",length(A)),rep("B",length(B)),rep("C",length(C))))
)

# 函数计算.................................................................................
test.1 <- kruskal.test(score ~ method, data = dat)
# 结果如下，p-value < a，表明拒绝H0，认为不同方法下的成绩有差异
#   Kruskal-Wallis rank sum test
# data:  score by method
# Kruskal-Wallis chi-squared = 7.9097, df = 2, p-value = 0.01916
# 手工计算.................................................................................
# 数据
A <- c(70,75,80,65,78,73,76,81)
B <- c(88,85,90,87,91,80,79)
C <- c(60,63,58,66,62,90,88)
df <- 3-1
nA <- length(A)
nB <- length(B)
nC <- length(C)
N <- nA+nB+nC
ABC <- c(A,B,C)
RABC <- rank(ABC)
RA <- RABC[1:nA] %>% sum()
RB <- RABC[(nA+1):(nA+nB)] %>% sum()
RC <- RABC[(nA+nB+1):(nA+nB+nC)] %>% sum()
ties <- table(ABC)                                   #记录ties
H <- 12/N/(N+1)*(RA^2/nA+RB^2/nB+RC^2/nC)-3*(N+1)
Hcor <- H/(1-sum(ties^3-ties)/(N^3-N))               #Hcor=7.909679

# 计算双侧检验p值
p <- 1-pchisq(Hcor, df)              #p=0.01916174

# 计算拒绝域临界值
thes <-  qchisq(1-0.05, df)          #thes=5.991465

# 结果表明：p<a，Z>thres，因此拒绝H0，认为不同方法下的成绩有差异

5.8 Kolmogorov-Smirnov检验

又叫KS检验

目的

一种基于累积分布函数（CDF）差异的检验方法，主要用于：

• 检验单样本：是否来自某个特定的理论分布（如正态分布、指数分布等）

• 检验两个独立样本：否来自同一分布

注意：

• 数据分布参数由样本估计（如正态分布的均值和方差），需使用Lilliefors修正

• 对离散数据或高维数据，可考虑卡方检验或Anderson-Darling检验

适用

• 检验类型：非参数检验（比较分布函数，无参数假设）。
• 数据类型：数值数据（虽然也有对离散变量的扩展方法）
• 数据体量：小样本或大样本均可（但大样本对微小差异更敏感，易拒绝H0）。

计算公式

单样本统计量：\(D_m=sup_x|F_m(x)-F(x)|\)

两样本统计量：\(D_{m,n}=sup_x|F_{1,m}(x)-F_{2,n}(x)|\)

其中：

1. 原数据应先进行升序排列：\(x_1,x_2,x_3,...,x_n\)升序排列为\(x_{(1)},x_{(2)},x_{(3)},...,x_{(n)}\)

2. 样本检验分布函数ECDF：\(F_m(x)=\frac{样本中小于等于x的数据点的数量}{m}\)

3. 参考分布函数CDF：\(F(x)\)

结果判定

H0：数据服从给定分布（单样本）或H0：两数据分布相同（双样本）

通常靠查表计算拒绝域临界值、蒙特卡洛模拟或大样本下的近似公式：\(p\approx2exp(-2nD^2)\)

案例

R语言中，使用stats::ks.test函数执行KS检验，其中参数如下：

• x：第一组数据（数值向量），用于造经验分布函数
• y：第二组数据（双样本检验时提供），或者是字符串，如”pnorm”、“pexp” 等（用于单样本检验时）
• alternative：备择假设的方向，默认”two.sided”，可选”less”或”greater”，经验分布大于理论分布的关系
• exact：是否使用精确 p 值计算，TRUE对双样本、样本较小、且没有重复值时有效，FALSE：使用渐进方法（大样本近似）；NULL（默认），自动判断是否可以使用精确法
• simulate.p.value：是否使用模拟（bootstrap）方式来计算p值，适用于样本小、重复值多、exact不可用时
• B：当simulate.p.value=TRUE 时，指定模拟的重复次数，默认是 B = 2000

案例1，单样本检验：检验样本是否服从特定参数的正态分布

# 设定显著性水平α=0.05
a <- 0.05

# 生成样本数据 
set.seed(123)
sample_data <- rnorm(100, mean = 0, sd = 1) + rt(100, 3)

# 执行KS检验（H0: 样本服从均值为3方差为1的正态分布，注意，如果不指定，则默认mean=0sd=1即标准正态分布）
test.1 <- ks.test(sample_data, "pnorm", mean = 3, sd = 1)
# 结果如下：证明p-value < a，说明拒绝原假设，认为数据不服从均值为3方差为1的正态分布
#   Asymptotic one-sample Kolmogorov-Smirnov test
# data:  sample_data
# D = 0.81903, p-value < 2.2e-16
# alternative hypothesis: two-sided

案例2，两样本检验：判断两个样本是否来自同一个分布

# 设定显著性水平α=0.05
a <- 0.05

# 生成两个样本
set.seed(123)
sample1 <- rnorm(100, mean = 0, sd = 1)
sample2 <- rnorm(100, mean = 1, sd = 1)   

# 执行双样本KS检验
test.2 <- ks.test(sample1, sample2)
# 结果如下：证明p-value < a，说明拒绝原假设，认为两数据不来自于同一主体（两数据分布不同）
#   Asymptotic two-sample Kolmogorov-Smirnov test
# data:  sample1 and sample2
# D = 0.37, p-value = 2.267e-06
# alternative hypothesis: two-sided

5.9 Lilliefors检验

又叫LI检验，是Kolmogorov-Smirnov（KS）检验的改进版本

目的

基于经验分布函数（EDF）与假设的理论分布函数：

• 检验某个样本：否来自某个理论分布（尤其是正态分布），评估数据对分布的偏离程度

适用

• 检验类型：非参数检验（改进的KS检验，自动估计正态分布参数）
• 数据类型：数值数据（连续型）
• 数据体量：小样本或大样本均可（比KS检验更适用于正态性检验）
• 与SW检验相比，LI检验对非正态性更敏感，尤其在尾部偏离时表现更好，但Shapiro-Wilk检验在小样本中功效更高
• 与K-S检验相比，而LI检验检分布时，不需要制定理论分布的参数（如均值方差），直接使用样本估计的均值和方差，因此更适用于参数未知的情况

计算公式

1. 估计样本均值和标准差：\(\hat \mu = \frac{1}{n}\sum_{i=1}^n{x_i}\)，\(\hat\sigma=\sqrt{\frac{1}{n}\sum_{i=1}^{n}{(x_i-\hat\mu)^2}}\)

2. 样本标准化：\(z_i=\frac{x_i-\hat\mu}{\hat\sigma}\)

3. 经验分布函数（EDF排序后数据）：\(F_n(z_{(i)})=\frac{i}{n}\)

4. 理论分布函数（CDF）：\(\Phi(z_i)=\int_{-\infty}^{z_i}{\frac{1}{\sqrt{2\pi}}e^{-t^2/2}dt}\)

5. 统计量D：\(D=max(|F_n(z_i)-\Phi(z_i)|, |F_n(z_i)-\Phi(z_{i-1})|)\)

结果判定

H0：数据服从正态分布 vs. H1：数据不服从正态分布

通常靠查表或者使用软件计算

案例

在R语言中，通过nortest::lillie.test执行LI检验，该函数只有一个参数x，代表被测数据

# 设定显著性水平α=0.05
a <- 0.05

# 加载包
library(nortest)

# 生成样本数据 
set.seed(123)
sample_data <- rnorm(100, mean = 31, sd = 1) 

# 执行LI检验（H0: 样本服从正态分布）
test.1 <- lillie.test(sample_data)
# 结果如下：证明p-value > a，说明不拒绝原假设，认为数据服从正态分布
#   Lilliefors (Kolmogorov-Smirnov) normality test
# data:  sample_data
# D = 0.058097, p-value = 0.5575

5.10 Anderson-Darling检验

又叫AD检验，是Kolmogorov-Smirnov（KS）检验的改进版本

目的

• 检验某个样本：是否来自某个特定分布（常用正态分布）

• 本质是比较样本的经验分布函数（EDF）与理论分布函数（CDF）之间的差异，但对两端（尾部）的偏离赋予更高的权重。

适用

• 检验类型：非参数检验

• 数据类型：数值数据（连续型）

• 数据体量：适用于中小样本，也可用于大样本

计算公式

1. 原数据应先进行升序排列：\(x_1,x_2,x_3,...,x_n\)升序排列为\(x_{(1)},x_{(2)},x_{(3)},...,x_{(n)}\)

2. 统计量\(A^2=-n-\frac{1}{n}\sum_{i=1}^n{[(2i-1)(lnF(x_{(i)})+ln(1-F(x_{(n+1-i)})))]}\)

3. 其中\(F(x)\)为理论CDF函数

结果判定

H0：样本来自指定的理论分布（例如正态分布）

通常靠查表或者使用软件计算

案例

在R语言中，通过nortest::ad.test执行AD检验，该函数只用于检验数据是否符合正态分布，且只有一个参数x，代表被测数据

# 设定显著性水平α=0.05
a <- 0.05

# 加载包
library(nortest)

# 生成样本数据 
set.seed(123)
sample_data <- rnorm(100, mean = 31, sd = 1) 

# 执行AD检验（H0: 样本服从正态分布）
test.1 <- ad.test(sample_data)
# 结果如下：证明p-value > a，说明不拒绝原假设，认为数据服从正态分布
#   Anderson-Darling normality test
# data:  sample_data
# A = 0.182, p-value = 0.9104

5.11 Shapiro-Wilk检验

又叫SW检验

目的

通过比较样本数据与理想正态数据的分位数相关性：

• 检验某个样本：否服从正态分布

适用

• 检验类型：参数检验（专用于正态性检验，依赖样本参数估计）
• 数据类型：数值数据（连续型）
• 数据体量：小样本n<50 时最优，大样本可能过于敏感（若数据离散或含大量重复值，检验功效可能降低）
• 分布对称性：对非对称分布（如偏态分布）的检测效果优于其他正态性检验（如K-S检验）

计算公式

统计量\(W=\frac{(\sum_{i=1}^n{a_ix_{(i)}})^2}{\sum_{i=1}^n{(x_i-\bar x})^2}\)

其中：\(x_{(i)}\)，样本数据排序后的第个数据；\(a_i\)，有正态分布期望值生成的系数

结果判定

H0：数据服从正态分布 vs. H1：数据不服从正态分布

W值越接近1，样本越可能服从正态分布；显著小于1则偏离正态

案例

R语言中，使用stats::shapiro.test执行SW检验，只有一个参数x（可以添加均值方差），要求是x中不能有缺失值、长度在3到5000之间

# 设定显著性水平α=0.05
a <- 0.05

# 生成样本数据 
set.seed(123)
sample_data <- rnorm(100, mean = 31, sd = 1) 

# 执行SW检验（H0: 样本服从正态分布）
test.1 <- shapiro.test(sample_data)
# 结果如下：证明p-value > a，说明不拒绝原假设，认为数据服从正态分布
#   Shapiro-Wilk normality test
# data:  sample_data
# W = 0.99626, p-value = 0.9954

5.12 方差分析

英文名Analysis of Variance，ANOVA

目的

通过方差分解，检验多组数据的：组间均值是否相同

具体来说，ANOVA将总体方差分解为组间均方（Mean Square Between, MSB）和组内均方（Mean Square Within, MSW）。其中，组间均方MSB衡量各组均值相对整体均值的离散程度，而组内均方MSW描述组内数据相对组内均值的离散程度（很多情况下，组内差异也被称为残差，Residuals）。这样，如果MSB很大，说明各组均值差异显著，而如果MSW很大，说明组内本身分布不均匀。而ANOVA就是在判断，当前观察到的“组间偏差”，是否远大于可以归因于“随机组内波动”的程度？换句话说，你们之间的差异，是因为你们本来就不同，还是只是偶然的观测波动。

注意：

• ANOVA的核心是F检验，但其前提、后续涉及多个检验方法，这里重点是在案例中演示R语言进行ANOVA的详细流程
• ANOVA的本质是线性回归的特例，即当自变量全是分类变量时，回归模型等价于方差分析
• 可以先做回归再做方差分析，例如rstatix::anova_test的例子anova_test(lm(yield ~ block + N*P*K, npk))，类似的，可利用回归分析+ANOVA，处理混合变量类型（连续+分类）、非平衡设计（各组样本量不等）、复杂模型结构（随机效应、非线性项）等问题

此外，还有很多扩展的方差分析，这里不做一一介绍：

1. 重复测量方差分析：让同一组受试者在不同时间点或条件下进行多次测量时使用
2. 混合设计方差分析：结合了“重复测量方差分析”和“多因素方差分析”的特点，适用于有重复测量因素和独立因素的混合设计
3. 方差分析与加权法：为不同组别赋予不同权重的方差分析，常用于组别大小不一致的情况
4. 随机效应方差分析：与固定效应模型不同，随机效应模型假设因子水平是随机选取的，且试图分析这些随机因素对因变量的影响
5. 多重响应变量的方差分析：MANOVA是ANOVA的扩展，允许同时分析多个相关的因变量，通过多元正态分布来评估组间的差异
6. 路径分析与结构方程模型中的方差分析：在结构方程模型（SEM）中，方差分析的思想被用于分析因变量与自变量之间的关系，尤其是在多因子模型中
7. 协方差分析（ANCOVA）：基础的ANOVA自变量都是分类型，如果额外加入数值型自变量作为协变量（只是限制因素，并不参与结果分析），则为协方差分析

适用

• 检验类型：参数检验，要求正态性、方差齐性
• 数据类型：因变量为数值型，自变量为分类型（2组或多组）
• 数据体量：建议各组样本量 ≥ 5（多因素时还包括交互组），总体样本量 ≥ 20

注意，常规ANOVA，以及不满足条件时的下替代方案，如下表所示：

方法	正态?	方差齐次?	多因素?	交互?
1.ANOVA	✅	✅	✅	✅
2.Welch’s ANOVA	✅	❌	❌	❌
3.非参数检验	❌	❌	❌	❌
4.非参数多因素扩展	❌	❌	✅	✅
5.稳健单因素ANOVA	❌	❌	❌	❌
6.稳健多因素ANOVA	❌	❌	✅	✅

其中：

• 单因素Welch’s ANOVA，可用于方差非齐次、未知、甚至齐次（此时结果与普通ANOVA相接近）
• 多因素时，每个因子的不同水平组合都可能拥有不同的方差，且主效应、交互效应不再独立，异方差会扭曲交互结构
• 除方差分析或非参数检验，还可以用广义最小二乘法（nlme::gls + varIdent）或者文件线性模型（lm + sandwich::vcovHC）来进行非正态、非齐次下的多因素、交互效应对比

计算公式（单因素方差分析）

• 统计量：\(F=\frac{MSB}{MSW}=\frac{SSB/(k-1)}{SSW/(N-k)}\)

• 组间均方（Mean Square Between, MSB）：\(MSB=\frac{\sum^{k}_{i=1}{n_i(\bar X_i-\bar X)^2}}{k-1}\)

• 组内均方（Mean Square Within, MSW）：\(MSW=\frac{\sum_{i=1}^{k}{\sum_{j=1}^{n_i}{(X_{ij}-\bar X_i)^2}}}{N-k}\)

• 自由度：组间（分子）\(df_1=k-1\)，组内（分母）\(df_2=N-k\)

• \(k\)，组数：\(n_i\)，第i组的样本量；\(N=\sum n_i\)，总样本量；\(\bar X_i\)，第i组的样本均值，\(\bar X\)，总体均值

注意：

• 计算MSB时，乘以\(n_i\)是为了反映每组在总体中的代表性或权重，即组的样本量越大，它的均值对整体的贡献越大，偏离总体均值所造成的变异也越大。

结果判定

统计量F服从自由度为\(df=(k-1,N-k)\)的F分布。设定显著性水平为α后，结果判定步骤如下：

1. 方差分析只有双侧检验，即H0: \(\mu_1=\mu_2=...=\mu_k\) vs. H1：\(\mu_i\)不全相等

2. \(p=2 \times min(P(X \geq F), P(X \leq F))\)，拒绝域为\(F \leq F_{\alpha /2}(k-1,N-k)\) 或 \(F \geq F_{1-\alpha/2}(k-1,N-k)\)

3. 如果\(p<\alpha\)或F落在拒绝域内，则拒绝H0

案例

案例方法按如下表格：

方法	正态?	方差齐次?	多因素?	交互?
案例1.ANOVA	✅	✅	✅	✅
案例2.Welch’s ANOVA	✅	❌	❌	❌
案例3.非参数检验	❌	❌	❌	❌
案例4.非参数多因素扩展	❌	❌	✅	✅
案例5.稳健单因素ANOVA	❌	❌	❌	❌
案例5.稳健多因素ANOVA	❌	❌	✅	✅

总结：ANOVA流程与R函数：

1. 正态性检验

   • Shapiro-Wilk检验（推荐）：satats::shapiro.test
   • Lilliefors检验： nortest::lillie.test
   • Anderson-Darling检验： nortest::ad.test

2. 方差齐次性检验

   • Bartlett检验，要求数据正态分布： stats::bartlett.test
   • Levene检验，适用于非正态数据： rstatix::levene_test或car::leveneTest
   • Fligner-Killeen检验，基于中位数，不依赖正态性：stats::fligner.test

3. ANOVA（或扩展的组间差异显著性判别方法）：

   • base包，支持单因素、多因素，结果为aov型： stats::aov
   • rstatix包，tidyverse风格，结果为df型，还支持方差修正、协方差、重复测量等： rstatix::anova_test
   • Welch’s ANOVA，数据正态、方差不齐次、仅单因素： stats::oneway.test
   • Kruskal-Wallis检验，非参数，仅单因素： stats::kruskal.test
   • Scheirer Ray Hare检验，非参数，双因素： rcompanion::scheirerRayHare
   • 稳健ANOVA，小样本、非正态、方差非齐次、一/二/三因素：WRS2::t1way、WRS2::t2way和WRS2::t3way

4. 事后检验（即多重比较，在组间差异显著时，进一步判断哪些组之间的差异显著）

   • TukeyHSD检验（又叫Tukey-Kramer检验），常规ANOVA，支持单因素或多因素（推荐）： stats::TukeyHSD
   • LSD-t检验，常规ANOVA，仅单因素： PMCMRplus::lsdTest
   • Dunnett-t检验，常规ANOVA，仅单因素，只比较对照组： PMCMRplus::dunnettTest
   • SNK-q检验，常规ANOVA，仅单因素： PMCMRplus::snkTest
   • Games-Howell检验，Welch ANOVA：rstatix::games_howell_test（推荐）或PMCMRplus::gamesHowellTest
   • Dunn检验，非参数，仅单因素（组数≤2报错）或主效应： FSA::dunnTest
   • Conover-Iman检验，非参数，仅单因素或主效应： PMCMRplus::kwAllPairsConoverTest
   • DSCF(Dwass–Steel–Critchlow–Fligner)检验，非参数，仅单因素：PMCMRplus::dscfAllPairsTest
   • 对齐秩变换，非参数，仅*型交互效应： ARTool::art+ARTool::art.con
   • 稳健ANOVA后的线性对比，仅单因素： WRS2::lincon
   • 注意非常规、多因素时，只能用Dunn检验或Conover-Iman检验指定的主效应+用对齐秩变换检验交互效应

5. 可视化

   • 基础绘图，可自定义选项多： library(ggplot2)
   • ggplot2扩展，显著性标记更方便： library(ggpubr)

注意：

1. 部分案例中用到了效应量\(\eta^2\)分析，其数值大小的解释方法

   效应量\(\eta^2\)范围	效应大小	直观解释
   0.01 ≤ \(\eta^2\) < 0.06	小效应	效应存在但可能无实际意义
   0.06 ≤ \(\eta^2\) < 0.14	中等效应	理论或实践中值得关注的效应
   \(\eta^2\) ≥ 0.14	大效应	明显的重要效应

案例1.1. 单因素ANOVA

R语言中，使用stats::aov执行基础ANOVA，其中参数包括：

• formula：模型形式，因变量_{自变量，例如y}A+B+A:B+C或（等价于）y~A*B+C，其中
  • A:B，表示只有A与B的交互效应，不含主效应
  • A*B，表示包含A主效应、B主效应、 A:B交互效应
• data：df型数据，其中包含formula中所用变量
• projections：是否保留模型投影矩阵，通常用于进一步的线性模型分析，默认F
• qr：是否保留QR分解结果（用于数值稳定性和计算效率，比如解正规方程），默认T
• contrasts：控制分类变量的编码方式，默认使用 options(“contrasts”) 中的设置

# 设定显著性水平α=0.05
a <- 0.05

# 生成四组数据
set.seed(123)
groupA <- rnorm(30, mean = 50, sd = 5)
groupB <- rnorm(30, mean = 55, sd = 5)
groupC <- rnorm(30, mean = 60, sd = 5)
groupD <- rnorm(30, mean = 64, sd = 5)

# 数据框整理
dat <- data.frame(
  value = c(groupA, groupB, groupC, groupD),
  group = factor(rep(c("A", "B", "C", "D"), each = 30))
)

# 单因素ANOVA..............................................................................
# 1. 正态性检验（每组），
# 使用Shapiro-Wilk检验、Lilliefors检验、Anderson-Darling检验
normtest1 <- tapply(dat$value, dat$group, stats::shapiro.test) 
normtest2 <- by(dat$value, dat$group, nortest::lillie.test)  
normtest3 <- by(dat$value, dat$group, nortest::ad.test)  
# 结果如下，所有组p>a，表明全不拒绝H0，都服从正态分布
# > normtest1
# $A
#   Shapiro-Wilk normality test
# data:  X[[i]]
# W = 0.97894, p-value = 0.7966
# $B
#   Shapiro-Wilk normality test
# data:  X[[i]]
# W = 0.98662, p-value = 0.9614
# $C
#   Shapiro-Wilk normality test
# data:  X[[i]]
# W = 0.98085, p-value = 0.8478
# $D
#   Shapiro-Wilk normality test
# data:  X[[i]]
# W = 0.97106, p-value = 0.5685
# > normtest2
# $A
#   Lilliefors (Kolmogorov-Smirnov) normality test
# data:  X[[i]]
# D = 0.091059, p-value = 0.7587
# $B
#   Lilliefors (Kolmogorov-Smirnov) normality test
# data:  X[[i]]
# D = 0.097828, p-value = 0.6559
# $C
#   Lilliefors (Kolmogorov-Smirnov) normality test
# data:  X[[i]]
# D = 0.079714, p-value = 0.8968
# $D
#   Lilliefors (Kolmogorov-Smirnov) normality test
# data:  X[[i]]
# D = 0.10226, p-value = 0.5864
# > normtest3
# dat$group: A
#   Anderson-Darling normality test
# data:  dd[x, ]
# A = 0.19951, p-value = 0.873
# -------------------------------------------------------------------------- 
# dat$group: B
#   Anderson-Darling normality test
# data:  dd[x, ]
# A = 0.20978, p-value = 0.8471
# -------------------------------------------------------------------------- 
# dat$group: C
#   Anderson-Darling normality test
# data:  dd[x, ]
# A = 0.22672, p-value = 0.7984
# -------------------------------------------------------------------------- 
# dat$group: D
#   Anderson-Darling normality test
# data:  dd[x, ]
# A = 0.31171, p-value = 0.5322

# 2. 方差齐性检验，使用Bartlett检验或Levene检验或Fligner-Killeen检验
homogvar1 <- stats::bartlett.test(value ~ group, data = dat)#Bartlett检验，要求数据正态分布
homogvar2 <- rstatix::levene_test(value ~ group, data = dat)#Levene检验，适用于非正态数据
homogvar3 <- car::leveneTest(value ~ group, data = dat)     #同上
homogvar4 <- stats::fligner.test(value ~ group, data = dat) #Fligner-Killeen检验，基于中位数，不依赖正态
# 结果如下，所有方法p>a，表明不能拒绝H0，认为满足方差齐次性
# > homogvar1 
#   Bartlett test of homogeneity of variances
# data:  value by group
# Bartlett's K-squared = 0.82355, df = 3, p-value = 0.8438
# > homogvar2
# # A tibble: 1 × 4
#     df1   df2 statistic     p
#   <int> <int>     <dbl> <dbl>
# 1     3   116     0.418 0.741
# > homogvar3
# Levene's Test for Homogeneity of Variance (center = median)
#        Df F value Pr(>F)
# group   3  0.4176 0.7407
#       116 
# > homogvar4
#   Fligner-Killeen test of homogeneity of variances
# data:  value by group
# Fligner-Killeen:med chi-squared = 1.5412, df = 3, p-value = 0.6728

# 3. 满足前置条件，进行单因素ANOVA
model <- aov(value ~ group, data = dat)
# 结果如下，p<a，表明拒绝原假设，认为各组均值不相等
# > summary(model)
#              Df Sum Sq Mean Sq F value Pr(>F)    
# group         3   3166  1055.5   52.14 <2e-16 ***
# Residuals   116   2348    20.2                   
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 对结果的解释：
# + 自由度表明，组数为4，样本总数为120（3=4-1,116=120-4）
# + Residuals表示残差，在ANOVA中，指组内的差异
# + 组间：平方和Sum Sq为3166.455，除以自由度DF，等于Mean Sq1055.5
# + 组内：平方和为2348.143，自由度116，均方为20.2
# + F统计量为52.14，p值小于2e-16，说明在0.05显著性水平下，拒绝原假设，认为各组均值不全相等
# + 至此可以判断至少有一组均值不等，但具体是哪一组，需进一步用多重比较进行检查

# 4. 事后检验（即多重比较），其作用是在ANOVA拒绝H0的情况下，检验哪些组之间存在显著差异
# 可选以下几种方法
# + TukeyHSD检验 (Tukey's Honest Significant Difference Test) （推荐）
# + LSD-t检验 (Least Significant Difference Test)
# + Dunnett-t检验
# + SNK-q检验 (Student-Newman-Keuls Test)
PHT_result1 <- stats::TukeyHSD(model)               #TukeyHSD检验
PHT_result2 <- PMCMRplus::lsdTest(model)            #LSD-t检验
PHT_result3 <- PMCMRplus::dunnettTest(model)        #Dunnett-t检验（只比较对照组）
PHT_result4 <- PMCMRplus::snkTest(model)            #SNK-q检验
# 结果如下，表明四组两两之间均值都不相同
# > PHT_result1
#   Tukey multiple comparisons of means
#     95% family-wise confidence level
# Fit: aov(formula = value ~ group, data = dat)
# $group
#          diff        lwr       upr     p adj
# B-A  6.127210  3.0990917  9.155329 0.0000037
# C-A 10.357621  7.3295020 13.385739 0.0000000
# D-A 13.766074 10.7379554 16.794193 0.0000000
# C-B  4.230410  1.2022915  7.258529 0.0022695
# D-B  7.638864  4.6107449 10.666982 0.0000000
# D-C  3.408453  0.3803346  6.436572 0.0207466
# 对结果的解释
# + TukeyHSD的H0是，第i组和第j组的均值相等
# + 它计算服从q分布（studentized range distribution）的studentized range统计量
# + TukeyHSD检验，控制了全局错误率，适合多组事后比较，相比之下如果两两进行t检验，容易引起多重比较问题使得整体第I类错误率（假阳性）上升
# > PHT_result2
#   Pairwise comparisons using Least Significant Difference Test
# data: value by group
#   A       B       C      
# B 6.3e-07 -       -      
# C 8.2e-15 0.00041 -      
# D < 2e-16 1.5e-09 0.00403
# P value adjustment method: none
# alternative hypothesis: two.sided
# > PHT_result3
#   A      
# B 1.5e-06
# C 1.5e-14
# D < 2e-16
# > PHT_result4
#   A       B       C      
# B 6.3e-07 -       -      
# C 5.6e-14 0.00041 -      
# D 9.2e-15 4.4e-09 0.00403

# 5. 结果可视化，使用ggpubr包简化绘图
# 加载包
library(ggplot2)
library(ggpubr)
vars <- unique(dat$group) %>% as.vector()      
varscomb <- combn(vars,2,simplify=F)
ggviolin(dat, x = "group", y = "value", fill = "group",
         palette = c("#00AFBB", "#E7B800", "#FC4E07", "#6EA8FD"), alpha=0.5,
         add = "boxplot", add.params = list(fill = "white"))+
  # 添加显著性标记，方法t.test、wilcox.test、anova等，comparisons时不能用anova
  stat_compare_means(
    comparisons = varscomb,
    method = "wilcox.test",                                     
    symnum.args = list(cutpoints = c(0, 0.001, 0.01, 0.05, Inf),
                       symbols = c("***", "**", "*", "NS")),
    label = "p.signif")+
  # 添加全局P值标记，没有comparisons是全局的
  stat_compare_means(label.y = 30, method = "anova")+               
  theme(legend.position = "none")

案例1.2. 多因素ANOVA

R语言中，还可以使用rstatix::anova_test执行方差分析，该函数更符合tidyverse编程风格，可执行

• 独立测量ANOVA：受试者间设计
• 重复测量ANOVA：受试者内设计
• 混合ANOVA：受试者内部和受试者之间的混合设计，也称为分裂图方差分析
• ANCOVA：协方差分析。

rstatix::anova_test中的参数如下：

1. data：data.frame形式的数据集

2. formula：公式，与stats::aov中的用法相同

3. dv：字符型，因变量（Dependent Variable），即要分析的连续型变量。在formula中已经可以指定，通常可以不用这个参数

4. wid：字符型，在重复测量设计中，表示每个观测的受试者ID，即每个观测的唯一标识

5. between：字符型，在重复测量设计中，用于指定组间因素

6. within：字符型，在重复测量设计中，用于指定组内因素（即被试内设计，如每个受试者接受不同的处理）

7. covariate：字符型，用于指定协变量，用来进行协方差分析（ANCOVA）

8. type：整型或字符型，用于指定ANOVA方差类型，默认值NULL，根据数据自动选择。其中有3种类型：

   • Type I（顺序方差分析，通常用于平衡设计）
   • Type II（平衡与不平衡设计均适用，常用于比较不同的主效应）
   • Type III（最常用，适用于不平衡设计，考虑所有主效应和交互效应）

9. effect.size：字符型，指定效应大小的计算方法，默认ges，即广义效应大小，可选其他计算方法，如eta, omega等

10. error：字符型，用于指定误差类型，例如用来调整误差项，常用于多重比较。比如：选择 “sphericity” 来检验球形性

11. white.adjust：逻辑值，默认FALSE，如果TRUE，将应用White’s修正来调整误差项，以适应异方差性问题

12. observed：逻辑值，是否显示观察值，默认是FALSE，如果TRUE，则会显示一些统计值，如均值等

13. detailed：逻辑值，是否显示详细输出，默认是FALSE，如果TRUE会提供更多的输出内容，如每个层级的均值和标准误等

# 加载包
library(data.table)
library(rstatix)

# 设定显著性水平α=0.05
a <- 0.05

# 数据来自datasets
dat <- ToothGrowth 
dat$dose <- factor(dat$dose)
dat$supp <- factor(dat$supp)

# 1. 正态性检验，ANOVA实际上要求误差项服从正态分布，因此这里使用多元回归+Shapiro-Wilk检验，考虑了所有预测变量，比分组检验对样本量要求更少，也更高效
model <- lm(len ~ supp * dose, data = ToothGrowth)
normtest <- stats::shapiro.test(residuals(model))
# 结果如下，表明p>a，不拒绝原H0，认为符合正态分布
#   Shapiro-Wilk normality test
# data:  residuals(model)
# W = 0.98162, p-value = 0.5008

# 2. 方差齐次性检验，使用Bartlett检验或Levene检验
hv_test1  <- stats::bartlett.test(len ~ interaction(supp, dose), data = dat)
hv_test2  <- car::leveneTest(len ~ supp * dose, data = dat)
# 结果如下，表明p>a，不拒绝原H0，认为符合方差齐次性要求
# > hv_test1
#   Bartlett test of homogeneity of variances
# data:  len by interaction(supp, dose)
# Bartlett's K-squared = 6.9273, df = 5, p-value = 0.2261
# > hv_test2
# Levene's Test for Homogeneity of Variance (center = median)
#       Df F value Pr(>F)
# group  5  1.7086 0.1484
#       54       

# 3. 满足前置条件，进行多因素ANOVA
model <- rstatix::anova_test(dat, len ~ supp * dose, detailed = T) %>%
  rstatix::get_anova_table() %>%
  setDT() %>% 
  .[, R2:=SSn/(SSn + SSd)]         #计算各效应的R²
# 结果如下，p<a，表明拒绝原假设，认为各组均值不相等
# > model
#       Effect      SSn     SSd   DFn   DFd      F        p  p<.05   ges        R2
#       <char>    <num>   <num> <num> <num>  <num>    <num> <char> <num>     <num>
# 1:      supp  205.350 712.106     1    54 15.572 2.31e-04      * 0.224 0.2238254
# 2:      dose 2426.434 712.106     2    54 92.000 4.05e-18      * 0.773 0.7731092
# 3: supp:dose  108.319 712.106     2    54  4.107 2.20e-02      * 0.132 0.1320279
# 对结果的解释：
# + SSn组间平方和，SSd组内（残差）平方和，DFn分子（组间）自由度，DFd分母（组内/残差）自由度，ges广义η2
# + 自由度表明，supp、dose各有2组，自由度2-1，supp:dose有6组，自由度(2-1)*(3-1)，残差自由度60-(1+2+2)=54
# + 三个效应的p值都小于0.05，说明拒绝原假设，其中的组间均值都存在显著差异
# + ges与R2的值相同，代表了各组自变量对于因变量的解释程度，如supp中，VC和OJ对牙齿长度的影响差异显著，解释了22.4%的变异
 
# 4. 事后检验（即多重比较）
# 注意，前面提到的四个，除了TukeyHSD适用于多因素，其余的如LSD-t检验、Dunnett-t检验、SNK-q检验都不好用
PHT_result1 <- stats::TukeyHSD(aov(len ~ supp * dose, data=dat))             #TukeyHSD检验
# 结果如下，结果表明，supp和dose中的各组之间均值都有显著差异，交互效应中，有四组不显著
#   Tukey multiple comparisons of means
#     95% family-wise confidence level
# Fit: aov(formula = len ~ supp * dose, data = dat)
# $supp
#       diff       lwr       upr     p adj
# VC-OJ -3.7 -5.579828 -1.820172 0.0002312
# $dose
#         diff       lwr       upr   p adj
# 1-0.5  9.130  6.362488 11.897512 0.0e+00
# 2-0.5 15.495 12.727488 18.262512 0.0e+00
# 2-1    6.365  3.597488  9.132512 2.7e-06
# $`supp:dose`
#                diff        lwr        upr     p adj
# VC:0.5-OJ:0.5 -5.25 -10.048124 -0.4518762 0.0242521
# OJ:1-OJ:0.5    9.47   4.671876 14.2681238 0.0000046
# VC:1-OJ:0.5    3.54  -1.258124  8.3381238 0.2640208
# OJ:2-OJ:0.5   12.83   8.031876 17.6281238 0.0000000
# VC:2-OJ:0.5   12.91   8.111876 17.7081238 0.0000000
# OJ:1-VC:0.5   14.72   9.921876 19.5181238 0.0000000
# VC:1-VC:0.5    8.79   3.991876 13.5881238 0.0000210
# OJ:2-VC:0.5   18.08  13.281876 22.8781238 0.0000000
# VC:2-VC:0.5   18.16  13.361876 22.9581238 0.0000000
# VC:1-OJ:1     -5.93 -10.728124 -1.1318762 0.0073930
# OJ:2-OJ:1      3.36  -1.438124  8.1581238 0.3187361
# VC:2-OJ:1      3.44  -1.358124  8.2381238 0.2936430
# OJ:2-VC:1      9.29   4.491876 14.0881238 0.0000069
# VC:2-VC:1      9.37   4.571876 14.1681238 0.0000058
# VC:2-OJ:2      0.08  -4.718124  4.8781238 1.0000000

# 5. 结果可视化，使用ggpubr包简化绘图
# 注意stat_compare_means只能标记与x轴变量，这里只能用分组标记supp，但是没有线只有显著性，高度也只能统一
# 加载包
library(ggplot2)
library(ggpubr)
vars <- unique(dat$supp) %>% as.vector()          
varscomb <- combn(vars,2,simplify=F)
ggviolin(dat, x="dose", y="len", fill="supp", alpha=0.5,
         add="boxplot", add.params=list(alpha=0, show.legend=F)) + 
  stat_compare_means(
           aes(group = supp),                     #指定比较的是supp
           method = "wilcox.test",                #或t.test
           label = "p.signif",                    #显示 * ** ***
           label.y = max(dat$len) * 1.2,          #比dose标签更高
           symnum.args = list(
             cutpoints = c(0, 0.001, 0.01, 0.05, Inf),
             symbols = c("***", "**", "*", "NS")
           )
         )+
  # 添加全局P值标记，没有comparisons是全局的
  stat_compare_means(label.x=2, label.y = 0,method = "anova", hjust=0)+               
  theme(legend.position = "none")

案例2. Welch’s ANOVA

当数据满足正态性、不满足方差齐次性，并且想做单因素ANOVA模型时，可用Welch’s ANOVA方法

R语言中，使用stats::oneway.test执行基础Welch’s ANOVA，其中参数包括：

1. formula：公式，格式为 因变量 ~ 分组变量

2. data：数据框，指定包含公式中变量的数据框。

3. subset：子集，选择数据的子集进行分析，例如subset = Sepal.Length > 4.5 #只分析长度大于4.5的样本

4. na.action：缺失值处理，指定如何处理缺失值，默认为na.omit（删除缺失值），还可选na.fail（如果存在缺失值，直接报错）

5. var.equal：方差齐性假设，是否假设各组方差相等，默认为FALSE，即使用Welch’s ANOVA，可选T，经典ANOVA

# 设定显著性水平α=0.05
a <- 0.05

# 示例数据
dat <- iris[c("Sepal.Length","Species")]

# 1.正态性检验
norm_test <- by(dat$Sepal.Length, dat$Species, stats::shapiro.test)
# 结果如下，p>a，不拒绝H0，认为每组数据都服从正态分布
# > norm_test$setosa$p.value = 0.4595132
# > norm_test$versicolor$p.value = 0.464737
# > norm_test$virginica$p.value = 0.2583147

# 2. 方差齐次性检验，使用Bartlett检验或Levene检验
hv_test1 <- stats::bartlett.test(Sepal.Length ~ Species, data = dat)
hv_test2 <- car::leveneTest(Sepal.Length ~ Species, data = dat)
# 结果如下：
# hv_test1结果表明，p<a，拒绝H0，认为方差非齐次
#   Bartlett test of homogeneity of variances
# data:  Sepal.Length by Species
# Bartlett's K-squared = 16.006, df = 2, p-value = 0.0003345
# hv_test2结果表明，p<a，拒绝H0，认为方差非齐次
# Levene's Test for Homogeneity of Variance (center = median)
#        Df F value   Pr(>F)   
# group   2  6.3527 0.002259 **
#       147                    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# 3. 未满足方差齐次性，使用Welch’s ANOVA
anova_result <- stats::oneway.test(Sepal.Length ~ Species, data = dat, var.equal = FALSE)
# 结果如下：p<a，说明拒绝H0，认为各组的均值不全相等
#   One-way analysis of means (not assuming equal variances)
# data:  Sepal.Length and Species
# F = 138.91, num df = 2.000, denom df = 92.211, p-value < 2.2e-16

# 注意anova_test可以处理异方差，但是，其作用是通过White校正（异方差一致性协方差矩阵，HCCM）调整标准误，使得在方差不齐时仍能获得有效的p值，结果如下。可以看出，其与Welch’s ANOVA并不相同。
# anova_result2 <- anova_test(dat, Sepal.Length ~ Species, white.adjust = T)
# > anova_result2
# ANOVA Table (type II tests)
#    Effect DFn DFd       F        p p<.05
# 1 Species   2 147 137.114 2.49e-34     *

# 4. 事后检验（即多重比较）
# 函数TukeyHSD仅适用于aov的结果，这里使用针对方差非齐次的Games-Howell检验
# 常用rstatix::games_howell_test和PMCMRplus::gamesHowellTest，推荐前者
ght_rest1 <- rstatix::games_howell_test(dat, Sepal.Length ~ Species)
ght_rest2 <- PMCMRplus::gamesHowellTest(Sepal.Length ~ Species, data = dat)
# 结果如下，表明四组两两之间均值都不相同
# > ght_rest1
# A tibble: 3 × 8
#   .y.          group1     group2     estimate conf.low conf.high    p.adj p.adj.signif
# * <chr>        <chr>      <chr>         <dbl>    <dbl>     <dbl>    <dbl> <chr>       
# 1 Sepal.Length setosa     versicolor    0.93     0.719     1.14  2.86e-10 ****        
# 2 Sepal.Length setosa     virginica     1.58     1.34      1.83  0        ****        
# 3 Sepal.Length versicolor virginica     0.652    0.376     0.928 5.58e- 7 ****  
# > ght_rest2
#   Pairwise comparisons using Games-Howell test
# data: Sepal.Length by Species
#            setosa  versicolor
# versicolor 2.9e-10 -         
# virginica  < 2e-16 5.6e-07   
# P value adjustment method: none
# alternative hypothesis: two.sided
 
# 5. 结果可视化，使用ggpubr包简化绘图，注意ggpubr中的ANOVA是基础的，与Welch’s ANOVA不同
# 加载包
library(ggplot2)
library(ggpubr)
vars <- unique(dat$Species) %>% as.vector()            
varscomb <- combn(vars,2,simplify=F)
ggboxplot(dat, x = "Species", y = "Sepal.Length", fill = "Species",
         palette = c("#00AFBB", "#E7B800", "#FC4E07", "#6EA8FD"), alpha=0.5)+
  # 添加显著性标记，方法t.test、wilcox.test、anova等，comparisons时不能用anova
  stat_compare_means(
    comparisons = varscomb,
    method = "wilcox.test",                                     
    symnum.args = list(cutpoints = c(0, 0.001, 0.01, 0.05, Inf),
                       symbols = c("***", "**", "*", "NS")),
    label = "p.signif")+
  # 添加全局P值标记，没有comparisons是全局的
  stat_compare_means(label.y = 10, method = "anova")+               
  theme(legend.position = "none")

案例3. 非参数检验

当数据不满足正态性、方差齐次性，并且想做单因素组间比较时，可使用非参数检验方法

这里主要指Kruskal-Wallis检验，可检验多个（≥2）独立样本：是否来自相同的总体，或者说这些样本在中位数或总体位置上是否相同。由于该方法是非参数检验，不要求正态性和方差齐次性，因此在这些情况下，可以作为单因素ANOVA的替代。

R语言中，使用stats::kruskal.test函数执行Kruskal-Wallis检验，其H0为各组分布（中位数）全相同。

# 设定显著性水平α=0.05
a <- 0.05

# 示例数据
dat <- iris[c("Sepal.Length","Species")]

# 1. Kruskal-Wallis检验
krus_rest <- stats::kruskal.test(Sepal.Length ~ Species, data = dat)
# 结果如下，p-value < a，表明拒绝H0，认为组别之间中位数（或分布）不相同
#   Kruskal-Wallis rank sum test
# data:  Sepal.Length by Species
# Kruskal-Wallis chi-squared = 96.937, df = 2, p-value < 2.2e-16

# 2. 事后检验（多重比较）
# 可以使用非参数方法Dunn检验、Conover-Iman检验或DSCF(Dwass–Steel–Critchlow–Fligner)检验
comp_rest1 <- FSA::dunnTest(Sepal.Length ~ Species, data = iris, method = "bonferroni")
comp_rest2 <- PMCMRplus::kwAllPairsConoverTest(iris$Sepal.Length, iris$Species, p.adjust.method = "holm")
comp_rest3 <- PMCMRplus::dscfAllPairsTest(Sepal.Length ~ Species, data = iris)
# 结果如下，p<a，表明三组之间都有显著差异
# > comp_rest1
# Dunn (1964) Kruskal-Wallis multiple comparison
#   p-values adjusted with the Bonferroni method.
#               Comparison         Z      P.unadj        P.adj
# 1    setosa - versicolor -6.106326 1.019504e-09 3.058513e-09
# 2     setosa - virginica -9.741785 2.000099e-22 6.000296e-22
# 3 versicolor - virginica -3.635459 2.774866e-04 8.324597e-04
# > comp_rest2
#   Pairwise comparisons using Conover's all-pairs test
# data: iris$Sepal.Length and iris$Species
#            setosa  versicolor
# versicolor < 2e-16 -         
# virginica  < 2e-16 8.6e-09   
# P value adjustment method: holm
# > comp_rest3
#   Pairwise comparisons using Dwass-Steele-Critchlow-Fligner all-pairs test
# data: Sepal.Length by Species
#            setosa  versicolor
# versicolor 2.6e-13 -         
# virginica  3.1e-14 1.7e-06   
# P value adjustment method: single-step

# 3. 结果可视化，注意这里，整体显著性用kruskal.test标记，组间还是用的wilcox.test
# 加载包
library(ggplot2)
library(ggpubr)
vars <- unique(dat$Species) %>% as.vector()           
varscomb <- combn(vars,2,simplify=F)
ggboxplot(dat, x = "Species", y = "Sepal.Length", fill = "Species",
         palette = c("#00AFBB", "#E7B800", "#FC4E07", "#6EA8FD"), alpha=0.5)+
  # 添加显著性标记，方法t.test、wilcox.test、anova等，comparisons时不能用anova
  stat_compare_means(
    comparisons = varscomb,
    method = "wilcox.test",                                     
    symnum.args = list(cutpoints = c(0, 0.001, 0.01, 0.05, Inf),
                       symbols = c("***", "**", "*", "NS")),
    label = "p.signif")+
  # 添加全局P值标记，没有comparisons是全局的
  stat_compare_means(label.y = 10, method = "kruskal.test")+               
  theme(legend.position = "none")

案例4. 非参数多因素扩展

当数据不满足正态性、方差齐次性，并且想做多因素组间比较时，可采用非参数多因素扩展

这里主要指Scheirer Ray Hare检验，属于非参数方法，不依赖数据满足正态性或方差齐性假设，可用于多组比较，包括多因素设计和交互效应分析，是Kruskal-Wallis检验的多因素扩展

R语言中，使用rcompanion::scheirerRayHare函数执行Scheirer Ray Hare检验，其参数和含义如下：

1. formula：公式

2. data：数据

3. y：向量，因变量（响应变量），若未用formula指定，需直接传入

4. x1, x2：分类向量（因子或字符），若未用formula指定，需直接传入

5. type：整数1或2，用于指定平方和（SS）的计算方法：

   • type=1，序贯型平方和（Sequential SS），按模型中的变量顺序计算（先主效应后交互效应）

   • type=2，部分型平方和（Partial SS，默认），控制其他变量的影响后计算（更常用）

   • 若无交互效应或变量间独立，1或2结果相似，若存在交互效应，推荐2

6. tie.correct：逻辑值，是否对数据中的结（ties，即相同值的秩）进行校正默认TRUE

7. ss：逻辑值，是否在结果中返回平方和（SS）和均方（MS），默认TRUE

8. verbose：逻辑值，是否打印详细结果（如检验统计量、p 值），默认TRUE

# 设定显著性水平α=0.05
a <- 0.05

# 数据来自datasets
dat <- ToothGrowth 
dat$dose <- factor(dat$dose)
dat$supp <- factor(dat$supp)

# 1. 执行Scheirer Ray Hare检验
model <- rcompanion::scheirerRayHare(len ~ supp * dose, dat, type = 2, tie.correct = T, verbose = T)

# 结果如下，结果表明，主效应dose的p<a，supp和交互效应p>a，说明dose组间差异显著，而supp与交互效应不显著
# > model
#           Df  Sum Sq      H p.value
# supp       1  1050.0  3.445 0.06343
# dose       2 12394.4 40.669 0.00000
# supp:dose  2   515.5  1.692 0.42923
# Residuals 54  4021.1    
# 注意上述结果，与相同数据、公式下的多因素ANOVA结果不同（案例1.2），主要体现在supp和交互效应是否显著上：
# 多因素ANOVA
# model <- anova_test(dat, len ~ supp * dose, detailed = T) %>%
#   get_anova_table() %>%
#   setDT() %>% 
#   .[, R2:=SSn/(SSn + SSd)]         #计算各效应的R²
# 结果如下，p<a，表明拒绝原假设，认为各组均值不相等
# > model
#       Effect      SSn     SSd   DFn   DFd      F        p  p<.05   ges        R2
#       <char>    <num>   <num> <num> <num>  <num>    <num> <char> <num>     <num>
# 1:      supp  205.350 712.106     1    54 15.572 2.31e-04      * 0.224 0.2238254
# 2:      dose 2426.434 712.106     2    54 92.000 4.05e-18      * 0.773 0.7731092
# 3: supp:dose  108.319 712.106     2    54  4.107 2.20e-02      * 0.132 0.1320279

# 2. 事后多重比较，判断哪些组别之间均值不等
# 可以使用非参数方法Dunn检验或Conover-Iman检验对某个主效应进行检验，
# 或使用Aligned Rank Transform检验交互效应
# 对supp主效应的Dunn检验
# res1 <- FSA::dunnTest(len ~ supp, data = dat, method = "holm")             #这里分组小于等于2，报错
res1 <- PMCMRplus::kwAllPairsConoverTest(dat$len, dat$supp, method = "BH")
# 对dose主效应的Conover-Iman检验
res2 <- PMCMRplus::kwAllPairsConoverTest(dat$len, dat$dose, method = "BH")
# 对supp:dose交互效应的对齐秩变换检验
res3 <- ARTool::art(len ~ supp * dose, data = dat) %>% ARTool::art.con("supp:dose") 
# 结果如下：
# > res1
#    OJ   
# VC 0.063
# > res2
#   0.5     1      
# 1 1.6e-07 -      
# 2 1.1e-11 2.0e-05
# > res3
#  contrast        estimate   SE df t.ratio p.value
#  OJ,0.5 - OJ,1     -22.25 3.86 54  -5.766  <.0001
#  OJ,0.5 - OJ,2     -30.65 3.86 54  -7.942  <.0001
#  OJ,0.5 - VC,0.5     9.60 3.86 54   2.488  0.1460
#  OJ,0.5 - VC,1      -7.40 3.86 54  -1.918  0.4030
#  OJ,0.5 - VC,2     -30.00 3.86 54  -7.774  <.0001
#  OJ,1 - OJ,2        -8.40 3.86 54  -2.177  0.2655
#  OJ,1 - VC,0.5      31.85 3.86 54   8.253  <.0001
#  OJ,1 - VC,1        14.85 3.86 54   3.848  0.0041
#  OJ,1 - VC,2        -7.75 3.86 54  -2.008  0.3513
#  OJ,2 - VC,0.5      40.25 3.86 54  10.430  <.0001
#  OJ,2 - VC,1        23.25 3.86 54   6.025  <.0001
#  OJ,2 - VC,2         0.65 3.86 54   0.168  1.0000
#  VC,0.5 - VC,1     -17.00 3.86 54  -4.405  0.0007
#  VC,0.5 - VC,2     -39.60 3.86 54 -10.261  <.0001
#  VC,1 - VC,2       -22.60 3.86 54  -5.856  <.0001
# P value adjustment: tukey method for comparing a family of 6 estimates 

# 3. 结果可视化，与1.2相同，这里省略

案例5. 稳健ANOVA

当数据为小样本、不满足正态性、方差齐次性，并且想做单因素或多因素组间比较时，可使用稳健ANOVA方法

这里主要指WRS2中的均值修剪方差分析，它是一种广义的Welch方法。它通过修剪均值进行单因素方差分析，可减少极端值对分析结果的影响，适用于数据含有离群点或方差不均的情况。

在R语言WRS2包中，应用t1way、t2way和t3way函数用于均值截断方差分析，它们分别用于单因素、双因素和三因素情况。其中参数如下：

1. formula: 公式

2. data: 数据框

3. tr：数值型，修剪比例，默认0.2，表示去除数据的上下各20%的部分，用来计算修剪均值

4. alpha：数值型，显著性水平，默认值是0.05，用于判断结果是否显著

5. nboot：数值型，启用自助法（bootstrap）时的模拟次数，默认是100次，用于计算置信区间或p值。

单因素案例

# 设定显著性水平α=0.05
a <- 0.05

# 数据来自WRS2
set.seed(123)
dat <- WRS2::viagra

# 1. 均值截断方差分析（单因素）
tanova_result <- WRS2::t1way(libido ~ dose, data = dat)
# 结果如下，p>a，说明不能拒绝H0，认为各种用药剂量对没有影响
# Show in New Window
# Call:
# t1way(formula = libido ~ dose, data = viagra)
# Test statistic: F = 3 
# Degrees of freedom 1: 2 
# Degrees of freedom 2: 4 
# p-value: 0.16 
# Explanatory measure of effect size: 0.79 
# Bootstrap CI: [0.57; 1.57]
# 注意这里效应量为0.79，即使结果不显著，效应也比较明显

# 2. 事后多重比较，判断哪些组别之间均值不等
after_result <- WRS2::lincon(libido ~ dose, data = dat)
# 结果如下，
# > after_result
# Call:
# lincon(formula = libido ~ dose, data = dat)
#                  psihat ci.lower ci.upper p.value
# placebo vs. low      -1 -5.31858  3.31858 0.43533
# placebo vs. high     -3 -7.31858  1.31858 0.18051
# low vs. high         -2 -6.31858  2.31858 0.31660
# 结果中，psihat是估计效应值，指每对比较的效应大小，表示两个组之间的均值差异
# 这里-1表示placebo比low低1个单位，类似的placebo比high低3个单位low比high低2个单位

# 3. 结果可视化，与1.1相同，这里省略

多因素案例

# 设定显著性水平α=0.05
a <- 0.05

# 数据来自WRS2
set.seed(123)
dat <- WRS2::movie

# 1. 均值截断方差分析（多因素）
tanova_result <- WRS2::t3way(aggressive ~ degree*gender*type, data = dat)
# 结果如下，只有gender:type的p<a，其余p>a，说明大部分组不能拒绝H0，大部分组间均值差异不显著
# > tanova_result
# Call:
# t3way(formula = aggressive ~ degree * gender * type, data = dat)
#                        value p.value
# degree             3.3089497   0.085
# gender             1.3845004   0.260
# type               1.7102122   0.207
# degree:gender      3.6013168   0.073
# degree:type        0.3918961   0.540
# gender:type        8.7595078   0.008
# degree:gender:type 3.2614250   0.087

# 2. 事后多重比较，判断哪些组别之间均值不等
# 正确方法，同案例4
# 可以使用非参数方法Dunn检验或Conover-Iman检验对某个主效应进行检验，
# 或使用Aligned Rank Transform检验交互效应
# 对degree主效应的Dunn检验
# res1 <- FSA::dunnTest(aggressive ~ degree, data = dat, method = "holm")    #这里分组小于等于2，报错
res1 <- PMCMRplus::kwAllPairsConoverTest(dat$aggressive, dat$degree, method = "BH")
# 对gender主效应的Conover-Iman检验
res2 <- PMCMRplus::kwAllPairsConoverTest(dat$aggressive, dat$gender, method = "BH")
# 对degree:gender交互效应的对齐秩变换检验
res3 <- ARTool::art(aggressive ~ degree * gender, data = dat) %>% ARTool::art.con("degree:gender") 
# 结果如下：
# > res1
#   Pairwise comparisons using Conover's all-pairs test
# data: dat$aggressive and dat$degree
#        no degree
# degree 0.17     
# P value adjustment method: single-step
# > res2
#   Pairwise comparisons using Conover's all-pairs test
# data: dat$aggressive and dat$gender
#        male
# female 0.71
# P value adjustment method: single-step
# > res3
#  contrast                          estimate   SE df t.ratio p.value
#  degree,female - degree,male          12.22 6.56 64   1.863  0.2541
#  degree,female - no degree,female     17.59 6.75 64   2.606  0.0540
#  degree,female - no degree,male        8.97 6.56 64   1.368  0.5240
#  degree,male - no degree,female        5.37 6.56 64   0.819  0.8454
#  degree,male - no degree,male         -3.25 6.36 64  -0.511  0.9563
#  no degree,female - no degree,male    -8.62 6.56 64  -1.314  0.5574
# P value adjustment: tukey method for comparing a family of 4 estimates 

# 错误方法：
# after_result <- WRS2::lincon(aggressive ~ degree*gender*type, data = dat)
# 注意，多因素情况，不要用WRS2::lincon，原因是
# 即便可以用A:B或A*B公式传入，但其运行的还是A的单因素分析，示例数据和结果如下
# 此外，单因素type如果只有2组，则结果不会有如下这么详细
# > WRS2::lincon(value ~ type1*type2, data = dat)
# Call:
# WRS2::lincon(formula = value ~ type1 * type2, data = dat)
#           psihat ci.lower ci.upper p.value
# A vs. B -0.19613 -1.02741  0.63515 0.71106
# A vs. C  0.14024 -0.89231  1.17279 0.71106
# B vs. C  0.33637 -0.75783  1.43058 0.71106
# > WRS2::lincon(value ~ type1:type2, data = dat)
# Call:
# WRS2::lincon(formula = value ~ type1:type2, data = dat)
#           psihat ci.lower ci.upper p.value
# A vs. B -0.19613 -1.02741  0.63515 0.71106
# A vs. C  0.14024 -0.89231  1.17279 0.71106
# B vs. C  0.33637 -0.75783  1.43058 0.71106
# > WRS2::lincon(value ~ type1, data = dat)
# Call:
# WRS2::lincon(formula = value ~ type1, data = dat)
#           psihat ci.lower ci.upper p.value
# A vs. B -0.19613 -1.02741  0.63515 0.71106
# A vs. C  0.14024 -0.89231  1.17279 0.71106
# B vs. C  0.33637 -0.75783  1.43058 0.71106

# 3. 结果可视化，与1.2相同，这里省略

案例6. 显著性标记绘图总结

总结了显著性标记的几种方法

注意

1. bar，errorbar和text的position参数要相同

2. ggplot2中，调整errorbar图层至bar图层前面，可实现显示上半个errorbar，在ggpubr包中，可使用error.plot = "upper_errorbar"显示半个

3. 对于条形图，ggplot()+geom_bar(stat="identity")与ggplot()+geom_col()效果相同

4. 不想让aes继承，可以用这个参数inherit.aes = FALSE

5. ggplot2中只有mean_sdl，其中为2*sd，而ggpubr中有mean_sd

6. 标记图层继承多行数据行时，geom_text可能会警告，解决办法是用annotate或者是用inherit.aes=F

7. 用ggpubr+rstatix实现检验和标记，此外，也可以用ggsignif::geom_signif包

8. Kruskal-Wallis检验，可识别多组独立样本的中位数是否有差别（至少有一组的中位数与其他组不同），相当于非参数版的单因素方差分析

9. Wilcoxon秩和检验，又名Mann-Whitney-U检验，可识别两组独立样本的中位数或分布差异，相当于非参数版的独立样本t检验

1. 单因素标记

（1）分布图（boxplot或violin）+内部检验+显著性成对标记（星号或p值）

此种情况，使用使用全部数据，无需统计变换，检验和标记也用ggpubr内部方法，相对简单，但可调节选项稍微少一些，虽然ggpubr绘图与ggplot2可以无缝衔接

# 加载包
library(ggplot2)
library(ggpubr)

# 设定显著性水平α=0.05
a <- 0.05

# 生成四组数据
set.seed(123)
groupA <- rnorm(30, mean = 50, sd = 5)
groupB <- rnorm(30, mean = 55, sd = 5)
groupC <- rnorm(30, mean = 60, sd = 5)
groupD <- rnorm(30, mean = 64, sd = 5)

# 数据框整理
dat <- data.frame(
  value = c(groupA, groupB, groupC, groupD),
  group = factor(rep(c("A", "B", "C", "D"), each = 30))
)
 
# 绘图
vars <- unique(dat$group) %>% as.vector()      
varscomb <- combn(vars,2,simplify=F)
ggviolin(dat, x = "group", y = "value", fill = "group",
         palette = c("#00AFBB", "#E7B800", "#FC4E07", "#6EA8FD"), alpha=0.5,
         add = "boxplot", add.params = list(fill = "white"))+
  # 添加标记，有comparisons是成对的，方法包括t.test、wilcox.test（成对时不能用anova这种全局的）
  stat_compare_means(
    comparisons = varscomb,
    method = "wilcox.test",                                     
    symnum.args = list(cutpoints = c(0, 0.001, 0.01, 0.05, Inf),
                       symbols = c("***", "**", "*", "NS")),
    label = "p.signif")+
  # 添加标记，没有comparisons是全局的，标记整体显著性，方法包括anova等（整体是不能用t.test这种成对的，会重叠）
  stat_compare_means(label.y = 30, method = "anova")+               
  theme(legend.position = "none")

（2）均值图（barplot或point）+误差（se或sd）+外部检验+显著性成对标记（星号或p值）

此种情况，使用使用全部数据+统计变换（此处用ggpubr包免变换），区别与上面，这里用外部检验，稍微放宽了对ggpubr的依赖，可调参数稍多一些。

# 加载包
library(ggplot2)
library(ggpubr)
# library(rstatix) 
library(magrittr)

# 设定显著性水平α=0.05
a <- 0.05

# 生成四组数据
set.seed(123)
groupA <- rnorm(30, mean = 50, sd = 5)
groupB <- rnorm(30, mean = 55, sd = 5)
groupC <- rnorm(30, mean = 60, sd = 5)
groupD <- rnorm(30, mean = 64, sd = 5)

# 数据框整理
dat <- data.frame(
  value = c(groupA, groupB, groupC, groupD),
  group = factor(rep(c("A", "B", "C", "D"), each = 30))
)
 
# 统计检验（Wilcoxon秩和检验，Bonferroni校正）
stat.test <- dat %>%
  rstatix::wilcox_test(value~group, p.adjust.method = "bonferroni") %>%
  rstatix::add_xy_position(x = "group", data = dat, step.increase = 0.3)     #添加作图属性

# 修改显著性符号（等效于上面的symnum.args）
stat.test <- stat.test %>%
  dplyr::mutate(
    p.signif_manu = dplyr::case_when(
      p.adj < 0.001 ~ "***",
      p.adj < 0.01  ~ "**",
      p.adj < 0.05  ~ "*",
      TRUE          ~ "NS"  # 不显著
    )
  )

# 计算全局检验（ANOVA或Kruskal-Wallis）
global_p <- aov(Sepal.Length ~ Species, data = iris) %>% summary() %>% {.[[1]]["Pr(>F)"][1, 1]}
# global_p <- kruskal.test(Sepal.Length ~ Species, data = iris)$p.value

# 将p值格式化为文本（带方法名称）
formatted_p <- paste0("ANOVA, p", 
                      ifelse(grepl("<", format.pval(global_p)), " ", " = "), 
                      format.pval(global_p))

# 绘图
ggbarplot(dat, x = "group", y = "value", fill = "group", 
          palette = c("#00AFBB", "#E7B800", "#FC4E07", "#6EA8FD"), alpha=0.5,
          add = "mean_se", error.plot = "upper_errorbar")+                     #ggplot没有upper_errorbar方法
  # 添加全局P值标记
  annotate("text", x=0.5, y=90, label=formatted_p, hjust=0)+
  # 添加配对显著性标记
  stat_pvalue_manual(
    stat.test,
    y.position = stat.test$y.position-20,  #整体下移20
    label = "p.signif_manu",               #label所在变量
    tip.length = 0.03,                     #括号（“连接线”）的两边长度（0=无横线）
    bracket.size = 0.5,                    #括号的线粗
    label.size = 5,                        #显著性标签（如 "***"）字体大小
    inherit.aes = FALSE                    #不继承主图的 aes 映射（自己指定的变量）
  )+
  theme_classic()+
  theme(legend.position = "none")

（3）均值图（barplot或point）+误差（se或sd）+外部检验+显著性整体标记（字母）

此种情况，使用使用全部数据+统计变换（此处用ggplot2包变换），这里使用ggplot2包，可调参数最多

# 加载包
library(ggplot2)
library(magrittr)
library(multcomp)
library(emmeans) 

# 设定显著性水平α=0.05
a <- 0.05

# 生成四组数据
set.seed(123)
groupA <- rnorm(30, mean = 50, sd = 5)
groupB <- rnorm(30, mean = 55, sd = 5)
groupC <- rnorm(30, mean = 60, sd = 5)
groupD <- rnorm(30, mean = 64, sd = 5)

# 数据框整理
dat <- data.frame(
  value = c(groupA, groupB, groupC, groupD),
  group = factor(rep(c("A", "B", "C", "D"), each = 30))
)
 
# 字母显著性创建方法1，用multcomp::glht的Tukey方法，结果只是一个向量，如需绘图还得整理（不推荐）
# cld_result1 <- aov(value ~ group, data = dat) %>% 
#   multcomp::glht(linfct = mcp(group = "Tukey")) %>% 
#   multcomp::cld(decreasing = T)          #这里参数Letters = LETTERS不起作用

# 字母显著性创建方法2，用emmeans::emmeans的线性假设，结果可能与aov的Tukey不同，但是df型，可以直接用
cld_result2 <- aov(value ~ group, data = dat) %>% 
  emmeans::emmeans(~ group) %>% 
  # graphics::pairs(adjust = "tukey") %>%    #显式应用Tukey，不要加这句，否则就是成对的字母显著性标记
  multcomp::cld(Letters = letters, decreasing = T) 
cld_result2$.group <- stringr::str_trim(cld_result2$.group)  #标记两侧有空格
 
# 计算全局检验（ANOVA或Kruskal-Wallis）
global_p <- aov(Sepal.Length ~ Species, data = iris) %>% summary() %>% {.[[1]]["Pr(>F)"][1, 1]}
# global_p <- kruskal.test(Sepal.Length ~ Species, data = iris)$p.value

# 将p值格式化为文本（带方法名称）
formatted_p <- paste0("ANOVA, p", 
                      ifelse(grepl("<", format.pval(global_p)), " ", " = "), 
                      format.pval(global_p))

# 绘图
ggplot(dat, aes(group, value, fill=group))+
  stat_summary(fun.data = mean_sd, geom = "errorbar", width=0.3)+       #ggplot2只能放在col下面以隐藏下半errorbar
  stat_summary(fun = mean, geom = "col", show.legend = F)+
  annotate("text", x=0.5, y=75, label=formatted_p, hjust=0)+
  geom_text(data = cld_result2, aes(y=emmean+10,label=.group))+
  theme_classic()

2. 双因素标记

（1）分布图（boxplot或violin）/均值误差图（barplot或point+se或sd）+外部检验+显著性分组标记（星号或p值）

此种情况，自定义了一个标记函数，对分布或变换数据无要求，使用ggplot2，标记每个x轴对应组别内部的，关于分类变量之间的显著性

# 加载包
library(magrittr)
library(data.table)
# library(rstatix)
library(ggplot2)
library(multcomp)    
library(ggsci)
library(stringr)
library(magrittr)

# 设定显著性水平α=0.05
a <- 0.05

# 数据来自datasets
dat <- ToothGrowth 
dat$dose <- factor(dat$dose)
dat$supp <- factor(dat$supp)

# 定义一个函数，生成双因素方差结果图的显著性成对标记图层
# 成对标记为，x轴方向+class分类变量，在doge形式上，标记每个x类别对应的class组间对比
pairs_sig_mark <- function(data, yn, xn, cn, xlength=0.4, ylengths=0.02, yheighs=0.1, theighs=0.02, is_pvalue=F,
                           p_dig=3, linewidth=0.5, size=4, ...){
  # data：绘图数据（在...中传递）
  # yn,xn,cn：分别是y轴、x轴、分类的变量所在列号，每个变量一个
  # xlength：横向标记线长度，绝对值
  # ylengths：竖向标记线（横向两端）长度乘数（长度为数据y范围*yscale），0-1之间
  # yheighs：横向标记线的高度位置尺度（高度为数据y最大值+y范围*yvscale），0-1之间
  # theighs：文字标记的高度位置尺度（高度为yvscale+范围*tvscale），0-1之间
  # is_pvalue：显著性标记p值还是星号
  # p_dig：p值小数位数，低于此精度p值显示小于号
  
  # 获取、定制数据基本信息 
  setDT(data)
  ocolname <- names(data)             
  sigdata <- data[, .SD, .SDcols = c(yn, xn, cn)] %>% 
    set_colnames(c("yy","xx","cc"))
  xx_level <- unique(sigdata$xx)                      #x轴上分类个个数
  yrange <- range(sigdata$yy) %>% diff() %>% abs()    #y轴范围尺度，全局数据
  
  # 每组显著性，x轴变量的每组下，对class变量单独做单因素方差分析，提取结果
  p_value <- rep(0, length(xx_level))
  for(i in seq_along(xx_level)){
    p_value[i] <- sigdata[xx==xx_level[i], ] %>% rstatix::anova_test(yy~cc, white.adjust = T) %>% .[,"p"]
  }
  if(isTRUE(is_pvalue)) {                              #返回p值还是星号
    signal_text <- ifelse(p_value<10^(-p_dig), paste0("p<1e-",p_dig), sprintf(paste0("p=%.",p_dig,"f"),p_value))
  } else {
    signal_text <- ifelse(p_value <= 0.001, "***",
                          ifelse(p_value <= 0.01,  "**",
                                 ifelse(p_value <= 0.05,  "*", "NS")))
  }
   
  # 显著性标线数据，横向部分
  sigdata_x <- copy(sigdata) %>% 
    .[,.(yymax=max(yy)+yrange*yheighs), by=c("xx")] %>%                  #指定横线y坐标
    .[, xcen:=1:length(xx_level)] %>%                                    #指定横线中心x坐标
    .[, c("xleft", "xright"):=list(xcen-xlength/2, xcen+xlength/2)]      #由中心向两端延伸出端点x坐标
  
  # 显著性标线数据，纵向部分
  sigdata_y <- copy(sigdata_x) %>% 
    melt(id=c("xx", "yymax", "xcen"), value.name = "xendpoint") %>%      #指定竖线x坐标、上端y坐标（即yymax）
    .[, yymin:=yymax-yrange*ylengths]                                    #指定竖线下端y坐标
  
  # 显著性文字
  sigdata_text <- copy(sigdata_x) %>% 
    .[,ycen:=yymax+yrange*theighs] %>%                                   #指定显著性文字y坐标
    .[,p_l:=signal_text]                                                 #指定显著性文字内容
  
  # 图层输出
  geom_list <- list(
    geom_linerange(data=sigdata_x, aes(y=yymax, xmin=xleft, xmax=xright), show.legend=F, inherit.aes=F, linewidth=linewidth, ...),
    geom_linerange(data=sigdata_y, aes(x=xendpoint, ymin=yymin, ymax=yymax), show.legend=F, inherit.aes=F, linewidth=linewidth, ...),
    geom_text(data=sigdata_text, aes(xcen, ycen, label=p_l), show.legend=F, inherit.aes=F, hjust=0.5, vjust=0, size=size, ...)
  )
  return(geom_list)
}

# 绘图
ggplot(data=dat, aes(dose, len, fill=supp))+
  # geom_violin(trim=F)+
  stat_summary(fun.data = ggpubr::mean_sd, geom="errorbar", width=0.2, position=position_dodge(width=0.6))+           #ggplot2中只有mean_sdl，其中为2*sd
  stat_summary(fun=mean, geom="col", width=0.6, position=position_dodge(width=0.6))+
  pairs_sig_mark(data=dat, yn=1,xn=3,cn=2, is_pvalue=T,linewidth=0.5,size=3)+
  theme_classic()

（2）均值和误差（barplot或point+se或sd）+外部检验+显著性整体标记（字母）

此种情况，使用使用全部数据+统计变换（此处用ggplot2包变换），这里使用ggplot2包，可调参数最多

# 加载包
library(magrittr)
library(data.table)
# library(rstatix)       #这个包会和其他包常用函数有冲突，如select、filter
library(ggplot2)
library(multcomp)    
library(ggsci)
library(stringr)
library(magrittr)

# 设定显著性水平α=0.05
a <- 0.05

# 数据来自datasets
dat <- ToothGrowth 
dat$dose <- factor(dat$dose)
dat$supp <- factor(dat$supp)

# 字母显著性创建方法2，用emmeans::emmeans的线性假设，结果可能与aov的Tukey不同，但是df型，可以直接用
cld_result_temp <- aov(len ~ supp*dose, data = dat) %>%          #注意这里需要考虑交互，否则结果不严谨
  emmeans::emmeans(~ supp|dose) %>% 
  multcomp::cld(Letters = letters, decreasing = T) %>% 
  setDT() %>% 
  .[,.group:=stringr::str_trim(.group)]                         #标记两侧有空格

cld_result <- copy(dat) %>% 
  setDT() %>% 
  .[, .(se=sd(len)/sqrt(.N)), by=c("supp","dose")] %>% 
  .[cld_result_temp, on=c("supp","dose")]

# 计算全局检验（ANOVA）
global_p <- aov(Sepal.Length ~ Species, data = iris) %>% summary() %>% {.[[1]]["Pr(>F)"][1, 1]} 

# 将p值格式化为文本（带方法名称）
formatted_p <- paste0("ANOVA, p", 
                      ifelse(grepl("<", format.pval(global_p)), " ", " = "), 
                      format.pval(global_p))

# 绘图
pd <- position_dodge(width=0.6)
ggplot(dat, aes(dose, len, fill=supp))+
  stat_summary(fun.data = mean_se, geom = "errorbar", width=0.3, position = pd)+       #ggplot2只能放在col下面以隐藏下半errorbar
  stat_summary(fun = mean, geom = "col", width=0.5, position = pd, show.legend = F)+
  annotate("text", x=0.5, y=30, label=formatted_p, hjust=0)+
  geom_text(data = cld_result, aes(x=dose, y=emmean+se+1, label=.group), position = pd, hjust=0.5, vjust=0)+
  theme_classic()

3. 自定义标记

早前写的，自定义函数版本，事先将数据汇总，并计算好显著性，然后画在图中，为了向bar图添加显著性标记，自定义了一个函数。

# 加载所需包
library(magrittr)
library(data.table)
library(ggplot2)
library(multcomp)   #计算avo的pvalue
library(ggsci)
# library(ggsignif) #当画箱形图时，可用这个包

# 定义一个函数，生成均值方差和显著性标记线、标签（此处用的是方差分析）
data_material <- function(data, xdrift=0.15, ydrift=0.04, equleg=T){
  # xdrift和ydrift为xy方向偏移量，其中xdrift应为bar宽度width的一半，equleg设置标记的腿是否等长
  # 提取计算列，备份原列名，然后赋予新列名，便于后续计算  
  con_data <- data[,1:3]
  ori_name <- colnames(con_data)
  names(con_data) <- c("xaxis", "compare", "value")
  
  #先按因子类型排序con_data
  setorder(con_data, xaxis, compare)         
  
  # 生成均值、方差列
  data_msd <- con_data[, .(m=mean(value), sd=sd(value)), by=c("xaxis", "compare")]
  
  # 方差分析，生成p值显著性标签
  p.frame <- data.table(xaxis=unique(con_data$xaxis), p.value=NA_real_, p.label=NA_character_)
  for(i in p.frame$xaxis){
    tres <- aov(value~compare, con_data[xaxis==i,])
    tp <- glht(tres, linfct = mcp(compare = "Tukey")) %>% summary()
    p.frame[xaxis == i, p.value:=tp$test$pvalues] 
    p.frame[xaxis == i, p.label:=ifelse(p.value<0.01, "***", ifelse(p.value<0.05, "**", ifelse(p.value<0.1, "*", "NS")))] 
  }
  
  # 生成标记线图层
  data_msd[, msd:=m+sd]                      #添加一列，为mean+sd  
  compare_lv1 <- levels(data_msd$compare)[1]
  compare_lv2 <- levels(data_msd$compare)[2]
  xcen <- seq_along(unique(data_msd$xaxis))  #x的中心位置，为fa1的变量个数
  xoff <- xdrift                             #x起止位置偏移量，为bar宽度的一半
  m_plus_sd <- data_msd[, max(msd), by=c("xaxis")]$V1     
  unitdrift <- ydrift*max(m_plus_sd)           #全图中最大的mean+sd乘以偏移比，为单位偏移量
  ycen <- m_plus_sd+unitdrift                  #y中心位置，为m+sd加单位偏移量
  yleft <- data_msd[compare==compare_lv1, msd]+unitdrift*0.5 #左脚最低位置，为m+sd加0.5个单位偏移量
  yrigh <- data_msd[compare==compare_lv2, msd]+unitdrift*0.5 #右脚最低位置，为m+sd加0.5个单位偏移量
  if(isTRUE(equleg)) yleft <- yrigh <- pmax(yleft,yrigh)     #如果相同，则取最大的那个
  
  # 生成标记线和标签数据，p值默认两位小数
  sig_seg <- data.table(x0=c(xcen-xoff, xcen-xoff, xcen+xoff), x1=c(xcen-xoff, xcen+xoff, xcen+xoff),
                        y0=c(yleft, ycen, ycen), y1=c(ycen, ycen, yrigh))
  sig_lab <- data.table(x=xcen, y=ycen+unitdrift*0.5, p.value=sprintf("%.2f", p.frame$p.value), p.label=p.frame$p.label)
  
  # 恢复数据列名
  names(data_msd) <- c(ori_name[1:2], "m", "sd","msd")
  names(p.frame) <- c(ori_name[1], "p.value", "p.label")
  
  # 输出结果
  res <- list(data.msd=data_msd, sig.segment=sig_seg, sig.label=sig_lab)
  return(res)
}

# 数据，第一列是x轴，第二列是对比类（只有两类），第三列是数据
# 其中，第一、二列必须是因子型，顺序是已定义的
data <- data.table(method = c('B', 'B', 'B', 'A', 'A', 'A', 'C', 'C', 'C', 'D', 'D', 'D', 
                              'E', 'E', 'E', 'F', 'F', 'F', 'G', 'G', 'G', 'H', 'H', 'H', 
                              'I', 'I', 'I', 'J', 'J', 'J', 'B', 'B', 'B', 'A', 'A', 'A', 
                              'C', 'C', 'C', 'D', 'D', 'D', 'E', 'E', 'E', 'F', 'F', 'F',
                              'G', 'G', 'G', 'H', 'H', 'H', 'I', 'I', 'I', 'J', 'J', 'J'), 
                   model = c('B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B',
                             'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B',
                             'B', 'B', 'B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A',
                             'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A',
                             'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A'),
                   MAE = c(9.04, 1.20, 4.80, 16.24, 7.80, 7.63, 15.38, 7.33, 5.10, 10.14,
                           6.77, 11.20, 2.49, 5.20, 2.33, 15.78, 7.00, 14.63, 8.18, 2.73,
                           4.13, 16.18, 7.87, 13.97, 14.51, 4.93, 9.87, 14.44, 5.93, 9.30,
                           7.60, 5.24, 6.94, 19.08, 7.59, 1.69, 16.95, 7.39, 5.09, 14.59,
                           6.03, 10.08, 8.66, 10.31, 2.71, 23.00, 10.28, 18.10, 15.75, 7.13,
                           6.50, 20.84, 7.90, 15.17, 21.33, 8.25, 17.28, 20.63, 8.07, 15.99))
method_lev <- c("A","B","C","D","E","F","G","H","I","J")
model_lev <- c("A", "B")
data[, method:=factor(method, levels = method_lev)]
data[, model:=factor(model, levels = model_lev)]

# 生成均值方差、和显著性数据
res <- data_material(data)
data_msd <- res$data.msd
sig_seg <- res$sig.segment
sig_lab <- res$sig.label

# 画图
ggplot(data_msd, aes(method, m, fill = model))+
  geom_errorbar(aes(ymin= m-sd, ymax=m+sd), position = position_dodge(0.6), width = 0.3)+
  geom_bar(stat = "identity", position = position_dodge(), width = 0.6)+
  scale_fill_npg(labels = c("Model_A", "Model_B"))+
  geom_segment(data = sig_seg, aes(x=x0, xend=x1, y=y0, yend=y1, fill=NULL), color="red")+
  geom_text(data=sig_lab, aes(x=x, y=y, label=p.value, fill=NULL), color="red", vjust=0, family = "serif")+    #可选，显示p值
  # geom_text(data=sig_lab, aes(x=x, y=y, label=p.label, fill=NULL), color="red", vjust=0, family = "serif")+  #可选，显示显著性符号
  scale_y_continuous(limits = c(0, 26), expand = c(0,0))+
  guides(fill = guide_legend(title = NULL))+
  labs(x = "Method", y = "MAE")+
  theme_bw()+
  theme(aspect.ratio = 1/1.8,
        axis.text = element_text(size = 10, family = "serif"),
        axis.title = element_text(size = 13, family = "serif"),
        text = element_text(family = "serif"),
        panel.grid = element_blank(),
        legend.text = element_text(size = 10),
        legend.position = c(0.12, 0.89),
        plot.title = element_text(hjust = 0.5))

dfsd

5.8 效应量计算

检验类型	是否推荐计算效应量	常用效应量
t 检验	✅ 是	Cohen’s d
Wilcoxon/Mann–Whitney	✅ 是	r（基于Z的转换）
卡方检验	✅ 是	Cramér’s V

2. 常见效应量类型

1. 均值差异类（适用于t检验） Cohen’s d 独立样本t检验 标准化均值差（单位：标准差） Hedges’ g 对Cohen’s d的小样本校正 小样本独立t检验 更保守的估计 Glass’ Δ 使用对照组标准差 实验组vs对照组 避免实验组方差影响

配对d 配对样本t检验 标准化配对差值

2. 关联强度类（适用于相关/回归分析） 类型 公式 适用场景 Pearson’s r 相关系数 r r 连续变量相关性分析 R² 回归模型解释方差比例 线性回归

目的

适用

• 检验类型：

• 数据类型：

• 数据体量：

计算公式

结果判定

统计量z服从标准正态分布。设定显著性水平为α后，检验结果判定步骤如下：

1. 若H0为=，则双侧检验，\(p=2P(X \geq |z|)\)，拒绝域为\(|z| \geq z_{1-\alpha/2}\) （由于对称，只取右侧）

2. 若H0为≤，则右侧检验，\(p=P(X \geq z)\)，拒绝域为\(z \geq z_{1-\alpha}\)

3. 若H0为≥，则左侧检验，\(p=P(X \leq z)\)，拒绝域为\(z \leq z_{\alpha}\)

案例

R语言中，使用stats::shapiro.test执行SW检验：

判定方法处给出检验的H0

总结处汇总一下，各种检验的R语言函数

📊 错误类型对照表

真实情况  判决	不拒绝 \(H_0\)	拒绝 \(H_0\)
\(H_0\) 为真	✅ 正确	❌ 第一类错误（Type I）
\(H_0\) 为假	❌ 第二类错误（Type II）	✅ 正确

总结：

类型	名称	含义	控制方式
Type I	第一类错误	错误拒绝了正确的 \(H_0\)	控制显著性水平 \(\alpha\)
Type II	第二类错误	错误接受了错误的 \(H_0\)	提高样本量 / 效能