TokyoR#111_ANOVA

第111回Tokyo.Rで喋った際のスライド資料です。

Transcript

1. #111 @kilometer00 2024.02.24 分散分析の基礎
2. Who！？ Who?
3. Who！？ ・ @kilometer ・特任教員 (Ph.D. Eng.) ・神経科学 ・⾏動計算論 ・データ可視化 ・R： ~ 15 years ・近況：つかまり⽴ち！？
4. 宣伝!!（書籍の翻訳に参加しました。） 絶賛販売中！
5. #111 @kilometer00 2024.02.24 分散分析の基礎
6. library(tidyverse) library(palmerpenguins) df <- penguins %>% na.omit() データの準備
7. パイプ演算⼦ pipe X %>% f() f(X) X %>% f(y) f(X, y) X %>% f %>% g g(f(X)) X %>% f(y, .) f(y, X) 直前の項を直後の関数の第⼀引数（or 指定引数）に取る {magrittr}パッケージ
8. # A tibble: 5 × 13 variable n min max median q1 q3 iqr mad mean sd se ci <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> 1 bill_le… 333 32.1 59.6 44.5 39.5 48.6 9.1 6.97 44.0 5.47 0.3 0.59 2 bill_de… 333 13.1 21.5 17.3 15.6 18.7 3.1 2.22 17.2 1.97 0.108 0.212 3 flipper… 333 172 231 197 190 213 23 16.3 201. 14.0 0.768 1.51 4 body_ma… 333 2700 6300 4050 3550 4775 1225 890. 4207. 805. 44.1 86.8 5 year 333 2007 2009 2008 2007 2009 2 1.48 2008. 0.813 0.045 0.088 df %>% rstatix::get_summary_stats() データの準備
9. # A tibble: 15 × 14 species variable n min max median q1 q3 iqr mad mean sd s <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl 1 Adelie bill_len… 146 32.1 46 38.8 36.7 40.8 4.05 2.96 38.8 2.66 0.22 2 Adelie bill_dep… 146 15.5 21.5 18.4 17.5 19 1.5 1.19 18.3 1.22 0.10 3 Adelie flipper_… 146 172 210 190 186 195 9 7.41 190. 6.52 0.54 4 Adelie body_mas… 146 2850 4775 3700 3362. 4000 638. 445. 3706. 459. 38.0 5 Adelie year 146 2007 2009 2008 2007 2009 2 1.48 2008. 0.812 0.06 6 Chinstrap bill_len… 68 40.9 58 49.6 46.4 51.1 4.72 3.63 48.8 3.34 0.40 7 Chinstrap bill_dep… 68 16.4 20.8 18.4 17.5 19.4 1.9 1.41 18.4 1.14 0.13 8 Chinstrap flipper_… 68 178 212 196 191 201 10 7.41 196. 7.13 0.86 9 Chinstrap body_mas… 68 2700 4800 3700 3488. 3950 462. 371. 3733. 384. 46.6 10 Chinstrap year 68 2007 2009 2008 2007 2009 2 1.48 2008. 0.863 0.10 11 Gentoo bill_len… 119 40.9 59.6 47.4 45.4 49.6 4.25 3.11 47.6 3.11 0.28 12 Gentoo bill_dep… 119 13.1 17.3 15 14.2 15.8 1.55 1.19 15.0 0.986 0.09 13 Gentoo flipper_… 119 203 231 216 212 222. 9.5 7.41 217. 6.58 0.60 14 Gentoo body_mas… 119 3950 6300 5050 4700 5500 800 593. 5092. 501. 46.0 15 Gentoo year 119 2007 2009 2008 2007 2009 2 1.48 2008. 0.789 0.07 df %>% group_by(species) %>% rstatix::get_summary_stats() データの準備
10. ggplot(data = df) + aes(x = body_mass_g, y = bill_length_mm) + geom_point() + geom_smooth(method = "lm", se = FALSE)
11. データ データから推定された回帰直線 残差 ! 𝑦 ε 𝑦 𝑥 𝑦 = ! 𝑎! + ! 𝑎" 𝑥 + ε ∑ ε# = 0, ∑ 𝑥$# ε# = 0 ! 𝑦 = 𝑦 − ε 線形回帰モデル
12. データ データから推定された回帰直線 真の回帰直線 誤差 残差 ! 𝑦 ε 𝑢 𝑦 𝑥 𝑌 𝑌 = α! + α" 𝑥 + 𝑢 𝑦 = ! 𝑎! + ! 𝑎" 𝑥 + ε 𝑢~𝑁(0, σ%) ∑ ε# = 0, ∑ 𝑥$# ε# = 0 ! 𝑦 = 𝑦 − ε 線形回帰モデル
13. 線形回帰モデル: 3つの変動 残差変動 RSS 全変動 TSS 回帰変動 ESS 𝑅𝑆𝑆 = $ !"# $ (𝑦! − ( 𝑦)% 𝑇𝑆𝑆 = $ !"# $ (𝑦! − + 𝑦)% 𝐸𝑆𝑆 = $ !"# $ (( 𝑦 − + 𝑦)% データ 推定値 データ 平均 推定値 平均 = -
14. 線形回帰モデル: 相関性の検定 帰無仮説 𝑯𝒐 : α𝟎 = % 𝒚, α𝟏 = 𝟎 これ
15. データ データから推定された回帰直線 真の回帰直線 誤差 残差 ! 𝑦 ε 𝑢 𝑦 𝑥 𝑌 𝑌 = α! + α" 𝑥 + 𝑢 𝑦 = ! 𝑎! + ! 𝑎" 𝑥 + ε 𝑢~𝑁(0, σ%) ∑ ε# = 0, ∑ 𝑥$# ε# = 0 ! 𝑦 = 𝑦 − ε 線形回帰モデル
16. 線形回帰モデル: 相関性の検定 帰無仮説が正しいとすると、 これ 𝑢~𝑁(0, σ!) 正規分布に従う
17. 𝜒!分布 標準正規分布𝒩(0,1)に従う独⽴な𝑘個の 確率変数𝑋" , 𝑋% , … , 𝑋$ について、 𝑍 = & 9:; < 𝑋9 = なる 𝑍 の確率分布を 𝜒% 分布と呼び、 𝑍~𝜒=(𝑘) と表す。(𝑘を⾃由度と呼ぶ)
18. 𝜒!分布 正規分布𝒩(𝑢, σ%)に従う独⽴な𝑘個の 確率変数𝑋" , 𝑋% , … , 𝑋$ について、 𝑍 = & 9:; < (𝑋9 −𝑢)= σ= なる𝑍は⾃由度𝑘 − 1の𝜒%分布に従う 𝑍~𝜒=(𝑘 − 1)
19. 𝑇𝑆𝑆~χ=(𝑛 − 1) 線形回帰モデル: 相関性の検定 帰無仮説が正しいとすると、 これ 𝑢~𝑁(0, σ!) 正規分布に従う
20. 線形回帰モデル: 相関性の検定 帰無仮説が正しいとすると、 𝑅𝑆𝑆~𝜒!(𝑛 − 𝑘 − 1) 𝑇𝑆𝑆~𝜒!(𝑛 − 1) 𝐸𝑆𝑆~𝜒!(𝑘)
21. 線形回帰モデル: 相関性の検定 帰無仮説が正しいとすると、 𝑓 = 𝐸𝑆𝑆/𝑘 𝑅𝑆𝑆/(𝑛 − 𝑘 − 1) 無相関(傾き=0)なら0 残差分散0なら0 𝒇~𝑭(𝒌, 𝒏 − 𝒌 − 𝟏)
22. F分布 𝑌 = 𝑋; /𝑚 𝑋= /𝑛 互いに独⽴な𝑋" ~𝜒%(𝑚)、𝑋% ~𝜒%(𝑛)について、 なる𝑌の確率分布をF分布と呼び、 𝑌~𝐹(𝑚, 𝑛) と表す。 (𝑚を第1⾃由度,𝑛を第2⾃由度と呼ぶ) カイ2乗分布
23. 線形回帰モデル: 相関性の検定 帰無仮説が正しいとすると、 𝒇~𝑭(𝒌, 𝒏 − 𝒌 − 𝟏) これ f[5,20]の 確率密度
24. これ 線形回帰モデル: 相関性の検定 帰無仮説が正しいとすると、 𝒇~𝑭(𝒌, 𝒏 − 𝒌 − 𝟏) データから 得られたf値 f[5,20]の 確率密度
25. これ 線形回帰モデル: 相関性の検定 帰無仮説が正しいとすると、 𝒇~𝑭(𝒌, 𝒏 − 𝒌 − 𝟏) f[5,20]の 確率密度 データから 得られたf値 無相関に近い より“甚だしい”相関
26. これ 線形回帰モデル: 相関性の検定 帰無仮説が正しいとすると、 𝒇~𝑭(𝒌, 𝒏 − 𝒌 − 𝟏) f[5,20]の 確率密度 データから 得られたf値 無相関に近い より“甚だしい”相関 上側の累積確率(p値)
27. 線形回帰モデル: 相関性の検定 帰無仮説が正しいとすると、 データから得られたf値よりも⼤きいf値が 得られる確率(p値)は⼗分に低い。 データに対し帰無仮説を正しいとする前提 が誤っていると考えられる。
28. 線形回帰モデル: 相関性の検定 帰無仮説 𝑯𝒐 : α𝟎 = % 𝒚, α𝟏 = 𝟎 これ 帰無仮説で現在のデータを説明するのは 不合理なので、仮説を棄却する。
29. 線形回帰モデル: 相関性の検定 帰無仮説 𝑯𝒐 : α𝟎 = % 𝒚, α𝟏 = 𝟎 これ↑が棄却されるとしたら、 何かしらの相関(α𝟏 ≠ 𝟎)が⽰唆される。 つまり変数𝒙によって変数𝒚の変動のうち 少なくとも⼀部が説明されると主張できる。
30. p値の解釈 検証される効果が存在しないと仮定した 上で、現在のデータおよびそれよりも極 端なデータが得られる確率の累積確率 濫⽤に注意が必要
31. 帰無仮説検定のよくある濫⽤ ❌ 相関が確かめられた 対⽴仮説は証明されていない。 まして、因果に⾔及するのは濫⽤。 ❌ p値がより低い⽅が効果が強い 正しい： 効果が強い場合にp値が低い 不成⽴： p値が低い場合に効果が強い ・p値はn数や偶然性に左右される ・甚だしさの仮定が違うから⽐較できない
32. ここが変だよ帰無仮説検定 帰無仮説は常に成り⽴たない。 → ピッタリゼロである確率はゼロ → どこまでもデータを増やせば 必ず棄却される。 差異を想定し、それが指定の⽔準で 検出される適切な数の標本を得よう。
34. 分散分析 ANOVA (analysis of variance) 𝑓 = 𝐸𝑆𝑆/𝑘 𝑅𝑆𝑆/(𝑛 − 𝑘 − 1) 帰無仮説の元で成⽴する分散(変動の⼆乗和)の⽐を 統計量(f値)として、帰無仮説を検定する。
35. 分散分析 ANOVA (analysis of variance) 説明変数が離散量の場合
36. ggplot(data = df) + aes(x = species, y = bill_length_mm, color = species) + geom_boxplot(outlier.color = NA) + geom_jitter(width = 0.1, alpha = 0.4) + theme(legend.position = "none") Box-plot
37. 分散分析 ANOVA (analysis of variance) 説明変数が離散量の場合
38. df %>% lm(bill_length_mm ~ species, data = .) %>% summary() #> Call: #> lm(formula = bill_length_mm ~ species, data = .) ... #> F-statistic: 397.3 on 2 and 330 DF, #> p-value: < 2.2e-16 分散分析 ANOVA (analysis of variance) 説明変数が離散量の場合
39. 分散分析 ANOVA (analysis of variance) 説明変数が離散量の場合 → 内部ではダミー変数が作られ線形回帰されている df %>% mutate(s1 = if_else(species == "Chinstrap", 1, 0), s2 = if_else(species == "Gentoo", 1, 0)) %>% lm(bill_length_mm ~ s1 + s2, data = .) %>% summary() ... #> F-statistic: 397.3 on 2 and 330 DF, #> p-value: < 2.2e-16
40. 分散分析 ANOVA (analysis of variance) 説明変数が離散量の場合 → 内部ではダミー変数が作られ線形回帰されている Gentoo [0, 0, 1] Chinstrap [0,1,0] Adelie [0,0,0] bill_length_mm
41. df %>% lm(bill_length_mm ~ species + sex, data = .) %>% summary() #> F-statistic: 502.8 on 3 and 329 DF, #> p-value: < 2.2e-16 分散分析 ANOVA (analysis of variance) 説明変数が2つの離散量の場合 𝑏𝑖𝑙𝑙"#$%&'(( ~ 𝑠𝑝𝑒𝑐𝑖𝑒𝑠 + 𝑠𝑒𝑥 Adelie Chinstrap Gentoo s1 0 1 0 s2 0 0 1 female male m 0 1 ダミー変数
42. df %>% lm(bill_length_mm ~ species + sex + species:sex, data = .) %>% summary() #> F-statistic: 305 on 5 and 327 DF, #> p-value: < 2.2e-16 分散分析 ANOVA (analysis of variance) 説明変数が2つの離散量に交互作⽤が加わる場合 𝑏𝑖𝑙𝑙"#$%&'(( ~ 𝑠𝑝𝑒𝑐𝑖𝑒𝑠 + 𝑠𝑒𝑥 + 𝑠𝑝𝑒𝑐𝑖𝑒𝑠: 𝑠𝑒𝑥 Adelie Chinstrap Gentoo 0 1 0 0 0 1 female male 0 1 Adelie Chinstrap Gentoo m1 0 1 0 m2 0 0 1
43. df %>% lm(bill_length_mm ~ species * sex, data = .) %>% summary() #> F-statistic: 305 on 5 and 327 DF, #> p-value: < 2.2e-16 分散分析 ANOVA (analysis of variance) 説明変数が2つの離散量に交互作⽤が加わる場合 𝑏𝑖𝑙𝑙"#$%&'(( ~ 𝑠𝑝𝑒𝑐𝑖𝑒𝑠 + 𝑠𝑒𝑥 + 𝑠𝑝𝑒𝑐𝑖𝑒𝑠: 𝑠𝑒𝑥 Adelie Chinstrap Gentoo 0 1 0 0 0 1 Adelie Chinstrap Gentoo m1 0 1 0 m2 0 0 1 female male 0 1
44. 分散分析 ANOVA (analysis of variance) 説明変数が離散量の場合 → 内部では離散量をダミー変数に変換したうえで 連続量の時と同じ計算を⾏っている。 → ただしANOVAの関⼼は、全体としての相関の有無 だけではない。 それぞれの変数の分散が⽬的変数の分散を 有意に説明しているかが知りたい。
45. df %>% lm(bill_length_mm ~ species * sex, data = .) %>% car::Anova() #> Anova Table (Type II tests) #> #> Response: bill_length_mm #> Sum Sq Df F value Pr(>F) #> species 6975.6 2 650.4786 <2e-16 *** #> sex 1135.7 1 211.8066 <2e-16 *** #> species:sex 24.5 2 2.2841 0.1035 #> Residuals 1753.3 327 分散分析 ANOVA (analysis of variance) それぞれの変数の分散が⽬的変数の分散を 有意に説明しているかが知りたい。 主効果 交互作⽤
46. df %>% rstatix::anova_test( bill_length_mm ~ species * sex ) ANOVA Table (type II tests) Effect DFn DFd F p p<.05 ges 1 species 2 327 650.479 1.06e-114 * 0.799 2 sex 1 327 211.807 2.42e-37 * 0.393 3 species:sex 2 327 2.284 1.03e-01 0.014 分散分析 ANOVA (analysis of variance) それぞれの変数の分散が⽬的変数の分散を 有意に説明しているかが知りたい。
47. 分散分析 ANOVA (analysis of variance) それぞれの変数の分散が⽬的変数の分散を 有意に説明しているかが知りたい。 種差：f[2, 327] = 650.4, p = 1.1e-144, ηg 2 = 0.799 性差：f[1, 327] = 211.8, p = 2.4e-37, ηg 2 = 0.393
48. df %>% rstatix::anova_test( bill_length_mm ~ species * sex ) ANOVA Table (type II tests) Effect DFn DFd F p p<.05 ges 1 species 2 327 650.479 1.06e-114 * 0.799 2 sex 1 327 211.807 2.42e-37 * 0.393 3 species:sex 2 327 2.284 1.03e-01 0.014 分散分析 ANOVA (analysis of variance) それぞれの変数の分散が⽬的変数の分散を 有意に説明しているかが知りたい。 効果量 (generalized eta squared)
49. 分散分析 ANOVA (analysis of variance) 効果量 effect size: 標準化された差異の⼤きさ η= = 𝐸𝑆𝑆 𝑇𝑆𝑆 ηB = = 𝐸𝑆𝑆 𝐸𝑆𝑆 + 𝑇𝑆𝑆 ηC = = 𝐸𝑆𝑆 σ𝐸𝑆𝑆 + ∑DEF 𝑆𝑆DEF + ∑G 𝐾 Olejnik & Algina, Psychological methods (2003) partial eta squared generalized eta squared eta squared 回帰される全ての変数の変動 全ての誤差項 測定要因なら0 操作要因なら1
50. 分散分析 ANOVA (analysis of variance) それぞれの変数の分散が⽬的変数の分散を 有意に説明しているかが知りたい。 種差：f[2, 327] = 650.4, p = 1.1e-144, ηg 2 = 0.799 性差：f[1, 327] = 211.8, p = 2.4e-37, ηg 2 = 0.393
51. 分散分析 ANOVA (analysis of variance) それぞれの変数の分散が⽬的変数の分散を 有意に説明しているかが知りたい。 種差：f[2, 327] = 650.4, p = 1.1e-144, ηg 2 = 0.799 性差：f[1, 327] = 211.8, p = 2.4e-37, ηg 2 = 0.393 個別具体的にどの群とどの群の間に 「有意な差」があるか知りたい
52. df %>% mutate( tag = str_c(species,":", sex) ) %>% rstatix::t_test( bill_length_mm ~ tag, p.adjust.method = "bonferroni" ) Post-hoc test t-testのBonferroni補正 (有意⽔準を調整)
53. df %>% rstatix::tukey_hsd( bill_length_mm ~ species * sex ) Post-hoc test Tukey Honest Significant Differences Test （正規性・等分散性を仮定）
54. install.packages("devtools") devtools::install_github( "PhDMeiwp/Steel.Dwass.test@master", force = TRUE) df_sd <- df %>% mutate(tag = str_c(species,"_", sex)) Steel.Dwass.test::Steel.Dwass( x = df_sd$bill_length_mm, group = df_sd$tag) Post-hoc test Steel-Dwass test （Tukeyのnon-para版）
55. 分散分析 ANOVA (analysis of variance) それぞれの変数の分散が⽬的変数の分散を 有意に説明しているかが知りたい。 種差：f[2, 327] = 650.4, p = 1.1e-144, ηg 2 = 0.799 性差：f[1, 327] = 211.8, p = 2.4e-37, ηg 2 = 0.393 * * *
56. 分散分析 ANOVA (analysis of variance) それぞれの変数の分散が⽬的変数の分散を 有意に説明しているかが知りたい。 種差：f[2, 327] = 650.4, p = 1.1e-144, ηg 2 = 0.799 性差：f[1, 327] = 211.8, p = 2.4e-37, ηg 2 = 0.393 * *
57. ちょっと待った、
58. bill_length_mmの変動には body_mass_gも寄与している。
59. 共分散分析 ANCOVA (analysis of covariance) df %>% rstatix::anova_test( bill_length_mm ~ body_mass_g + species * sex ) 連続量(共変量) 離散量 離散量 交互作⽤あり ANOVA Table (type II tests) Effect DFn DFd F p p<.05 ges 1 body_mass_g 1 326 23.462 1.97e-06 * 0.067 2 species 2 326 468.046 1.50e-96 * 0.742 3 sex 1 326 49.702 1.07e-11 * 0.132 4 species:sex 2 326 4.061 1.80e-02 * 0.024 body_mass_gの影響を考慮した群間⽐較
60. body_mass_g + species * sex このモデルを可視化する。 fit <- df %>% lm(bill_length_mm ~ body_mass_g + species * sex, data = .) df_pred <- df %>% mutate(pred = predict(fit)) ggplot(df_pred) + aes(x = body_mass_g, y = bill_length_mm, color = species) + geom_point(aes(shape = sex)) + geom_line(aes(y = pred, linetype = sex))
61. body_mass_g + species * sex このモデルを可視化する。
62. body_mass_g + species * sex → 傾きは全群共通, 切⽚は各群異なる
63. body_mass_g + species * sex → 傾きは全群共通, 切⽚は各群異なる → これは「良いモデル」だろうか？
64. モデル選択 (model selection) 定量的な基準に基づいて「良いモデル」を 選び出そう。 ⾚池情報量基準 (AIC) 同⼀データに対して回帰された複数のモデルについて、 AICがより⼩さいモデルは、新規データに 対する予測精度がより⾼いと⾔える。
65. モデル選択 (model selection) models <- list( function(df){ lm(y ~ x, data = df)}, function(df){ lm(y ~ x + species, data = df)}, function(df){ lm(y ~ x + sex, data = df)}, function(df){ lm(y ~ x + species + sex, data = df)}, function(df){ lm(y ~ x + species * sex, data = df)} ) df_fit <- df %>% mutate(y = bill_length_mm, x = body_mass_g) %>% group_nest() %>% mutate(models = list(models)) %>% unnest(models) %>% mutate(fit = map2(models, data, ~.x(.y))) %>% mutate(AIC = map_dbl(fit, AIC)) %>% mutate(dAIC = AIC - min(AIC)) xとyの相関に対して、傾き・切⽚にspecies, sexが影響を与える様々なモデル
66. > df_fit # A tibble: 20 × 5 data models fit AIC dAIC <list> <list> <list> <dbl> <dbl> 1 <tibble [333 × 10]> <fn> <lm> 1939. 448. 2 <tibble [333 × 10]> <fn> <lm> 1539. 48.4 3 <tibble [333 × 10]> <fn> <lm> 1936. 445. 4 <tibble [333 × 10]> <fn> <lm> 1495. 4.19 5 <tibble [333 × 10]> <fn> <lm> 1491. 0 models <- list( function(df){ lm(y ~ x, data = df)}, function(df){ lm(y ~ x + species, data = df)}, function(df){ lm(y ~ x + sex, data = df)}, function(df){ lm(y ~ x + species + sex, data = df)}, function(df){ lm(y ~ x + species * sex, data = df)} )
67. body_mass_g + species * sex → 無事に「良いモデル」だった。
68. モデル選択 (model selection) 線形回帰モデルの範囲であれば MASS::stepAIC()関数が便利。 lm(bill_length_mm ~ body_mass_g * sex * species, data = df) %>% MASS::stepAIC() %>% .$anova ... Final Model: bill_length_mm ~ body_mass_g + sex + species + sex:species sex * species 全ての交互作⽤を仮定
69. 解析の流れ 1. 興味のある⽬的変数の分布を可視化 2. 影響を与える変動要因を考慮して モデル選択 3. 最良モデルを可視化 4. 仮説に従って統計的検定
70. Enjoy!!