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

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

3. データを整理する

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

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

実験生物学でよく使う基本的な統計解析

• 基礎知識と心構え
• 選び方
• 実践

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

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

この実習で扱わないこと

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

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

• 母集団の性質を知りたいが実測は不可能
• そこでランダムに抽出した標本から母集団の性質を推定する
  • 実験で得られた値など（ランダムでなければ意味がない）
• 集団にはばらつきがあるので毎回異なる値が得られる
• →ばらつきを考慮して母集団の性質を推定する

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

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

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

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

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

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

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

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

• 統計モデリング（重回帰分析）

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

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

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

• 例えばこんなとき…
• とある卒研生が「栄養豊富な餌で育ったハエは身体が大きくなる」ことを主張しようとしてデータを集めた。
• これまでのデータによると、通常の餌で育てたハエの体長は平均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\)値は検定統計量に基づいて求める

• 検定統計量とは…
  • 仮説検定で用いられる統計量
  • （めっちゃざっくり言うと、とある前提に基づいて標本の性質をまとめた仮説検定の基準になる値）
  • 標本から計算された値を「検定統計量の実現値」と呼ぶ
  • データに適した検定統計量を使う必要がある
    • \(Z\)統計量、\(t\)統計量、\(F\)統計量、カイ二乗統計量など
• 帰無仮説のもとでの検定統計量の確率分布を元に、「検定統計量の実現値が偶然得られる確率（=\(p\)値）」を計算する
• \(p\)値が有意水準（\(\alpha=0.05\)）より低ければ帰無仮説を棄却する

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

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

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

• とある卒研生（p.6）の例だと…母集団は正規分布に従い、平均（=3.0）と 標準偏差（=1.0）がわかっている
• 平均と標準偏差が既知の正規分布の場合、\(Z\)統計量が使える
  • \(Z = \frac{\bar{X} - \mu}{\sqrt{\sigma^2/n}}\)…平均値の差を標準誤差\({\sqrt{\sigma^2/n}}\)で標準化
    • \(\bar{X}\)…標本の平均 = 3.74
    • \(\mu\)…母集団の平均 = 3.0
    • \(\sigma\)…母集団の標準偏差 = 1.0
    • \(n\)…サンプルサイズ = 10
• 帰無仮説のもとでは\(Z\)は標準正規分布（平均0, 標準偏差１の正規分布）に従う
• 今回の場合、\(Z = \frac{3.74 - 3}{\sqrt{1^2/10}}\) = 2.340085

• 帰無仮説のもとで\(Z\)が従う標準正規分布では、\(Z\)が灰色の範囲（-1.96 < \(Z\) < 1.96）外の値をとる確率は5%（=有意水準\(\alpha\)）

• 今回の場合、\(Z=2.340085\) ←範囲外！
• 帰無仮説のもとで\(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\)値はサンプルサイズが大きいと有意になりやすい
• \(p\)値が小さいからといって差が大きいとは限らない
• 差の大きさを示すには効果量（effect size）
  • easystatsパッケージ群が便利かも（日本語解説）

守るべき前提条件によって方法は変わる

penguinsで練習しよう

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

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

## # A tibble: 6 x 8
##   species island bill_length_mm bill_depth_mm flipper_length_… body_mass_g sex  
##   <fct>   <fct>           <dbl>         <dbl>            <int>       <int> <fct>
## 1 Adelie  Torge…           39.1          18.7              181        3750 male 
## 2 Adelie  Torge…           39.5          17.4              186        3800 fema…
## 3 Adelie  Torge…           40.3          18                195        3250 fema…
## 4 Adelie  Torge…           NA            NA                 NA          NA <NA> 
## 5 Adelie  Torge…           36.7          19.3              193        3450 fema…
## 6 Adelie  Torge…           39.3          20.6              190        3650 male 
## # … with 1 more variable: year <int>

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

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

[練習問題]

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

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

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

## # A tibble: 6 x 8
##   species island bill_length_mm bill_depth_mm flipper_length_… body_mass_g sex  
##   <fct>   <fct>           <dbl>         <dbl>            <int>       <int> <fct>
## 1 Adelie  Torge…           39.1          18.7              181        3750 male 
## 2 Adelie  Torge…           39.5          17.4              186        3800 fema…
## 3 Adelie  Torge…           40.3          18                195        3250 fema…
## 4 Adelie  Torge…           NA            NA                 NA          NA <NA> 
## 5 Adelie  Torge…           36.7          19.3              193        3450 fema…
## 6 Adelie  Torge…           39.3          20.6              190        3650 male 
## # … with 1 more variable: year <int>

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

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

