担当 : 田中冬彦 016 年 4 月 19 日 @ 統計モデリング 統計モデリング 第二回配布資料 文献 : A. J. Dobson and A. G. Barnett: An Introduction to Generalized Linear Models 3rd ed., CRC Press. 配布資料の PDF は以下からも DL できます. 短縮 URL http://tinyurl.com/lxb7kb8
第二回 ( 今回 ) 線形モデル 今後の予定 第三回一般化線形モデル Location 第四回ベイズ統計 ( 導入 ) Google map から転載 画像は You tube, 新唐人日本 011 年 5 月 5 日付ニュース より * ニュースの内容の信ぴょう性には言及しません ( たとえば, 日本住血吸虫の場合, 経口感染でないことが実験で示されています ) ガンバ大阪ホームページより http://www.gamba-osaka.net/stadium/ 第五回ベイズファクター Google map から転載
復習もかねて ~ 相関分析の初歩
アンケート 番から 学習したいこと ( おもに統計 3 研究室以外 ) 複数のデータの関係性を見つける方法を知りたい 多変量解析, グラフィカルモデリング, 因果推論 etc. ( 上級の話題 ) 復習もかねて, 相関分析の話を少し ( 阪大学部 1 年統計学 B-/C- より )
標本相関係数 変量データで因果関係が不明瞭な場合 二つの項目の関連を見るには? 例 英語の点数 (X) と統計の点数 (Y) (100 点満点 ) 男女カップルにおける男の身長 (X) と女の身長 (Y) (X, Y) 88, 90 45, 78 56, 100 (X, Y) 188, 190 145, 178 156, 00... まずは散布図を作って, 標本相関係数を求める
標本相関係数 散布図 変量データで片方を横軸にもう片方を縦軸にとってプロットしたもの 標本相関係数 変量データの線形的な関係 (y=ax+b といった関係 ) をはかる指標 r = xy ただし, ここで s S x xy S y 1 r 1 xy s x = S xy =
散布図と標本相関係数の例 r= r= r= Y 44 46 48 50 5 54 56 r xy = 0.00357 Y 44 46 48 50 5 54 56 58 r xy = 0.658 Y 45 50 55 r xy = 0.79-4 - 0 4 6 X -6-4 - 0 4 6 X -4-0 4 6 X
考えてみよう r xy = s S x xy S y =? r= r= Y 0.0 0. 0.4 0.6 0.8 1.0 Y 0.0 0. 0.4 0.6 0.8 1.0-1.0-0.5 0.0 0.5 1.0 X -1.0-0.5 0.0 0.5 1.0 X a. 0.109 b. 0.89 c. -0.009 d. -0.74 e. -0.37 f. 0.44
まとめ : 標本相関係数 標本相関係数 変量データの をはかる指標 r = xy s S x xy S y 1 r 1 xy 注意点 ( 統計学 B-1/C-1 の内容 ) y = ax + b できない!! といった関係 ( 非線形という ) をはかることは r の値のみでなく必ず散布図を目で確認することが大事 多変量データの散布図は第三回で出てきます
復習もかねて ~ 分析前の作業
本講義のスタンスデータは所与とする. データの下ごしらえ 本来 : 仮説立案とその検証を目的とする場合, データの収集段階で統計の専門家が入って収集方法 ( データ数など ) を検討すべき 分析前の作業 分析前に生データに手を加えることが多い データの出典に関する情報, 収集方法, 分野固有の専門知識 が必要 そのまま調理したらダメ! ( アク抜き必須 )
考えてみよう : データの下ごしらえ 1. 実際に生データを分析したことがある人は 分析前に データの一部を削ったことがあるかもしれない その理由を思い出してみる. 分析経験がない人は 自由に想像してみること 正解を求めているわけではない!
分析前の作業 分析データの精錬作業 データの中身の確認 ( 数値チェック ) プロット, 図示 ( 視覚的に見ておかしくないか?) 加工 削除 ( 慎重に!)
1. 線形モデル
1.1 単回帰と統計モデル
回帰分析 (B-/C- 資料より ) 例題 : みずほの部屋探し O 大学新入生のみずほさんは賃貸情報をネットで検索. 以下のようなデータを得ました. 豊中キャンパス近くの賃貸物件 (1K) 最寄り駅からの距離 ( 徒歩 ): 3 5 6 10 11 17 一カ月の賃料 ( 万円 ): 8 7.3 6. 4 4. 3.5 傾向をみるため, 横軸に距離, 縦軸に賃料をとりプロット ( 点を打つ )
データのプロット Kaiki 10^4 YEN 0 4 6 8 10 R プログラム例 x <- c(3, 5,6, 10, 11, 17); y <- c(8, 7.3, 6., 4, 4., 3.5); plot(x,y, pch=18, col=, xlim=c(0, 0), ylim=c(0, 10), main="kaiki", xlab="min Walk", ylab="10^4 YEN"); abline(h=0, lty=, col="gray"); # hori line abline(v=0, lty=, col="gray"); # vert line 0 5 10 15 0 Min Walk
線形モデルの導入 線形モデル ( 説明変数 目的変数が各々 1 変数 ; 線形回帰モデル ) Y = α + βx + ε i i i i =1,,,6 ε ~ N(0, σ ) i x: 最寄駅からの距離 ( 分 : 徒歩換算 ), y: 一か月の家賃 ( 万円 ) モデルのパラメータ, α, β; σ パラメータを最尤推定法によって推定し以下の直線を引いてみる Y = αˆ + βx ˆ
推定したパラメータでの直線 R プログラム例 ( 回帰分析 ) x <- c(3, 5,6, 10, 11, 17); y <- c(8, 7.3, 6., 4, 4., 3.5); res <- lm(y~x); ahat <- res$coefficients[1]; bhat <- res$coefficients[]; R プログラム例 ( 回帰直線 ) plot(x,y, pch=18, col=, xlim=c(0, 0), ylim=c(0, 10), main="kaiki", xlab="min Walk", ylab="10^4 YEN"); abline(h=0, lty=, col="gray"); # hori line abline(v=0, lty=, col="gray"); # vert line abline(a=ahat, b=bhat); 10^4 YEN 0 4 6 8 10 Kaiki 0 5 10 15 0 Min Walk
考えてみよう 下図は最小二乗法を用いて求めた直線を引いたもの. 最尤推定法に基づくものと同じだろうか, 違うだろうか Kaiki 最小二乗法 min α, β 6 { yi ( α + βx i )} i= 1 10^4 YEN 0 4 6 8 10 0 5 10 15 0 Min Walk
考えてみよう 重要! 最小二乗法 最尤推定法 min α, β 6 { yi ( α + βx i )} i= 1 Y = α + β + ε i x i i ε ~ i N(0, σ ) Kaiki 10^4 YEN 0 4 6 8 10 0 5 10 15 0 Min Walk
統計モデリングの概念的な式 統計モデリングの立場 Y = f (x) + ε ( 真の ) 関係式確率的な変動 ε 部分に確率分布を導入 ( 仮定 ) データが少ないことによる効果 ( 統計誤差 ) や実験で不可避の測定誤差を踏まえた議論が可能になる!
統計モデリングの概念的な式 注意すべきこと Y = f (x) + ε ε 1. の項は加法的とは限らない. 正しい関係式があるとは限らないし, 理解しやすい近似式の方が有効なこともある 3. 独立で同一な正規分布を使う理由は計算上の都合 ( 正規分布以外を仮定, 異分散や相関をもたせたモデル = 上級の話題 ) ε 4. 部分に 正しい確率分布 はないと考えてよい.
1. 複数の線形モデル
Birthweight vs Gestational Age データ例 ( 余計なものは取り除いてある ) 胎内にいた期間 [ 週 ], 出生時の体重 [g], 男児 (b)/ 女児 (g) 男児, 女児ともに標本サイズは 1 ずつ 1 40 968 b 38 795 b 3 40 3163 b 4 35 95 b 5 36 65 b 6 37 847 b 7 41 39 b 8 40 3473 b 9 37 68 b 10 38 3176 b 11 40 341 b 1 38 975 b 13 40 3317 g 14 36 79 g 15 40 935 g 16 38 754 g 17 4 310 g 18 39 817 g 19 40 316 g 0 37 539 g 1 36 41 g 38 991 g 3 39 875 g 4 40 331 g
データのプロット 見てわかること Weight 400 600 800 3000 300 3400 b g Chap. 35 36 37 38 39 40 41 4 定量的な確認 : ρ = 0. 744 たとえば相関係数の計算 Age
線形モデルの導入 線形モデル ( 説明変数 目的変数が各々 1 変数 ) E[ Y i ] = µ = α + βx i Y = α + βx + ε i i i i i =1,,,4 ε ~ N(0, σ ) i モデルのパラメータ, x: 胎内にいた期間 ( 週 ), y: 出生時の体重 α, β; σ パラメータを最尤推定法によって推定し以下の直線を引いてみる Y = αˆ + βx ˆ
回帰直線を引いてみる Chap. R プログラム例 データ > dat AGE WEI TYPE 1 40 968 b 38 795 b 3 39 875 g 4 40 331 g 線形回帰 Weight 400 600 800 3000 300 3400 b g > dat.res <- lm(wei~ AGE, data= dat); 35 36 37 38 39 40 41 4 回帰直線 y = ˆ α + βx ˆ ˆ α = 1484, ˆ β = 115. 5 Age
男児 女児に分けて推定すると データを男児 (j=1), 女児 (j=) に分けて分析 Chap. 線形モデル Y ji = α + β x + ε j 回帰直線 y = ˆ α + ˆ β x j j j ji α = 169, ˆ β ˆ1 1 = = 14, ˆ β ˆ = α ji 11 130 ε ji ~ N(0, σ ) Weight 400 600 800 3000 300 3400 b g 35 36 37 38 39 40 41 4 Age 分析の課題 線形回帰はよさそうだが, 男児と女児で分けて考えた方がいいのか?
仮説検定によるモデル選択 (1/) 種類の線形モデルで仮説検定 ( 一般的な形 ) j =1,..., J; k = 1,..., K H0: 傾きは等しい線形モデル (0) Yjk = j + βx jk α + ε jk ε jk ~ N(0, σ ) H1: 傾きが異なる線形モデル (1) Yjk = j + β j α x + ε jk jk ε jk ~ N(0, σ ) 検定統計量 : J S S S 0 1 ~ 1 = ; K = 1 J ( K ) J 1 F J 1, J ( K ) S 0 (0) (1) : = Y jk ˆ α j ˆ βx S = ˆ 1 α j j, k jk : Y ˆ jk β x j, k j jk
仮説検定によるモデル選択 (/) 種類の線形モデルで仮説検定 H0: 傾きは等しい線形モデル VS H1: 傾きが異なる線形モデル ( 有意水準 0.05) 仮説 H0 が正しいとすると, 検定統計量の値は 0 から 4.35 程度におさまるはず 実際に ( がんばって ) 計算すると S0 S1 J ( K ) = 0.19 << S J 1 1 4.35 S 0 (0) (1) : = Y jk ˆ α j ˆ βx S = ˆ 1 α j j, k jk : Y ˆ jk β x j, k j jk J = ; K = 1 (H0 は棄却されないため ) 傾きが異なるとは言えない
一般のモデル選択 ( ネイマン ピアソン流の ) 仮説検定 途中の計算がとても大変 初学者にはかなりわかりにくい 使用上の制約が多い 赤池情報量規準 (Akaike Information Criterion) を用いたモデル選択 AIC = L( θ ˆ) + p L (θ ) θˆ 尤度関数 最尤推定での推定値 p パラメータ数 尤度関数の最大値 L ( ˆ) θ = max L( θ ) ( 最尤推定は尤度関数を最大化 ) とパ θ ラメータ数からAICを計算 ; 相対的にAICの小さいモデルを選ぶ 本講義は統計モデルの紹介が主, AIC などは機械的に使用してよい.