スライド 1

2016年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

今後の予定

Google map から転載 Location

第二回（今回） 線形モデル

第三回 一般化線形モデル

*ニュースの内容の信ぴょう性には言及しません (たとえば, 日本住血吸虫の場合, 経口感染でないことが実験で示されています。) 画像はYou tube, 【新唐人日本2011年5月5日付ニュース】より

第四回 ベイズ統計（導入）

ガンバ大阪ホームページより http://www2.gamba-osaka.net/stadium/ Google map から転載

第五回 ベイズファクター

復習もかねて

アンケート２番から

復習もかねて

, 相関分析の話を少し (阪大 学部１年 統計学 B-2/C-2 より)

学習したいこと （おもに統計３研究室以外）

複数のデータの関係性を見つける方法

を知りたい

多変量解析, グラフィカルモデリング,

因果推論 etc. （上級の話題）

標本相関係数

２変量データで因果関係が不明瞭な場合

二つの項目の関連を見るには？

英語の点数(X)と統計の点数(Y) (100点満点)

男女カップルにおける男の身長(X)と女の身長(Y)

例

→ まずは

散布図

を作って,

標本相関係数

を求める

(X, Y)

88, 90

45, 78

56, 100

(X, Y)

188, 190

145, 178

156, 200

....

標本相関係数

散布図

２変量データで 片方を横軸にもう片方を縦軸にとってプロットしたもの

標本相関係数

２変量データの線形的 な関係(y=ax+b といった関係)をはかる指標

2 2 y x xy xy

S

s

S

r

=

1

1

≤

≤

−

r

_xy

=

2

x

s

ただし, ここで

=

xy

S

散布図と標本相関係数の例

-4 -2 0 2 4 6 44 46 48 50 52 54 56 r= X Y 00357 . 0 − = xy r -6 -4 -2 0 2 4 6 44 46 48 50 52 54 56 58 r= X Y 658 . 0 = xy r -4 -2 0 2 4 6 45 50 55 r= X Y 792 . 0 − = xy r

考えてみよう

-1.0 -0.5 0.0 0.5 1.0 0 .0 0 .2 0 .4 0 .6 0 .8 1 .0 r= X Y ? 2 2 = = y x xy xy S s S r -1.0 -0.5 0.0 0.5 1.0 0 .0 0 .2 0 .4 0 .6 0 .8 1 .0 r= X Y

a. 0.109 b. 0.89 c. -0.009

d. -0.74 e. -0.37 f. 0.44

まとめ：標本相関係数

標本相関係数

２変量データの

をはかる指標

2 2 y x xy xy

S

s

S

r

=

1

1

≤

≤

−

r

_xy

といった関係（非線形という）をはかることは

できない！！

→ r の値のみでなく

必ず散布図を目で確認する

_{ことが大事}

b

ax

y

=

2

+

注意点 (統計学B-1/C-1 の内容)

多変量データの散布図は第三回で出てきます

復習もかねて

～分析前の作業

データの下ごしらえ

分析前に生データに手を加えることが多い

データの出典に関する情報

, 収集方法, 分野固有の専門知識

が必要

本講義のスタンス

データは所与とする

.

本来： 仮説立案とその検証を目的とする場合, データの収集段階で 統計の専門家が入って収集方法(データ数など)を検討すべき

分析前の作業

そのまま調理したらダメ！ （アク抜き必須）

考えてみよう：データの下ごしらえ

1. 実際に生データを分析したことがある人は、分析前に、データの 一部を削ったことがあるかもしれない。その理由を思い出してみる

２．分析経験がない人は、自由に想像してみること

分析前の作業

分析データの精錬作業

データの中身の確認

（数値チェック）

プロット, 図示

（視覚的に見ておかしくないか？）

加工・削除

（慎重に！）

回帰分析（B-2/C-2資料より）

O大学 新入生のみずほさんは賃貸情報を

ネットで検索. 以下のようなデータを得ました.

例題： みずほの部屋探し

→ 傾向をみるため, 横軸に距離, 縦軸に賃料をとりプロット(点を打つ)

最寄り駅からの距離 (徒歩)： 3 5 6 10 11 17

一カ月の賃料 (万円)：

8 7.3 6.2 4 4.2 3.5

豊中キャンパス近くの賃貸物件（１K）

データのプロット

