心理学者のための統計モデリング 清水裕士 関西学院大学社会学部.

Presentation on theme: "心理学者のための統計モデリング 清水裕士 関西学院大学社会学部."— Presentation transcript:

1 心理学者のための統計モデリング 清水裕士 関西学院大学社会学部

2 自己紹介 清水裕士 専門 趣味 Web 関西学院大学社会学部 社会心理学・グループダイナミックス 心理統計 フリーの統計ソフトHADの開発

3 導入

4 本講演のテーマ GLMM 階層ベイズ Generalized Linear Mixed Model 一般化線形混合モデル
階層モデルをベイズで解く 階層モデルとは？ ベイズで解くメリットは？

5 GLMM 線形モデル GLMM = G + M + LM LM: Linear Model 正規分布を仮定した，線形モデル
重回帰，分散分析，共分散分析などを含む GLMM = G + M + LM GLM：Generalized Linear Model　一般化する話 LMM：Linear Mixed Model　混ぜる話

6 なぜGLMM?

7 「みどり本」 久保(2012)　データ解析のための統計モデリング入門

8

9

10 なぜ心理学でGLMM? 分散分析の限界 GLMM 正規分布するデータしか扱えない 変量効果を組み合わせると十分な推定ができない
仮説検定の枠組みでしか理解ができない GLMM 正規分布から解放される 変量効果を複数組み合わせられる 統計モデリングの視点を導入できる

11 このスライドの概要 LMの話：LM Gの話：GLM 回帰分析のおさらい 心理統計学と統計モデリングの違い 一般化線形モデルの話

12 続きのスライドの概要 Mの話：LMM GLMMと階層ベイズ 線形混合モデルの話 階層モデルの話 GとMが両方入ることの問題
HLMはここにはいる 階層モデルの話 GLMMと階層ベイズ GとMが両方入ることの問題 スマートな解決法としての階層ベイズ 階層ベイズをRで実行する

13 サンプルデータ なぜか2014年度のプロ野球データ 使う変数は・・・ プロ野球データフリークからダウンロード 野手199名の成績と翌年の年俸
野手199名の成績と翌年の年俸 GLMMを扱うにはちょうどいいデータ ここからダウンロード可 ただし，Macユーザーは文字コードの変換が必要 使う変数は・・・ 年俸，HR，HIT，打数あたりです。

14 プロ野球の打者成績と年俸

15 Rにデータを読み込む baseball.csvを読み込む わからない人向け 読み込み方は自由
データのオブジェクト名はdatにしておいてください わからない人向け Rの作業ディレクトリと同じフォルダに保存 >dat <- read.csv("baseball.csv") と書く dat <- read.csv("baseball.csv",fileEncoding = "CP932") Macの人は上のコードを実行してください

16 おさらい 線形モデル

17 回帰分析 ある変数で，別の変数を予測する 回帰分析 たとえば，野球選手の野手の年俸を予測する
ホームラン1本打ったら，年俸はいくら増えるだろうか？ 年俸を決める目安を知っておくと，次の年俸を決めるのに便利。 回帰分析

18 回帰係数 回帰式とは Y＝a + b X + e あるいは，Y’ = a + b X (ただし，Y’=Y – e)
Y’は予測値，aが切片，bが傾き，eは誤差 野球選手の年俸の例で言えば・・・ 年俸　=　a +　b ホームラン数 + e bをホームラン数の年俸に対する回帰係数という HRが1本増えたときの，平均的な年俸の増額を意味する

19 散布図

20 HRが一本増えたときに予想される年俸の平均増加額
散布図 HRが0本のときに予想される年俸 HRが一本増えたときに予想される年俸の平均増加額 回帰係数b 切片a 予測値Y’

21 回帰係数の計算方法 回帰係数の算出 Y ＝ a + b X + e 要は，回帰係数は，Xが1単位増加したときの，Yの増加量
eの二乗和を最も小さくするようにbを決定する 最小二乗法という方法で計算・・・説明は省略 要は，回帰係数は，Xが1単位増加したときの，Yの増加量 回帰係数b　=　共分散Sxy ／ Xの分散 切片a　=　Yの平均値－回帰係数×Xの平均値

