第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  xi 記述統計学の立場関数形f（X)をうまく決めて、 Yˆi  f ( X i ) Yをうまく記述したい（説明できなかった部分の）残差･誤差残差を小さくしたい ˆ i  yi  yi (２乗和で評価するので) 残差の二乗和を小さくしたい X 最小二乗法 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 – 未知数で微分すると S E  2 ( yi  a  bxi )  0 a i 1 S E  2 ( yi  a  bxi ) xi  0 b • 正規方程式上式を整理して y i  a n  b xi  0  xi yi  a xi  b xi  0 2 – それらを解いて aˆ   xi  y x x y n x  ( x ) 2 i 2 i i i 2 i i n  xi yi   xi  yi ˆ b 2 n xi  ( xi ) 2 5 最小二乗法（重回帰モデル） yi  a  bwi  czi • 線形モデル • 残差平方和 – n n S E    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 • 残差平方和 S E    i 2   t  ( y  X )t ( y X ) i 1 – 未知数で微分すると  2 X ( y  Xˆ )  0 • 正規方程式 t t t 1 t ˆ ˆ X y  X X より、   ( X X ) ( X y) t 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 yˆi 残差･誤差 xi 回帰平方和平均値周りのバラツキ（全平方和） ST   ( yi  y ) i X S R   ( yˆi  y ) 2 i 2 残差平方和 S E   i2   ( yi  yˆi ) 2 i i  i  yi  yˆi 決定係数 R2  SR S  1 E ST ST 説明できた平方和の割合 11 決定係数平均値周りのバラツキ（全平方和） ST   ( yi  y ) 2 i 回帰平方和 S R   ( yˆi  y ) 2 i 残差平方和 S E   i2   ( yi  yˆi ) 2 i i S S 決定係数 R squared 説明変数を増やせ R2  R  1 E ST ST 説明できた平方和の割合ば1に近づく本当はあまり意味のない変数を増やしてしまう危険性あり S E /(n  k ) n 1 ~2 自由度調整済み決定係数 R  1  1 (1  R 2 ) ST /(n  1) nk 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 ( X i ) ~ 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 ( X i ) ~ N(0, ) 2 残差の実現値の同時発生確率（正規分布確率密度の積）が最大になるように、関数 f (X ) を決める。 Y L 1 e ( i /  ) 2 2 n n 2 L  log L   log(2 )  n log    i /  2 2 i 1 i 1, n 2 この最大化のために最終項の残差二乗和 n n 2 S E    i  ( yi  f ( X i ))2 を最小化するように、 X i 1 i 1 関数 f (X ) を定める必要がある（記述統計学で使われてきた）最小二乗法は、（推測統計学の考え方に従う）最尤推定法でもある。 16 AIC 赤池の情報量基準 • AIC =－２（最大対数尤度)＋２パラメータ数 • パラメータ数＝ｋ＋1(回帰係数と残差の分散) • 最大対数尤度＝  n log( 2ˆ 2 )  n 2 2 であるから， 2 ˆ AIC  n log(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 t  t X ( X t X ) 1 ]   2 ( X t X ) 1 X t X ( X t X ) 1   2 ( X t X ) 1 • j番目のパラメータの分散：  V jj j番目の対角要素 • 誤差の母分散 2 はわからないから、残差平方和 2 から求めた分散で代用 ˆ ( y  y ) SE  i i 2 2 ˆ  nk  nk 20 • 最小二乗推定量の確率分布（２）正規分布に従う変数の「ずれ」の評価 – 母分散がわからないので、サンプル（残差平方和）からの推定値を用いる→t分布 • t統計量の作成 ˆ 2  tj  SE  nk 2 ˆ ( y  y )  i i nk ˆ 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  1e2 x log y  log 1   2 x y  1  2 x  3 x 2 z  0   2 x y  1  2 x  3w 24