1. 課題を見つける／仮説を立てる

2. 実験や観察をしてデータを集める

3. データを整理する

4. データを解析して仮説を検証する

• 仮説の検証にはエビデンスが必要
• でも生物の表現型には ばらつきがある
• 実験や観察で見られた差や効果は偶然の産物じゃないの？
  • これに答えるのが統計解析
  • 推論や予測などにも使うよ

統計解析の基本的な考え方

• 基礎知識と心構え
• 統計解析法の選び方

Rのみ（課金なし）でできる統計解析の実践

• さまざまな統計解析用パッケージ
• より発展的な解析の足がかりへ

この実習で扱わないこと

• 理論背景や難しい計算
• 発展的な解析（ベイズ統計とか）

1. まず、自分が 何を示したいのかを整理する

   • 「A群とB群は異なる」とか「Aが増えるとBも増える」とか

2. 自分のデータの性質を知り、適切な統計手法を探す

   • 量的（連続/離散）か質的か、分布や分散など
     • 帰無仮説や統計モデルを設定する

3. 統計手法を実践する方法を探す

   • 手計算（大変）
   • 統計ソフトウェア（簡単だが有料、選択肢少ない）
   • Rのパッケージ（無料、選択肢が多い）

手近にできるからと言って不適切な統計方法を使うのはダメ！

統計解析…標本に基づいて、母集団の性質を推定する

• 母集団の性質を知りたいが実測は不可能
• ランダム抽出した標本（実測値など）から母集団を 推定
• 標本の ばらつきを考慮する必要→統計解析の出番！

• 応答変数（Y）（または目的変数, 従属変数）…調べたい値
• 説明変数（X）（または独立変数）…Yのばらつきを説明する値

• 質的変数（ex. A, B, C…）
  • 名義尺度（電話番号、背番号、バスの系統番号など）
    • 順序もないし加減などの演算もできない
  • 順序尺度（レベルわけ、レースの着順など）
    • 順序による比較ができるが加減などの演算はできない
• 量的変数（ex. -2, -1.3, 0, 0.5, 0.7, 1.1…）
  • 間隔尺度（日付、摂氏温度や華氏温度など）
    • 目盛が等間隔で間隔に意味がある
    • 加減の演算にも意味がある
  • 比例尺度（質量、長さ、エネルギー、絶対温度など）
    • 間隔と比率に意味がある
    • 比にも、乗除の演算にも意味がある

応答変数（Y）も説明変数（X）も量的な場合…

• モデリング
  • 一般化線形モデル（直線回帰／ポアソン回帰など）
  • ベイズ統計モデル…さらに進んだ人に
• （相関分析）…素朴すぎて限界がある

応答変数（Y）は質的、説明変数（X）は量的な場合…

• 割合データは質的変数を変形したもの

• モデリング

  • 一般化線形モデル（ロジスティック回帰）
    • 応答変数のカテゴリ数が２個のみ（３以上は難しい）

応答変数（Y）は量的、説明変数（X）は質的な場合…

• モデリング…量的と質的が混じる場合はこっち
  • 一般化線形モデル（説明変数をダミー変数として扱う）
• 仮説検定
  • ２群（ t検定, Wilcoxon検定, Brunner-Munzel検定）
  • 多群（ ANOVA, ART-ANOVA）

応答変数（Y）も説明変数（X）も質的な場合…

• 仮説検定
  • 割合の検定（独立性の検定）
    • χ（カイ）二乗検定
    • フィッシャーの正確確率検定

1. まず、自分が何を示したいのかを整理する

   • 「A群とB群は異なる」とか「Aが増えるとBも増える」とか

2. 自分のデータの性質を知り、使えそうな統計手法を探す

3. 統計手法を実践する方法を探す

   • まず仮説検定かモデリングかを決める
   • それぞれのなかから適切な方法を探す
     • 仮説検定←こっちから実践してみましょう
     • モデリング（一般化線形モデル）

• 例えばこんなとき…
• あなたは「良い餌で育ったハエは大きくなる」ことを主張しようとしてデータを集めた。
• これまでのデータによると、普通の餌で育てたハエの体長は平均3.0mm, 標準偏差1.0mmの正規分布に従う。
• あなたが良い餌で育てたハエ10匹の体長は、3.0, 3.0, 3.2, 3.4, 3.6, 4.0, 4.2, 4.2, 4.3, 4.5 mmであった。
• この結果から、あなたは良い餌で育ったハエは大きくなると結論づけた。
• ところが教授はこれに対して 「たまたま大きいハエがとれたんじゃないの？」とコメントした。
• さて、この教授をどのように納得させたらよいだろうか？

「こんな標本が偶然得られる可能性は低いことを示す」

帰無仮説 （普通の餌と変わらない母集団から偶然得られた）

対立仮説 （偶然ではない＝異なる母集団から得られた）

※帰無仮説と対立仮説は排他的

1. 帰無仮説のもとで 今回得られた標本が偶然生じる確率（\(p\)値, \(p\)-value）を計算する
2. \(p\)値が有意水準（\(\alpha\)）より低い場合、帰無仮説を棄却する
   • 伝統的に\(p\) < 0.05で棄却（偶然起こるのは1/20回以下）
3. その結果、 対立仮説を採用する（=「有意差がある」）

\(p\)値は検定統計量に基づいて求める

• 検定統計量とは…
  • 仮説検定で用いられる統計量（いろいろある）の総称
  • めっちゃざっくり言うと：とある前提に基づいて標本の性質をまとめた 仮説検定の基準になる値
  • 標本から計算された値を 検定統計量の実現値と呼ぶ
• 検定統計量の実現値と同じか、さらに極端な値が偶然得られる確率（=\(p\)値）を計算する
  • 帰無仮説のもとでの検定統計量の確率分布を元に計算する
• \(p\)値が有意水準（\(\alpha=0.05\)）より低ければ帰無仮説を棄却する

• 確率分布とは…ある確率変数がどんな確率でどんな値をとるかを表した分布
• 確率変数とは…実際に結果が得られるまでどんな値が得られるかが決まっていない変数
  • 同じ手続きでデータを得ても異なる値が得られるもの = 多くの実験／観測データはこれに当てはまる