22 Rで回帰分析 lm()を使う >fit.lm <- lm(salary ~ HR, dat)
>summary(fit.lm) HR1本打つと，675万年俸が増える HRが0本でも3200万の年俸はもらえる

23 予測直線も作ってみる plot()で散布図を作ってから，線を書く >plot(dat$HR,dat$salary)
>abline(fit.lm)

24 年俸はどうみても正規分布じゃない salaryのヒストグラム >hist(dat$salary)

25 年収って，対数正規分布らしい

26 OK，対数変換してやろうじゃない 年俸の対数の分布 >hist(log(dat$salary)) それっぽくなった！

27 対数年俸とHRの回帰モデル 対数年俸を目的変数に回帰 >fit.lm2 <- lm(log(salary) ~ HR, dat)
>summary(fit.lm2) まぁどっちにしても有意やな（

28 心理統計としての回帰分析 従属変数が正規分布に従ってるかチェック 有意な説明変数を探す 有意な標準化偏回帰係数を報告する
してようがしてまいが，結局は重回帰にかける 場合によっては対数変換して重回帰 有意な説明変数を探す ステップワイズ法を使って有意な変数だけを選ぶこともある 有意な標準化偏回帰係数を報告する おおよその効果量として使える これは効果が大きいね，小さいね，っていう話をする

29 心理統計としての回帰分析 実験と分散分析 回帰分析も基本的には同じ使われ方 効果があるかどうかに興味がある
実験室内の行動を予測してもあまり意味がない 行動に違いがあるかどうかだけが知りたい 回帰分析も基本的には同じ使われ方 変数に効果があるかどうかに興味がある 尺度値を予測してもあまり意味がない 尺度値が変数に応じて変わるかどうかが知りたい

30 回帰分析は予測にも使えるよ HRが10本の人の年俸 3200万＋675万×10 = 9975万 予測区間を使えばもっと便利

31 心理統計から統計モデリングへ GLMMを理解するために モデリングが心理学で必要かどうか
まずは，回帰分析を例に統計モデリングの世界を理解する必要がある 心理統計とは，哲学が若干違う点に注意 モデリングが心理学で必要かどうか この講演の最後にまた心理統計との関連を議論 明日あたり，飲みながら議論しましょう

32 統計モデリング

33 統計モデリング データ発生メカニズムを知る 確率分布のパラメータを推測する データが何らかの分布から発生している
その発生メカニズムを「確率分布」を使って表現 データの予測に興味がある 確率分布のパラメータを推測する 確率分布の形を決める値のこと 正規分布なら，平均値と分散 正規分布以外の確率分布も扱う

34 真のモデル≒データ発生メカニズム 直接は知りえない 直接観測できる データの発生 真の分布 データ

35 真の分布≒乱数生成装置 分布に従って乱数を発生させる装置 現象は乱数生成装置から出た乱数 知りたいこと
平均50，標準偏差10の正規分布の場合 50を中心に，30～70ぐらいの値がでてくる 現象は乱数生成装置から出た乱数 ホームランが30本の人，20本の人・・・ 知りたいこと 真の分布の種類と，そのパラメータ

36 例：平均値を推定する 平均値の推定も統計モデリング Y ~ N(μ, σ) Yが平均μ，標準偏差σの正規分布に従う
正規分布という確率モデルの，平均値というパラメータを推定する いわゆる算術平均は，正規分布の平均値パラメータの最尤推定量

37 平均値を推定する 年俸を無理やり正規分布に当てはめて推定 >fit.lm0 <- lm(salary~1,dat)
>summary(fit.lm0) 平均=74.66，標準偏差=87.95 平均値の信頼区間　62.36 ~ 86.95

38 真の分布と統計モデリング 真の分布は完全にはわからない 統計モデリングの目的 真の分布が何かは直接にはわからない
われわれが知ってる確率分布とも限らない でも，その分布から生成されたデータのいくつかは手に入っている 統計モデリングの目的 限られたデータから，最も真の分布に近いと思われる確率分布モデル＝予測分布を作成する

39 予測分布 仮定した確率モデルから予測される分布 平均=74.66，標準偏差=87.95の正規分布
95%予測区間=　 ~ 予測値が負！

