TokyoR#111_ANOVA - Speaker Deck

Slide 1

Slide 1 text

#111 @kilometer00 2024.02.24 分散分析の基礎

Slide 2

Slide 2 text

Who！？ Who?

Slide 3

Slide 3 text

Who！？・ @kilometer ・特任教員 (Ph.D. Eng.) ・神経科学・⾏動計算論・データ可視化・R： ~ 15 years ・近況：つかまり⽴ち！？

Slide 4

Slide 4 text

宣伝!!（書籍の翻訳に参加しました。）絶賛販売中！

Slide 5

Slide 5 text

#111 @kilometer00 2024.02.24 分散分析の基礎

Slide 6

Slide 6 text

library(tidyverse) library(palmerpenguins) df <- penguins %>% na.omit() データの準備

Slide 7

Slide 7 text

パイプ演算⼦ 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}パッケージ

Slide 8

Slide 8 text

# A tibble: 5 × 13 variable n min max median q1 q3 iqr mad mean sd se ci 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() データの準備

Slide 9

Slide 9 text

# A tibble: 15 × 14 species variable n min max median q1 q3 iqr mad mean sd s % group_by(species) %>% rstatix::get_summary_stats() データの準備

Slide 10

Slide 10 text

ggplot(data = df) + aes(x = body_mass_g, y = bill_length_mm) + geom_point() + geom_smooth(method = "lm", se = FALSE)

Slide 11

Slide 11 text

データデータから推定された回帰直線残差 ! 𝑦 ε 𝑦 𝑥 𝑦 = ! 𝑎! + ! 𝑎" 𝑥 + ε ∑ ε# = 0, ∑ 𝑥$# ε# = 0 ! 𝑦 = 𝑦 − ε 線形回帰モデル

Slide 12

Slide 12 text

データデータから推定された回帰直線真の回帰直線誤差残差 ! 𝑦 ε 𝑢 𝑦 𝑥 𝑌 𝑌 = α! + α" 𝑥 + 𝑢 𝑦 = ! 𝑎! + ! 𝑎" 𝑥 + ε 𝑢~𝑁(0, σ%) ∑ ε# = 0, ∑ 𝑥$# ε# = 0 ! 𝑦 = 𝑦 − ε 線形回帰モデル

Slide 13

Slide 13 text

線形回帰モデル: 3つの変動残差変動 RSS 全変動 TSS 回帰変動 ESS 𝑅𝑆𝑆 = $ !"# $ (𝑦! − ( 𝑦)% 𝑇𝑆𝑆 = $ !"# $ (𝑦! − + 𝑦)% 𝐸𝑆𝑆 = $ !"# $ (( 𝑦 − + 𝑦)% データ推定値データ平均推定値平均 = -

Slide 14

Slide 14 text

線形回帰モデル: 相関性の検定帰無仮説 𝑯𝒐 : α𝟎 = % 𝒚, α𝟏 = 𝟎 これ

Slide 15

Slide 15 text

Slide 16

Slide 16 text

線形回帰モデル: 相関性の検定帰無仮説が正しいとすると、これ 𝑢~𝑁(0, σ!) 正規分布に従う

Slide 17

Slide 17 text

𝜒!分布標準正規分布𝒩(0,1)に従う独⽴な𝑘個の確率変数𝑋" , 𝑋% , … , 𝑋$ について、 𝑍 = & 9:; < 𝑋9 = なる 𝑍 の確率分布を 𝜒% 分布と呼び、 𝑍~𝜒=(𝑘) と表す。(𝑘を⾃由度と呼ぶ)

Slide 18

Slide 18 text

𝜒!分布正規分布𝒩(𝑢, σ%)に従う独⽴な𝑘個の確率変数𝑋" , 𝑋% , … , 𝑋$ について、 𝑍 = & 9:; < (𝑋9 −𝑢)= σ= なる𝑍は⾃由度𝑘 − 1の𝜒%分布に従う 𝑍~𝜒=(𝑘 − 1)

Slide 19

Slide 19 text

𝑇𝑆𝑆~χ=(𝑛 − 1) 線形回帰モデル: 相関性の検定帰無仮説が正しいとすると、これ 𝑢~𝑁(0, σ!) 正規分布に従う

Slide 20

Slide 20 text

