※資料作成は岩嵜航さん(東北大学)にご協力いただきました。

(左右キーで進みます!)

Return to HOME

統計解析をしよう!

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

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

3. データを整理する

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


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

広くて深い統計解析の世界

  • データの荒海を泳ぎ切ってもどこにも 「究極の真実」 などありはしないのだ。
  • 統計学はその時その場かぎりでの 「最良の結論」 を導く便法にすぎないのだ。
    (三中信宏)

この実習で扱うのは…

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

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

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

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

この実習で扱わないこと

  • 理論背景や難しい計算
  • 発展的な解析(ベイズ統計とか)
Return to HOME

どのように統計解析をするのか

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

    • 「A群とB群は異なる」とか「Aが増えるとBも増える」とか
  2. 自分のデータの性質を知り、適切な統計手法を探す

    • 量的(連続/離散)か質的か、分布や分散など
      • 帰無仮説や統計モデルを設定する
  3. 統計手法を実践する方法を探す

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

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

Return to HOME

統計の基本(1/3)母集団と標本

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

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

統計の基本(2/3)応答変数と説明変数

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

Return to HOME

統計の基本(3/3)どんな変数なのか

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

選び方:Step1
統計モデリングか統計的仮説検定か

  • △…傾き = 1の直線的な関係じゃない限り、素朴すぎる
  • XとYの関係が複雑な場合は実験系を見直すか、ベイズ統計などを用いる
Return to HOME

まず仮説検定かモデリングかを選ぼう

統計モデリング

  • 応答変数(Y)を説明変数(X)の式にあてはめる
  • (= 関係性を説明する)

仮説検定

  • 各標本群が異なる母集団から得られたかを調べる
  • (= 差があるかを検定する)

Return to HOME

モデリングか仮説検定か(1/4)

応答変数(Y)説明変数(X)量的な場合…

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

モデリングか仮説検定か(2/4)

応答変数(Y)質的説明変数(X)量的な場合…

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

  • モデリング

    • 一般化線形モデル(ロジスティック回帰)
      • 応答変数のカテゴリ数が2個のみ(3以上は難しい)
Return to HOME

モデリングか仮説検定か(3/4)

応答変数(Y)量的説明変数(X)質的な場合…

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

モデリングか仮説検定か(4/4)

応答変数(Y)説明変数(X)質的な場合…

  • 仮説検定
    • 割合の検定(独立性の検定)
      • χ(カイ)二乗検定
      • フィッシャーの正確確率検定
Return to HOME

どのように統計解析をするのか

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

    • 「A群とB群は異なる」とか「Aが増えるとBも増える」とか
  2. 自分のデータの性質を知り、使えそうな統計手法を探す

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

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

統計的仮説検定を使う場面とは

  • 例えばこんなとき…
  • あなたは「良い餌で育ったハエは大きくなる」ことを主張しようとしてデータを集めた。
  • これまでのデータによると、普通の餌で育てたハエの体長は平均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であった。
  • この結果から、あなたは良い餌で育ったハエは大きくなると結論づけた。
  • ところが教授はこれに対して 「たまたま大きいハエがとれたんじゃないの?」とコメントした。
  • さて、この教授をどのように納得させたらよいだろうか?
Return to HOME

仮説検定の考え方

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

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

対立仮説 (偶然ではない=異なる母集団から得られた)

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

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

\(p\)値をどのように求めるか(1/2)

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

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

確率分布って何?

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

\(p\)値をどのように求めるか(2/2)

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

簡単な例(p.16)で見てみよう(1/2)

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

簡単な例(p.16)で見てみよう(2/2)

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

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

検定統計量もいろいろある

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

  • \(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\)に意味はない

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

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

仮説検定のよくある間違い

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

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

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

  • サンプルサイズが大きいと精度が上がる →僅かな差でも有意差が検出される
  • \(p\)値が小さいからといって差が大きいとは限らない
  • 差の大きさを示すのは効果量(effect size)
Return to HOME

選び方:Step2−1 仮説検定法(1/2)

1要因2群の検定 & 割合の検定

Return to HOME

選び方:Step2−1 仮説検定法(2/2)

多重比較(3群以上) & 2要因の検定

Return to HOME