40 真の分布と統計モデリング 直接は知りえない 直接観測できる 真の分布 データ 統計モデリング

41 どのモデルが真の分布に近いか どのようなモデルが「よい」のか 予測誤差が小さいモデル≒よいモデル？ モデルの良さには，いろんな基準がある
真の分布とモデルの近さをどう決めるか 予測誤差が小さいモデル≒よいモデル？ 将来をよりよく予測できるモデルがよい 真の分布が生成するデータをあてることができるようなモデルを選びたい

42 平均値を違う確率モデルで推定 対数正規分布で推定 平均値 = 72.50 平均値の95%信頼区間 = 63.15~ 83.23
95％予測区間 = 6.36~ 点推定値はそれほど変わらないが，信頼区間，予測区間は大きく異なる

43 正規分布と対数正規分布 どちらが真の分布に近い？ パラメータはそれ単体ではなく予測に使う 平均値の推定値はあまりかわらない
74.66と72.50 信頼区間は少し違うが，まぁあまりかわらない じゃあどちらでもいい？ パラメータはそれ単体ではなく予測に使う パラメータは分布の形を決めるために推定 データの分布と予測分布の近さを評価する

44 データのあてはまりを考える 今回のデータとモデルの近さ 尤度＝データとモデルの近さ
モデルを正しいとした場合の，データが得られる確率で評価することができる これを尤度と呼ぶ 尤度＝データとモデルの近さ 尤度を大きくするモデルは，データを最もうまく説明できているモデルといえる -2*対数尤度は，逆にモデルの予測とデータの誤差を意味している

45 正規分布と対数正規分布 対数尤度を比較 対数正規分布のほうが対数尤度が大きい 正規分布 対数正規分布 対数尤度 = -1172.74
-2*対数尤度 = 対数尤度 = -2*対数尤度 =

46 今度はモデルを複雑にしてみる また正規分布に戻って・・・ 今度は回帰分析

47 パラメータを増やすと尤度は増える 回帰直線と予測区間 説明変数が増えると，予測区間は狭くなる →予測精度が上がっている = 尤度も高くなる

48 データへのあてはまりをよくする モデルの複雑さとデータのあてはまり オーバーフィッティング
モデルを複雑にすればするほど，データとの当てはまりのよさ（尤度）は高くなる 複雑なモデル＝パラメータが多いモデル オーバーフィッティング しかし，それはデータ生成の確率的な変動を無視して，過剰にあてはまりを高くしてしまう これをオーバーフィッティングという 真の分布の代わりにデータによる推定値を用いることのバイアスは，オーバーフィッティングによって生じる

49 真のモデルが一次式のデータ 破線は予測区間 真のモデルは y ~ N( x, 0.8)

50 2次式で推定

51 3次式

52 4次式

53 データとの当てはまり 次数が上がるほど，データとの当てはまりはよい

54 将来のデータを予測したい 今回のデータ「だけ」に当てはまっても仕方ない モデルの複雑さをむやみに高くしないこと
どのデータにも，サンプリングのバイアスがある たまたまこういうデータである，という可能性 今回のデータだけに当てはまっても，将来の予測が上手くいくとは限らない モデルの複雑さをむやみに高くしないこと パラメータが多いと尤度は高くなる パラメータが多すぎると予測誤差は大きくなる ほどほどのバランスが必要

55 将来のデータへの予測誤差 モデルは固定して，同じ真の分布から発生した，別の50個のデータに当てはめたときの‐2*対数尤度の平均

56 でも実際は真のモデルはわからない たった一つのデータから予測誤差が知りたい できた そんなことできるわけ・・・ 情報量規準
モデルの尤度にパラメータ数の関数を罰則として加えることで，たった一つのデータから平均的な予測誤差の期待値を推定できる

57 平均的な予測誤差を評価する 赤池情報量規準（AIC） AICの特徴 モデルが最尤法で推定できて，サンプルサイズが十分大きいとき，
この値は，平均予測誤差の近似値となる AICの特徴 AICは，データへのあてはまりの良さ（尤度）と，パラメータ数で決まる 対数尤度が高いほど，パラメータ数が少ないほど平均予測誤差は小さくなる

