# ボンフェローニ補正
(1 - (1 - 0.05/2)^2) * 100[1] 4.9375Rで統計的仮説検定:共分散分析
Rでデータサイエンス
共分散分析
- 「傾きに差は無い」との帰無仮説が棄却されない場合に、「切片に差は無い」との帰無仮説の検定に進む。
- 2回の検定の各有意水準は2.5%とする。
サンプルデータ
- 傾きは同一、切片は異なるサンプルデータを作成する。
library(dplyr)
library(ggplot2)
library(ggforce)
seed <- 20230915
set.seed(seed = seed)
n <- 20
x1 <- rnorm(n = n)
x2 <- rnorm(n = n)
x3 <- rnorm(n = n)
b1 <- 2
b2 <- 2
b3 <- 2
a1 <- 1
a2 <- 2
a3 <- 3
y1 <- b1 * x1 + a1 + rnorm(n)
y2 <- b2 * x2 + a2 + rnorm(n)
y3 <- b3 * x3 + a3 + rnorm(n)
df_sample <- data.frame(group = rep(c("A", "B", "C"), each = n), y = c(y1, y2, y3), x = c(x1, x2, x3))
df_sample$group <- df_sample$group %>%
factor(levels = c("A", "B", "C"))
glimpse(df_sample)Rows: 60
Columns: 3
$ group <fct> A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, B, B…
$ y <dbl> -0.870158643, 4.040351797, 1.436650434, 0.007324155, 3.715461872…
$ x <dbl> -0.850201613, 1.882664975, 0.009093009, 0.002701531, 0.167838678…point_color <- (scales::hue_pal())(3)
# scales::show_col(point_color)
g1 <- ggplot(data = df_sample, mapping = aes(x = x, y = y, col = group, shape = group)) + geom_point(size = 3) + scale_color_manual(values = point_color) + theme_minimal() + geom_mark_ellipse(mapping = aes(fill = group), alpha = 0.1, expand = unit(0.5, "mm"), show.legend = F) + geom_hline(yintercept = 0) + geom_vline(xintercept = 0)
g1 + geom_smooth(method = "lm", se = F, show.legend = F)
傾きの平行性の検定
\[ \begin{align} y_1=a_1+b_1x_1\\ y_2=a_2+b_2x_2\\ y_3=a_3+b_3x_3\\ \end{align} \] であるとき、
\[b_1=b_2=b_3\] を帰無仮説とする。
各群ごとの平均値
(df_mean <- data.frame(group = c("A", "B", "C"), x = c(mean(x1), mean(x2), mean(x3)), y = c(mean(y1), mean(y2), mean(y3)))) group x y
1 A 0.1999866 1.418299
2 B 0.1458300 2.488350
3 C -0.1749465 2.673628# (参考) 各群の平均値を追記した散布図
# サークルで囲んだ3つのポイントが追記した各群の平均値(以降同様)
g2 <- g1 + geom_point(mapping = aes(x = x, y = y, col = group), data = df_mean, size = 10, shape = 21, stroke = 2, show.legend = F) + geom_point(mapping = aes(x = x, y = y, col = group), data = df_mean, size = 3, show.legend = F)
g2 + geom_smooth(method = "lm", se = F, show.legend = F)
傾きを共通とする回帰直線
# 群間差の影響(切片)を取り除くために各群の平均値を(0,0)にする。
x1_corr <- x1 - mean(x1)
x2_corr <- x2 - mean(x2)
x3_corr <- x3 - mean(x3)
y1_corr <- y1 - mean(y1)
y2_corr <- y2 - mean(y2)
y3_corr <- y3 - mean(y3)
df_sample_corr <- data.frame(group = rep(c("A", "B", "C"), each = n), y = c(y1_corr, y2_corr, y3_corr), x = c(x1_corr, x2_corr, x3_corr))
df_sample_corr$group <- df_sample_corr$group %>%
factor(levels = c("A", "B", "C"))
glimpse(df_sample)Rows: 60
Columns: 3
$ group <fct> A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, B, B…
$ y <dbl> -0.870158643, 4.040351797, 1.436650434, 0.007324155, 3.715461872…
$ x <dbl> -0.850201613, 1.882664975, 0.009093009, 0.002701531, 0.167838678…# 補正により各群の平均値は(0,0)となる。
# xの平均値
aggregate(x = df_sample_corr$x, by = list(df_sample_corr$group), mean) Group.1 x
1 A -0.000000000000000006938894
2 B 0.000000000000000026388126
3 C 0.000000000000000001393200# yの平均値
aggregate(x = df_sample_corr$y, by = list(df_sample_corr$group), mean) Group.1 x
1 A 0.00000000000000007775898
2 B -0.00000000000000009989839
3 C -0.00000000000000019978593# 平均値を補正した散布図
g3 <- ggplot(data = df_sample_corr, mapping = aes(x = x, y = y, col = group, shape = group)) + geom_point(size = 3) + scale_color_manual(values = point_color) + theme_minimal() + geom_mark_ellipse(mapping = aes(fill = group), alpha = 0.1, expand = unit(0.5, "mm"), show.legend = F) + geom_hline(yintercept = 0) + geom_vline(xintercept = 0)
g3 + geom_smooth(method = "lm", se = F, show.legend = F)
# 各群の回帰直線は全て原点(0,0)を通過する。
# 黒色実践が傾きを共通とした回帰直線
g3 + geom_smooth(mapping = aes(x = x, y = y), data = df_sample_corr, method = "lm", se = F, show.legend = F, inherit.aes = F, col = "black")
共通の傾きの算出
# 傾きを共通とする回帰直線の係数
(result <- lm(y ~ x, data = df_sample_corr) %>%
summary())
# 切片は0となる。
result$coefficients[2, ]
# 以下では手作業で同傾きを算出する。
Call:
lm(formula = y ~ x, data = df_sample_corr)
Residuals:
Min 1Q Median 3Q Max
-2.56347 -0.58252 -0.06709 0.63856 2.35940
Coefficients:
Estimate Std. Error t value
(Intercept) -0.0000000000000002293 0.1340720271235040506 0.00
x 1.9359218490185201578 0.1345446550661526242 14.39
Pr(>|t|)
(Intercept) 1
x <0.0000000000000002 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.039 on 58 degrees of freedom
Multiple R-squared: 0.7812, Adjusted R-squared: 0.7774
F-statistic: 207 on 1 and 58 DF, p-value: < 0.00000000000000022
Estimate Std. Error
1.935921849018520157770240075 0.134544655066152624156572415
t value Pr(>|t|)
14.388693836010581605933111859 0.000000000000000000008582622 # 偏差平方和
Sxx1 <- sum((x1 - mean(x1))^2)
Sxx2 <- sum((x2 - mean(x2))^2)
Sxx3 <- sum((x3 - mean(x3))^2)
list(Sxx1 = Sxx1, Sxx2 = Sxx2, Sxx3 = Sxx3)$Sxx1
[1] 14.93417
$Sxx2
[1] 25.33319
$Sxx3
[1] 19.31184# 偏差積和
Sxy1 <- sum((x1 - mean(x1)) * (y1 - mean(y1)))
Sxy2 <- sum((x2 - mean(x2)) * (y2 - mean(y2)))
Sxy3 <- sum((x3 - mean(x3)) * (y3 - mean(y3)))
list(Sxy1 = Sxy1, Sxy2 = Sxy2, Sxy3 = Sxy3)$Sxy1
[1] 26.19793
$Sxy2
[1] 51.25051
$Sxy3
[1] 37.89225# 共通の傾き
(b_hat <- (Sxy1 + Sxy2 + Sxy3)/(Sxx1 + Sxx2 + Sxx3))[1] 1.935922# 共通の傾きの下での切片
a1_hat <- mean(y1) - b_hat * mean(x1)
a2_hat <- mean(y2) - b_hat * mean(x2)
a3_hat <- mean(y3) - b_hat * mean(x3)
list(a1_hat = a1_hat, a2_hat = a2_hat, a3_hat = a3_hat)
# 切片の違いがグループ間の違いを示す。$a1_hat
[1] 1.03114
$a2_hat
[1] 2.206035
$a3_hat
[1] 3.012311# (参考)傾き共通としたときの各群の回帰直線
g2 + stat_function(fun = function(x) x * b_hat + a1_hat, color = point_color[1], size = 1, show.legend = F) + stat_function(fun = function(x) x * b_hat + a2_hat, color = point_color[2], size = 1, show.legend = F) + stat_function(fun = function(x) x * b_hat + a3_hat, color = point_color[3], size = 1, show.legend = F) + geom_hline(yintercept = 0) + geom_vline(xintercept = 0)
# 各群の回帰直線はそれぞれの群の平均値を通過する。
F検定
# 共通の傾きの場合の残差平方和
S11 <- sum((y1 - (x1 * b_hat + a1_hat))^2)
S12 <- sum((y2 - (x2 * b_hat + a2_hat))^2)
S13 <- sum((y3 - (x3 * b_hat + a3_hat))^2)
S1 <- S11 + S12 + S13
list(S1 = S1, S11 = S11, S12 = S12, S13 = S13)$S1
[1] 62.55407
$S11
[1] 17.42738
$S12
[1] 22.79732
$S13
[1] 22.32937# 各群ごとの残差平方和(各群ごとの傾き)
S21 <- sum((lm(y1 ~ x1)$resi)^2)
S22 <- sum((lm(y2 ~ x2)$resi)^2)
S23 <- sum((lm(y3 ~ x3)$resi)^2)
S2 <- S21 + S22 + S23
list(S2 = S2, S21 = S21, S22 = S22, S23 = S23)$S2
[1] 61.85545
$S21
[1] 16.93436
$S22
[1] 22.60498
$S23
[1] 22.31611- 傾き共通の回帰式の残差平方和(S1)と個別の回帰式の残差平方和(S2)の差が個別の回帰式の残差平方和(S2)より有意に小さい場合、全ての回帰式の傾きは等しい(平行である)との帰無仮説は棄却されない。
# 回帰式の自由度
# 共通式
phi_b1 <- 1 + 3 # 共通の傾き1つと個別の切片3つ。
# 個別式
phi_b2 <- (1 + 1) * 3 # 3群の回帰式それぞれに傾き1つと切片1つ。
list(phi_b1 = phi_b1, phi_b2 = phi_b2)$phi_b1
[1] 4
$phi_b2
[1] 6# 残差の自由度
# 共通式
phi_e1 <- nrow(df_sample) - phi_b1
# 個別式
phi_e2 <- nrow(df_sample) - phi_b2
list(phi_e1 = phi_e1, phi_e2 = phi_e2)$phi_e1
[1] 56
$phi_e2
[1] 54# 残差の差の平方平均
MSA <- (S1 - S2)/(phi_e1 - phi_e2)
# 個別回帰式の残差の平方平均
MSe <- S2/phi_e2
list(MSA = MSA, MSe = MSe)$MSA
[1] 0.3493118
$MSe
[1] 1.145471# 分散比(F値)
(f_value <- MSA/MSe)[1] 0.3049503# p値
pf(q = f_value, df1 = phi_e1 - phi_e2, df2 = phi_e2, lower.tail = FALSE)[1] 0.7384211- 「3群個別の回帰式の傾きに差は無い」との帰無仮説は棄却されない。
切片の差の検定
\[ \begin{align} y_1=a_1+bx_1\\ y_2=a_2+bx_2\\ y_3=a_3+bx_3\\ \end{align} \] であるとき、
\[a_1=a_2=a_3\] を帰無仮説とする。
- 3つの回帰式の傾きと切片いずれにも差が無いのであれば、3群全てのデータ(以降、全データ)をまとめた回帰の残差と群毎の回帰の残差に差がないことになる。
全データの回帰
(common_result <- lm(y ~ x, data = df_sample))
Call:
lm(formula = y ~ x, data = df_sample)
Coefficients:
(Intercept) x
2.09 1.82 F検定
# 全データの残差平方和
S31 <- sum((predict(object = common_result, newdata = data.frame(x = x1)) - y1)^2)
S32 <- sum((predict(object = common_result, newdata = data.frame(x = x2)) - y2)^2)
S33 <- sum((predict(object = common_result, newdata = data.frame(x = x3)) - y3)^2)
S3 <- S31 + S32 + S33
list(S3 = S3, S31 = S31, S32 = S32, S33 = S33)$S3
[1] 101.4339
$S31
[1] 38.44129
$S32
[1] 24.00485
$S33
[1] 38.98772# 全データの回帰式の自由度
phi_b3 <- 1 + 1 # 共通回帰式 y=bx+a の傾き1つと切片1つ# 全データの残差の自由度
(phi_e3 <- nrow(df_sample) - phi_b3)[1] 58# 平方平均
# 傾き共通回帰式の残差
MSe1 <- S1/phi_e1
# 全データの回帰式による残差の増加
MSA3 <- (S3 - S1)/(phi_e3 - phi_e1)
list(MSe1 = MSe1, MSA3 = MSA3)$MSe1
[1] 1.117037
$MSA3
[1] 19.43989# 分散比(F値)
(f_value <- MSA3/MSe1)[1] 17.40309# p値
pf(q = f_value, df1 = phi_e3 - phi_e1, df2 = phi_e1, lower.tail = FALSE)[1] 0.000001324446- 「切片に差は無い」との帰無仮説は棄却される。
参考引用資料
- https://chemstat.hatenablog.com/entry/2020/09/12/180346
- https://chemstat.hatenablog.com/entry/2020/10/04/075824
- https://chemstat.hatenablog.com/entry/2020/10/04/074929
- https://scientist-press.com/wp-content/uploads/2019/07/seminar3.pdf
- http://www.snap-tck.com/room04/c01/stat/stat13/stat1302.html
- https://stackoverflow.com/questions/8197559/emulate-ggplot2-default-color-palette
最終更新
Sys.time()[1] "2024-04-12 05:37:03 JST"R、Quarto、Package
R.Version()$version.string[1] "R version 4.3.3 (2024-02-29 ucrt)"quarto::quarto_version()[1] '1.4.542'packageVersion(pkg = "tidyverse")[1] '2.0.0'