いよいよ…仮説検定をしてみよう!

仮説検定をしてみよう!

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

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>
  • アデリーペンギンの体重に性差があるか?
Return to HOME

仮説検定をしてみよう!

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

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

[練習問題]

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

仮説検定をしてみよう!

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

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>
Return to HOME

仮説検定をしてみよう!

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

gp = ggplot(data = penguins_adelie) +
  geom_point(aes(x = sex, y = body_mass_g), 
             position = position_jitter(width = 0.05))
gp
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

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

仮説検定をしてみよう!

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

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>
Return to HOME

仮説検定をしてみよう!

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

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

  • いい感じ!さて…
Return to HOME

仮説検定をしてみよう!

  • Y(体重)は量的、X(性)は質的で1種類、群数は2(=オス or メス)
  • データの 正規性 等分散性によって使える方法は変わる
Return to HOME

正規性:標本が正規分布に従うかどうか

正規分布とは

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

等分散性:分散の程度が違うかどうか

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

仮説検定をしてみよう!

仮説検定の進め方

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

仮説検定をしてみよう!

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で直線上に乗っているような気がする
Return to HOME

仮説検定をしてみよう!

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

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

仮説検定をしてみよう!

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
  • 正規性がないとは言えない→「正規性あり」と判断
Return to HOME

仮説検定をしてみよう!

2. 等分散性があるか

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

仮説検定をしてみよう!

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なので「等分散性なし」と判断
Return to HOME

正規性のある2群の比較

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

  • 正規性はあり
  • 等分散性はなし
  • Studentのt検定(正規性と等分散性が前提条件)は使えない
  • → Welchのt検定(正規性のみが前提条件)を使う
Return to HOME

正規性のある2群の比較

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

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

正規性のある2群の比較

アデリーペンギンの体重に性差があるか?(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
  • 帰無仮説「平均が等しい」は棄却→有意差あり
Return to HOME

正規性のある2群の比較

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

gp = ggplot(data = penguins_adelie) +
  geom_point(aes(x = sex, y = body_mass_g), size = 0.8,
             position = position_jitter(width = 0.07), alpha = 0.5) +
  annotate("text", x = 1.5, y = 5200, label = "***", size = 5, color = "gray30") +
  annotate("segment", x = 1.2, xend = 1.8, y = 5000, yend = 5000, color = "gray30") +
  coord_cartesian(ylim = c(0, 6000))
gp

ようやくグラフが完成!

Return to HOME

統計解析の結果を出力する

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

統計解析の結果を出力する(1/2)

文章として出力する(効果量も教えてくれる)

  • ※1行目の#を消してインストール
#install.packages("easystats")
library(easystats)
report(t.test(data = penguins_adelie, body_mass_g ~ sex)) 
## Effect sizes were labelled following Cohen's (1988) recommendations.
## 
## The Welch Two Sample t-test testing the difference of body_mass_g by sex (mean
## in group female = 3368.84, mean in group male = 4043.49) suggests that the
## effect is negative, statistically significant, and large (difference = -674.66,
## 95% CI [-776.30, -573.01], t(135.69) = -13.13, p < .001; Cohen's d = -2.25, 95%
## CI [-2.68, -1.82])
Return to HOME

統計解析の結果を出力する(2/2)

CSVとして書き出す(フォルダにCSVができるはず)

report(t.test(data = penguins_adelie, body_mass_g ~ sex)) %>%
  as.data.frame %>%                        #データフレームの形にして
  write_csv(file = "result_ttest.csv") %>% #csvに書き出し
  print()                                  #表示
## Warning: Unable to retrieve data from htest object.
##   Returning an approximate effect size using t_to_d().
## Welch Two Sample t-test
## 
## Parameter   | Group | Mean_Group1 | Mean_Group2 | Difference |             95% CI | t(135.69) |      p |     d |          d  CI
## -------------------------------------------------------------------------------------------------------------------------------
## body_mass_g |   sex |     3368.84 |     4043.49 |    -674.66 | [-776.30, -573.01] |    -13.13 | < .001 | -2.25 | [-2.68, -1.82]
## 
## Alternative hypothesis: two.sided
Return to HOME

正規性のある2群の比較のいろいろ

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)
Return to HOME