• 確率の総和は１となる
• さまざまな典型的な確率分布が知られる

1. 母集団に関する帰無仮説と対立仮説を設定する
2. どの検定統計量（=検定方法）を使うかを選ぶ
   • 検定の目的と標本の性質によって変わる
     • \(Z\)統計量、\(t\)統計量、\(F\)統計量、カイ二乗統計量など
3. 有意水準\(\alpha\)を決める（伝統的には\(\alpha = 0.05\)）
4. 標本から 検定統計量の実現値を求める
5. 帰無仮説のもとでの 検定統計量の確率分布と 検定統計量の実現値を比較し、 検定統計量の実現値が帰無仮説のもとで偶然得られる確率（=\(p\)値）を計算する
   
6. \(p\)値が有意水準より低ければ帰無仮説を棄却する
   

• 4と5はRが自動的にやってくれます

• 帰無仮説：平均（=3.0）と標準偏差（=1.0）が正規分布に従う母集団から標本が得られた
• この場合、検定統計量として\(Z\)統計量が使える
  • \(Z = \frac{\bar{X} - \mu}{\sqrt{\sigma^2/n}}\)…平均値の差を標準誤差\({\sqrt{\sigma^2/n}}\)で標準化した値
  • \(Z\)は標準正規分布（平均0, 標準偏差１の正規分布）に従う
• \(Z\)（検定統計量）の実現値を計算する（\(\frac{3.74 - 3}{\sqrt{1^2/10}}\) = 2.340085）
  • \(\bar{X}\)…標本の平均 = 3.74
  • \(\mu\)…母集団の平均 = 3.0
  • \(\sigma\)…母集団の標準偏差 = 1.0
  • \(n\)…サンプルサイズ = 10

• 標準正規分布では、\(Z\)は95%の割合で灰色の範囲（-1.96 < \(Z\) < 1.96）の値をとる

• 今回の場合、\(Z=2.340085\) ←95%の範囲外！
• \(Z=2.340085\)になる確率（=\(p\)値）は0.01927935
• → 有意水準\(\alpha=0.05\)より低い
• → 帰無仮説を棄却し、対立仮説を採用する

標本の性質や検定の目的によって使う検定統計量は異なる

• \(Z\)統計量：\(\frac{\bar{X} - \mu}{\sqrt{\sigma^2/n}}\)…平均値の差を標準誤差\({\sqrt{\sigma^2/n}}\)で標準化

  • 母集団が平均と分散が既知の正規分布であるとき
  • 帰無仮説のもとでは\(Z\)は標準正規分布に従う

• \(t\)統計量：\(\frac{\bar{X} - \mu}{\sqrt{\hat{\sigma}^2/n}}\)…母分散\(\sigma^2\)の代わりに不偏分散\(\hat{\sigma}^2\)を使う

  • 母集団が平均が既知、分散が未知の正規分布であるとき
  • 帰無仮説のもとでは\(t\)は自由度\(df = n-1\)の\(t\)分布に従う

• 他にも\(F\)統計量やカイ二乗統計量など
  

• 不適切な統計量（=方法）を用いると誤った\(p\)値が得られる

統計解析法にはそれぞれ 守るべき前提条件がある

• データに合った\(p\)値を求めるために必要
• 適切な方法を使わないと間違った\(p\)値が得られてしまう

有意差を求めて 間違った解析法を使うのは不正／不誠実

• 間違った方法によって得られた\(p < 0.05\)に意味はない

自分のデータに合った解析を選べるようになろう！

• より良い実験方法を模索するのと一緒で、より良い解析法を常に検討していく姿勢が大切！

\(p\)>0.05 →「差がない」と結論するのは間違い

• \(p\)>0.05 → 帰無仮説が棄却できない ≠ 差がない
• 仮説検定の立場からは どちらとも言えない
• 態度を保留するのが正解

\(p\)値が小さい →「差が大きい」と結論するのも間違い

• サンプルサイズが大きいと精度が上がる →僅かな差でも有意差が検出される
• \(p\)値が小さいからといって差が大きいとは限らない
• 差の大きさを示すのは効果量（effect size）
  • easystatsパッケージ群が便利かも（日本語解説）

penguinsで練習しよう

【復習】データを読み込んでざっくり眺める

library(palmerpenguins)               #penguinsを読み込み
library(tidyverse)                    #tidyverseを読み込み
head(penguins)                        #penguinsの最初を見せて！

## # A tibble: 6 × 8
##   species island    bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
##   <fct>   <fct>              <dbl>         <dbl>             <int>       <int>
## 1 Adelie  Torgersen           39.1          18.7               181        3750
## 2 Adelie  Torgersen           39.5          17.4               186        3800
## 3 Adelie  Torgersen           40.3          18                 195        3250
## 4 Adelie  Torgersen           NA            NA                  NA          NA
## 5 Adelie  Torgersen           36.7          19.3               193        3450
## 6 Adelie  Torgersen           39.3          20.6               190        3650
## # ℹ 2 more variables: sex <fct>, year <int>

• アデリーペンギンの体重に性差があるか？

アデリーペンギンの体重に性差があるか？

統計解析をする前に必ずデータをプロットすること！！

[練習問題]

• まずはデータをプロットしよう
  • 元データはpenguins
  • dplyrのfilterでspecies == “Adelie”だけ抽出
  • ggplot2でプロット
    • X軸…sex
    • Y軸…body_mass_g
      Return to HOME

アデリーペンギンの体重に性差があるか？

penguins_adelie = penguins %>%
  filter(species == "Adelie")
head(penguins_adelie)

## # A tibble: 6 × 8
##   species island    bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
##   <fct>   <fct>              <dbl>         <dbl>             <int>       <int>
## 1 Adelie  Torgersen           39.1          18.7               181        3750
## 2 Adelie  Torgersen           39.5          17.4               186        3800
## 3 Adelie  Torgersen           40.3          18                 195        3250
## 4 Adelie  Torgersen           NA            NA                  NA          NA
## 5 Adelie  Torgersen           36.7          19.3               193        3450
## 6 Adelie  Torgersen           39.3          20.6               190        3650
## # ℹ 2 more variables: sex <fct>, year <int>