58 情報量規準最小のモデルを選ぶ AICで予測誤差が小さいモデルを選ぶ AICの仮定
パラメータ（データではない）が正規分布になる サンプルサイズが大きい 真の分布をモデルが含んでいる

59 AICの比較 1次の情報量規準が最も小さい

60 いろんな情報量規準 cAIC WAIC BIC AICの有限サンプルサイズバージョン AICのベイズ版 より拡張的な指標
真の分布が正規分布の場合のみ利用可能 WAIC AICのベイズ版　より拡張的な指標 真の分布がモデルに含まれてなくても使える パラメータが正規分布で近似できなくても使える BIC モデルが正しい確率が高いモデルを選ぶ

61 統計モデリングとしての回帰分析

62 正規分布を使って予測 年俸をHRで予測したい

63 数式で表現してみる 回帰式 統計モデリング的な表現 Y = α + β X + e ただし，e ~ N(0, σ) 𝑌 = α + β X
残差は，平均0，標準偏差σの正規分布に従う 統計モデリング的な表現 𝑌 = α + β X 𝑌 はYの予測値 Y ~ N( 𝑌 , σ) Yは平均 𝑌 ，標準偏差σの正規分布に従う ただし， 𝑌 はXの値によって変化する条件付きのパラメータ

64 回帰分析の仮定 データは独立に生成されている データは同じ（分散の）分布から生成されている データは正規分布から生成されている
独立性（正確には無相関の仮定） データ内に局所的な相関がない データは同じ（分散の）分布から生成されている 同一性（均一分散の仮定） データはすべて同じ（分散の）分布に従っている データは正規分布から生成されている 正規性

65 正規分布を使って予測 予測直線を中心とした正規分布に従う かなり無理がある・・ 全然同一じゃない・・ 予測外れまくり・・・

66 予測分布からデータ生成 真の分布の近似 N( 𝑌 , σ)に推定したパラメータを入れたものが予測分布
>predict <- fit.lm$ >sigma <- (sum(fit.lm$residuals^2) / 199) ^ 0.5 >y <- rnorm(199,predict,sigma) >hist(y) 予測分布から生成したデータ 本来の年俸のデータ

67 統計モデリングと回帰分析 説明変数の有意性はそれほど重視しない 変数変換をすると情報量規準を比較できない
モデル全体でよい予測ができるかどうか もちろん有意性検定は，AICを向上させるのに役立つ場合はある 変数変換をすると情報量規準を比較できない 対数変換をするのではなく，対数正規分布で直接モデリングするほうがよい

68 一般化線形モデル

69 一般化線形モデルとは 線形モデル 一般化線形モデル データが正規分布から生成されたと仮定するモデル
データが指数分布族の確率分布から生成されたと仮定するモデル

70 指数分布族 シンプルな構造を持った確率分布 離散分布 連続分布 いろいろわかりやすい性質を持っている
二項分布，ポアソン分布，負の二項分布など 連続分布 正規分布，対数正規分布，ガンマ分布など

71 対数正規分布

72 対数正規分布とは 対数をとると正規分布になる分布 年俸を対数正規分布に当てはめる 大きいものほどより大きくなるような性質のもの 年収データ
ネットワークサイズ 年俸を対数正規分布に当てはめる Rでは，直接対数正規分布に当てはめる関数はたぶんない 自分で作ればある

73 パラメータの推定は・・・ 対数変換してlm()で求めればよい
>fit.lm.log <- lm(log(salary)~HR,dat) >summary(fit.lm.log) 回帰係数は対数変換された値なので，直接は解釈が難しい

74 予測曲線と予測区間 対数正規分布を当てはめたモデル 正規分布よりはうまいこと当てはまっている

75 予測分布 対数正規分布を使って乱数を生成 predict <- fit.lm.log$fitted.values
sigma <- (sum(fit.lm.log$residuals^2) / 199) y <- rlnorm(199,predict,sigma) hist(y) 予測分布から生成したデータ 本来の年俸のデータ

76 ベルヌーイ分布

77 2値データの予測 年俸が1億をこえるかをHRで予測したい 1億未満と1億以上を区別する変数

