Transcript 回帰分析
第5回 回帰モデル •回帰モデル –単回帰モデル、重回帰モデル –非線形式の回帰分析 •最小二乗法（記述統計） –回帰式の適合度のチェック（分散分析） •最小二乗推定量の性質（推測統計） –パラメータ推定量の不偏性 –パラメータ推定量の確率分布とｔ検定 1 回帰モデル • いくつかの変数間に相関関係が存在 • ある変数の値を、別の変数を用いて説明 従属変数、目的変数 説明式を作成 独立変数、説明変数 x 被説明変数 変数Y,実現値yi 推計値yi=f(xi) 変数X,実現値xi i Y Y yi yˆi Yˆi f ( X i ) Yˆi f ( X i , Zi ) Z xi X X 2 記述統計学と推測統計学 記述統計学 多数データの 数学的要約 ･記述 母集団の データ 推測統計学 (仮想的) 母集団 「標本統計量」 無作為 抽出 標本集団 のデータ 少数データの 数学的要約 ･記述 確率的推測･記述 ｢母集団の確率分布･母数｣ ｢標本統計量の確率分布」 3 目的変数を[記述]する 従属変数、目的変数 説明式を作成 独立変数、説明変数 被説明変数 変数Y,実現値yi 推計値yi=f(xi) 変数X,実現値xi Y yi yˆi 記述統計学の立場 関数形f（X)をうまく決めて、 Yˆi f ( X i ) Yをうまく記述したい （説明できなかった部分の） 残差･誤差 残差を小さくしたい i i i (２乗和で評価するので) 残差の二乗和を小さくしたい X y yˆ xi 最小二乗法 min i2 ( yi yˆi )2 ( yi f ( X i ))2 i i i 4 最小二乗法（単回帰モデル） • 線形モデル yi a bxi n n 2 2 S ( y a bx ) • 残差平方和 E i i i i 1 – 未知数で微分すると SE 2( yi a bxi ) 0 a i 1 SE 2( yi a bxi ) xi 0 b • 正規方程式 上式を整理して y a n b x 0 x y a x b x 0 i i i i i 2 i – それらを解いて 2 x y x xy ˆa i 2i i 2 i i n xi ( xi ) n xi yi xi yi ˆ b n xi 2 ( xi )2 5 最小二乗法（重回帰モデル） yi a bwi czi • 線形モデル • 残差平方和 – n n SE i ( yi a bwi czi )2 2 i 1 未知数で微分すると i 1 SE a 2( yi a bwi czi ) 0 SE b 2( yi a bwi czi )wi 0 SE c 2( yi a bwi czi ) zi 0 • 正規方程式 上式を整理して y a n b w c z 0 w y a w b w c w z 0 z y a z b z w c z 0 i i i 2 i i i i i i i i i i i 2 i – 単回帰モデルと同様の連立一次方程式が出来る 6 最小二乗法（重回帰モデル）行列表示 yˆi a bwi czi • 線形モデル yˆ X y X • 残差 n 2 • 残差平方和 SE i t ( y X )t ( y X ) i 1 – 未知数で微分すると 2X ( y Xˆ ) 0 t • 正規方程式 t t t 1 t ˆ ˆ X y X X より、 ( X X ) ( X y) 7 行列を用いた重回帰式の計算例 lm(formula = 強さ ~ 薬剤A + 薬剤B) Coefficients: (Intercept) 薬剤A 薬剤B 39.728 6.433 1.750 8 Rによる計算例 Call: #データフレームを作成する lm(formula = 強さ ~ 薬剤A) concrete <- data.frame(Call: + 強度 =c(40.0,45.8,53.0,42.3,46.7,54.1,44.2,47.1,58.0), lm(formula Residuals: = 強さ ~ 薬剤A + 薬剤B) + 混和A =c(0,1,2,0,1,2,0,1,2), Min 1Q Median 3Q Max Residuals: -2.1111 -1.3444 -0.8111 0.8222 3.6556 + 混和B =c(0,0,0,1,1,1,2,2,2) Min 1Q Median 3Q Max ) -2.5611 -0.3611 0.2722 0.8222 1.9056 Coefficients: #データフレームから変数に代入(付置)する Estimate Std. Error t value Pr(>|t|) 強さ <-concrete$強度 Coefficients: (Intercept) 41.478 1.128 36.761 2.86e-09 *** Estimate Std. Error 薬剤A <- concrete$混和A 薬剤A 6.433 0.874 t value 7.361 Pr(>|t|) 0.000154 *** (Intercept) 39.7278 1.0076 39.426 1.78e-08 *** --薬剤B <- concrete$混和B 薬剤A 0.6171 4.56e-05 Signif. codes:6.4333 0 ‘***’ 0.001 ‘**’ 10.426 0.01 ‘*’ 0.05 ‘.’*** 0.1 ‘ ’ #単回帰分析の結果をkekka1に入れ，その概要を出力 薬剤B 1.7500 0.6171 2.836 0.0297 * 1 kekka1 <- lm(強さ~薬剤A) --summary(kekka1) Signif. ‘***’ 0.001 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ Residualcodes: standard0 error: 2.141‘**’ on 0.01 7 degrees of freedom 1Multiple R-squared: 0.8856, #重回帰分析の結果をkekka2に入れ，その概要を出力 Adjusted R-squared: 0.8692 F-statistic: 54.18 on 1 and 7 DF, p-value: 0.0001545 kekka2 <- lm(強さ~薬剤A+薬剤B) Residual standard error: 1.511 on 6 degrees of freedom summary(kekka2) Multiple R-squared: 0.9511, Adjusted R-squared: 0.9348 9 F-statistic: 58.37 on 2 and 6 DF, p-value: 0.0001168 重共線性の問題 • 多数の説明変数の間に相関がある場合 – 目的変数への効果を一意に分離できない – 係数の推計値が安定しない（直感に反する符 号を取るなど） すべての観測値が ほぼ一直線上にある Y この直線を含むような平面で あれば、どの式を使っても当て はまりにはほとんど差はない Z X 直線上にない場所のYの予測 値には大きな差がでる 10 目的変数をどの程度[記述]出来たか？ Xiによる説明式がない場合 yiの推計値 Y として、 平均値 y を y 使うしかない Yˆi f ( X i ) なる説明式がある場合 回帰（Xiで説明 できた y から のずれ） Y y yˆi yˆi y 残差･誤差 xi 回帰平方和 平均値周りの バラツキ （全平方和） ST ( yi y ) i X SR ( yˆi y)2 i 2 残差平方和 SE i2 ( yi yˆi )2 i i i yi yˆi 決定係数 R2 SR S 1 E ST ST 説明できた 平方和の割合 11 決定係数 平均値周りの バラツキ （全平方和） ST ( yi y )2 i 回帰平方和 SR ( yˆi y)2 i 残差平方和 SE i2 ( yi yˆi )2 i i S S 決定係数 R squared R2 R 1 E 説明変数を増やせ ST ST 説明できた平方和の割合 ば1に近づく 本当はあまり意味のない変数を増やしてしまう危険性あり SE /(n k ) n 1 ~2 自由度調整済み決定係数 R 1 1 (1 R2 ) ST /(n 1) nk Adjusted R squared Residual standard error: 1.511 on 6 degrees of freedom Multiple R-squared: 0.9511, Adjusted R-squared: 0.9348 12 回帰による記述（説明）力の検定 – 回帰平方和が統計的に大きな意味を持っている か？ – 分散分析表を作り、F検定を行って判断する。 • 帰無仮説：回帰平方和は誤差平方和と同程度の大きさ （回帰式は、誤差に比べて大きな説明力はない） • 対立仮説：回帰平方和は誤差平方和より大きい （回帰式によって誤差よりもかなり大きい部分が説明でき た） 13 回帰による記述（説明）力の検定例 • 帰無仮説：回帰平方和は誤差平方和と同程度の大きさ （回帰式は、誤差に比べて大きな説明力はない）→棄却 • 対立仮説：回帰平方和は誤差平方和より大きい （回帰式により、誤差よりもかなり大きい部分が説明できた） Multiple R-squared: 0.9511, Adjusted R-squared: 0.9348 F-statistic: 58.37 on 2 and 6 DF, p-value: 0.0001168 14 推測統計学の導入 誤差(残差)が確率分布 (正規分布)をもつと仮定 i yi f ( Xi ) ~ N(0, ) 2 もともとの観測値は、真の値を中心に確率的に ばらついていたと考える yi f ( X i ) i Y Y yi yˆi Yˆi f ( X i ) 残差･誤差 i yi yˆi xi X X 関数形 f (X ) を決めると、誤差の実現値 1 ,, n が決まる。 それが正規分布のもとで発生する確率をチェックできる。 15 最尤推定法としての最小二乗法 誤差(残差)が確率分布 (正規分布)をもつと仮定 i yi f ( Xi ) ~ N(0, ) 2 残差の実現値の同時発生確率（正規分布確率密 度の積）が最大になるように、関数 f (X ) を決める。 1 ( / ) L e Y 2 i 1,n 2 2 i n n L log L log(2 ) n log i 2 / 2 2 i 1 この最大化のために最終項の残差二乗和 n n 2 SE i ( yi f ( X i ))2 を最小化するように、 X i 1 i 1 関数 f (X ) を定める必要がある （記述統計学で使われてきた）最小二乗法は、 （推測統計学の考え方に従う）最尤推定法でもある。 16 AIC 赤池の情報量基準 • AIC =－２（最大対数尤度)＋２パラメータ数 • パラメータ数＝ｋ＋1(回帰係数と残差の分散) • 最大対数尤度＝ n log(2ˆ 2 ) n 2 2 であるから， AIC n log(2ˆ 2 ) n 2(k 1) この値が小さいほうが,無駄の ない効率的なモデル式ができ ていることを意味している. > AIC(kekka1) [1] 42.98059 > AIC(kekka2) [1] 37.32718 17 パラメータ推定量の確率分布 • サンプルの観測値は、真の回帰式から誤差の分 だけ、上下にずれている。 • 誤差を含む観測値から計算したパラメータの推 定値も、真の値からずれている。 Y Y 切片の 推定値 真の 切片は０ 傾きの推定値 真の 傾きは０ X X 各サンプルの誤差が平均0分散 2 の正規分布に従うこと を用い、パラメータの推定量 ˆk が従う分布を計算すれば、 18 真の値からかけ離れている確率をチェックできる。 最小二乗推定量の確率分布 • 最小二乗推定量の行列表示 ˆ ( X t X )1 X t y y X • 観測値の誤差構造 ˆ ( X t X )1 X t ( X ) ( X t X )1 ( X t X ) ( X t X )1 X t ( X t X )1 X t – ここで、 ( X t X )1 X t は確率変数ではない定数。 – これらより、 • 不偏性 – 誤差は正と負を同程度にとるから、E[ε]=0 E[ˆ ] ( X t X )1 X t E[ ] ( X t X )1 X t 0 – 多数のサンプルを取りながら繰り返して推定す ると、その平均値は真の値に一致 19 最小二乗推定量の確率分布 （２） • 最小二乗推定量の確率分布 ˆ ( X t X )1 X t • 観測値の誤差が、分散 2 の正規分布に従うなら、 それに定数項をかけて足したものは正規分布に従う。 – パラメータの分散共分散行列は、 E[(ˆ )(ˆ )t ] E[( X t X )1 X tt X ( X t X )1 ] 2 ( X t X )1 X t X ( X t X )1 2 ( X t X )1 2 • j番目のパラメータの分散： V jj j番目の対角要素 • 誤差の母分散 2 はわからないから、残差平方和 2 から求めた分散で代用 ˆ ( y y ) S ˆ 2 E i i nk nk 20 • 最小二乗推定量の確率分布 （２） 正規分布に従う変数の「ずれ」の評価 – 母分散がわからないので、 サンプル（残差平方和）からの 推定値を用いる→t分布 • t統計量の作成 2 ˆ ( y y ) S ˆ 2 E i i nk nk tj ˆ j j – 自由度n-kのt分布に従う ˆ 2V jj – 通常、帰無仮説として j 0 （ｊ番目の要因の影響はない)とおけば、 ˆ 2V jj を用いて検定できる。 検定統計量 t j ˆ j Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 39.7278 1.0076 39.426 1.78e-08 *** 薬剤A 6.4333 0.6171 10.426 4.56e-05 *** 薬剤B 1.7500 0.6171 2.836 0.0297 * 21 回帰式のチェック plot(強さ,薬剤A) plot(kekka2$fitted.values,強さ,asp=1) あるいは， plot(kekka2)によって,4種類の結果診断グラ フが出力される． 22 AICに基づくステップワイズ法に よる変数選択 • すべての変数を含むモデルをlm()で推定 • その結果が入ったオブジェクトに，step()関 数を適用する > step(kekka2) Start: AIC=9.79 強さ ~ 薬剤A + 薬剤B Df Sum of Sq RSS AIC <none> 13.707 9.786 - 薬剤B 1 18.375 32.082 15.440 - 薬剤A 1 248.327 262.034 34.341 Call: lm(formula = 強さ ~ 薬剤A + 薬剤B) Coefficients: (Intercept) 薬剤A 薬剤B 39.728 6.433 1.750 23 非線形式の回帰分析 • 非線形の関係式でも、置き換えによって線 形化できる場合がある y 1e2 x log y log 1 2 x y 1 2 x 3 x 2 z 0 2 x y 1 2 x 3w 24