アデリーペンギンの体重に性差があるか？

penguins_adelie = penguins %>%
  filter(species == "Adelie", !is.na(sex)) #NAを除く
head(penguins_adelie)

## # A tibble: 6 × 8
##   species island    bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
##   <fct>   <fct>              <dbl>         <dbl>             <int>       <int>
## 1 Adelie  Torgersen           39.1          18.7               181        3750
## 2 Adelie  Torgersen           39.5          17.4               186        3800
## 3 Adelie  Torgersen           40.3          18                 195        3250
## 4 Adelie  Torgersen           36.7          19.3               193        3450
## 5 Adelie  Torgersen           39.3          20.6               190        3650
## 6 Adelie  Torgersen           38.9          17.8               181        3625
## # ℹ 2 more variables: sex <fct>, year <int>

アデリーペンギンの体重に性差があるか？

gp = ggplot(data = penguins_adelie) +
  geom_point(aes(x = sex, y = body_mass_g), 
             position = position_jitter(width = 0.05))
gp

• いい感じ！さて…

正規分布とは

• 左右対象の釣鐘型の分布
• 現実世界でしばしば見られる
  • ある一つの数値を目標とした作業で生じる偶然的な誤差
  • 性別と年齢を固定したときの身長の分布平均と分散
• 平均と分散（または標準偏差）がわかればどのような分布になるか一意に決まる（めちゃ便利）

• この場合、オスのほうがメスより分散が大きい
• 実際はFテストやBartlettテストの結果で判断する

仮説検定の進め方

1. データをプロットして眺める
2. 各標本に正規性があるか調べる
   • Q-Q plot
   • Shapiro-Wilkテスト
3. 等分散性はあるか調べる
   • Fテスト
4. 2, 3の結果を受けて検定方法を決め、検定する

1. 各標本に正規性があるか（Q-Q plot）

gp = ggplot(data = penguins_adelie) +
       stat_qq(aes(sample = body_mass_g, color = sex)) + 
       facet_wrap( ~ sex, scales = "free")
gp

• 正規Q-Q plotで直線上に乗っているような気がする

1. 各標本に正規性があるか

• Shapiro-Wilkテスト
  • shapiro.test(【調べたい標本】)
    • 【調べたい標本】はベクトルc(a, b, c…)で入れる
    • 【data.frame】$【列名】でも指定できる
  • 帰無仮説：正規性がある
  • p<0.05で正規性がないことを示す
  • （正規性があることは言えない）
  • １標本の値を入れる（２標本を混ぜない）
    • 平均の異なる正規分布を２つ混ぜたら正規分布ではなくなる

1. 各標本に正規性があるか（Shapiro-Wilkテスト）

sw_test_results = penguins_adelie %>%   #penguins_adelieを使う
  group_by(sex) %>%                     #sexでグループ
  summarise(statistic = shapiro.test(body_mass_g)$statistic, 
        #body_mass_gにshapiro.testを行い、そのstatistic値をstatisticという列に
            p.value = shapiro.test(body_mass_g)$p.value,    #同様
            .groups     = "drop") #グループを解除
head(sw_test_results)

## # A tibble: 2 × 3
##   sex    statistic p.value
##   <fct>      <dbl>   <dbl>
## 1 female     0.977   0.199
## 2 male       0.983   0.416

• 正規性がないとは言えない→「正規性あり」と判断

2. 等分散性があるか

• Fテスト
  • ベクトルとして
    • var.test(【標本A】, 【標本B】)
      • 標本はベクトルc(a, b, c…)で入れる
  • データフレームとして
    • var.test(data = 【data.frame】, 【調べたい値の列名】~ 【グループ分けの列名】)
  • 帰無仮説：等分散性がある
  • p<0.05で等分散性がないことを示す

2. 等分散性があるか（Fテスト）

var.test(data = penguins_adelie, body_mass_g ~ sex)

## 
##  F test to compare two variances
## 
## data:  body_mass_g by sex
## F = 0.60331, num df = 72, denom df = 72, p-value = 0.03356
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.3787186 0.9611059
## sample estimates:
## ratio of variances 
##          0.6033147

• p<0.05なので「等分散性なし」と判断

アデリーペンギンの体重に性差があるか？

• 正規性はあり
• 等分散性はなし
• Studentのt検定（正規性と等分散性が前提条件）は使えない
• → Welchのt検定（正規性のみが前提条件）を使う
  • 選び方2-1を参照

正規性あり等分散性なしの２群の検定

• Welchのt検定
  • ベクトルとして
    • t.test(【標本A】, 【標本B】)
      
  • データフレームとして
    • t.test(data = 【data.frame】, 【調べたい値の列名】~ 【グループ分けの列名】)
      
  • 帰無仮説：平均が等しい
    
  • p<0.05で平均が異なる母集団に由来することを示す

アデリーペンギンの体重に性差があるか？(Welchのt検定)

t.test(data = penguins_adelie, body_mass_g ~ sex)

## 
##  Welch Two Sample t-test
## 
## data:  body_mass_g by sex
## t = -13.126, df = 135.69, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group female and group male is not equal to 0
## 95 percent confidence interval:
##  -776.3012 -573.0139
## sample estimates:
## mean in group female   mean in group male 
##             3368.836             4043.493

• 帰無仮説「平均が等しい」は棄却→有意差あり

• easystatsパッケージを使う
  • パッケージをインストールしてから…
  • report(【統計解析】)
    • このままだと文章型で出力される
  • 後ろに%>% as.data.frameとすると表型になる
  • さらに後ろに%>% write_csvで書き出せる

var.equal = の設定で「等分散性あり」の検定に変えられる

• 等分散性なし→Welchのt検定
  • t.test(data, Y ~ X, var.equal = FALSE)…デフォルト
    • 何も設定しなかったらこっち
• 等分散性あり→Studentのt検定
  • t.test(data, Y ~ X, var.equal = TRUE)
  • 等分散性が未知ならF検定なしでWelchのt検定をすべきという考え方がメジャーになりつつある

paired = TRUEの設定で「対応あり」に変えられる