78 1億を超える確率を推定 二項分布を使う n回の試行で，失敗するか成功するかのどちらかが生起するような現象についての分布
二値データの場合，試行が１回だけの場合の二項分布を考える ベルヌーイ分布ともいう

79 1億を超える確率を推定 確率pをHRで予測 しかし，確率の値は０～１の範囲しかとらない 予測値が０～１の範囲に収まる必要がある
そこで，予測範囲を-∞～∞の範囲に拡張するために，ロジット変換をした確率を予測する logit(pi) = log(pi/(1-pi))= b0 + b1 Xi1 + b2 Xi2 ・・・ + ei これをロジットリンク関数と呼ぶ またベルヌーイ分布＋ロジットリンクのモデルを特にロジスティック回帰分析と呼ぶ

80 ロジット変換とは

81 ホームランで１億を超えるかを予測

82 Rで推定 glm()を使う >dat$ichioku <- ifelse(dat$salary<=100,0,1)
>fit.logistic <- glm(ichioku~HR,dat,family=binomial) >summary(fit.logistic) HR1本打つと，ロジット得点が0.156点上がる・・・ 解釈しにくい

83 オッズ比 HR１本打ったときに，１億を超える確率が何倍になるか オッズ比 = exp(回帰係数)
exp(fit.logistic$coefficients) HR1本打つと，1億を超える確率が1.17倍になる

84 二項分布

85 二項分布 N回中，r回成功する確率を推定 打率を直接予測しちゃダメ？ 比率，割合のデータ たとえば打率を推定したい場合は・・・
安打を成功，打数を試行数とするモデリングをする 打率を直接予測しちゃダメ？ 選手によって打数が異なる 点推定値は変わらなくても区間推定の結果が異なる より打数が多い選手の打率の推定精度は高くなり，打数が少ない選手は推定精度が低くなるはず 今回は，リーグで打率に差があるかを検討

86 基本はロジスティック回帰と同じ 確率は0～1の範囲をとる Rではglm()を使うが・・・ よって，ロジットリンクでモデリング
glmではbinomialを指定するとロジットになる Rではglm()を使うが・・・ 従属変数として成功数と失敗数を指定 >fit.bin <- glm(cbind(HIT,ATbats-HIT)~league-1,dat,family=binomial) >summary(fit.bin)

87 各リーグの打率の推定値 それぞれの係数をロジスティック変換 psychパッケージのlogistic()を使うと楽
>logistic(fit.bin$coefficients[1]) >logistic(fit.bin$coefficients[2])

88 ただし，この分析は後に述べるか分散の問題を含んでいるので真に受けてはいけない
打率に差はある？ 素直にリーグの回帰係数を見る >fit.bin2 <- glm(cbind(HIT,ATbats-HIT)~league,dat,family=binomial) >summary(fit.bin2) 一応，差はあるらしい・・・ ただし，この分析は後に述べるか分散の問題を含んでいるので真に受けてはいけない

89 ポアソン分布

90 ホームランのデータ 明らかに正規分布ではない

91 打席数で回帰 上手くフィットしていない・・・？ 打席数が増えるほど，HRの分散が大きくなっている・・・？

92 カウントデータの分析 ある事象が起きた回数のデータ ポアソン分布 特に，事象の起こる確率が引く場合 データは0を含む自然数
0側に歪んだデータ ポアソン分布 起こる確率が低い事象が起きた回数についての分布 ホームランなどは，それに該当

93 ポアソン分布の性質 値は非負の整数 平均と分散が等しい 平均が大きいほど，分散も大きい パラメータは平均値であるλ λ=3.5 λ=10

94 ポアソン回帰分析 データがポアソン分布に従うと仮定 説明変数がの上昇に従い，ポアソン分布のパラメータλがどれほど上昇するかを予測する
打席数が200のときの分布と，打席数が400のときの分布に違い 打席数が200増えると，ホームランは異なるポアソン分布に従う 平均が増えるほど分散が大きくなることをモデル化することができる

95 ポアソン回帰分析 ポアソン分布は非負 予測値が負になってしまうと，分布の範囲からはずれてしまう
そこで，ポアソン分布のパラメータであるλを対数変換した値を予測することで解決する Log(λi) = b0 + b1 Xi1 + b2 Xi2 ・・・ + ei これを対数リンク関数と呼ぶ

