计算框架
因子A有\(A_1,A_2,...,A_k\) 共k个水平
总变异total sum of squares
\[
SS_T=\sum_{j=1}^{k}\sum_{i=1}^{n_j}X_{ij}^2-\frac{(\sum_{j=1}^{k}\sum_{i=1}^{n_j}X_{ij})^2}{n}=SS_{组间}+SS_{组内}
\]
组间变异 between groups sum of squares \[
SS_{组间}=\sum_{j=1}^{k}\frac{(\sum_{i=1}^{n_j}X_{ij})^2}{n_j}-\frac{(\sum_{j=1}^{k}\sum_{i=1}^{n_j}X_{ij})^2}{n}
\]
自由度\(\nu=n-1,\nu_{组间}=k-1,\nu_{组内}=\sum_{j=1}^{k}(n_j-1)=n-k\)
Between groups mean square \(MS_{组间}=\frac{SS_{组间}}{k-1}\)
Within groups mean square \(MS_{组内}=\frac{SS_{组内}}{n-k}\)
\[H_0:\mu_1=\mu_2=...=\mu_k\]
\[
\frac{SS_T}{\sigma^2}\sim \chi^2(\nu),\nu=n-1
\] \[
\frac{SS_{组内}}{\sigma^2}\sim \chi^2(\nu),\nu=n-k
\] 因此，
\[
\frac{SS_{组间}}{\sigma^2}=\frac{SS_T}{\sigma^2}-\frac{SS_{组内}}{\sigma^2}\ \ \ \sim \chi^2(\nu),\nu=k-1
\]
检验统计量
\[
F=\frac{\frac{SS_{组间}}{(k-1)\sigma^2}}{\frac{SS_{组内}}{(n-k)\sigma^2}}=\frac{\frac{SS_{组间}}{k-1}}{\frac{SS_{组内}}{n-k}}=\frac{MS_{组间}}{MS_{组内}}\ \ \ \sim \chi^2(\nu),\nu=k-1
\]
数据来源
2017.临床研究中的统计分析和图形表达实例详解 第2版
Show the code # 宽格式数据框（和原表格一致） library ( tidyverse ) df_wide <- data.frame ( normal = c ( 332.96 , 297.64 , 312.57 , 295.47 , 284.25 , 307.97 , 292.12 , 244.61 , 261.46 , 286.46 , 322.49 , 282.42 ) ,
middle = c ( 253.21 , 235.87 , 269.30 , 258.90 , 254.39 , 200.87 , 227.79 , 237.05 , 216.85 , 238.03 , 238.19 , 243.49 ) ,
high = c ( 232.55 , 217.71 , 216.15 , 220.72 , 219.46 , 247.47 , 280.75 , 196.01 , 208.24 , 198.41 , 240.35 , 219.56 )
) df_long <- df_wide |> pivot_longer ( cols = everything ( ) , names_to = "level" ,
values_to = "value" )
R实现
stats::aov()
Show the code df_long $ level <- factor ( df_long $ level ) df_aov <- aov ( value ~ level ,data = df_long ) anova ( df_aov ) #> Analysis of Variance Table #> #> Response: value #> Df Sum Sq Mean Sq F value Pr(>F) #> level 2 31292 15646 31.355 2.343e-08 *** #> Residuals 33 16467 499 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
stats::lm()
Show the code lm_aov <- lm ( formula = value ~ level , data = df_long ) lm_aov #> #> Call: #> lm(formula = value ~ level, data = df_long) #> #> Coefficients: #> (Intercept) levelmiddle levelnormal #> 224.78 14.71 68.59 anova ( lm_aov ) #> Analysis of Variance Table #> #> Response: value #> Df Sum Sq Mean Sq F value Pr(>F) #> level 2 31292 15646 31.355 2.343e-08 *** #> Residuals 33 16467 499 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 model.matrix ( ~ level , data = df_long ) #> (Intercept) levelmiddle levelnormal #> 1 1 0 1 #> 2 1 1 0 #> 3 1 0 0 #> 4 1 0 1 #> 5 1 1 0 #> 6 1 0 0 #> 7 1 0 1 #> 8 1 1 0 #> 9 1 0 0 #> 10 1 0 1 #> 11 1 1 0 #> 12 1 0 0 #> 13 1 0 1 #> 14 1 1 0 #> 15 1 0 0 #> 16 1 0 1 #> 17 1 1 0 #> 18 1 0 0 #> 19 1 0 1 #> 20 1 1 0 #> 21 1 0 0 #> 22 1 0 1 #> 23 1 1 0 #> 24 1 0 0 #> 25 1 0 1 #> 26 1 1 0 #> 27 1 0 0 #> 28 1 0 1 #> 29 1 1 0 #> 30 1 0 0 #> 31 1 0 1 #> 32 1 1 0 #> 33 1 0 0 #> 34 1 0 1 #> 35 1 1 0 #> 36 1 0 0 #> attr(,"assign") #> [1] 0 1 1 #> attr(,"contrasts") #> attr(,"contrasts")$level #> [1] "contr.treatment" lm ( formula = value ~ 0 + level , data = df_long ) #> #> Call: #> lm(formula = value ~ 0 + level, data = df_long) #> #> Coefficients: #> levelhigh levelmiddle levelnormal #> 224.8 239.5 293.4 model.matrix ( ~ 0 + level , data = df_long ) #> levelhigh levelmiddle levelnormal #> 1 0 0 1 #> 2 0 1 0 #> 3 1 0 0 #> 4 0 0 1 #> 5 0 1 0 #> 6 1 0 0 #> 7 0 0 1 #> 8 0 1 0 #> 9 1 0 0 #> 10 0 0 1 #> 11 0 1 0 #> 12 1 0 0 #> 13 0 0 1 #> 14 0 1 0 #> 15 1 0 0 #> 16 0 0 1 #> 17 0 1 0 #> 18 1 0 0 #> 19 0 0 1 #> 20 0 1 0 #> 21 1 0 0 #> 22 0 0 1 #> 23 0 1 0 #> 24 1 0 0 #> 25 0 0 1 #> 26 0 1 0 #> 27 1 0 0 #> 28 0 0 1 #> 29 0 1 0 #> 30 1 0 0 #> 31 0 0 1 #> 32 0 1 0 #> 33 1 0 0 #> 34 0 0 1 #> 35 0 1 0 #> 36 1 0 0 #> attr(,"assign") #> [1] 1 1 1 #> attr(,"contrasts") #> attr(,"contrasts")$level #> [1] "contr.treatment"
ez::ezANOVA()
Show the code df_long <- df_long |> rowid_to_column ( var = "id" ) |> relocate ( id ,.before = 1 ) |> mutate ( id= factor ( id ) )
ez :: ezANOVA ( data = df_long , dv = value ,
wid = id ,
between = level ,
type = 3 ,
detailed = T )
#> $ANOVA #> Effect DFn DFd SSn SSd F p p<.05 #> 1 (Intercept) 1 33 2296103.8 16466.87 4601.44760 5.097646e-37 * #> 2 level 2 33 31291.8 16466.87 31.35476 2.342720e-08 * #> ges #> 1 0.9928794 #> 2 0.6552067 #> #> $`Levene's Test for Homogeneity of Variance` #> DFn DFd SSn SSd F p p<.05 #> 1 2 33 135.1174 7847.544 0.2840937 0.7545184
前提假设
https://www.statmethods.net/stats/rdiagnostics.html
方差齐性
http://www.cookbook-r.com/Statistical_analysis/Homogeneity_of_variance/
Bartlett’s test
数据满足正态性
比较每一组方差的加权算术均值和几何均值。
\[
H_0:\sigma_1^2=\sigma_2^2=...=\sigma_k^2
\] 当样本量\(n_j\) ≥5时，检验统计量(各组样本量相等)
\[
B=\frac{(n-1)[kln\bar S^2-\sum_{j=1}^{k}lnS_j^2]}{1+\frac{k+1}{3k(n-1)}} \sim\ \chi^2(\nu)\ ,\nu=k-1
\] 其中n是每一组的样本量，\(S_j^2\) 是某一组的样本方差，\(\bar S^2\) 是所有k个组样本方差的平均值。
当各组样本量不等时，
\[
B=\frac{\sum_{j=1}^{k}(n_j-1)ln\frac{\bar S^2}{S_j^2}}{1+\frac{1}{3(k-1)}(\sum_{j=1}^{k}\frac {1}{n_j-1}-\frac{1}{\sum_{j=1}^{k}(n_j-1)})} \sim\ \chi^2(\nu)\ ,\nu=k-1
\] 其中\(h_j\) 是某一组的样本量，\(\bar S^2=(\sum_{j=1}^{k}(n_j-1)S_j^2)/(\sum_{j=1}^{k}(n_j-1))\) 是所有k个组样本方差的加权平均值。
Show the code bartlett.test ( value ~ level ,data = df_long ) #> #> Bartlett test of homogeneity of variances #> #> data: value by level #> Bartlett's K-squared = 0.83651, df = 2, p-value = 0.6582
Levene’s test
数据不满足正态性
k个随机样本是独立的
随机变量X是连续的
Show the code car :: leveneTest ( value ~ level ,data = df_long ,center = mean ) # 同SPSS #> Levene's Test for Homogeneity of Variance (center = mean) #> Df F value Pr(>F) #> group 2 0.3195 0.7287 #> 33
Levene 变换：
\[
Z_{ij}=|X_{ij}-\bar X_{.\ j}|(i=1,2...,n_j;\ j=1,2,...,k)
\]
检验统计量（基于变换后的F检验）
\[
W=\frac{MS_{组间}}{MS_{组内}}=\frac{\sum_jn_j(\bar Z_{.j}-\bar Z_{..})^2/(k-1)}{\sum_j\sum_i(Z_{ij}-\bar Z_{.j})^2/(n-k)} \sim F(\nu_1,\nu_2) \ \ \ \ \nu_1=k-1,\nu_2=n-k
\]
事后比较（Post hoc）
成对比较的数量 \(N=\frac{k!}{2!(k-2)!},k≥3\) ，导致犯第Ⅰ类错误的概率迅速增加，\([1-(1-\alpha)^N]\) 。
Tukey’s test
Tukey’s test 也被称为Tukey’s honestly significant difference (Tukey’s HSD) test。
k个均值从大到小排列；
均值最大的组依次与均值最小，第二小，……，第二大比较；
均值第二大的组以同样的方式比较；
以此类推
在各组样本量相等的情况下，如果在两个均值之间未发现显著差异，则推断这两个均值所包含的任何均值之间不存在显著差异，并且不再检验所包含均值之间的差异。
studentized range statistic \(q=\frac{\bar X_{max}-\bar X_{min}}{S_{\bar X_{max}-\bar X_{min}}}\) ，其中\(S_{\bar X_{max}-\bar X_{min}}=\sqrt{\frac{MS_{组内}}{n}}\) ，\(n\) 是每一个治疗组的样本量。
如果各组样本量不等,则\(S_{\bar X_{max}-\bar X_{min}}=\sqrt{\frac{MS_{组内}}{2}(\frac{1}{n_i}+\frac{1}{n_j})}\)
检验统计量
\[
HSD=q_{(k,\nu_{组内}),1-\alpha} \times S_{\bar X_{max}-\bar X_{min}},\nu_{组内}=k(n_j-1)
\]
对于任意i,j且\(\bar X_i＞\bar X_j\) ，如果\(\bar X_i-\bar X_j＞HSD\) ，那么拒绝\(H_0\) ，说明这两组存在显著差异。
\[
H_0:\mu_i=\mu_j(i≠j)
\]
Show the code posthoc_tukey <- TukeyHSD ( df_aov ) posthoc_tukey #> Tukey multiple comparisons of means #> 95% family-wise confidence level #> #> Fit: aov(formula = value ~ level, data = df_long) #> #> $level #> diff lwr upr p adj #> middle-high 14.71333 -7.664142 37.09081 0.2542971 #> normal-high 68.58667 46.209192 90.96414 0.0000000 #> normal-middle 53.87333 31.495858 76.25081 0.0000037 library ( multcomp ) posthoc_tukey_glht <- glht ( df_aov , linfct = mcp ( level = "Tukey" ) ) summary ( posthoc_tukey_glht ) #> #> Simultaneous Tests for General Linear Hypotheses #> #> Multiple Comparisons of Means: Tukey Contrasts #> #> #> Fit: aov(formula = value ~ level, data = df_long) #> #> Linear Hypotheses: #> Estimate Std. Error t value Pr(>|t|) #> middle - high == 0 14.71 9.12 1.613 0.254 #> normal - high == 0 68.59 9.12 7.521 <0.001 *** #> normal - middle == 0 53.87 9.12 5.907 <0.001 *** #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #> (Adjusted p values reported -- single-step method)
LSD-t test
least significant difference t-test
挑选任意感兴趣的两组进行比较
\[
H_0:\mu_i=\mu_j(i≠j)
\]
\[
LSD-t=\frac{\bar X_i-\bar X_j}{\sqrt{MS_{组内}(\frac{1}{n_i}+\frac{1}{n_j})}} \sim t(\nu) \ ,\ \nu=\nu_{组内}=n-k,a=k
\]
Show the code library ( agricolae ) # LSD 检验 lsd_result <- LSD.test ( df_aov , # 方差分析模型
"level" , # 分组变量
p.adj = "none" , # LSD 不校正（真实 LSD 定义）
console = TRUE # 直接输出结果
) #> #> Study: df_aov ~ "level" #> #> LSD t Test for value #> #> Mean Square Error: 498.996 #> #> level, means and individual ( 95 %) CI #> #> value std r se LCL UCL Min Max Q25 #> high 224.7817 23.24461 12 6.448488 211.6621 237.9012 196.01 280.75 214.1725 #> middle 239.4950 18.72159 12 6.448488 226.3755 252.6145 200.87 269.30 233.8500 #> normal 293.3683 24.62068 12 6.448488 280.2488 306.4879 244.61 332.96 283.7925 #> Q50 Q75 #> high 219.510 234.500 #> middle 238.110 253.505 #> normal 293.795 309.120 #> #> Alpha: 0.05 ; DF Error: 33 #> Critical Value of t: 2.034515 #> #> least Significant Difference: 18.55384 #> #> Treatments with the same letter are not significantly different. #> #> value groups #> normal 293.3683 a #> middle 239.4950 b #> high 224.7817 b # LSD 值 = 18.55 （两组均值差值 >18.55 即为显著） # 正常钙组 与 中剂量组：字母不同（a vs b）→ 差异极显著（P<0.05） # 正常钙组 与 高剂量组：字母不同（a vs b）→ 差异极显著（P<0.05） # 中剂量组 与 高剂量组：字母相同（b vs b）→ 差异不显著（P>0.05）
Dunnett’s test
Dunnett’s test 也称为q’ -test,是两独立样本t-test的一种修正。Dunnett’s test 假设数据符合正态分布，并且各组的方差相等。 控制对照组（C）与其他每个实验组（T）比较。
\(H_0:\mu_C=\mu_T\)
\[
q'=\frac{\bar X_T-\bar X_C}{\sqrt{MS_{组内}(\frac{1}{n_T}+\frac{1}{n_C})}} \sim q'(\nu,a) \ \ \nu=\nu_{组内},a=k
\]
临界值 \(q'_{(a,\nu_E),1-\alpha/2}\)
Show the code library ( multcomp ) dunnett_result <- glht ( df_aov , linfct = mcp ( level = c ( "middle - normal = 0" , "high - normal = 0" ) ) ) # 查看 Dunnett's test 结果 summary ( dunnett_result ) #> #> Simultaneous Tests for General Linear Hypotheses #> #> Multiple Comparisons of Means: User-defined Contrasts #> #> #> Fit: aov(formula = value ~ level, data = df_long) #> #> Linear Hypotheses: #> Estimate Std. Error t value Pr(>|t|) #> middle - normal == 0 -53.87 9.12 -5.907 2.51e-06 *** #> high - normal == 0 -68.59 9.12 -7.521 2.38e-08 *** #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #> (Adjusted p values reported -- single-step method)
Welch’s ANOVA Test
Show the code welch_aov <- oneway.test ( value ~ level , data = df_long ,var.equal = F ) welch_aov #> #> One-way analysis of means (not assuming equal variances) #> #> data: value and level #> F = 26.721, num df = 2.000, denom df = 21.665, p-value = 1.417e-06 library ( rstatix ) df_long |> welch_anova_test ( value ~ level ) #> # A tibble: 1 × 7 #> .y. n statistic DFn DFd p method #> * <chr> <int> <dbl> <dbl> <dbl> <dbl> <chr> #> 1 value 36 26.7 2 21.7 0.00000142 Welch ANOVA # 方差不齐时，两两 Welch t 检验（不合并方差，输出p值） pairwise.t.test ( df_long $ value ,
df_long $ level ,
p.adjust.method = "none" , # 不校正 = LSD风格
pool.sd = FALSE # 不使用合并方差，即两两Welch t
) #> #> Pairwise comparisons using t tests with non-pooled SD #> #> data: df_long$value and df_long$level #> #> high middle #> middle 0.1 - #> normal 4.9e-07 6.0e-06 #> #> P value adjustment method: none
数据变换
当方差分析的正态性假设或方差齐性假设不为真时，通常使用(1)数据变换方法；(2)非参数检验方法 比较均值差异
平方根变换
当每个水平组的方差与均值成比例，尤其是样本来自泊松分布
\[
Y=\sqrt{X}
\]
当数据中有零或非常小的值时，
\[
Y=\sqrt{X+a} \ \ \ \ a=0.5或0.1
\]
对数变换
当数据方差不齐且每个水平组标准差与均值成比例时
\[
Y=\log{X} \ \ \ \ base=e或10
\]
当数据中有零或负值时，
\[
Y=\log{(X+a)} \ \ \ \ a为实数，使得X+a>0
\]
反正弦平方根变换
率，服从二项分布\(B(n,\pi)\)
\[
Y=\arcsin {\sqrt{\pi}} \ \ \ \
\]