• t.test(data, Y ~ X, pared = TRUE)
• 式(y ~ x)での設定ができないバグ？（R4.4.1, 2024/7/3）

くちばしの形は雌雄で異なるか？

• 細長さは(bill_length_mm)÷(bill_depth_mm)で指標が作れそう
• こういう割り算データは正規性がないと思ったほうがいい
• （平均値や中央値が0や1に近いときの分布を考えてみよう）

くちばしの細長さは雌雄で異なるか？

[練習問題]

• アデリーペンギンのくちばしの細長さを計算しよう
• 計算できたらプロットしてみよう

くちばしの細長さは雌雄で異なるか？

[練習問題]

penguins_adelie = penguins %>%
  filter(species == "Adelie", !is.na(sex)) %>% #NAを除く
  mutate(ratio = bill_length_mm/bill_depth_mm)
head(penguins_adelie)

## # A tibble: 6 × 9
##   species island    bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
##   <fct>   <fct>              <dbl>         <dbl>             <int>       <int>
## 1 Adelie  Torgersen           39.1          18.7               181        3750
## 2 Adelie  Torgersen           39.5          17.4               186        3800
## 3 Adelie  Torgersen           40.3          18                 195        3250
## 4 Adelie  Torgersen           36.7          19.3               193        3450
## 5 Adelie  Torgersen           39.3          20.6               190        3650
## 6 Adelie  Torgersen           38.9          17.8               181        3625
## # ℹ 3 more variables: sex <fct>, year <int>, ratio <dbl>

くちばしの細長さは雌雄で異なるか？

[練習問題]

gp = ggplot(data = penguins_adelie) +
  geom_point(aes(x = sex, y = ratio), 
             position = position_jitter(width = 0.05))
gp

くちばしの細長さは雌雄で異なるか？

1. 各標本に正規性はあるか（Q-Q plot）

gp = ggplot(data = penguins_adelie) +
       stat_qq(aes(sample = ratio, color = sex)) + 
       facet_wrap( ~ sex, scales = "free")
gp

くちばしの細長さは雌雄で異なるか？

1. 各標本に正規性はあるか（Shapiro-Wilkテスト）

sw_test_results = penguins_adelie %>%   #penguins_adelieを使う
  group_by(sex) %>%                     #sexでグループ
  summarise(statistic = shapiro.test(ratio)$statistic, 
            p.value = shapiro.test(ratio)$p.value, 
            .groups     = "drop") 
head(sw_test_results)

## # A tibble: 2 × 3
##   sex    statistic p.value
##   <fct>      <dbl>   <dbl>
## 1 female     0.982   0.397
## 2 male       0.984   0.480

くちばしの細長さは雌雄で異なるか？

2. 等分散性があるか（Fテスト）

var.test(data = penguins_adelie, ratio ~ sex)

## 
##  F test to compare two variances
## 
## data:  ratio by sex
## F = 0.78018, num df = 72, denom df = 72, p-value = 0.2946
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.4897412 1.2428573
## sample estimates:
## ratio of variances 
##          0.7801785

くちばしの細長さは雌雄で異なるか？

• 正規性はなし
• 等分散性はあり
• 対応はない
• → Wilcoxonの順位和検定（U検定と同じ）を行う
  • 選び方２−１（２群の検定）を参照

正規性なし等分散性ありの２群の比較

• Wilcoxonの順位和検定
  • デフォルトのwilcox.test()ではタイ（同値）があった場合に正確なp値が出ないのでcoinパッケージを使う
    • ベクトルとして
      • wilcox_test(【標本A】, 【標本B】, distribution = “exact”)
        
    • データフレームとして
      • wilcox_test(data = 【data.frame】, 【調べたい値の列名】~ 【グループ分けの列名】, distribution = “exact”)
        
  • 帰無仮説：２群間に差がない
    
  • p<0.05で異なる母集団に由来することを示す

くちばしの細長さは雌雄で異なるか？

まずインストールする（１行目の#を消して走らせる）

#install.packages("coin")
library(coin)             # 読み込み
wilcox_test(data = penguins_adelie, ratio ~ sex, distribution = "exact")

## 
##  Exact Wilcoxon-Mann-Whitney Test
## 
## data:  ratio by sex (female, male)
## Z = -0.29941, p-value = 0.7661
## alternative hypothesis: true mu is not equal to 0

• 帰無仮説「２群間に差がない」は棄却されない→有意差なし

対応なし→Wilcoxonの順位和検定

• wilcox_test(data, Y ~ X, distribution = “exact”)

対応あり→Wilcoxonの符号順位検定

• wilcoxsign_test(data, Y ~ X, distribution = “exact”)

• Brunner-Munzel検定
  • brunnermunzelパッケージを使う
  • ベクトルとして
    • brunnermunzel.test(【標本A】, 【標本B】)
      
  • データフレームとして
    • brunnermunzel.test(data = 【data.frame】, 【調べたい値の列名】~ 【グループ分けの列名】)
      
  • 帰無仮説：２群間に差がない
  • p<0.05で異なる母集団に由来することを示す
    
  • 「サンプル数が少ない」というエラーが出たらbrunnermunzel.permutation.test

さっきのデータで使ってみよう

まずインストールする（１行目の#を消して走らせる）

#install.packages("brunnermunzel")
library(brunnermunzel)             # 読み込み
brunnermunzel.test(data = penguins_adelie, ratio ~ sex)

## 
##  Brunner-Munzel Test
## 
## data:  ratio by sex
## Brunner-Munzel Test Statistic = 0.29701, df = 143.58, p-value = 0.7669
## 95 percent confidence interval:
##  0.4188189 0.6098919
## sample estimates:
## P(X<Y)+.5*P(X=Y) 
##        0.5143554

1. まずはプロット
2. 正規性, 等分散性, 対応あり/なしを確認
   • Q-Q plot, Shapiro-Wilkテスト
   • F検定
3. 検定方法を決める
   • Studentのt検定（対応あり／なし）
   • Welchのt検定（対応あり／なし）
   • Wilcoxonの順位和検定
   • Wilcoxonの符号順位検定
   • Brunner-Munzel検定（対応あり／なし）

２群の比較を単に繰り返すのはダメ