## Warning: Removed 1 rows containing missing values (geom_point).

• オスのほうが大きそう…！でもNAとかあるなぁ…

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

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

## # A tibble: 6 x 8
##   species island bill_length_mm bill_depth_mm flipper_length_… body_mass_g sex  
##   <fct>   <fct>           <dbl>         <dbl>            <int>       <int> <fct>
## 1 Adelie  Torge…           39.1          18.7              181        3750 male 
## 2 Adelie  Torge…           39.5          17.4              186        3800 fema…
## 3 Adelie  Torge…           40.3          18                195        3250 fema…
## 4 Adelie  Torge…           36.7          19.3              193        3450 fema…
## 5 Adelie  Torge…           39.3          20.6              190        3650 male 
## 6 Adelie  Torge…           38.9          17.8              181        3625 fema…
## # … with 1 more variable: 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 x 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なので「等分散性なし」と判断

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

• 正規性はあり
• 等分散性はなし
• 普通のt検定（正規性と等分散性が前提条件）は使えない
• → Welchのt検定（正規性のみが前提条件）を使う
  • フローチャートを参照

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

• 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 = の設定で対応ありの検定に変えられる

対応なし

• t.test(data, Y ~ X, paired = FALSE)…デフォルト

対応あり

• t.test(data, Y ~ X, paired = TRUE)
  • 同じ個体から複数回とったデータを比較したいとき

［練習問題］

• 対応なし、等分散性ありの検定をしてみよう
• 対応あり、等分散性なしの検定をしてみよう

［練習問題］

• 対応なし、等分散性ありの検定をしてみよう

t.test(data = penguins_adelie, body_mass_g ~ sex, 
       paired = FALSE, var.equal = TRUE) # paired = FALSEは省略可能

## 
##  Two Sample t-test
## 
## data:  body_mass_g by sex
## t = -13.126, df = 144, 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.2484 -573.0666
## sample estimates:
## mean in group female   mean in group male 
##             3368.836             4043.493

［練習問題］

• 対応あり、等分散性なしの検定をしてみよう

t.test(data = penguins_adelie, body_mass_g ~ sex, 
       paired = TRUE, var.equal = FALSE) # var.equal = FALSEは省略可能

## 
##  Paired t-test
## 
## data:  body_mass_g by sex
## t = -14.264, df = 72, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -768.944 -580.371
## sample estimates:
## mean of the differences 
##               -674.6575

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

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

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

[練習問題]

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

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

[練習問題]

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

## # A tibble: 6 x 9
##   species island bill_length_mm bill_depth_mm flipper_length_… body_mass_g sex  
##   <fct>   <fct>           <dbl>         <dbl>            <int>       <int> <fct>
## 1 Adelie  Torge…           39.1          18.7              181        3750 male 
## 2 Adelie  Torge…           39.5          17.4              186        3800 fema…
## 3 Adelie  Torge…           40.3          18                195        3250 fema…
## 4 Adelie  Torge…           36.7          19.3              193        3450 fema…
## 5 Adelie  Torge…           39.3          20.6              190        3650 male 
## 6 Adelie  Torge…           38.9          17.8              181        3625 fema…
## # … with 2 more variables: 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 x 3
##   sex    statistic p.value
##   <fct>      <dbl>   <dbl>
## 1 female     0.975 0.164  
## 2 male       0.954 0.00962

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

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

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

## 
##  F test to compare two variances
## 
## data:  ratio by sex
## F = 0.73443, num df = 72, denom df = 72, p-value = 0.1928
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.4610261 1.1699847
## sample estimates:
## ratio of variances 
##          0.7344341

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

• 正規性はなし
• 等分散性はあり
• 対応はない
• → 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検定

• 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.3901081 0.5811811
## sample estimates:
## P(X<Y)+.5*P(X=Y) 
##        0.4856446

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 x 8
##   species island bill_length_mm bill_depth_mm flipper_length_… body_mass_g sex  
##   <fct>   <fct>           <dbl>         <dbl>            <int>       <int> <fct>
## 1 Adelie  Torge…           39.1          18.7              181        3750 male 
## 2 Adelie  Torge…           39.5          17.4              186        3800 fema…
## 3 Adelie  Torge…           40.3          18                195        3250 fema…
## 4 Adelie  Torge…           NA            NA                 NA          NA <NA> 
## 5 Adelie  Torge…           36.7          19.3              193        3450 fema…
## 6 Adelie  Torge…           39.3          20.6              190        3650 male 
## # … with 1 more variable: year <int>

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

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

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

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