線形回帰モデル: 相関性の検定帰無仮説が正しいとすると、 𝑅𝑆𝑆~𝜒!(𝑛 − 𝑘 − 1) 𝑇𝑆𝑆~𝜒!(𝑛 − 1) 𝐸𝑆𝑆~𝜒!(𝑘)

Slide 21

Slide 21 text

線形回帰モデル: 相関性の検定帰無仮説が正しいとすると、 𝑓 = 𝐸𝑆𝑆/𝑘 𝑅𝑆𝑆/(𝑛 − 𝑘 − 1) 無相関(傾き=0)なら0 残差分散0なら0 𝒇~𝑭(𝒌, 𝒏 − 𝒌 − 𝟏)

Slide 22

Slide 22 text

F分布 𝑌 = 𝑋; /𝑚 𝑋= /𝑛 互いに独⽴な𝑋" ~𝜒%(𝑚)、𝑋% ~𝜒%(𝑛)について、なる𝑌の確率分布をF分布と呼び、 𝑌~𝐹(𝑚, 𝑛) と表す。 (𝑚を第1⾃由度,𝑛を第2⾃由度と呼ぶ) カイ2乗分布

Slide 23

Slide 23 text

線形回帰モデル: 相関性の検定帰無仮説が正しいとすると、 𝒇~𝑭(𝒌, 𝒏 − 𝒌 − 𝟏) これ f[5,20]の確率密度

Slide 24

Slide 24 text

これ線形回帰モデル: 相関性の検定帰無仮説が正しいとすると、 𝒇~𝑭(𝒌, 𝒏 − 𝒌 − 𝟏) データから得られたf値 f[5,20]の確率密度

Slide 25

Slide 25 text

これ線形回帰モデル: 相関性の検定帰無仮説が正しいとすると、 𝒇~𝑭(𝒌, 𝒏 − 𝒌 − 𝟏) f[5,20]の確率密度データから得られたf値無相関に近いより“甚だしい”相関

Slide 26

Slide 26 text

これ線形回帰モデル: 相関性の検定帰無仮説が正しいとすると、 𝒇~𝑭(𝒌, 𝒏 − 𝒌 − 𝟏) f[5,20]の確率密度データから得られたf値無相関に近いより“甚だしい”相関上側の累積確率(p値)

Slide 27

Slide 27 text

線形回帰モデル: 相関性の検定帰無仮説が正しいとすると、データから得られたf値よりも⼤きいf値が得られる確率(p値)は⼗分に低い。データに対し帰無仮説を正しいとする前提が誤っていると考えられる。

Slide 28

Slide 28 text

線形回帰モデル: 相関性の検定帰無仮説 𝑯𝒐 : α𝟎 = % 𝒚, α𝟏 = 𝟎 これ帰無仮説で現在のデータを説明するのは不合理なので、仮説を棄却する。

Slide 29

Slide 29 text

線形回帰モデル: 相関性の検定帰無仮説 𝑯𝒐 : α𝟎 = % 𝒚, α𝟏 = 𝟎 これ↑が棄却されるとしたら、何かしらの相関(α𝟏 ≠ 𝟎)が⽰唆される。つまり変数𝒙によって変数𝒚の変動のうち少なくとも⼀部が説明されると主張できる。

Slide 30

Slide 30 text

p値の解釈検証される効果が存在しないと仮定した上で、現在のデータおよびそれよりも極端なデータが得られる確率の累積確率濫⽤に注意が必要

Slide 31

Slide 31 text

帰無仮説検定のよくある濫⽤ ❌ 相関が確かめられた対⽴仮説は証明されていない。まして、因果に⾔及するのは濫⽤。 ❌ p値がより低い⽅が効果が強い正しい：効果が強い場合にp値が低い不成⽴： p値が低い場合に効果が強い・p値はn数や偶然性に左右される・甚だしさの仮定が違うから⽐較できない

Slide 32

Slide 32 text

ここが変だよ帰無仮説検定帰無仮説は常に成り⽴たない。 → ピッタリゼロである確率はゼロ → どこまでもデータを増やせば必ず棄却される。差異を想定し、それが指定の⽔準で検出される適切な数の標本を得よう。

Slide 33

Slide 33 text

No content

Slide 34

Slide 34 text

分散分析 ANOVA (analysis of variance) 𝑓 = 𝐸𝑆𝑆/𝑘 𝑅𝑆𝑆/(𝑛 − 𝑘 − 1) 帰無仮説の元で成⽴する分散(変動の⼆乗和)の⽐を統計量(f値)として、帰無仮説を検定する。