• 仮説検定の考え方を思い出そう

1. 同じ母集団由来であった場合に今回得られた標本の差が生じる確率が低い（<5%）→帰無仮説を棄却する
2. その結果、対立仮説を採用する

どこかで有意差が出る確率はどんどん高まる

• 100回検定すれば5回は偶然有意差が出る 阪大医学・老年・腎臓内科学HP

有意差があった群のみ検定したことにする？

→もちろん 不正！！

多重比較法を使うか、有意水準を補正する

• 今回は多重比較を紹介します
• 有意水準の補正の方法はこちらなど

【復習】読み込んでデータをざっくり眺める

library(palmerpenguins)               #penguinsを読み込み
library(tidyverse)                    #tidyverseを読み込み
head(penguins)                        #penguinsの最初を見せて！

## # A tibble: 6 × 8
##   species island    bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
##   <fct>   <fct>              <dbl>         <dbl>             <int>       <int>
## 1 Adelie  Torgersen           39.1          18.7               181        3750
## 2 Adelie  Torgersen           39.5          17.4               186        3800
## 3 Adelie  Torgersen           40.3          18                 195        3250
## 4 Adelie  Torgersen           NA            NA                  NA          NA
## 5 Adelie  Torgersen           36.7          19.3               193        3450
## 6 Adelie  Torgersen           39.3          20.6               190        3650
## # ℹ 2 more variables: sex <fct>, year <int>

• ペンギンの体重は種ごとに異なるか？

ペンギンの体重は種ごとに異なるか？

• 体重に性差があるのでメスに限定して調べよう
• まずはデータをプロットしよう

penguins_female = penguins %>%
  filter(sex == "female")
head(penguins_female)

## # A tibble: 6 × 8
##   species island    bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
##   <fct>   <fct>              <dbl>         <dbl>             <int>       <int>
## 1 Adelie  Torgersen           39.5          17.4               186        3800
## 2 Adelie  Torgersen           40.3          18                 195        3250
## 3 Adelie  Torgersen           36.7          19.3               193        3450
## 4 Adelie  Torgersen           38.9          17.8               181        3625
## 5 Adelie  Torgersen           41.1          17.6               182        3200
## 6 Adelie  Torgersen           36.6          17.8               185        3700
## # ℹ 2 more variables: sex <fct>, year <int>

ペンギンの体重は種ごとに異なるか？

• まずはデータをプロットしよう

gp = ggplot(data = penguins_female) +
  geom_point(aes(x = species, y = body_mass_g), 
             position = position_jitter(width = 0.05))
gp

• 差がありそうな雰囲気…！

３群以上の比較を行うときも進め方は同じ

1. データをプロットして眺める
2. 各標本に正規性があるか調べる
   • Q-Q plot
   • Shapiro-Wilkテスト
3. 等分散性はあるか調べる
   • Bartlettテスト ←Fテストの多群版
4. 2, 3の結果を受けて検定方法を決め、検定する

1. 各標本に正規性はあるか（Q-Q plot）

gp = ggplot(data = penguins_female) +
  stat_qq(aes(sample = body_mass_g, color = species)) + 
  facet_wrap( ~ species, scales = "free")
gp

• 正規Q-Q plotで直線上に乗っているような気がする

1. 各標本に正規性はあるか（Shapiro-Wilkテスト）

sw_test_results = penguins_female %>%   #penguins_adelieを使う
  group_by(species) %>%                     #speciesでグループ
  summarise(statistic = shapiro.test(body_mass_g)$statistic, 
        #body_mass_gにshapiro.testを行い、そのstatistic値をstatisticという列に
            p.value = shapiro.test(body_mass_g)$p.value,    #同様
            .groups     = "drop") #グループを解除
head(sw_test_results)

## # A tibble: 3 × 3
##   species   statistic p.value
##   <fct>         <dbl>   <dbl>
## 1 Adelie        0.977   0.199
## 2 Chinstrap     0.963   0.306
## 3 Gentoo        0.981   0.511

• 「正規性あり」と判断

2. 等分散性があるか

• Bartlettテスト
  • ベクトルとして
    • bartlett.test(x = list(A, B, C))
  • データフレームとして
    • bartlett.test(data = 【data.frame】, 【調べたい値の列名】~ 【グループ分けの列名】)
  • 帰無仮説：等分散性がある
  • p<0.05で等分散性がないことを示す

2. 等分散性があるか（Bartlettテスト）

bartlett.test(data = penguins_female, body_mass_g ~ species)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  body_mass_g by species
## Bartlett's K-squared = 0.19828, df = 2, p-value = 0.9056

• 「等分散性あり」と判断

ペンギンの体重は種ごとに異なるか？

• 正規性あり
• 等分散性あり
• ３群以上
• → 一元分散分析（1-way ANOVA）を使う
  • 選び方２−１（多重比較）を参照

• Type I：要因が１つ（一元分散分析）のときに使う
  • aov()やanova()で行える
• Type II, III：要因が２つ以上のときに使う
  • aov()やanova()ではTypeを選択できない（Type Iのみ）
  • carパッケージのAnova()関数を用いる
  • Type IIIを使うときにはoptions(contrasts = c(“contr.sum”, “contr.sum”))でオプションを設定しておく
  • 詳しくはこちらなど
• → 今回はType Iを使う

正規性あり等分散性ありの多群比較

• 一元分散分析（ANOVA）
  • summary(aov(data = 【data.frame】, 【調べたい値の列名】~ 【グループ分けの列名】))
    
  • 帰無仮説：全ての群が同一の母集団に由来する（=全ての群に差がない）
  • p<0.05で異なる母集団に由来することを示す

さて、ペンギンの体重は種ごとに異なるのか？

summary(aov(data = penguins_female, body_mass_g ~ species))

##              Df   Sum Sq  Mean Sq F value Pr(>F)    
## species       2 60350016 30175008   393.2 <2e-16 ***
## Residuals   162 12430757    76733                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• 有意差はある
• 帰無仮説（群間に差がない）は棄却された
• どのペアに差があるのか？
  • post-hocテストを行う

正規性あり等分散性ありの多群比較のpost-hocテスト