正規性のない2群の比較

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

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

正規性のない2群の比較

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

[練習問題]

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

正規性のない2群の比較

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

[練習問題]

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>
Return to HOME

正規性のない2群の比較

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

[練習問題]

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

Return to HOME

正規性のない2群の比較

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

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

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

Return to HOME

正規性のない2群の比較

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

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
Return to HOME

正規性のない2群の比較

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

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
Return to HOME

正規性のない2群の比較

正規性のない2群の比較

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

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

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

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

まずインストールする(1行目の#を消して走らせる)

#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
  • 帰無仮説「2群間に差がない」は棄却されない→有意差なし
Return to HOME

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

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

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

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

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

正規性なし、等分散性なしの2群の比較

Return to HOME

正規性なし、等分散性なしの2群の比較

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

正規性なし、等分散性なしの2群の比較

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

まずインストールする(1行目の#を消して走らせる)

#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
Return to HOME

【まとめ】1要因2群間の検定法

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

多群(3群以上)の比較をしよう

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

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

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

Return to HOME

多群(3群以上)の比較をしよう

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

→もちろん 不正!!

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

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

多群(3群以上)の比較をしよう

選び方2−1:多重比較 & 2要因の検定

Return to HOME

多群(3群以上)の比較をしよう

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

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>
  • ペンギンの体重は種ごとに異なるか?
Return to HOME

多群(3群以上)の比較をしよう

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

  • 体重に性差があるのでメスに限定して調べよう
  • まずはデータをプロットしよう
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>
Return to HOME

多群(3群以上)の比較をしよう

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

  • まずはデータをプロットしよう
gp = ggplot(data = penguins_female) +
  geom_point(aes(x = species, y = body_mass_g), 
             position = position_jitter(width = 0.05))
gp

  • 差がありそうな雰囲気…!
Return to HOME

多群(3群以上)の比較をしよう

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

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

多群(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で直線上に乗っているような気がする
Return to HOME

多群(3群以上)の比較をしよう

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
  • 「正規性あり」と判断
Return to HOME

多群(3群以上)の比較をしよう

2. 等分散性があるか

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

多群(3群以上)の比較をしよう

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
  • 「等分散性あり」と判断
Return to HOME

多群(3群以上)の比較をしよう

分散分析(ANOVA)には種類がある

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

多群(3群以上)の比較をしよう

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

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

多群(3群以上)の比較をしよう

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

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テストを行う
Return to HOME

多群(3群以上)の比較をしよう

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

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

多群(3群以上)の比較をしよう

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

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
Return to HOME

多群(3群以上)の比較をしよう

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

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 or values outside the scale range
## (`geom_point()`).

できあがり!

Return to HOME

多群(3群以上)の比較をしよう

正規性や等分散性がない場合は…

  • 整列ランク変換(Aligned Rank Transform, ART)ANOVA
    • ARToolパッケージを使う
    • データを整列ランク変換
      • model = art(data = 【data.frame】, 【調べたい値の列名】~ 【グループ分けの列名】)
      • summary(model) → 全ての数値が0になっていることを確認
    • 分散分析を実行
      • anova(model)
    • ANOVAが有意だったらposthocを実行
      • art.con(model, “【グループ分けの列名】”)
    • p<0.05で異なる母集団に由来することを示す
    • 詳しい使い方や2-wayはこちら
  • 古典的にはKruscal-Wallis検定のあとにpost-hocとしてSteel-Dwassを行っていた
  • 最近は2要因を扱えるART ANOVAが主流になりつつある
Return to HOME

多群(3群以上)の比較をしよう

ペンギンの体重は種ごとに異なるか?(ART ANOVA)

  • データを整列ランク変換
#install.packages("ARTool") #インストールする(#を外す)
library(ARTool) #おまじない
model = art(data = penguins_female, body_mass_g ~ species)
summary(model)
## Aligned Rank Transform of Factorial Model
## 
## Call:
## art(formula = body_mass_g ~ species, data = penguins_female)
## 
## Column sums of aligned responses (should all be ~0):
## species 
##       0 
## 
## F values of ANOVAs on aligned responses not of interest (should all be ~0):
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 
  • 全ての値が0であることを確認する
Return to HOME

多群(3群以上)の比較をしよう

ペンギンの体重は種ごとに異なるか?(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を行う
Return to HOME

多群(3群以上)の比較をしよう

ペンギンの体重は種ごとに異なるか?(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

【まとめ】1要因多群間の検定法

  1. まずはプロット
  2. 正規性, 等分散性, 対応あり/なしを確認
    • Q-Q plot, Shapiro-Wilkテスト
    • Bartlettテスト
  3. 検定方法を決める

2要因の場合は…

  • 2-way ANOVA(Type IIかIII)や2-way ART-ANOVAを使う
Return to HOME

割合の検定をしてみよう!

1要因2群の検定 & 割合の検定

割合の検定をしてみよう!

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

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
Return to HOME

割合の検定をしてみよう!

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

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

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

割合の検定をしてみよう!(1/2)

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

割合の検定をしてみよう!(2/2)

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

割合の検定をしてみよう!

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

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
Return to HOME

割合の検定をしてみよう!

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

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

  • 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
Return to HOME

割合の検定をしてみよう!

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

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

  • 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
Return to HOME

【まとめ】割合の検定

どのように統計解析をするのか

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

    • 「A群とB群は異なる」とか「Aが増えるとBも増える」とか
  2. 自分のデータの性質を知り、適切な統計手法を探す

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

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

選び方:Step1
統計モデリングか統計的仮説検定か

  • △…傾き = 1の直線的な関係じゃない限り、素朴すぎる
  • XとYの関係が複雑な場合は実験系を見直すか、ベイズ統計などを用いる
Return to HOME

選び方:Step2−2統計モデリング

統計モデリング

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

Return to HOME

一般化線形モデル(GLM)とは(1/2)

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

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

一般化線形モデル(GLM)とは(2/2)

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

Coded by Heavy Watal

Return to HOME

GLMにはglm関数を使う(1/4)

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

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

Return to HOME

GLMにはglm関数を使う(2/4)

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

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

GLMにはglm関数を使う(3/4)

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

  • 正規分布(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
Return to HOME

GLMにはglm関数を使う(4/4)

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

  • そのまま(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が使える
Return to HOME

GLMを実践してみよう!

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>
  • 体重の大きいペンギンは翼長も長い?
Return to HOME

GLMを実践してみよう!

体重の大きいアデリーは翼長も長い?(まずはプロット)

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

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

GLMを実践してみよう!

体重の大きいアデリーは翼長も長い?(GLM)

glm_result = glm(flipper_length_mm ~ body_mass_g,        #線形予測子
                 family = gaussian(link = "identity"),   #誤差分布とリンク関数
                 data = penguins_adelie)                 #データ
summary(glm_result)                  #結果を格納したglm_resultの詳細を表示する
## 
## Call:
## glm(formula = flipper_length_mm ~ body_mass_g, family = gaussian(link = "identity"), 
##     data = penguins_adelie)
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.656e+02  3.918e+00   42.27  < 2e-16 ***
## body_mass_g 6.610e-03  1.049e-03    6.30  3.4e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 33.57457)
## 
##     Null deviance: 6167.5  on 145  degrees of freedom
## Residual deviance: 4834.7  on 144  degrees of freedom
## AIC: 931.33
## 
## Number of Fisher Scoring iterations: 2
Return to HOME

GLMを実践してみよう!

結果の見方

  • 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\)に回帰できる
Return to HOME

GLMを実践してみよう!

結果をcsvで出力する(easystatsのreport)

report(glm_result) %>%
  as.data.frame() %>%                     #データフレーム型にする
  write_csv(file = "result_glm.csv") %>%  # CSVで出力
  print()                                 #ここでも表示
## Parameter   | Coefficient |           95% CI | t(144) |      p | Std. Coef. | Std. Coef. 95% CI |    Fit
## --------------------------------------------------------------------------------------------------------
## (Intercept) |      165.60 | [157.92, 173.28] |  42.27 | < .001 |  -1.93e-15 |     [-0.14, 0.14] |       
## body mass g |    6.61e-03 | [  0.00,   0.01] |   6.30 | < .001 |       0.46 |     [ 0.32, 0.61] |       
##             |             |                  |        |        |            |                   |       
## AIC         |             |                  |        |        |            |                   | 931.33
## AICc        |             |                  |        |        |            |                   | 931.50
## BIC         |             |                  |        |        |            |                   | 940.28
## R2          |             |                  |        |        |            |                   |   0.22
## Sigma       |             |                  |        |        |            |                   |   5.79
Return to HOME

GLMを実践してみよう!

結果をプロットに重ねる

  • ggeffectsパッケージを使う(日本語解説
  • glm結果から、予測値を計算してくれる
  • インストールしてから(1行目の#を消して走らせる)、おまじないをして使う
#install.packages("ggeffects")
library(ggeffects) #予測値をglm_predictという変数に格納
glm_predict = ggpredict(glm_result, terms = "body_mass_g") #termsにはX軸の列名を入れる
glm_predict %>% as_tibble() #中身はこんなデータフレーム
## # A tibble: 21 × 6
##        x predicted std.error conf.low conf.high group
##    <int>     <dbl>     <dbl>    <dbl>     <dbl> <fct>
##  1  2800      184.     1.06      182.      186. 1    
##  2  2900      185.     0.972     183.      187. 1    
##  3  3000      185.     0.883     184.      187. 1    
##  4  3100      186.     0.797     185.      188. 1    
##  5  3200      187.     0.716     185.      188. 1    
##  6  3300      187.     0.642     186.      189. 1    
##  7  3400      188.     0.577     187.      189. 1    
##  8  3500      189.     0.526     188.      190. 1    
##  9  3600      189.     0.492     188.      190. 1    
## 10  3700      190.     0.480     189.      191. 1    
## # ℹ 11 more rows
Return to HOME

GLMを実践してみよう!

結果をプロットに重ねる(ggeffects)

  • 生データ(penguins_adelie)とglm予測値(回帰線)(glm_predict)を重ねる
plot(glm_predict, show_data = T)

Return to HOME

GLMを実践してみよう!

結果をプロットに重ねる(ggplot2)

  • 生データ(penguins_adelie)とglm予測値(回帰線)(glm_predict)を重ねる
gp = ggplot() +
  geom_point(data = penguins_adelie, aes(x = body_mass_g, y = flipper_length_mm), alpha = 0.5)  +
  geom_line(data = glm_predict, aes(x = x, y = predicted)) +
  geom_ribbon(data = glm_predict, aes(x = x, ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  annotate("text", x = 4200, y = 210, 
           label = "(Y ~ 0.006677X + 165.2, p = 3.4e-09)", size = 4)
gp

Return to HOME

GLMを実践してみよう!

残差が正規分布していたかを確認

Coded by Heavy Watal

shapiro.test(glm_result$residuals)  #残差の正規性をShapiro検定
## 
##  Shapiro-Wilk normality test
## 
## data:  glm_result$residuals
## W = 0.99179, p-value = 0.562
  • してたっぽいですね!
Return to HOME

直線回帰すればいいってもんじゃない

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

直線回帰すればいいってもんじゃない

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

Coded by Heavy Watal

Return to HOME

よく使うモデル1:ポアソン回帰

  • 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

Return to HOME

ポアソン分布とは

  • 平均\(\lambda\)で単位時間(空間)あたりに発生する事象の回数。
  • e.g., 1時間あたりの交通事故数、メッシュ区画内の生物個体数
## 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!}\)

Coded by Heavy Watal

Return to HOME

ポアソン回帰でGLMしよう

植物のサイズが種子数に影響するか?

  • 久保先生の種子数データを解析しよう
  • 上限のないカウントデータ(離散値)はポアソン回帰が良いかも
seed = 
  read_csv("https://kuboweb.github.io/-kubo/stat/iwanamibook/fig/poisson/data3a.csv") %>%
  rename(no.seed = y, plant_size = x, treatment = f)
head(seed)
## # A tibble: 6 × 3
##   no.seed plant_size treatment
##     <dbl>      <dbl> <chr>    
## 1       6       8.31 C        
## 2       6       9.44 C        
## 3       6       9.5  C        
## 4      12       9.07 C        
## 5      10      10.2  C        
## 6       4       8.32 C
Return to HOME

ポアソン回帰でGLMしよう

植物のサイズが種子数に影響するか?(まずはプロット)

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

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

ポアソン回帰でGLMしよう

植物のサイズが種子数に影響するか?(GLM)

glm_result = glm(no.seed ~ plant_size, 
                 family = poisson(link = "log"), #ポアソン回帰の設定
                 data = seed)
summary(glm_result)
## 
## Call:
## glm(formula = no.seed ~ plant_size, family = poisson(link = "log"), 
##     data = seed)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.29172    0.36369   3.552 0.000383 ***
## plant_size   0.07566    0.03560   2.125 0.033580 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 89.507  on 99  degrees of freedom
## Residual deviance: 84.993  on 98  degrees of freedom
## AIC: 474.77
## 
## Number of Fisher Scoring iterations: 4
Return to HOME

ポアソン回帰でGLMしよう

植物のサイズが種子数に影響するか?(GLM)

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

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

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

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

Return to HOME

ポアソン回帰でGLMしよう

植物のサイズが種子数に影響するか?(ggeffects)

glm_predict = ggpredict(glm_result, terms = "plant_size")
plot(glm_predict, show_data = T)
## Data points may overlap. Use the `jitter` argument to add some amount of
##   random variation to the location of data points and avoid overplotting.

Return to HOME

ポアソン回帰でGLMしよう

植物のサイズが種子数に影響するか?(ggplot2)

gp = ggplot() +
  geom_point(data = seed, aes(x = plant_size , y = no.seed)) +
  geom_line(data = glm_predict, aes(x = x , y = predicted)) +
  geom_ribbon(data = glm_predict, aes(x = x , ymin = conf.low, ymax = conf.high), alpha = 0.2) 
gp

Return to HOME

よく使うモデル2:ロジスティック回帰

  • 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

Return to HOME

二項分布とは

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

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

Coded by Heavy Watal

Return to HOME

ロジスティック回帰でGLMしよう

植物サイズは種子生存率に影響するか?(対データ)

  • 久保先生の種子の生存率データを解析しよう
  • 上限のあるカウントデータ(離散値)はロジスティック回帰が良いかも
seed2= 
  read_csv("https://kuboweb.github.io/-kubo/stat/iwanamibook/fig/binomial/data4a.csv") %>%
  rename(observed_seed = N, survived_seed = y, plant_size = x, treatment = f)
head(seed2)
## # A tibble: 6 × 4
##   observed_seed survived_seed plant_size treatment
##           <dbl>         <dbl>      <dbl> <chr>    
## 1             8             1       9.76 C        
## 2             8             6      10.5  C        
## 3             8             5      10.8  C        
## 4             8             6      10.9  C        
## 5             8             1       9.37 C        
## 6             8             1       8.81 C
Return to HOME

ロジスティック回帰でGLMしよう

植物サイズは種子生存率に影響するか?(まずプロット)

gp = ggplot(data = seed2) +
  geom_point(aes(x = plant_size , y = survived_seed, color = treatment))
gp

  • サイズが大きいほうが生存率が高いような…?
  • treatmentでも違いがありそう
Return to HOME

ロジスティック回帰でGLMしよう

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

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

ロジスティック回帰でGLMしよう

植物サイズは種子生存率に影響するか?(データ整形)

  • 2値データ(成功数と失敗数)を示す2列を作る
  • 今回は生存した種子数(survived_seed)と死んだ種子数(dead_seed)
seed3 = seed2 %>% 
  mutate(dead_seed = observed_seed-survived_seed)
head(seed3)
## # A tibble: 6 × 5
##   observed_seed survived_seed plant_size treatment dead_seed
##           <dbl>         <dbl>      <dbl> <chr>         <dbl>
## 1             8             1       9.76 C                 7
## 2             8             6      10.5  C                 2
## 3             8             5      10.8  C                 3
## 4             8             6      10.9  C                 2
## 5             8             1       9.37 C                 7
## 6             8             1       8.81 C                 7
Return to HOME

ロジスティック回帰でGLMしよう

植物サイズは種子生存率に影響するか?(GLM)

glm_result = glm(cbind(survived_seed, dead_seed) ~ plant_size, 
                 family = binomial(link = "logit"), #ロジスティック回帰の設定
                 data = seed3)
summary(glm_result)
## 
## Call:
## glm(formula = cbind(survived_seed, dead_seed) ~ plant_size, family = binomial(link = "logit"), 
##     data = seed3)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -13.7785     1.0773  -12.79   <2e-16 ***
## plant_size    1.4626     0.1107   13.21   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 499.23  on 99  degrees of freedom
## Residual deviance: 217.17  on 98  degrees of freedom
## AIC: 364.35
## 
## Number of Fisher Scoring iterations: 5
Return to HOME

ロジスティック回帰でGLMしよう

植物サイズは種子生存率に影響するか?(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)の効果も考慮に入れてみよう!
  • 重回帰(説明変数が2つ:植物サイズ + 処理)
Return to HOME

GLMで重回帰をしよう

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

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

GLMで重回帰をしよう

植物サイズと処理は種子生存率に影響するか?(GLM)

glm_result = glm(cbind(survived_seed, dead_seed) ~ plant_size + treatment, #重回帰
                 family = binomial(link = "logit"), #ロジスティック回帰の設定
                 data = seed3)
summary(glm_result)
## 
## Call:
## glm(formula = cbind(survived_seed, dead_seed) ~ plant_size + 
##     treatment, family = binomial(link = "logit"), data = seed3)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -19.5361     1.4138  -13.82   <2e-16 ***
## plant_size    1.9524     0.1389   14.06   <2e-16 ***
## treatmentT    2.0215     0.2313    8.74   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 499.23  on 99  degrees of freedom
## Residual deviance: 123.03  on 97  degrees of freedom
## AIC: 272.21
## 
## Number of Fisher Scoring iterations: 5
Return to HOME

GLMで重回帰をしよう

植物サイズと処理は種子生存率に影響するか?(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)}}\)に回帰できる
Return to HOME

GLMで重回帰をしよう

植物サイズと処理は種子生存率に影響するか?(ggeffects)

glm_predict = ggpredict(glm_result, terms = c("plant_size", "treatment"))
glm_predict %>% as_tibble()
## # A tibble: 22 × 6
##        x predicted std.error conf.low conf.high group
##    <dbl>     <dbl>     <dbl>    <dbl>     <dbl> <fct>
##  1   7.5   0.00744     0.389  0.00349    0.0158 C    
##  2   7.5   0.0536      0.276  0.0319     0.0887 T    
##  3   8     0.0195      0.325  0.0104     0.0362 C    
##  4   8     0.131       0.220  0.0890     0.188  T    
##  5   8.5   0.0502      0.263  0.0306     0.0812 C    
##  6   8.5   0.285       0.172  0.222      0.358  T    
##  7   9     0.123       0.206  0.0857     0.173  C    
##  8   9     0.514       0.144  0.444      0.584  T    
##  9   9.5   0.271       0.159  0.214      0.337  C    
## 10   9.5   0.738       0.146  0.679      0.789  T    
## # ℹ 12 more rows
Return to HOME

GLMで重回帰をしよう

植物サイズと処理は種子生存率に影響するか?(ggeffects)

plot(glm_predict, show_data = T)
## Raw data not available.

Return to HOME

GLMで重回帰をしよう

植物サイズと処理は種子生存率に影響するか?(ggplot2)

gp = ggplot(data = seed2) +
  geom_point(aes(x = plant_size , y = survived_seed/observed_seed, color = treatment)) +
  geom_line(data = glm_predict, aes(x = x , y = predicted, color = group)) +
  geom_ribbon(data = glm_predict, aes(x = x , ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2) 
gp

  • 植物体が大きいと種子の生存率が有意に良くなる
  • treatmentをすると種子の生存率が有意に良くなる
Return to HOME

ロジスティック回帰でGLMしよう2

母親の体重は低体重出生に影響するか(2値データ)

#install.packages("MASS")  #まずインストール
library(MASS)              #まず読み込み
head(birthwt, n = 2)       #MASSに入っているbirthwtというサンプルデータを使う
##    low age lwt race smoke ptl ht ui ftv  bwt
## 85   0  19 182    2     0   0  0  1   0 2523
## 86   0  33 155    3     0   0  0  0   3 2551
  • low:新生児の体重が2.5kg未満か否か(2.5kg未満が1)
  • age:母親の年齢
  • lwt:母親の体重(ポンド単位)
  • smoke:母親の喫煙有無(喫煙有りが1)
  • ptl:母親の早産経験の有無(経験有りが1)
  • ht:母親の高血圧の有無(有りが1)
  • ui:母親の子宮神経過敏の有(有りが1)
Return to HOME

ロジスティック回帰でGLMしよう2

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

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

  • 体重(lwt)が多いほうが低体重出生ではない(y = 0)傾向…?
Return to HOME

ロジスティック回帰でGLMしよう2

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

glm_result = glm(low ~ lwt,
                 family = binomial(link = "logit"), #ロジスティック回帰の設定
                 data = birthwt)
summary(glm_result)
## 
## Call:
## glm(formula = low ~ lwt, family = binomial(link = "logit"), data = birthwt)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  0.99831    0.78529   1.271   0.2036  
## lwt         -0.01406    0.00617  -2.279   0.0227 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 228.69  on 187  degrees of freedom
## AIC: 232.69
## 
## Number of Fisher Scoring iterations: 4
Return to HOME

ロジスティック回帰でGLMしよう2

母親の体重は低体重出生に影響するか(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)}}\)に回帰できる
Return to HOME

ロジスティック回帰でGLMしよう2

母親の体重は低体重出生に影響するか(ggeffects)

glm_predict = ggpredict(glm_result, terms = "lwt[all]")
plot(glm_predict, show_data = T)
## Data points may overlap. Use the `jitter` argument to add some amount of
##   random variation to the location of data points and avoid overplotting.

Return to HOME

ロジスティック回帰でGLMしよう2

母親の体重は低体重出生に影響するか(ggplot2)

gp = ggplot() +
  geom_point(data = birthwt, aes(x = lwt , y = low)) +
  geom_line(data = glm_predict, aes(x = x , y = predicted)) +
  geom_ribbon(data = glm_predict, aes(x = x , ymin = conf.low, ymax = conf.high), alpha = 0.2) 
gp

  • 母親の体重(lwt)が多いと有意に低体重出生率が低下する
  • 回帰曲線は低体重出生の割合を示すことに注意
Return to HOME

いろいろな分布 + リンク関数が選べる

  • 正規分布(gaussian)…確率変数の和、平均値

    • \(-∞\leqq Y < ∞\)(下限なし, 上限なし), 連続値
    • 選べるリンク関数:identity
  • 二項分布(binomial)…成功率p、試行回数nのうちの成功回数

    • \(0\leqq Y \leqq max\)(下限あり, 上限あり), 離散値
    • 選べるリンク関数:logit, probit, cloglog
  • ポアソン分布(poisson)…単位時間あたり平均λ回起こる事象の発生回数

    • \(0\leqq Y \leqq ∞\), (下限あり, 上限なし), 離散値
    • 選べるリンク関数:log, indentity, sqrt
  • ガンマ分布(Gamma)…ポアソン過程でk回起こるまでの待ち時間

    • \(0\leqq Y \leqq ∞\), (下限あり, 上限なし), 連続値
    • 選べるリンク関数:inverse, identity, log
  • 他にinverse.gaussian, quasi, quasibinomial, quasipoisson.

  • それぞれの分布やリンク関数の形を調べて選ぼう

Return to HOME

どのモデルを選んだら良いの?

“All models are wrong, but some are useful”

  • 唯一無二の完璧なモデルは存在しない
  • 目的に応じて少しでも「マシな」モデルを選ぶ
  1. メカニズム的に納得できるもの
    • ポアソン過程のカウントデータならポアソン分布
      • 過分散(overdispersion)があったら負の二項分布など
    • n回中k回のように割合的なカウントデータなら二項分布
  2. データを可視化してみて、それっぽい形・性質のもの
    • 左右対称のひと山ならとりあえず正規分布
    • 負の値をとらないならガンマ分布
    • 直線的か、指数関数的か、頭打ちかなど
Return to HOME

【まとめ】統計モデリング

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

自分のデータを統計解析しよう

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