n ˆb = i=1 (yi − y¯ + ˆb¯ x)xi n i=1 n x2i ˆb = i=1 n i=1 n i=1 n x2i ˆb = x2i ˆb = n ☞ x2i n i=1 n i=1 n x2i をかける両辺に i=1 (yi − y¯ + ˆb¯ x)xi (yi − y¯)xi + ˆb¯ x n n xi xi = n¯ xより ☞ i=1 i=1 (yi − y¯)xi + n(¯ x)2ˆb ☞ ˆbでまとめると i=1 n x2i − n(¯ x)2 ˆb = i=1 (yi − y¯)xi i=1 ここで、左辺には補題 2 の (1) 式を, 右辺には補題 2 の (2) 式を用いますと n (yi − y¯)(xi − x ¯) となります. よって i=1 i=1 n ˆb = n (xi − x ¯)2 ˆb = i=1 (yi − y¯)(xi − x ¯) n (xi − x ¯)2 i=1 これより結果が示せました. 4.6 多項式回帰 (散布図をうまく表す曲線を考えてみよう) ここでは前の章で説明した回帰直線では、予測がうまくいかない場合について説明します. まず以下の散布図を考えて見ましょう. y O x 図 4.6 この散布図をうまく表す直線をあてはめてみるとどうなるでしょうか? 40 y (x10 , y10 ) y10 yˆ x10 O x 図 4.7 大雑把な傾向を表すにはこれでいいのかもしれませんが，例えば図 4.7 のような，x10 が与えられたときの予測値 yˆ = ˆb0 + ˆb1 x10 と実際の値 y10 を比較してみますと直線だとあてはまりがイマイチですね. そこで即ち 3 次式を使って散布図をうまく表すような曲線を作ってみることにします. y y = b0 + b1 x + b2 x2 + b3 x3 d1 = y1 − b0 − b1 x1 −b2 x2 − b3 x3 x O 対のデータ (x1 , y1 ), (x2 , y2 ), ..., (xn , yn ) に対して, 偏差平方和 n 2 yi − (b0 + b1 xi + b2 x2i + b3 x3i ) D(b0 , b1 , b2 , b3 ) = i=1 を最小にするような b0 = ˆb0 , b1 = ˆb1 , b2 = ˆb2 , b3 = ˆb3 から y = ˆb0 + ˆb1 x + ˆb2 x2 + ˆb3 x3 という曲線を作ろうということです. そうすると上の図のような曲線を得ることができ, この式のことを多項式回帰と呼びます. ここで実際に ˆb0 , ˆb1 , ˆb2 , ˆb3 の値は Excel などの統計ソフトが計算してくれますので心配はいりません. 4.7 多項式回帰の次数を増やすとどうなるの? ここでは一般の多項式回帰のモデルを考えてみましょう. Y = b0 + b1 x + b2 x2 + · · · + bk xk + (誤差) (4.8) ここで以下のデータを考えます. これはオリンピックにおける日本が取った金メダルの数を表したものです．ここで 1980 年のモスクワオリンピックは政治的な理由で参加していません． 41 番号 (x) 年種目数金メダル数金メダルの比率 (Y ) 1 2 3 4 5 6 1952 1956 151 153 1 4 0.006623 0.026144 1960 1964 1968 1972 150 167 182 205 4 16 11 13 0.026667 0.095808 0.06044 0.063415 7 1976 1980 1984 1988 1992 1996 2000 198 204 224 260 284 271 300 9 0.045455 10 4 3 3 5 表 4.4 0.044643 0.015385 0.010563 0.01107 0.016667 8 9 10 11 12 13 ここで，Y の値を金メダルの比率 (これは取得した金メダル数を種目数で割ったものです) としまして，説明変数 x はその年を表す番号（例えば 1984 ならば 9 という感じ）とします．この x と Y の組のデータに対して散布図を描き，それに対して多項式回帰を描いてみると以下のようになります． 42 y = −0.0021x + 0.0498 R2 = 0.0976 y = −0.0014x2 + 0.0173x + 0.0024 R2 = 0.517 図 4.8:1 次回帰図 4.9:2 次回帰 y = 0.0003x3 − 0.0077x2 +0.00534x − 0.0468 R2 = 0.7399 y = 2 × 10−5 x4 − 0.0002x3 −0.0026x2 + 0.0362x − 0.0309 R2 = 0.7484 図 4.10:3 次回帰図 4.11:4 次回帰 0.02 0.04 y 0.06 0.08 y = −6 × 10−6 x5 + 0.0002x4 −0.0031x3 + 0.013x2 + 0.0003x −0.0059 R2 = 0.7569 2 4 6 8 10 12 x 図 4.13:10 次回帰これはドイヒーな結果 · · · a 図 4.12:5 次回帰 a ドイヒーとはヒドイの業界用語. 深夜番組でしばしばみられる 43 つまり多項式の次数 (つまり k の値) を大きくすればするほど，回帰式が自由に動かすことができて，いろいろな散布図に対応できるようになるのがわかりますよね．このいろいろな形に対応できる度合いのことを自由度と言います．また単回帰モデルで説明した決定係数ですが, 多項式回帰モデルに対しても同様に作ることができます. 回帰曲線を y = ˆb0 + ˆb1 x + · · · + ˆbk xk とし, その x の部分に x1 , ..., xn を代入したときの予測値をそれぞれ yˆ1 , ..., yˆn とします. つまり yˆi = ˆb0 + ˆb1 xi +· · ·+ ˆbk xki (i = 1, ..., n) です. この時, 以下の決定係数 R2 を以下のように定義します. 1 n (yˆi − yˆ¯)2 n 予測値 y ˆ , ..., y ˆ の分散 1 n = i=1 R2 = y1 , ..., yn の分散 1 n (yi − y¯)2 n i=1 n 1 1 ここで y¯ = yi , y¯ˆ = n i=1 n わかりますよね. ✓ R2 の性質 (4.9) n yˆi です. これは単回帰モデルと同様の形であることが i=1 ✏ • R2 0 =⇒ 回帰曲線は効果がない. • データが説明変数から何% 説明することができるかを表しています. ✒ ✑ 2 【学生】 R が高い方が良い回帰式なのでしょうか? 【先生】図 4.8∼ 図 4.12 を見てもらうとわかるように, 実は次数を増やすと自動的に R2 は大きくなってしまうのです. 【学生】フムフム. ということは次数が大きい方が良い回帰式なのでしょうか? 【先生】多項式回帰モデル (4.8) において多項式の次数 k を大きくしますと，図 9.13 のように回帰曲線があまりにグニャグニャしすぎて予測に使えなくなったり, 最小二乗推定量 ˆb0 , ˆb1 ,..., ˆbk が数学的に求められなくなったりするのです. (これをモデルが不安定になると言います) これらをまとめますと以下のようになります. ✓ ✏ 多項式の次数を大きくする． =⇒ 多項式の次数を小さくする． =⇒ ✒ 自由度は増えるが, モデルが不安定になる自由度は減るが, モデルは安定 ✑ 【学生】ではどのような次数で折り合いをつければよいのでしょうか? 【先生】これを調整するための方法としまして自由度調整済み決定係数というものがあります．これは以下の形で与えられます． ✓ 自由度調整済み決定係数 ¯2 = 1 − R n−1 (1 − R2 ) n−k−1 ✏ (4.10) ここで k は多項式の次数，n はデータ数，R2 は決定係数．ただし上の 1 − n−1 ¯ 2 = 0 とします． (1 − R2 ) が負の値をとるときは R n−k−1 ✒ ✑ 【例 4.3】 10 個の組のデータ (x1 , y1 ), (x2 , y2 ), ..., (x10 , y10 ) が与えられた時，単回帰モデル Y = b0 + b1 x + の最小二乗推定量を求めたときに R2 = 0.8 であったとしましょ 44 う．このときには ¯2 = 1 − R 10 − 1 9 (1 − 0.8) = 1 − · 0.2 10 − 1 − 1 8 0.78 (4.11) です．一方で y = b0 + b1 x + b2 x2 という式で最小二乗推定量を求めたときに R2 = 0.81 であったとしましょう．その時 ¯2 = 1 − R 10 − 1 9 (1 − 0.81) = 1 − · 0.19 10 − 2 − 1 7 0.76 となります．【先生】実際にはないことですが, 例えば n = 15, R2 = 0.8 に固定して, k の値を 1 から 4 まで変化させてみると k ¯ R2 1 2 3 4 0.7846 0.7667 0.7455 0.72 ¯ 2 が小となります. つまり R2 がほとんど変わらない場合は, 次数 k を大きくすると R さくなることがわかります. 一方で n = 15, k = 2 に固定して, R2 の値を 0.6 から 0.9 まで変化させてみると R2 ¯2 R 0.6 0.7 0.8 0.9 0.5333 0.65 0.7667 0.8833 ¯ 2 も大きくなります. これをまとめますと, となり, 決定係数 R2 が大きくなりますと R 自由度調整済み決定係数は，多項式の時数 k とモデルの当てはまりのよさ R2 を同時に考えているものだということがわかりますよね．ですから多項式の次数 k を増やして ¯ 2 は減少するということになも、当てはまりのよさ R2 があまり改善されない場合は R ります. ¯ 2 は直感的に作ったのではなく，理論的な妥当さを持っていますが、もちろんこの R 特に知らなくても差し支えありませんのでここでは省略します．さて実際に表 4.4 において与えられた日本のオリンピックにおけるデータにおける自由度調整済み決定係数を Excel で求めてみると次のようになります． R2 補正 R2 1 2 0.0976 0.517 0.0074 0.4097 3 4 0.7399 0.7484 0.6424 0.6046 5 0.7569 0.5543 多項式の次数これより 3 次の回帰式がもっともよいことがわかります．問 4.5 ある対のデータ (x1 , y1 ), ..., (x20 , y20 ) に対して 3 次の多項式を当てはめた所, 回帰式が y = 0.1 + 0.9x + 0.7x2 + 0.4x3 となり, R2 = 0.75 であった. このとき自由度調整済み決定係数を求めよ. 45 第5章 5.1 重回帰分析 (説明変数が 2 つ以上ある場合) 3 次元データ次に 3 組のデータ (x11 , x12 , y1 ), (x21 , x22 , y2 ), ...,(xn1 , xn2 , yn ) が与えられた場合を考えてみましょう．これは 3 次元データと呼びますが, 具体的には以下のようなものがあります. • アパートの家賃 Y , 広さ x1 , 駅からの距離 x2 • 肌の乾燥の度合い Y ，睡眠時間 x1 ，実際の年齢 x2 このような 3 次元データに対しても 2 次元データと同様に散布図を描くことができます. 例えば (1, 2, 3) というデータなら，以下の図ように点をとります． y (1, 2, 3) 3 2 x2 1 x1 このようにして点を描いていくと以下の図のように 3 次元空間における散布図を作ることができます． y x2 x1 5.2 b0 , b1 , b2 の最小二乗推定量 3 次元データ (x1 , x2 , Y ) が与えられたとき, x1 , x2 が Y にどのような影響を与えているのかを考えてみましょう. 46 説明変数応答変数 x1 y x2 これは説明変数が 2 個の場合だと考えることができます. この章では, 説明変数が 2 つある場合の回帰分析について述べます. 3 次元空間における平面は y = b0 + b1 x1 + b2 x2 と表せることが知られています．もしこれが平面の式になることが分からない場合は，割り切って平面の式はそのようなものだと思って下さい．ここで前の章の回帰直線の時と同様に，3 次元空間における散布図をうまく表すような平面をおくことを考えます．この時, y1 と b0 + b1 x11 + b2 x12 の差を d1 , y2 と b0 + b1 x21 + b2 x22 の差を d2 とし, 残りも同様にしますと d1 = y1 − (b0 + b1 x11 + b2 x12 ) d2 = y2 − (b0 + b1 x21 + b2 x22 ) ········· dn = yn − (b0 + b1 xn1 + b2 xn2 ) となります．これは以下の図のようになります． y (x11 , x21 , y1 ) d1 = y1 − (b0 + b1 x11 + b2 x12 ) x2 平面の式 y = b0 + b1 x1 + b2 x2 x1 これより d1 ,...,dn をそれぞれ 2 乗したものの和を最小にするような平面を求めることを考えましょう．つまり n n d2i = i=1 2 yi − (b0 + b1 x1i + b2 x2i ) i=1 を最小になるような b0 , b1 , b2 の値を求めようということです．これは Excel を用いればすぐに計算してくれますが, 一応公式がありますので，ここでは結果だけを示します．証明はこの後に書いておきますので興味がある方は見てくだ 47 さい．記号の準備 1 n Sx1 Y = 1 n Sx2 Y = Sx1 x2 = x ¯1 = 1 n 1 n n (x1i − x ¯1 )(yi − y¯), Sx1 x1 = i=1 n (x2i − x ¯2 )(yi − y¯), Sx2 x2 = i=1 n 1 n 1 n n (x1i − x ¯1 )2 , i=1 n (x2i − x ¯2 )2 , i=1 (x1i − x ¯1 )(x2i − x ¯2 ), y¯ = i=1 n x1i , x ¯2 = i=1 1 n 1 n n yi , i=1 n x2i i=1 とする. ✓ b0 , b1 , b2 の最小二乗推定量 ✏ b0 , b1 , b2 の推定量 ˆb0 , ˆb1 , ˆb2 は, 以下の形で与えられる. ˆb1 = Sx2 x2 Sx1 Y − Sx1 x2 Sx2 Y , Sx1 x1 Sx2 x2 − Sx21 x2 ˆb2 = Sx1 x1 Sx2 Y − Sx1 x2 Sx1 Y , Sx1 x1 Sx2 x2 − Sx21 x2 ˆb0 = y¯ − ˆb1 x ¯1 − ˆb2 x ¯2 . ✒ ✑ 一見複雑そうですが, ˆb1 , ˆb2 の分母は同じ形をしているのに注意すれば、少しは簡単に求められると思います．【生徒】結局, これは散布図をうまく表すような平面と考えてよいのでしょうか? 【先生】はい, その通りです. ちなみにこの平面 y = ˆb0 + ˆb1 x1 + ˆb2 x2 は回帰平面もしくは重回帰式と言い, 説明変数が 2 個以上あるモデルを重回帰モデルといいます. y x2 y = ˆb0 + ˆb1 x1 + ˆb2 x2 x1 【例 5.1】レポートの成績 x1 と期末試験の点数 x2 と学期末の成績 Y には, どのような関係があるのかを重回帰モデルで考えてみましょう． A君 B君 C君 D君 E君 F君平均値成績 (Y ) レポート (x1 ) 期末試験 (x2 ) 66(y1 ) 85(y2 ) 81(y3 ) 70(y4 ) 65(y5 ) 75(x11 ) 82(x12 ) 77(x13 ) 69(x14 ) 76(x15 ) 61(x21 ) 92(x22 ) 87(x23 ) 65(x24 ) 50(x25 ) 65(y6 ) 77(x16 ) 56(x26 ) 72(¯ y) 76(¯ x1 ) 68.5(¯ x2 ) 48 表 5.1 これより x1 , x2 ，Y を 2 乗したものと x1 Y などをまず計算しますと以下の表のようになります． No. Y2 x21 x22 x1 Y x2 Y x1 x 2 1 2 3 4 4356 7225 6561 4900 4225 5625 6724 5929 4761 5776 3721 8464 7569 4225 2500 4950 6970 6237 4830 4940 4026 7820 7047 4550 3250 4575 7544 6699 4485 3800 4225 5929 3136 5005 3640 4312 5248.7 5790.7 4935.8 5488.7 5055.5 5235.8 5 6 平均値ここで y¯ = 72, x ¯1 = 76, x ¯2 = 68.5 より Sx1 x1 = Sx2 x2 = Sx1 x2 1 6 1 6 1 = 6 Sx 1 Y = Sx 2 Y = 1 6 1 6 6 x21i − x ¯21 = 5790.7 − (76)2 = 14.7, i=1 6 x22i − x ¯22 = 4935.8 − (68.5)2 = 243.55, i=1 6 x1i x2i − x ¯1 · x ¯2 = 5235.8 − 76 · 68.5 = 29.8, i=1 6 x1i yi − x ¯1 · y¯ = 5488.7 − 72 · 76 = 16.7, i=1 6 x2i yi − x ¯2 · y¯ = 5055.5 − 68.5 · 72 = 123.5, i=1 Sx1 x1 Sx2 x2 − Sx21 x2 = 14.7 · 243.55 − (29.8)2 = 2692.145. よって ˆb1 = 243.55 · 16.7 − 29.8 · 123.5 0.1437, 2692.145 ˆb2 = 14.7 · 123.5 − 29.8 · 16.7 0.4895, 2692.145 ˆb0 = 72 − 0.1437 · 76 − 0.4895 · 68.5 = 27.54805. これより, 重回帰式 y = 27.5481 + 0.1437x1 + 0.4895x2 が得ることができました．問 5.1 1 以下はあるクラスから無作為に選んだ 10 人の 1 学期から 3 学期までの評価とする. 1, 2 学期を説明変数 x1 , x2 , 3 学期を応答変数 (Y ) とし, 回帰係数を求めよ. 1 学期 2 学期 3 学期 7 7 8 8 7 7 8 8 8 8 8 8 8 7 7 8 7 7 7 7 6 7 6 7 8 8 6 8 6 7 49 5.3 予測 4.3 章の単回帰分析で述べたように、重回帰分析における回帰式を用いると予測することができます. 【例 5.2】【例 5.1】においてレポート x1 と期末試験 x2 と学期末の成績 Y には成績 = 27.5481 + 0.1437 × レポート + 0.4895 × 期末試験という関係があることがわかりました. これより例えば、金剛地君がレポート 65 点, 期末試験 70 点だったとしましょう. そのとき学期末の成績は 27.5481 + 0.1437 × 65 + 0.4895 × 70 = 71.1536 となり、約 71 点と予測されます. 問 5.2 問 5.1 のデータのクラスにおいて, 1 学期は 6, 2 学期は 7 の評価の生徒がいた. この生徒の 3 学期の成績を予測せよ. 5.4 当てはまりの良さをみる (決定係数) 説明変数が 2 つ場合も説明変数が 1 つの場合と同様にして定義できます．即ち 1 n n 予測値yˆ1 , ..., yˆn の分散 R2 = = i=1 y1 , ..., yn の分散 1 n n i=1 (yˆi − y¯ˆ)2 (5.1) (yi − y¯)2 となります．ここで単回帰モデルと違うのは yˆi = ˆb0 + ˆb1 x1i + ˆb2 x2i となっている点で n 1 す. また y¯ ˆ= yˆi となります．R2 の数値は 4 章と同様に以下のことが言えます. n i=1 ✓ R2 の性質 ✏ • R2 1 =⇒ 回帰直線の当てはまりがよい • R2 0 =⇒ 回帰直線の当てはまりが悪い • 重回帰モデルでデータを何% 説明することができるかを表す. ✒ 問 5.3 問 5.1 のデータに対する R2 を求めよ. 5.5 回帰係数の単位の変更による影響【先生】【例 5.1】において, 重回帰式は y = 27.5481 + 0.1437x1 + 0.4895x2 でした. ここで以下の例を考えてみましょう． A君 B君 C君 D君 E君 F君平均値学期末成績 (Y ) レポート (x1 ) 期末試験 (x2 ) 66(y1 ) 85(y2 ) 81(y3 ) 70(y4 ) 65(y5 ) 7.5(x11 ) 8.2(x12 ) 7.7(x13 ) 6.9(x14 ) 7.6(x15 ) 61(x21 ) 92(x22 ) 87(x23 ) 65(x24 ) 50(x25 ) 65(y6 ) 7.7(x16 ) 56(x26 ) 72(¯ y) 76(¯ x1 ) 6.85(¯ x2 ) 50 ✑ 表 5.2 これは Y と x2 は 100 点満点で成績を評価したのに対して，x1 は 10 点満点で評価しています．このとき，上のデータに対する重回帰式は y = 27.5481 + 1.437x1 + 0.4895x2 となります．一方で表 5.1 の重回帰モデルにおいて x1 の係数は ˆb2 = 0.1437 でした. 実は以下のことが知られています. ✓ ✏ 最小二乗法で求められた係数 ˆb0 , ˆb1 , ˆb2 は Y , x1 , x2 の単位によって変わる!! ✒ ✑ 1 に x1 における 100 点満点の点数を 10 点満点に単位を変えると，即ち x1 の単位を 10 すると，x1 の係数 ˆb1 の値は 10 倍された値になります．同様にして，x1 の単位を a 倍 1 すると，x1 の推定された係数 ˆb1 の値は倍された値になります． a 5.6 係数の標準誤差【学生】最小二乗法で推定された係数は推定ですから誤差が出てきますよね. その誤差を表す方法はあるのですか? 【先生】あります。それは標準誤差と呼ばれているものです. 例えば ˆb1 の標準誤差とは何かということを大雑把に説明します．少々強引ですが，表 5.1 の 6 人に加えて, その 6 人と同じレポート点, 期末試験の点数を持つ組を 9 組集め, 合計 10 組とします. そしてこの 10 組に対して学期末の成績 Y を調べて, 以下のような表を作ります．（現実には難しいでしょうけれど · · · ）表1 Y x1 x2 表 10 Y x1 x2 A君 B君 C君 D君 E君 F君 66 85 81 70 66 65 75 82 77 69 76 77 61 92 87 • • • 65 50 56 M君 N君 O君 P君 Q君 R君 67 84 82 69 66 64 75 82 77 69 76 77 61 92 87 65 50 56 x1 と x2 が同じであったとしても, 小テストや出席など個人差があるので Y の値は変わる! これらの 10 個の表に対して x1 の係数 b1 に対する最小二乗推定量をそれぞれ計算すると, 以下のようになります. ˆb1 , 1 ˆb2 , 1 ..., ˆb10 1 ここで ˆb10 1 であれば上の表 10 における b1 の推定値を表しています. この時，ˆb1 , ˆb2 ,...,ˆb10 は全て同じ値でなく各表のデータに依存して少しずつ異なる値 1 1 1 をとりますので, 標準偏差を求めると以下のようになります. Sˆb1 = ここで ˆb1 = 1 10 10 1 10 10 (ˆbi1 − ˆb1 )2 (5.2) i=1 ˆbi . 1 i=1 そうすると，これは最小二乗推定量 ˆb1 のばらつきを表していると考えられます. 実はこの Sˆ を ˆb1 の標準誤差として考えればよいということです. ただし厳密には，上 b1 51 のような 10 個の表はなく，前の章で与えられたような 1 個の表のデータだけから Sˆb1 に代わるものを与える必要があります. これは理論的な方法を用いると，たった 1 つの表から b1 の標準誤差 Sˆb1 に代わるもの σˆb1 を求めることができます. i ただしここでは Excel が計算をしくれるので, 細かい点は省略します. 5.7 t値実際には Y と x1 と x2 において Y = 1 + 0.5x1 + (誤差) (5.3) という関係があるモデルから 10 個のデータ (x11 , x21 , y1 ), (x12 , x22 , y2 ), ..., (x1 10 , x2 10 , y10 ) (5.4) が与えられているとしましょう．つまり Y は x2 の影響を受けないということです．しかし私たちは実際には上の関係式が成り立つかどうかわかりませんので， Y = b0 + b1 x1 + b2 x2 + (誤差) というモデルから最小二乗推定量 ˆb0 , ˆb1 , ˆb2 を求めることになります．そうすると ˆb2 が 0 に近い値を取ります. ここで ˆb2 の 0 からのずれを測るために ✓ ✏ t値= 回帰直線の係数 − 0 係数の標準誤差 (標準偏差) ✒ ✑ ˆb2 − 0 と定義します．ii 例えば ˆb2 の t 値ならばとなります．ここで標準誤差で ˆb2 の標準誤差割るのは, 単位の影響を受けずに 0 からの離れ具合を測ろうとしています. これより ˆb2 の t 値が 0 から離れているならば説明変数 x2 はモデルに必要な変数であると考え, ˆb2 の t 値が 0 に近い値を取るのであれば, 説明変数 x2 はモデルに不必要であると判断します. 細かいことは省略しますが, t 値が 0 に近いかどうかの判定は以下のようになります. ✓ ✏ 説明変数が 2 個の場合の t 値の性質 • ˆbi の t 値 > t(n − 3) ⇐⇒ 説明変数 xi は必要 • ˆbi の t 値 ≤ t(n − 3) ⇐⇒ 説明変数 xi は不必要 ✒ ここで t(n) は判定値と呼ばれるもので, 以下の様に定められます. ✑ n 1 2 3 4 5 6 7 8 9 10 t(n) 12.71 4.3 3.18 2.78 2.57 2.45 2.36 2.31 2.26 2.23 n 12 14 16 18 20 30 40 60 120 t(n) 2.18 2.14 2.12 2.1 2.09 2.04 2.02 2 1.98 【生徒】この判定値はどこから出てきたのですか【先生】これは t 分布という確率分布から出てきたのですが, 細かい話は 10 章でします. 単純に統計処理をするだけであれば, t 値が 0 に近いかどうかの境界の値を与えて i 標準誤差の厳密な意味での定義は佐和 ii t [6], 田中, 脇本 [9] 値を利用した少し変わった解析法として上田 [2]p37 があります 52 るというぐらいの認識で大丈夫です. また実際には t 値よりも, 次に説明する p 値の方が使われます. 【例 5.3】表 5.1 のデータに対して Excel を用いたところ, 以下のような結果を得た. b0 b1 b2 推定値標準誤差 27.8104 0.139907 0.489878 16.318 0.23567 0.057829 このとき b2 の推定値 ˆb2 の標準誤差は 0.057829 なので, ˆb2 の t 値は 0.489878 0.057829 8.47 となります. 一方で表 5.1 のデータ数は 6 個なので n = 6 となります. これより判定値は t(n − 3) = t(6 − 3) = t(3) = 3.18 となります. ここで 8.47 > 3.18 から説明変数 x2 は必要であることがわかります. 問 5.4 【例 5.3】における b1 の t 値を求め, 説明変数 x1 が必要かどうか判定せよ. 5.8 p値実は t 値を使った場合の判定値はデータ数 n の影響を受けるので、少々読みづらい数値でして, これを改善するために p 値iii というものを紹介します. ここでは b1 の最小二乗推定量 ˆb1 に対する p 値を求め方を説明します．前にも述べたように，データ (5.4) に対して，Y = b0 + b1 x1 + b2 x2 + (誤差) というモデルから最小二乗推定量 ˆb0 , ˆb1 , ˆb2 を求めます. 実際には Y = 1 + 0.5x1 + (つまり b1 = 0.5, b2 = 0) という関係があるわけですから，推定された ˆb2 は 0 ではないが 0 に近い値をとります．推定値 ˆb2 の 0 からのずれを確率を使って表したのが p 値です．p 値は Excel で自動的に出力されて, 次の様に判断します. ✓ ✏ p 値の性質 • ˆbi の p 値 > 0.05 =⇒ 説明変数 xi は必要でない • ˆbi の p 値 ≤ 0.05 =⇒ 説明変数 xi は必要である． ✒ ✑ p 値は ˆb2 と b2 = 0 の離れ具合を単なる距離でなく，確率を使ってその離れ具合を測ったことに強調しておきます．つまり p 値が 0 に近いほど b2 = 0 である確率が高いと考えます. 5.9 自由度調整済み決定係数 (最適なモデルを求める) 表 5.1 のデータを Excel で回帰分析を行うと以下のような結果が得られます。係数標準誤差 t P-値下限 95% 上限 95% 切片 27.8104 16.318 1.704278 0.186877 -24.1208 79.7416 X値1 X値2 0.139907 0.489878 0.23567 0.057829 0.593659 8.471149 0.594526 0.003454 -0.6101 0.30584 0.889914 0.673916 iii この p 値は，仮説検定と呼ばれる統計手法が用いられてます．詳しくは佐和 [6] を参照のこと 53 ここで x1 の p 値が 0.05 よりかなり大きくなっているので、説明変数 x1 は必要でないということになります。ということは自然に考えれば学期末成績 Y と期末試験点 x2 だけから回帰分析を行えば良いと考えるでしょう。ではどのような基準で説明変数を選べばよいのでしょうか? 実はこれを調整するためには、以前紹介しました自由度調整済み決定係数を使います。 ✓ ✏ 自由度調整済み決定係数 ¯2 = 1 − R n−1 (1 − R2 ) n−k−1 ここで k は説明変数の数 (説明変数が 2 個あるとき,k = 2 とおく)，n はデータ数， n−1 ¯2 = 0 と R2 は決定係数．ただし 1 − (1 − R2 ) が負の値をとるときは R n−k−1 する． ✒ ✑ ここで当てはまりの良さを表す R2 は説明変数の数を増やせば高い値をとりますが、説明変数をたくさん使うと推定された係数が不安定になります. 一方で説明変数の数を少 ¯ 2 は当てはなくすると推定された係数は安定しますが、当てはまりが悪くなります. R ¯ 2 が高いほどよいモデルであり, 正確なまりの良さと説明変数の数を調整するもので R 予測が出来ます. さて Y と x2 から回帰分析を行ない、自由度調整済み決定係数 (Excel では補正 R2 と書かれている欄) の値を調べましたところ応答変数説明変数自由度調整済み決定係数 Y x1 , x2 0.952707 Y x2 0.96036 となりました。これより x1 は無視して Y と x2 からモデルを作ったほうが良いということがわかりました。 5.10 説明変数が k 個ある場合の t 値, p 値ここでは 5.1∼5.8 章の内容をさらに一般化させて以下のモデルを考えてみましょう. Y = b0 + b1 x1 + b2 x2 + · · · + bk xk + (誤差), (5.5) ここで b0 , b1 , ..., bk は未知の係数で，x1 , ..., xk は説明変数，は誤差となります. これより， n 個のデータ (x11 , x21 , .., xk1 , y1 ), (x12 , x22 , .., xk2 , y2 ),..., (x1n , x2n , ..., xkn , yn ) が与えられたとしたら説明変数が 2 個の場合と同じようにして yi と b0 +b1 x1i +b2 x2i +· · ·+bk xki (i = 1, 2, .., n) との差の 2 乗の和が最小になるように，即ち n 2 yi − (b0 + b1 x1i + b2 x2i + · · · + bk xki ) (5.6) i=1 が最小になるような b0 = ˆb0 , b1 = ˆb1 ,...,bk = ˆbk を求めるということになります. この計算も Excel を含む統計ソフトが計算してくれますのでそれに任せることにしましょう. また決定係数 R2 も同様に定義します. 1 n (yˆi − y¯ˆ)2 n i=1 予測値yˆ1 , ..., yˆn の分散 = R = . y1 , ..., yn の分散 1 n (yi − y¯)2 n i=1 2 54 (5.7) 1 ここで yˆi = ˆb0 + ˆb1 x1i + ˆb2 x2i + · · · + ˆbk xki , y¯ ˆ= n n yˆi となります．R2 の数値解釈 i=1 は 4 章と同様です. ✓ 説明変数が k 個の場合の自由度調整済み決定係数 ¯2 = 1 − R ✏ n−1 (1 − R2 ) n−k−1 ここで k は説明変数の数，n はデータ数，R2 は決定係数．ただし上の 1 − n−1 ¯ 2 = 0 とする． (1 − R2 ) が負の値をとるときは R n−k−1 ✒ ✑ ˆ ˆ また b1 , ..., bk に対する t 値も説明変数が 2 つの場合と同じですから ˆbi の t 値 = ˆbi − 0 ˆbi の標準誤差 , (i = 1, ..., k) となります. ✓ 説明変数が k 個の場合の t 値の性質 ✏ • ˆbi の t 値 > t(n − k − 1) ⇐⇒ 説明変数 xi は必要 • ˆbi の t 値 ≤ t(n − k − 1) ⇐⇒ 説明変数 xi は不必要 ✒ ✑ 説明変数が k 個の場合も同様にやはり t 値の大きな説明変数ほど, そのモデルに必要な変数ということが言えます. p 値も説明変数が 2 個の場合とまったく同じです. 5.11 説明変数の選択 5.9 章の自由度調整済み決定係数のところで話したように、重回帰分析においては全ての説明変数を取り込むとは限りません. この章では, どの説明変数を取り込めばよいのかということを説明したいと思います. 最初に総当たり法を紹介します. ここで x1 , x2 , x3 , x4 という説明変数が 4 個あったとします. これらの中から応答変数 Y と x1 と x2 の回帰モデルを用いて自由度調整済 ¯ 2 の値を調べます。また応答変数 Y と説明変数 x2 , x3 , x4 の回帰モデルをみ決定係数 R ¯ 2 の値を調べます. このような操作を、全ての組み合用いて自由度調整済み決定係数 R ¯ 2 を持つモデルを探すというわせの説明変数に対して回帰分析を行ない、最も大きな R のも一つの手です. この方法は最も確実な方法なのですが, とても大変な作業ですので, 一般的には以下の手法が用いられます。 ✓ ステップワイズ法 ✏ 全ての説明変数を入れた回帰式を考えて、そこから p 値の大きな説明変数から減 ¯ 2 が最も大きくなるような説明変数の組み合わせのもらしていく方法。その中で R のを採用する。 ✒ ✑ これは具体例を用いて説明します。これは東西線早稲田駅と JR 山手線高田馬場駅付近にあるアパートの物件のデータです。 55 最寄り駅徒歩エアコン広さ (m2 ) 敷金礼金管理費築家賃 0 0 0 0 1 1 1 1 0 0 1 1 0 0 0 1 6 6 2 8 15 12 12 7 1 8 11 4 12 12 5 10 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 9.08 9.91 8.26 11.55 13 14.4 14.4 17 25 19.04 18 18.22 30 17.5 17.5 17.95 1 1 1 1 1 2 2 2 1 2 0 2 1 1 1 2 2 1 1 1 2 2 2 2 2 2 0 2 2 2 1 1 0.1 0.2 0.2 0.15 0.2 0.15 0.15 0.3 0.3 0.2 0.1 0 0.3 0.3 0 0.2 32 35 55 25 35 30 30 16 19 20 23 45 19 19 12 17 2 2.7 2.8 2.9 4 4.75 4.85 5.7 5.9 6 6.1 6.2 6.3 6.4 6.5 6.5 1 0 0 0 1 7 9 9 8 9 1 1 1 1 1 21.45 19.87 19.87 23.14 16.14 2 2 2 2 2 2 2 2 2 1 0 0.3 0.3 0.3 0 21 14 14 13 14 6.5 6.8 6.8 6.8 6.9 0 0 0 6 9 9 1 1 1 16.6 21.53 21.53 2 2 2 2 2 2 0.2 0.3 0.3 11 14 14 6.9 7 7 1 1 10 8 1 1 18.63 17.3 2 2 1 2 0.2 0.2 12 20 7 7 1 1 1 8 11 9 1 1 1 21 18.39 16.14 2 2 2 2 1 1 0.2 0.3 0 16 4 14 7 7.2 7.2 1 1 1 7 9 9 1 1 1 21.75 19.87 17.6 1 2 2 1 2 2 0 0.2 0.1 16 16 19 7.2 7.2 7.2 1 0 1 9 7 7 1 1 1 21.45 23.2 21.05 2 2 1 1 2 1 0.2 0.3 0.2 10 18 11 7.4 7.5 7.5 1 1 10 10 1 1 20.25 22.5 2 2 2 2 0.3 0 7 13 7.5 7.5 1 9 1 21.12 2 2 0.2 16 7.6 0 1 12 8 1 1 24.75 20.7 1 2 2 2 0 0.2 25 5 7.6 7.7 0 1 0 1 1 0 0 1 1 1 1 1 1 0 0 0 8 5 4 11 14 6 6 9 13 8 8 3 11 4 4 10 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 27.02 18.69 21 24.84 28.36 19.83 19.83 19.87 24 38.41 24.75 31.44 41.81 30.86 30.86 39 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 0 2 2 2 2 1.5 2 2 2 2 2 2 0 9 7.7 0.2 0.2 0 0.2 0.2 0.2 0.2 0.3 0.2 0.3 0.2 0 0.5 0.5 0 2 9 8 5 4 4 3 13 29 9 8 7 6 6 47 7.8 7.8 7.9 7.9 8 8 8 8.2 8.5 8.7 8.9 9.3 11.9 11.9 12.5 ここで最寄駅の項目は 1 であれば高田馬場駅付近, 0 であれば早稲田駅付近を表し、エアコンの項目は 1 であればエアコン付, 0 であればエアコンなしを意味します。また築とはアパートが築何年経ったかを表しています。徒歩の単位は分とし、敷金、礼金、管理費、家賃の単位は万円とします。また最寄り駅、徒歩、エアコン、広さ、敷金、礼金、管理費、築をそれぞれ x1 , x2 , x3 , x4 , x5 , x6 , x7 , x8 とします。 56 回帰統計重相関 R 重決定 R2 補正 R2 0.898341 0.807017 0.774169 0.920375 56 標準誤差観測数係数標準誤差 t P-値切片 1.266663 0.907165 1.396288 0.169188 最寄り駅 -0.20454 -0.09117 2.149768 0.180686 0.938277 -0.44889 0.903207 -0.01202 0.305033 0.046945 0.693689 0.021513 0.330813 0.291348 1.125515 0.01422 -0.67054 -1.94211 3.099039 8.398795 2.836275 -1.54073 0.802484 -0.84542 0.505798 0.05813 0.003275 6.51E-11 0.006713 0.130089 0.426313 0.402161 徒歩 (分) エアコン広さ (m2 ) 敷金礼金管理費 (万円) 築そうすると最寄り駅の p 値は 0.05 より大きいことから, x1 を削除して、再度回帰分析を行うと以下のようになります。回帰統計重相関 R 重決定 R2 補正 R2 標準誤差観測数切片徒歩 (分) エアコン広さ (m2 ) 敷金礼金管理費 (万円) 築 0.897313 0.805171 0.776758 0.915083 56 係数標準誤差 t P-値 1.241857 -0.09809 2.026275 0.18252 0.886808 -0.39096 1.087251 -0.01218 0.901199 0.045535 0.664948 0.021216 0.319935 0.276646 1.085259 0.014137 1.378006 -2.15404 3.047267 8.602897 2.77184 -1.41322 1.001836 -0.8616 0.174592 0.036287 0.003748 2.74E-11 0.007909 0.164044 0.321447 0.39319 次に築 (x8 ) を削除すると以下のようになります. 回帰統計重相関 R 重決定 R2 補正 R2 0.895633 0.802158 0.777932 0.912674 56 標準誤差観測数 57 切片徒歩 (分) エアコン広さ (m2 ) 敷金礼金管理費 (万円) 係数標準誤差 t P-値 0.755784 -0.10251 2.244736 0.700926 0.045126 0.613091 1.078265 -2.27157 3.661344 0.286195 0.02754 0.000614 0.183492 0.975171 -0.45442 1.359934 0.02113 0.302254 0.265959 1.035355 8.683825 3.22633 -1.70861 1.313496 1.75E-11 0.002237 0.09385 0.195133 さらに管理費 (x7 ) を削除すると回帰統計重相関 R 重決定 R2 補正 R2 標準誤差観測数切片徒歩 (分) エアコン広さ (m2 ) 敷金礼金 0.891735 0.795192 0.774711 0.919269 56 係数標準誤差 t P-値 0.939161 -0.11019 2.259547 0.183727 0.994304 -0.40765 0.691846 0.045068 0.617417 0.021282 0.304085 0.265469 1.357471 -2.445 3.659679 8.632882 3.269825 -1.53557 0.180727 0.018051 0.000608 1.77E-11 0.001952 0.130947 ¯ 2 の値を比較しますととなります。ここで R 説明変数補正 R2 x1 , x2 , x3 , x4 , x5 , x6 , x7 , x8 x2 , x3 , x4 , x5 , x6 , x7 , x8 x2 , x3 , x4 , x5 , x6 , x7 x2 , x3 , x4 , x5 , x6 0.774169 0.776758 0.777932 0.774711 これより説明変数 x2 , x3 , x4 , x5 , x6 , x7 を用いた回帰モデルが最も妥当なものだといえます。よって家賃 =0.7558 − 0.1025 × 徒歩 (分) + 2.2447 × エアコン + 0.1835 × 広さ (m2 ) + 0.9752 × 敷金 − 0.4544 × 礼金 + 1.3599 × 管理費 (万円) となり, 管理費や早稲田駅の近くか高田馬場駅の近くかというのは家賃に影響を与えていないことがわかります。 5.12 ˆb0 と ˆb1 と ˆb2 の導きかた（数学的に）これも 4 章と同様にして求めます. まず偏差平方和を求め、それを D(b0 , b1 , b2 ) とする. n (yi − b0 − b1 x1i − b2 x2i )2 D(b0 , b1 , b2 ) = i=1 58 とします. ここで上の式は b0 = ˆb0 , b1 = ˆb1 , b2 = ˆb2 で最小値をとるとしましょう. ここで b1 = ˆb1 , b2 = ˆb2 を代入したもの D(b0 , ˆb1 , ˆb2 ) を b0 の関数と考えて、D(b0 , ˆb1 , ˆb2 ) を最小にする b0 = ˆb0 を求めると ˆb0 = y¯ − ˆb1 x¯1 − ˆb2 x¯2 (5.8) となります. 同様に D(ˆb0 , b1 , ˆb2 ) を b1 の関数と考えて、D(ˆb0 , b1 , ˆb2 ) を最小にする b1 = ˆb1 を求めると n 1 ˆb1 = n i=1 x21i n x1i yi − ˆb0 i=1 x1i − ˆb2 i=1 n x1i x2i (5.9) i=1 同様に D(ˆb0 , ˆb1 , b2 ) を b2 の関数と考えて、D(ˆb0 , ˆb1 , b2 ) を最小にする b2 = ˆb2 を求めると n 1 ˆb2 = n i=1 x22i n x2i yi − ˆb0 i=1 x2i − ˆb1 i=1 n (5.10) x1i x2i i=1 となります. ここで上の ˆb1 , ˆb2 に関する式は分母を払うと以下のような式を得ます n 0= i=1 n 0= i=1 x1i yi − ˆb0 x2i yi − ˆb0 n i=1 n x1i − ˆb1 x2i − ˆb1 i=1 n i=1 n x21i − ˆb2 n x1i x2i i=1 x1i x2i − ˆb2 i=1 n x22i (5.11) i=1 この式に上の ˆb0 = · · · の式を代入し、連立方程式を解くと次のような解となります. ˆb1 = Sx2 x2 Sx1 Y − Sx1 x2 Sx2 Y , Sx1 x1 Sx2 x2 − Sx21 x2 ˆb2 = Sx1 x1 Sx2 Y − Sx1 x2 Sx1 Y . Sx1 x1 Sx2 x2 − Sx21 x2 よって ˆb1 , ˆb2 も求めることができました. 59