• Tukeyの多重比較
  • TukeyHSD(aov(data = 【data.frame】, 【調べたい値の列名】~ 【グループ分けの列名】))
    
  • 帰無仮説：２群間に差がない
    
  • p<0.05で異なる母集団に由来することを示す

ペンギンの体重は種ごとに異なるか？

TukeyHSD(aov(data = penguins_female, body_mass_g ~ species))

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = body_mass_g ~ species, data = penguins_female)
## 
## $species
##                       diff        lwr       upr     p adj
## Chinstrap-Adelie  158.3703   22.32078  294.4197 0.0179471
## Gentoo-Adelie    1310.9058 1195.64908 1426.1624 0.0000000
## Gentoo-Chinstrap 1152.5355 1011.00620 1294.0648 0.0000000

ペンギンの体重は種ごとに異なるか？（ART ANOVA）

• 分散分析を実行

anova(model)

## Analysis of Variance of Aligned Rank Transformed Data
## 
## Table Type: Anova Table (Type III tests) 
## Model: No Repeated Measures (lm)
## Response: art(body_mass_g)
## 
##           Df Df.res F value     Pr(>F)    
## 1 species  2    162  191.44 < 2.22e-16 ***
## ---
## Signif. codes:   0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• 有意なのでpost-hocを行う

ペンギンの体重は種ごとに異なるか？（ART ANOVA）

• ANOVAが有意だったらposthocを実行

art.con(model, "species")

##  contrast           estimate   SE  df t.ratio p.value
##  Adelie - Chinstrap    -17.6 5.44 162  -3.240  0.0041
##  Adelie - Gentoo       -88.0 4.61 162 -19.104  <.0001
##  Chinstrap - Gentoo    -70.4 5.66 162 -12.443  <.0001
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

• できた！
  Return to HOME

1. まずはプロット
2. 正規性, 等分散性, 対応あり/なしを確認
   • Q-Q plot, Shapiro-Wilkテスト
   • Bartlettテスト
3. 検定方法を決める
   • ANOVA → Tukey
   • ART ANOVA → Tukey
   • 選び方２−１（多重比較）を参照

２要因の場合は…

• 2-way ANOVA（Type IIかIII）や2-way ART-ANOVAを使う

ペンギンの性比は種ごとに異なるか？

penguins_sexratio_sum = penguins %>% 
  group_by(species, sex) %>%            # speciesとsexごとにグループ分け
  summarise(length = length(species),   # N数を数える
            .groups     = "drop") %>%   # groupを解除
  filter(!is.na(sex)) %>%               # sex == NAを除去
  print()                               # 表示

## # A tibble: 6 × 3
##   species   sex    length
##   <fct>     <fct>   <int>
## 1 Adelie    female     73
## 2 Adelie    male       73
## 3 Chinstrap female     34
## 4 Chinstrap male       34
## 5 Gentoo    female     58
## 6 Gentoo    male       61

ペンギンの性比は種ごとに異なるか？

gp = ggplot(data = penguins_sexratio_sum, aes(x = species, y = length)) +
  geom_bar(aes(fill = sex), stat = "identity", position = "fill") +
  scale_y_continuous(labels = scales::percent) +
  coord_flip()                    # X軸とY軸を逆にする
gp

• まぁ違わなそうですけどね…

• カイ２乗検定
  • ベクトルとして
    • chisq.test(【ベクトルA】,【ベクトルB】)
      
    • A…【data.frame】$【調べたい値の列名】
    • B…【data.frame】$【グループ分けの列名】
  • 行列として
    • chisq.test(【行列】)
  • 帰無仮説：２群の比率に差がない
    
  • サンプル数が小さい場合はフィッシャーの正確確率検定
  • ３群以上ある場合は補正を行う
    • rcompanionパッケージのpairwiseNominalIndependence関数を使う

• サンプル数が小さい時
• フィッシャーの正確確率検定
  • fisher.test(【ベクトルA】,【ベクトルB】)
    • A…【data.frame】$【調べたい値の列名】
    • B…【data.frame】$【グループ分けの列名】
  • 帰無仮説：２群の比率に差がない

ペンギンの性比は種ごとに異なるか？

penguins_sexratio = penguins %>%
  filter(!is.na(sex))              # NAを除く
chisq.test(penguins_sexratio$sex, penguins_sexratio$species)

## 
##  Pearson's Chi-squared test
## 
## data:  penguins_sexratio$sex and penguins_sexratio$species
## X-squared = 0.048607, df = 2, p-value = 0.976

ペンギンの性比は種ごとに異なるか？

クロス集計表を使う場合…

• penguins_sexratio_sumを横長にしてクロス集計表を作る

penguins_sexratio_sum_cross = penguins_sexratio_sum %>%
  pivot_wider(names_from = sex, values_from = length)%>% #横長に
  tibble::column_to_rownames(var = "species") %>% #speciesを行名に
  as.matrix() #行列に変換
penguins_sexratio_sum_cross

##           female male
## Adelie        73   73
## Chinstrap     34   34
## Gentoo        58   61

ペンギンの性比は種ごとに異なるか？

クロス集計表を使う場合…

• penguins_sexratio_sumに対して検定する

chisq.test(penguins_sexratio_sum_cross)

## 
##  Pearson's Chi-squared test
## 
## data:  penguins_sexratio_sum_cross
## X-squared = 0.048607, df = 2, p-value = 0.976

• 3群以上の場合はこちら

1. まずはプロット
2. 検定方法を決める
   • カイ２乗検定
   • Fisherの正確確率検定
   • 選び方２−１（２要因）を参照
   • 3群以上の場合は有意水準の補正を行う

1. まず、自分が何を示したいのかを整理する

   • 「A群とB群は異なる」とか「Aが増えるとBも増える」とか

2. 自分のデータの性質を知り、適切な統計手法を探す

3. 統計手法を実践する方法を探す

   • まず仮説検定かモデリングかを決める
   • それぞれのなかから適切な方法を探す
     • 仮説検定
     • モデリング（一般化線形モデル）←やってみよう！

統計モデリング

• 応答変数（Y）を説明変数（X）の式にあてはめる（= 関係性を説明する）
• 質的変数にも使える（ANOVAなどと同等の解析が可能）
• 傾きや切片に関しては仮説検定（帰無仮説：推定量が0）を行う