## # A tibble: 6 x 8
##   species island bill_length_mm bill_depth_mm flipper_length_… body_mass_g sex  
##   <fct>   <fct>           <dbl>         <dbl>            <int>       <int> <fct>
## 1 Adelie  Torge…           39.5          17.4              186        3800 fema…
## 2 Adelie  Torge…           40.3          18                195        3250 fema…
## 3 Adelie  Torge…           36.7          19.3              193        3450 fema…
## 4 Adelie  Torge…           38.9          17.8              181        3625 fema…
## 5 Adelie  Torge…           41.1          17.6              182        3200 fema…
## 6 Adelie  Torge…           36.6          17.8              185        3700 fema…
## # … with 1 more variable: 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 x 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

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

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

• 正規性あり
• 等分散性あり
• ３群以上
• → 一元分散分析（ANOVA）を使う

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

• 一元分散分析（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

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

gp = ggplot(data = penguins) +
  geom_point(aes(x = species, y = body_mass_g), size = 0.5, 
             position = position_jitter(width = 0.05), alpha = 0.3) +
  annotate("text", x = c(1, 2, 3), y = 7000, label = c("a", "b", "c")) +
  coord_cartesian(ylim = c(0, 8000))
gp

## Warning: Removed 2 rows containing missing values (geom_point).

できあがり！

正規性なし、等分散性ありの場合は…

• Kruscal-Wallis検定
  • kruskal.test(data = 【data.frame】, 【調べたい値の列名】~ 【グループ分けの列名】)
    
  • 帰無仮説：全ての群に差がない
    
  • p<0.05で異なる母集団に由来することを示す
    
  • ただし近似のp値を計算するのでサンプルサイズが小さい場合は注意

post-hocは…

• Steel-Dwass検定
  • NSM3パッケージを使う
  • pSDCFlig(【調べたい値をベクトルで】, 【グループ分けをベクトルで】, method = 【方法】)
    • methodは３種類
      • “Exact”…正確
      • “Asymptotic”… 漸近
      • “Monte Carlo”…モンテカルロ法
        • 正確以外の場合はサンプルサイズに注意
  • 帰無仮説：２群間に差がない
    
  • p<0.05で異なる母集団に由来することを示す

ペンギンの体重は種ごとに異なるか？（Kruscal-Wallis）

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

## 
##  Kruskal-Wallis rank sum test
## 
## data:  body_mass_g by species
## Kruskal-Wallis chi-squared = 115.24, df = 2, p-value < 2.2e-16

1. まずはプロット
2. 正規性, 等分散性, 対応あり/なしを確認
   • Q-Q plot, Shapiro-Wilkテスト
   • Bartlettテスト
3. 検定方法を決める
   • ANOVA → Tukey
   • Kruscal-Wallis → Steel-Dwass
   • フローチャートを参照

２要因の場合は…

• Two-way ANOVAやFriedman検定など

統計モデリング

• 応答変数（Y）を説明変数（X）の式にあてはめる
• （= 関係性を説明する）
• 実は質的変数にも普通に使え、ANOVA等と同等の解析も可能

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

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

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

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

• glm(【線形予測子】, family = 【誤差構造】(link = 【リンク関数】), data = 【データフレーム】)
  • 例：glm(bill_length_mm ~ bill_depth_mm + species, 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 x 8
##   species island bill_length_mm bill_depth_mm flipper_length_… body_mass_g sex  
##   <fct>   <fct>           <dbl>         <dbl>            <int>       <int> <fct>
## 1 Adelie  Torge…           39.1          18.7              181        3750 male 
## 2 Adelie  Torge…           39.5          17.4              186        3800 fema…
## 3 Adelie  Torge…           40.3          18                195        3250 fema…
## 4 Adelie  Torge…           NA            NA                 NA          NA <NA> 
## 5 Adelie  Torge…           36.7          19.3              193        3450 fema…
## 6 Adelie  Torge…           39.3          20.6              190        3650 male 
## # … with 1 more variable: 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., １時間あたりの交通事故数、メッシュ区画内の生物個体数

\(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\)×種)
    • 説明変数同士は独立である必要あり（相関関係などがない）

植物サイズと処理は種子生存率に影響するか？（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. さらに自由なモデリング（ベイズ統計など）

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

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 x 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】$【グループ分けの列名】
  • 帰無仮説：２群の比率に差がない
    
  • サンプル数が小さい場合はフィッシャーの正確確率検定
• フィッシャーの正確確率検定
  • 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

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