0 5 10 15 20 0 2 4 6 8 10 Kaiki Min Walk 10^ 4 Y E N x <- c(3, 5,6, 10, 11, 17); y <- c(8, 7.3, 6.2, 4, 4.2, 3.5);

plot(x,y, pch=18, col=2, xlim=c(0, 20), ylim=c(0, 10), main="Kaiki", xlab="Min Walk", ylab="10^4 YEN"); abline(h=0, lty=2, col="gray"); # hori line

abline(v=0, lty=2, col="gray"); # vert line

6

,

,

2

,

1 

=

i

x: 最寄駅からの距離(分：徒歩換算), y:一か月の家賃 (万円)

線形モデルの導入

)

,

0

(

~

σ

2

ε

_i

N

i

i

i

x

Y

=

α

+

β

+

ε

モデルのパラメータ

,

α

,

β

;

σ

2

パラメータを最尤推定法によって推定し以下の直線を引いてみる

x

Y

=

α

ˆ

+

β

ˆ

線形モデル （説明変数・目的変数が各々１変数

; 線形回帰モデル）

推定したパラメータでの直線

0 5 10 15 20 0 2 4 6 8 10 Kaiki Min Walk 10^ 4 Y E N x <- c(3, 5,6, 10, 11, 17); y <- c(8, 7.3, 6.2, 4, 4.2, 3.5); res <- lm(y~x); ahat <- res$coefficients[1]; bhat <- res$coefficients[2]; R プログラム例 （回帰分析）

plot(x,y, pch=18, col=2, xlim=c(0, 20), ylim=c(0, 10), main="Kaiki", xlab="Min Walk",

ylab="10^4 YEN");

abline(h=0, lty=2, col="gray"); # hori line abline(v=0, lty=2, col="gray"); # vert line abline(a=ahat, b=bhat);

考えてみよう

最小二乗法

0 5 10 15 20 0 2 4 6 8 10 Kaiki Min Walk 10^ 4 Y E N

下図は最小二乗法を用いて求めた直線を引いたもの.

最尤推定法に基づくものと同じだろうか, 違うだろうか

{

}

∑

=

+

−

6 1 2 ,

(

)

min

i i i

x

y

α

β

β α

考えてみよう２【重要！】

最小二乗法

0 5 10 15 20 0 2 4 6 8 10 Kaiki Min Walk 10^ 4 Y E N

{

}

∑

=

+

−

6 1 2 ,

(

)

min

i i i

x

y

α

β

β α

最尤推定法

)

,

0

(

~

σ

2

ε

_i

N

i i i

x

Y

=

α

+

β

+

ε

統計モデリングの立場

部分に確率分布を導入 (仮定)

ε

+

=

f

(x

)

Y

（真の）関係式

_{確率的な変動}

ε

統計モデリングの概念的な式

データが少ないことによる効果（統計誤差）や実験で不可避の

測定誤差を踏まえた議論が可能になる！

注意すべきこと

１． の項は加法的とは限らない

２． 正しい関係式があるとは限らないし, 理解しやすい近似式の方が有

効なこともある

３． 独立で同一な正規分布を使う理由は計算上の都合

（正規分布以外を仮定, 異分散や相関をもたせたモデル＝上級の話題）

４． 部分に「正しい確率分布」はないと考えてよい.

ε

+

=

f

(x

)

Y

統計モデリングの

概念的な

_式

ε

ε

Birthweight vs Gestational Age

(余計なものは取り除いてある)

胎内にいた期間

[週], 出生時の体重 [g], 男児(b)/女児(g)

男児

, 女児 ともに標本サイズは １２ずつ

データ例

1 40 2968 b 2 38 2795 b 3 40 3163 b 4 35 2925 b 5 36 2625 b 6 37 2847 b 7 41 3292 b 8 40 3473 b 9 37 2628 b 10 38 3176 b 11 40 3421 b 12 38 2975 b 13 40 3317 g 14 36 2729 g 15 40 2935 g 16 38 2754 g 17 42 3210 g 18 39 2817 g 19 40 3126 g 20 37 2539 g 21 36 2412 g 22 38 2991 g 23 39 2875 g 24 40 3231 g

35 36 37 38 39 40 41 42 2400 2600 2800 3000 3200 3400 Chap. 2 Age W ei ght b g

定量的な確認：

たとえば相関係数の計算

ρ

=

0

.

744

データのプロット

見てわかること

24

,

,

2

,

1 

=

i

x: 胎内にいた期間(週), y:出生時の体重

i

i

i

x

Y

E

[

]

=

µ

=

α

+

β

)