• Generalized linear model
• 残差（観察値と推定値の差）を特定の確率分布とした線形モデル

• 直線以外にも指数関数orシグモイド的な回帰などができる
• ただしあまりぐにゃぐにゃせず、Xが増えるに従ってYは増えるか減るかどちらか

• 以下の３つの要素から成る
  • リンク関数…XとYが大体どんな関係かの回帰式
    • 例：直線、指数、シグモイド…
  • 線形予測子…回帰式のパラメータ（Y ~ \(\beta_0\) + \(\beta_1\)X)
  • 誤差構造…応答変数の予測値の周りのばらつきの分布

３つの要素を設定するだけ

• glm(線形予測子, family = 誤差構造(link = リンク関数), data = 【データフレーム】)
  • 例：glm(y ~ x1 + x2, family = gaussian(link = identity), data = penguins)

線形予測子…data.frameの列名を指定

• 単回帰【応答変数】 ~ 【説明変数】
  • 例）bill_length_mm ~ bill_depth_mm
    • 嘴の長さ ~ \(\beta_0\) + (\(\beta_1\)×嘴の幅)
• 重回帰（交互作用なし）【応答】 ~ 【説明１】+【説明２】
  • 例）bill_length_mm ~ bill_depth_mm + species
    • 嘴の長さ ~ \(\beta_0\) + (\(\beta_1\)×嘴の幅) + (\(\beta_2\)×種)
• 重回帰（交互作用あり）【応答】 ~ 【説明１】*【説明２】
  • 解釈が難しい場合は避ける
  • 例）bill_length_mm ~ bill_depth_mm * species
    • 嘴の長さ ~ \(\beta_0\) + (\(\beta_1\)×嘴の幅) + (\(\beta_2\)×種) + (\(\beta_3\)×嘴の幅と種の交互作用)

誤差構造…下記から選び、（）内の分布名を設定

• 正規分布（gaussian）…確率変数の和、平均値
  • \(-∞\leqq Y < ∞\)（下限なし, 上限なし）, 連続値
• 二項分布（binomial）…成功率p、試行回数nのうちの成功回数
  • \(0\leqq Y \leqq max\)（下限あり, 上限あり）, 離散値
• ポアソン分布（poisson）…単位時間あたり平均λ回起こる事象の発生回数
  • \(0\leqq Y \leqq ∞\), （下限あり, 上限なし）, 離散値
• ガンマ分布（Gamma）…ポアソン過程でk回起こるまでの待ち時間
  • \(0\leqq Y \leqq ∞\), （下限あり, 上限なし）, 連続値
• 他にinverse.gaussian, quasi, quasibinomial, quasipoisson

リンク関数…下記から選び、（）内の関数名を設定

• そのまま（identity）…直線回帰
  • \(Y = \beta_0 + \beta_1 X\)
  • 誤差構造：正規分布
• 対数（log）…指数回帰
  • \(log(Y) = \beta_0 + \beta_1 X \therefore Y = e^{\beta_0 + \beta_1 X}\)
  • 誤差構造：ポアソン分布など
• ロジット（logit）…ロジスティック回帰（シグモイド型）
  • \(logit(Y) = \beta_0 + \beta_1 X \therefore Y = \frac{1}{1+e^-(\beta_0 + \beta_1 X)}\)
  • 誤差構造：二項分布
• 他にprobit, cloglog, sqrt, 1/mu^2, inverseが使える

penguinsデータで練習しよう

library(palmerpenguins)               #penguinsを読み込み
library(tidyverse)                    #tidyverseを読み込み
head(penguins)                        #penguinsの最初を見せて！

## # A tibble: 6 × 8
##   species island    bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
##   <fct>   <fct>              <dbl>         <dbl>             <int>       <int>
## 1 Adelie  Torgersen           39.1          18.7               181        3750
## 2 Adelie  Torgersen           39.5          17.4               186        3800
## 3 Adelie  Torgersen           40.3          18                 195        3250
## 4 Adelie  Torgersen           NA            NA                  NA          NA
## 5 Adelie  Torgersen           36.7          19.3               193        3450
## 6 Adelie  Torgersen           39.3          20.6               190        3650
## # ℹ 2 more variables: sex <fct>, year <int>

• 体重の大きいペンギンは翼長も長い？

体重の大きいアデリーは翼長も長い？（まずはプロット）

gp = ggplot(data = penguins_adelie) +
  geom_point(aes(x = body_mass_g, y = flipper_length_mm)) 
gp

• 体重が大きければ翼長も長いような気がする
• 上昇は直線的な感じだし、翼長は正規分布でよさそう

結果の見方

• Coefficients: 回帰係数
  • Estimates: 推定値
  • Pr(>|t|): 帰無仮説（切片や係数が0である）に対するp値。
    • (Intercept)：切片\(\beta_0\)
    • body_mass_g：body_mass_gの係数\(\beta_1\)

今回の結果は…

• \(\beta_0\): 1.652e+02 = 165.2 (Pr(>|t|)< 2e-16)
• \(\beta_1\): 6.677e-03 = 0.006677 (Pr(>|t|)=3.4e-09)
• →翼長への体サイズの効果は有意である
• →\(Y = 0.006677X + 165.2\)に回帰できる

• Yの観察値は常に正の値なのに、回帰は負に入っている
• 縦軸は整数。Xが増えるとYのばらつきは増える。

• Yの観察値は常に正の値なのに、回帰は負に入っている
• 縦軸は整数。Xが増えるとYのばらつきは増える。
• データに合った統計モデルを使うとマシ。

• Yがポアソン分布に従う
  • 0以上の離散値（整数とか）で、指数関数的に増える
• 設定：family = poisson(link = “log”)
  • 誤差構造：poisson
  • リンク関数：\(logY = \beta_0 + \beta_1 X\)（\(\therefore Y = e^{(\beta_0 + \beta_1 X)}\)）

    Coded by Heavy Watal

• 平均\(\lambda\)で単位時間（空間）あたりに発生する事象の回数。
• e.g., １時間あたりの交通事故数、メッシュ区画内の生物個体数

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