Slide 35

Slide 35 text

分散分析 ANOVA (analysis of variance) 説明変数が離散量の場合

Slide 36

Slide 36 text

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

Slide 37

Slide 37 text

分散分析 ANOVA (analysis of variance) 説明変数が離散量の場合

Slide 38

Slide 38 text

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) 説明変数が離散量の場合

Slide 39

Slide 39 text

分散分析 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

Slide 40

Slide 40 text

分散分析 ANOVA (analysis of variance) 説明変数が離散量の場合 → 内部ではダミー変数が作られ線形回帰されている Gentoo [0, 0, 1] Chinstrap [0,1,0] Adelie [0,0,0] bill_length_mm

Slide 41

Slide 41 text

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 ダミー変数

Slide 42

Slide 42 text

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

Slide 43

Slide 43 text

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

Slide 44

Slide 44 text

分散分析 ANOVA (analysis of variance) 説明変数が離散量の場合 → 内部では離散量をダミー変数に変換したうえで連続量の時と同じ計算を⾏っている。 → ただしANOVAの関⼼は、全体としての相関の有無だけではない。それぞれの変数の分散が⽬的変数の分散を有意に説明しているかが知りたい。

Slide 45

Slide 45 text

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) それぞれの変数の分散が⽬的変数の分散を有意に説明しているかが知りたい。主効果交互作⽤

Slide 46

Slide 46 text

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) それぞれの変数の分散が⽬的変数の分散を有意に説明しているかが知りたい。

Slide 47

Slide 47 text

分散分析 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

Slide 48

Slide 48 text

Slide 49

Slide 49 text

分散分析 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

Slide 50

Slide 50 text

Slide 51

Slide 51 text

分散分析 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 個別具体的にどの群とどの群の間に「有意な差」があるか知りたい

Slide 52

Slide 52 text

df %>% mutate( tag = str_c(species,":", sex) ) %>% rstatix::t_test( bill_length_mm ~ tag, p.adjust.method = "bonferroni" ) Post-hoc test t-testのBonferroni補正 (有意⽔準を調整)

Slide 53

Slide 53 text

df %>% rstatix::tukey_hsd( bill_length_mm ~ species * sex ) Post-hoc test Tukey Honest Significant Differences Test （正規性・等分散性を仮定）

Slide 54

Slide 54 text

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版）

Slide 55

Slide 55 text

分散分析 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 * * *

Slide 56

Slide 56 text

分散分析 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 * *

Slide 57

Slide 57 text

ちょっと待った、

Slide 58

Slide 58 text

bill_length_mmの変動には body_mass_gも寄与している。

Slide 59

Slide 59 text

共分散分析 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の影響を考慮した群間⽐較

Slide 60

Slide 60 text

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))

Slide 61

Slide 61 text

body_mass_g + species * sex このモデルを可視化する。

Slide 62

Slide 62 text

body_mass_g + species * sex → 傾きは全群共通, 切⽚は各群異なる

Slide 63

Slide 63 text

body_mass_g + species * sex → 傾きは全群共通, 切⽚は各群異なる → これは「良いモデル」だろうか？

Slide 64

Slide 64 text

モデル選択 (model selection) 定量的な基準に基づいて「良いモデル」を選び出そう。⾚池情報量基準 (AIC) 同⼀データに対して回帰された複数のモデルについて、 AICがより⼩さいモデルは、新規データに対する予測精度がより⾼いと⾔える。

Slide 65

Slide 65 text

モデル選択 (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が影響を与える様々なモデル

Slide 66

Slide 66 text

> df_fit # A tibble: 20 × 5 data models fit AIC dAIC 1 1939. 448. 2 1539. 48.4 3 1936. 445. 4 1495. 4.19 5 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)} )

Slide 67

Slide 67 text

body_mass_g + species * sex → 無事に「良いモデル」だった。

Slide 68

Slide 68 text

モデル選択 (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 全ての交互作⽤を仮定

Slide 69

Slide 69 text

解析の流れ 1. 興味のある⽬的変数の分布を可視化 2. 影響を与える変動要因を考慮してモデル選択 3. 最良モデルを可視化 4. 仮説に従って統計的検定

Slide 70

Slide 70 text

Enjoy!!