96 Rでポアソン回帰分析 同じくglm()を使う
>fit.pois <- glm(HR~ATbase,dat,family=poisson) >summary(fit.pois) 係数が小さく見えるが，これは対数リンクを使っているため z値が25以上・・・すげぇ有意！（に見える）

97 ポアソン回帰

98 予測区間を表示すると・・・ 離散分布なので予測分布がガタガタ はずれまくり！

99 ポアソン回帰の解釈 1打席あたりHRが出る程度は，その人の打席数に依存する 打席数が多い人ほど，1打席あたりにHRを打つ割合が高くなっていく
ただし，ポアソン分布には重要な欠点が・・ パラメータがλひとつのため，過分散の問題

100 負の二項分布

101 ポアソン分布の過分散問題 離散分布はパラメータが1つのものが多い 分散は自動的に決まってしまう 二項分布・・・確率pのみ
ポアソン分布・・・平均λのみ 分散は自動的に決まってしまう 二項分布・・・分散=Np(1-p) ポアソン分布・・・分散=λ 真の分布がより分散が大きい場合，予測が大きく外れてしまう・・・過分散問題

102 過分散の対処法 同じタイプの分布で2パラメータのものを使う 変量効果を用いて過分散を推定 二項分布→ベータ二項分布
ポアソン分布→負の二項分布 変量効果を用いて過分散を推定 明日お話しします

103 負の二項分布についての詳細 こちらのスライドをどうぞ 負の二項分布

104 負の二項分布 ポアソン分布の過分散も推定 HRをポアソンと負の二項分布で推定 全然ダメ！ ポアソン分布 上手くいきそう！ 負の二項分布

105 Rで負の二項回帰分析 HRを打席数で予測 MASSパッケージのglm.nbを使う z値が13.5・・・ ポアソンの半分
library(MASS) fit.nb <- glm.nb(HR~ATbase,dat) summary(fit.nb) 過分散パラメータ z値が13.5・・・　ポアソンの半分 もちろん負の二項分布のほうが正しい

106 負の二項回帰分析 きれいに収まった

107 ポアソンと負の二項分布 AICを比較 ポアソン分布の使いどころ ポアソン分布：AIC=1461 負の二項分布：AIC=1033
負の二項分布のほうが圧倒的に小さい！ ポアソン分布の使いどころ 社会科学のデータでポアソンが当てはまる例は少ない気がする 個人差がでかいので 常に負の二項分布との比較を忘れないようにすること！ ポアソン分布で過分散が生じるとかなり危険なTypeⅠエラー！

108 GLMと心理統計学

109 GLMってあんまり使わない 2値データにロジスティック回帰 ・・・は比較的よく使うが・・・
ポアソンや負の二項分布，ガンマ，対数正規分布などはあまり使われない なぜか

110 有意性検定だけなら・・ あまり結果は変わらない 同じ有意なら別に正規分布でいいじゃん
1パラメータの分布は別にして，ある分布で有意だが別の分布では非有意になる，ということはあまりない もちろん，推定しているパラメータが違うので同じにはならない 同じ有意なら別に正規分布でいいじゃん という判断も多少はあるような気がする

111 実験計画法との関係 ANOVA：実験計画法との相性が良い GLM：実験計画法と相性が悪い 効果の解釈が容易 効果量を簡単に計算できる
グラフを簡単に書ける GLM：実験計画法と相性が悪い 効果の解釈がしづらい 効果量がまだうまく計算できない グラフも書きづらい 実験室の行動を予測しても意味がない

112 GLMは心理統計では不要？ もちろん，そうではないと思う 理由0：正規分布以外のデータはごろごろある 理由1：GLMの効果量もあるにはある
分析法を理由にとるデータに制限をかけるのは変 理由1：GLMの効果量もあるにはある 今後，発展されていくべきところ 理由2：効果量の区間推定をちゃんとする モデリングを軽視すると，区間推定に大きなバイアス 差の信頼区間から，効果量の信頼区間へ

113 続きます 心理学者のためのGLMM・階層ベイズ