,

0

(

~

σ

2

ε

_i

N

線形モデルの導入

i

i

i

x

Y

=

α

+

β

+

ε

モデルのパラメータ

,

α

,

β

;

σ

2

パラメータを最尤推定法によって推定し以下の直線を引いてみる

x

Y

=

α

ˆ

+

β

ˆ

線形モデル （説明変数・目的変数が各々１変数）

回帰直線を引いてみる

データ

> dat2

AGE WEI TYPE 1 40 2968 b 2 38 2795 b ・・・・ 23 39 2875 g 24 40 3231 g 線形回帰

> dat2.res <- lm(WEI~ AGE, data= dat2);

35 36 37 38 39 40 41 42 2400 2600 2800 3000 3200 3400 Chap. 2 Age W ei ght b g 回帰直線

x

y

=

α

ˆ

+

β

ˆ

α

ˆ

=

−

1484

,

β

ˆ

=

115

.

5

R プログラム例

データを男児(j=1), 女児(j=2)に分けて分析 線形モデル

)

,

0

(

~

σ

2

ε

_ji

N

ji ji j j ji

x

Y

=

α

+

β

+

ε

回帰直線

x

y

=

α

ˆ

_j

+

β

ˆ

_j

112

ˆ

,

1269

ˆ

₁

=

−

β

₁

=

α

35 36 37 38 39 40 41 42 2400 2600 2800 3000 3200 3400 Chap. 2 Age W ei ght b g

130

ˆ

,

2142

ˆ

₂

=

−

β

₂

=

α

線形回帰はよさそうだが

, 男児と女児で分けて考えた方がいいのか？

分析の課題

男児・女児に分けて推定すると

２種類の線形モデルで仮説検定（一般的な形） H1: 傾きが異なる線形モデル

)

,

0

(

~

σ

2

ε

_jk

N

jk jk j j jk

x

Y

=

α

(1)

+

β

+

ε

仮説検定によるモデル選択 (1/2)

H０: 傾きは等しい線形モデル jk jk j jk

x

Y

=

α

(0)

+

β

+

ε

ε

_~

₍

₀

_,

σ

2

₎

N

jk 検定統計量：

∑

−

−

=

k j jk j jk

x

Y

S

, 2 ) 0 ( 0

:

α

ˆ

β

ˆ

=

∑

−

−

k j jk j j jk

x

Y

S

, 2 ) 1 ( 1

:

α

ˆ

β

ˆ

) 2 ( , 1 1 1 0

~

1

)

2

(

− −

−

−

−

K J J

F

J

K

J

S

S

S

12 ; 2 = = K J K k J j =1,..., ; =1,...,

２種類の線形モデルで仮説検定 H1: 傾きが異なる線形モデル H０: 傾きは等しい線形モデル VS (有意水準 0.05)仮説H0が正しいとすると, 検定統計量の値は0から4.35程度に おさまるはず → 実際に（がんばって）計算すると・・・

35

.

4

19

.

0

1

)

2

(

1 1 0

=

<<

−

−

−

J

K

J

S

S

S

12 ; 2 = = K J

（

H0 は棄却されないため）傾きが異なるとは言えない

仮説検定によるモデル選択 (2/2)

∑

−

−

=

k j jk j jk

x

Y

S

, 2 ) 0 ( 0

:

α

ˆ

β

ˆ

=

∑

−

−

k j jk j j jk

x

Y

S

, 2 ) 1 ( 1

:

α

ˆ

β

ˆ

・途中の計算がとても大変・・・・ ・初学者にはかなりわかりにくい ・使用上の制約が多い

p

L

AIC

=

−

2

(

θ

ˆ

)

+

2

一般のモデル選択

（ネイマン・ピアソン流の）仮説検定

赤池情報量規準

(Akaike Information Criterion)

を用いたモデル選択

尤度関数の最大値 (最尤推定は尤度関数を最大化)とパ ラメータ数からAICを計算; 相対的にAICの小さいモデルを選ぶ p パラメータ数 ) ( max ) ˆ (

θ

θ

θ L L = ) (

θ

L 尤度関数

θ

ˆ 最尤推定での推定値

本講義は統計モデルの紹介が主

, AICなどは機械的に使用してよい.