\(Prob(X = k \mid \lambda) = \frac {\lambda^k e^{-\lambda}} {k!}\)

植物のサイズが種子数に影響するか？（まずはプロット）

gp = ggplot(data = seed) +
  geom_point(aes(x =plant_size , y = no.seed)) 
gp

• なんとなくサイズが大きいほうが種子数が多いような…？

植物のサイズが種子数に影響するか？（GLM）

• \(\beta_0\): 1.29172 (Pr(>|t|)=0.000383)

• \(\beta_1\): 0.07566 (Pr(>|t|)=0.033580)
  

• 種子数への植物サイズの効果は有意である

• \(Y = e^{(1.29172 + 0.07566 X)}\)に回帰できる

• Yが二項分布に従う
  • 0以上で上限のある離散値でシグモイド型に変化する
• 書き方：family = binomial(link = “logit”)
  • 誤差構造：binomial
  • リンク関数：\(logitY= \beta_0 + \beta_1 X\)（\(\therefore Y = \frac{1}{1+e^{-(\beta_0 + \beta_1 X)}}\)）

    Coded by Heavy Watal

• 確率\(p\)で当たるクジを\(n\)回引いてX回当たる確率。平均は\(np\)。

\({Prob}(X = k \mid n,~p) = \binom n k p^k (1 - p)^{n - k}\)

線形予測子の書き方は２種類がある

• 誤差構造とリンク関数は…
  • family = binomial(link = “logit”)
• 線形予測子（応答変数Yの指定の仕方）は…
  • A：成功数（Y1） vs 失敗数（Y2）などの整数の対データ
    • glm(cbind(Y1, Y2) ~ X)
      • 割り算した値は入れられない
  • B：成功（1） or 失敗（0）などの２値データ（Y）
    • glm(Y ~ X)
      • 0/1以外でも2値であればOK

植物サイズは種子生存率に影響するか？（GLM）

• \(\beta_0\): -13.7785 (Pr(>|t|)<2e-16)
• \(\beta_1\): 1.4626 (Pr(>|t|)<2e-16)
  • 種子生存率への植物サイズの効果は有意である
• \(Y = \frac{1}{1+e^{-(-13.7785 + 1.4626 X)}}\)に回帰できる

• 処理（treatment）の効果も考慮に入れてみよう！
• →重回帰（説明変数が２つ：植物サイズ + 処理）

重回帰=複数の説明変数で応答変数を説明したいときに使う

• 線形予測子
  • 【応答】 ~ 【説明１】+【説明２】…交互作用がないとき
  • 例）bill_length_mm ~ bill_depth_mm + species
    • 嘴の長さ ~ \(\beta_0\) + (\(\beta_1\)×嘴の幅) + (\(\beta_2\)×種)
    • 説明変数同士は独立である（相関関係などがない）必要あり
  • 【応答】 ~ 【説明１】*【説明２】…交互作用があるとき
    • 嘴の長さ ~ \(\beta_0\) + (\(\beta_1\)×嘴の幅) * (\(\beta_2\)×種)

植物サイズと処理は種子生存率に影響するか？（GLM）

• \(\beta_0\): -19.5361 (Pr(>|t|)<2e-16)
• \(\beta_1\): 1.9524 (Pr(>|t|)<2e-16)
  • plant sizeが1上がったときの効果
• \(\beta_2\): 2.0215 (Pr(>|t|)<2e-16)
  • treatmentがC→Tになったときの効果
• 植物サイズも処理も種子生存率への効果は有意である
• \(Y = \frac{1}{1+e^{-(-19.5361 + 1.9524 X_1 + 2.0215 X_2)}}\)に回帰できる

母親の体重は低体重出生に影響するか（まずプロット）

gp = ggplot(data = birthwt, aes(x = lwt , y = low)) +
  geom_point() 
gp

• 体重（lwt）が多いほうが低体重出生ではない（y = 0）傾向…？

母親の体重は低体重出生に影響するか（GLM）

• \(\beta_0\): 0.99831 (Pr(>|t|) = 0.2036)
• \(\beta_1\): -0.01406 (Pr(>|t|) = 0.0227)
• →母親の体重の、子の低体重出生への効果は有意である
• →母親の体重が重いと、子の低体重出生率は減少する
• \(Y = \frac{1}{1+e^{-(0.99831 -0.01406 X)}}\)に回帰できる

“All models are wrong, but some are useful”

• 唯一無二の完璧なモデルは存在しない
• 目的に応じて少しでも「マシな」モデルを選ぶ

1. メカニズム的に納得できるもの
   • ポアソン過程のカウントデータならポアソン分布
     • 過分散（overdispersion）があったら負の二項分布など
   • n回中k回のように割合的なカウントデータなら二項分布
2. データを可視化してみて、それっぽい形・性質のもの
   • 左右対称のひと山ならとりあえず正規分布
   • 負の値をとらないならガンマ分布
   • 直線的か、指数関数的か、頭打ちかなど

1. まずはプロット
2. メカニズムやプロットからモデルを選ぶ
   • リンク関数…XとYが大体どんな関係かの回帰式
     • 例：直線、指数、シグモイド…
   • 線形予測子…回帰式のパラメータ（Y ~ \(\beta_0\) + \(\beta_1\)X)
   • 誤差構造…応答変数の予測値の周りのばらつきの分布
     • 確率分布や回帰式についてもっと詳しく知りたい人は…
     • 統計モデリング概論（Heavy Watal）
3. GLMして結果を解釈する
4. どのモデルがいちばん”良い”かを客観的に評価するには…
   • モデル選択(AICなど)（今回は割愛）
5. さらに自由なモデリング（ベイズ統計など）

1. どんな統計解析ができるか考えてからデータを取ろう
   • どの統計解析を用いたら仮説を検証できるか？
     • 連続値？離散値？応答変数？説明変数？
     • 統計解析の前提条件もよく確認する
2. 統計解析に渡しやすい形でデータを取る
   • あるいは元データを変換する（tidyr）
3. まずはプロット！
   • 生データを見ることが一番大切
4. 統計解析する
5. 結果を出力し、解釈／考察する

• では自由課題をやってみよう！