回帰分析(後編)

統計モデル解析特論I/II : 第05回 (11/05/19)

　今回は、回帰分析の後半である。

　● 目次: 回帰分析 : 連続変量の予測
　　 1.　単回帰分析 : 予測等に使う、連続変量の関係 [第04回の資料へジャンプ]
　　 2.　「体重の大きい者を除外」して実行するには? [第04回の資料へジャンプ]
　　 3.　有効桁数に注意せよ : どこまでが「意味ある桁」か? [第04回の資料へジャンプ]
　　 4.　重回帰分析 : 2変量以上の説明する変量(説明変量)で 1変量(目的変量)を説明
　　 5.　特定グループでの解析
　　 6.　[要点] 解析する上での注意点
　　 7.　4つの尺度と回帰分析
　　 8.　身近な統計の話題から
　　 9.　回帰分析における変数選択、総当たり法

重回帰分析 : 2変量以上の説明する変量(説明変量)で 1変量(目的変量)を説明
- 説明変量が複数になる : 単 ===> 重
- 体重を身長と胸囲で説明したい。予測したい。
- [体重]=a[身長]+b[胸囲]+c : 回帰係数を求めたい。
- 単回帰とアイディアは同じ
  - 説明される変量(目的変量)と平行に残差を取る
  - 残差には正規分布が仮定されている
  - 残差の二乗和を最小にする(最小二乗法)
1. プログラム : les0501.sas
  /* Lesson 05-1 */ /* File Name = les0501.sas 11/05/19 */ options nocenter linesize=78 pagesize=30; proc printto log = 'Kougi19/les0501_log.txt' print = 'Kougi19/les0501_Results.txt' new; /* ods listing gpath='Kougi19/SAS_ODS99'; */ data gakusei; infile 'Kougi19/all19bu.csv' firstobs=2 truncover lrecl=80 dlm=',' dsd missover ; input sex $ shintyou taijyuu kyoui jitaku $ kodukai carryer $ tsuuwa; /* if sex^='M' & sex^='F' then delete; */ if shintyou=. | taijyuu=. | kyoui=. then delete; proc print data=gakusei(obs=10); run; proc reg data=gakusei; : 回帰分析 model taijyuu=shintyou kyoui; : 複数変量を指定 output out=outreg1 predicted=pred1 residual=resid1; : 結果項目の保存 run; : proc print data=outreg1(obs=15); run; : proc plot data=outreg1; : 散布図を描く where shintyou^=. and taijyuu^=. and kyoui^=.; : 解析に使ったデータのみ plot taijyuu*shintyou; : plot taijyuu*kyoui; : plot taijyuu*pred1; : 観測値と予測値 plot resid1*pred1 /vref=0; : 残差と予測値(残差解析) plot resid1*shintyou/vref=0; : 残差と説明変量(残差解析) plot resid1*kyoui /vref=0; : 残差と説明変量(残差解析) plot resid1*taijyuu /vref=0; : 残差と目的変量(残差解析) run; : : proc univariate data=outreg1 plot normal; : 残差を正規プロットして確かめる var resid1; : run; :
2. 出力結果 : les0501_Results.txt , les0501-output.pdf
  2019年11月 4日月曜日 15時24分08秒 2 REG プロシジャモデル : MODEL1 従属変数 : taijyuu 読み込んだオブザベーション数 122 使用されたオブザベーション数 122 分散分析要因自由度平方和平均平方 F 値 Pr > F Model 2 8150.86244 4075.43122 85.77 <.0001 Error 119 5654.68871 47.51839 Corrected Total 121 13806 Root MSE 6.89336 R2 乗 0.5904 従属変数の平均 59.06066 調整済み R2 乗 0.5835 変動係数 11.67166 2019年11月 4日月曜日 15時24分08秒 3 REG プロシジャモデル : MODEL1 従属変数 : taijyuu パラメータの推定パラメータ変数自由度推定値標準誤差 t 値 Pr > |t| Intercept 1 -104.22857 12.56735 -8.29 <.0001 shintyou 1 0.78259 0.07613 10.28 <.0001 kyoui 1 0.37198 0.08029 4.63 <.0001 2019年11月 4日月曜日 15時24分12秒 4 Obs sex shintyou taijyuu kyoui jitaku kodukai carryer tsuuwa pred1 resid1 1 F 146.7 41.0 85 J 10000 Vodafone 6000 42.1952 -1.1952 2 F 148.0 43.0 80 J 50000 DoCoMo 4000 41.3527 1.6473 3 F 150.0 46.0 86 40000 . 45.1497 0.8503 4 F 151.7 41.5 80 J 35000 . 44.2483 -2.7483 5 F 152.0 35.0 77 J 60000 DoCoMo 2000 43.3671 -8.3671 6 F 153.0 46.5 87 G 10000 . 47.8695 -1.3695 7 F 153.0 55.0 78 J 30000 . 44.5217 10.4783 8 F 154.4 44.0 75 J 9000 au 2000 44.5014 -0.5014 9 F 155.0 48.0 83 G 180000 . 47.9467 0.0533 10 F 156.0 42.0 85 J 0 DoCoMo 15000 49.4733 -7.4733 2019年11月 4日月曜日 15時24分12秒 5 Plot of taijyuu*shintyou. Legend: A = 1 obs, B = 2 obs, etc. 100 + A | | A | A | A | A A 75 + A A | A A BBA A A A A A | BB C A BAAA A A A taijyuu | A A A CA C AA B | B A AAAA D AA AA A AAA A | A A AA A BABB AAA A 50 + AC B CD AAB B A | A A A B A A BA A | AA A A | | A | 25 + --+------------+------------+------------+------------+------------+-- 140 150 160 170 180 190 shintyou 2019年11月 4日月曜日 15時24分12秒 6 Plot of taijyuu*kyoui. Legend: A = 1 obs, B = 2 obs, etc. 100 + A | | A | A | A | A A 75 + A A | A BE AB A | B A CAABAB B A taijyuu | A A AGA B A | BC ADCDB A | A A A AAA CD AB 50 + A B DD G B | B DAA BA | B B | | A | 25 + ---+-------------+-------------+-------------+-------------+-- 40 60 80 100 120 kyoui 2019年11月 4日月曜日 15時24分12秒 9 Plot of resid1*shintyou. Legend: A = 1 obs, B = 2 obs, etc. 40 + | A | | | A | A 残 20 + 差 | | A A | A A A AA | A A B ABA AA A | AA A A BA A BA AAAA 0 +----------AA--A---A-AAB-A--BA-----AAAA-AFA-A--AA-A-A-A--------------- | A B B DC ABB BA B B BAA A A A A A | A AA CB A A A A ABA A | A A A A | A A | -20 + --+------------+------------+------------+------------+------------+-- 140 150 160 170 180 190 shintyou 2019年11月 4日月曜日 15時24分12秒 10 Plot of resid1*kyoui. Legend: A = 1 obs, B = 2 obs, etc. 40 + | A | | | A | A 残 20 + 差 | | A A | A AA B | A AB A BD | A A AAC AAA B B A 0 +-----------------------A---C--CAABFBBD--BB-----A------------- | EBCAHAEEA B | A ABBA E AC A | A AB | A A | -20 + ---+-------------+-------------+-------------+-------------+-- 40 60 80 100 120 kyoui 2019年11月 4日月曜日 15時24分12秒 11 Plot of resid1*taijyuu. Legend: A = 1 obs, B = 2 obs, etc. 40 + | A | | | A | A 残 20 + 差 | | A A | A A A AA | A AA B CBA | A A A B CAACA A 0 +-----------------AAAAACB-B-B-ADBD-C--AA---------------------- | A CADCD AD BCBA B AA | A ABB B AAAC B A | A A B | A A | -20 + ---+-------------+-------------+-------------+-------------+-- 20 40 60 80 100 taijyuu
3. 結果の見方
  - 対象になったのは 122名。
  - 説明変量群が予測に役立っているか?
    - 回帰に役立っているか : Prob>F : 小さいと有意
    - 「役立っている」と言える : 0.01% だから 1% で有意
  - 決定係数 : R-Square ( 相関係数 : R )
    - 目的変量が説明変量でどの程度説明しているかの割合。
    - 1 に近いほど当てはまりが良いと言える。: 59.0%
    - 説明変量数が増えると大きくなるのが一般的。
  - 回帰係数 : Parameter Estimate
    - 回帰式: a=0.783, b=0.372, c=-104
  - ある特定の説明変量が予測に役立っているか?
    - 回帰係数の検定(帰無仮説:係数=0 か?) : Prob>|T| : 小さいと有意
    - 両方とも(身長も胸囲も)有意
    - 「各係数は 0ではない」と言える : 0.01% だから 1% で有意
  - 残差の性質 ===> 正規性 : 残差プロット、残差解析
    - 残差(予測誤差)は正規分布をしていると仮定してモデルが構築されている。
    - この仮定が覆ると、回帰分析として成立していないことになる。
    - 残差が正規分布をしているか確認する必要がある。
    - 均等に散らばっているか?
    - 傾向はないか? : もし傾向があると言うことになれば正規性の仮定が崩れている
    - 体重の大きい 3例程度が外れ値と考えられるか要確認 ===> [演習1](第5節)
    - ...
  - ...

特定グループでの解析

「男性のみ」と言う特定のグループに対して、同様の解析を行うには?

プログラム : les0502.sas

/* Lesson 05-2 */ /* File Name = les0502.sas 11/05/19 */ options nocenter linesize=78 pagesize=30; proc printto log = 'Kougi19/les0502_log.txt' print = 'Kougi19/les0502_Results.txt' new; /* ods listing gpath='Kougi19/SAS_ODS99'; */ data gakusei; infile 'Kougi19/all19bu.csv' firstobs=2 truncover lrecl=80 dlm=',' dsd missover ; input sex $ shintyou taijyuu kyoui jitaku $ kodukai carryer $ tsuuwa; if sex^='M' & sex^='F' then delete; : 性別不明は除外 if shintyou=. | taijyuu=. | kyoui=. then delete; : 欠損のあるデータは除外 proc print data=gakusei(obs=10); run; proc corr data=gakusei; : 相関係数 where sex='M'; : 男性について run; : : proc reg data=gakusei; : 回帰分析 model taijyuu=shintyou kyoui; : where sex='M'; : 男性について output out=outreg1 predicted=pred1 residual=resid1; : run; : proc print data=outreg1(obs=15); run; proc plot data=outreg1; where sex='M'; : 対象データについて plot taijyuu*shintyou; plot taijyuu*kyoui; plot taijyuu*pred1; plot resid1*(pred1 shintyou kyoui taijyuu)/vref=0; : まとめて記述 /* plot resid1*pred1 /vref=0; plot resid1*shintyou/vref=0; plot resid1*kyoui /vref=0; plot resid1*taijyuu /vref=0; */ run; proc univariate data=outreg1 plot normal; var resid1; run;

出力結果 : les0502_Results.txt , les0502-output.pdf
2019年11月 5日火曜日 08時41分38秒 2 CORR プロシジャ 5 変数 : shintyou taijyuu kyoui kodukai tsuuwa 単純統計量変数 N 平均標準偏差合計 shintyou 71 172.53662 5.93508 12250 taijyuu 71 64.72817 9.06513 4596 kyoui 71 88.09859 9.68527 6255 kodukai 67 56358 66472 3776000 tsuuwa 14 6632 4248 92850 単純統計量変数最小値最大値 shintyou 156.00000 185.00000 taijyuu 46.00000 100.00000 kyoui 46.00000 112.00000 kodukai 0 350000 tsuuwa 350.00000 15000 2019年11月 5日火曜日 08時41分38秒 3 CORR プロシジャ Pearson の相関係数 H0: Rho=0 に対する Prob > |r| オブザベーション数 shintyou taijyuu kyoui kodukai tsuuwa shintyou 1.00000 0.39968 0.15872 0.09516 0.11552 0.0006 0.1862 0.4437 0.6941 71 71 71 67 14 taijyuu 0.39968 1.00000 0.40227 0.11042 0.14591 0.0006 0.0005 0.3737 0.6187 71 71 71 67 14 kyoui 0.15872 0.40227 1.00000 -0.37945 -0.38661 0.1862 0.0005 0.0015 0.1721 71 71 71 67 14 kodukai 0.09516 0.11042 -0.37945 1.00000 0.53783 0.4437 0.3737 0.0015 0.0473 67 67 67 67 14 tsuuwa 0.11552 0.14591 -0.38661 0.53783 1.00000 0.6941 0.6187 0.1721 0.0473 14 14 14 14 14 2019年11月 5日火曜日 08時41分38秒 4 REG プロシジャモデル : MODEL1 従属変数 : taijyuu 読み込んだオブザベーション数 71 使用されたオブザベーション数 71 分散分析要因自由度平方和平均平方 F 値 Pr > F Model 2 1596.38065 798.19033 13.06 <.0001 Error 68 4155.98301 61.11740 Corrected Total 70 5752.36366 Root MSE 7.81776 R2 乗 0.2775 従属変数の平均 64.72817 調整済み R2 乗 0.2563 変動係数 12.07784 2019年11月 5日火曜日 08時41分38秒 5 REG プロシジャモデル : MODEL1 従属変数 : taijyuu パラメータの推定パラメータ変数自由度推定値標準誤差 t 値 Pr > |t| Intercept 1 -54.72134 27.50850 -1.99 0.0507 shintyou 1 0.52620 0.15946 3.30 0.0015 kyoui 1 0.32534 0.09772 3.33 0.0014 2019年11月 5日火曜日 08時41分41秒 6 Obs sex shintyou taijyuu kyoui jitaku kodukai carryer tsuuwa pred1 resid1 1 M 156 61 90 J 0 . 56.6453 4.3547 2 M 160 48 76 J 5000 No . 54.1954 -6.1954 3 M 162 62 100 J 0 . 63.0558 -1.0558 4 M 163 46 80 J 50000 . 57.0753 -11.0753 5 M 163 51 82 J 10000 . 57.7260 -6.7260 6 M 163 60 78 80000 . 56.4246 3.5754 7 M 165 55 90 J 10000 DoCoMo 7000 61.3810 -6.3810 8 M 166 55 85 100000 . 60.2806 -5.2806 9 M 166 60 88 G 150000 . 61.2566 -1.2566 10 M 167 60 94 J . . 63.7348 -3.7348 11 M 167 65 92 78000 . 63.0841 1.9159 12 M 167 65 96 J 30000 . 64.3854 0.6146 13 M 168 54 56 40000 . 51.8982 2.1018 14 M 168 58 88 G 15000 . 62.3089 -4.3089 15 M 168 62 80 G 120000 . 59.7063 2.2937 2019年11月 5日火曜日 08時41分41秒 7 Plot of taijyuu*shintyou. Legend: A = 1 obs, B = 2 obs, etc. 100 + A | | | A A | | A 80 + | A A A | A A taijyuu | A A A AA A A A A | B B A A B A A A A A A | A A C A B B B 60 + A A A A B A A A A A A A A | A A A A B A A A | A | A | A A | 40 + --+----------+----------+----------+----------+----------+----------+- 155 160 165 170 175 180 185 shintyou 2019年11月 5日火曜日 08時41分41秒 8 Plot of taijyuu*kyoui. Legend: A = 1 obs, B = 2 obs, etc. 100 + A | | | A A | | A 80 + | B A | AA taijyuu | AD AB A | A A CAAAAB B A | A A AEA B A 60 + BB AABAD A | A A A B CA | A | A | A A | 40 + ---+-------------+-------------+-------------+-------------+-- 40 60 80 100 120 kyoui 2019年11月 5日火曜日 08時41分41秒 11 Plot of resid1*shintyou. Legend: A = 1 obs, B = 2 obs, etc. 40 + | | A | | | A A 残 20 + 差 | A | | A A A | A A A | A A A C A A AA B A A 0 +----------------A--------A-A------B------A---A---A---A--A---A-------- | A B A D A C A B A A A A A | A A A A A A A A A B A | A A A | A | -20 + --+----------+----------+----------+----------+----------+----------+- 155 160 165 170 175 180 185 shintyou 2019年11月 5日火曜日 08時41分41秒 12 Plot of resid1*kyoui. Legend: A = 1 obs, B = 2 obs, etc. 40 + | | A | | | A A 残 20 + 差 | A | | A B | BA | A AB A B CA A A A 0 +-----------------------------------AAD-AAB--A---------------- | A A A ABACCAAAC A | A AAA B BD | A AA | A | -20 + ---+-------------+-------------+-------------+-------------+-- 40 60 80 100 120 kyoui 2019年11月 5日火曜日 08時41分41秒 13 Plot of resid1*taijyuu. Legend: A = 1 obs, B = 2 obs, etc. 40 + | | A | | | A A 残 20 + 差 | A | | A AA | AA A | A AAA C AA CA A 0 +-----------------------A-AB--BA-A-BA--------------------------------- | AA A EACC BA B | A A BA AAC B | A A A | A | -20 + --+----------+----------+----------+----------+----------+----------+- 40 50 60 70 80 90 100 taijyuu

結果の見方
- 単変量毎の相関が有意なのは、身長と体重、体重と胸囲の間。
- 対象になったのは 71名。
- 回帰に役立っているか : 役立っている : 0.01% だから 1% で有意
- 決定係数(R-square)は 27.8%
- 個々の説明変量が予測に役立っているか?
  - 係数がゼロか? : 定数項も身長も胸囲も有意(1% で有意)
- 残差の性質 ===> 正規性 : 残差プロット、残差解析
  - 均等に散らばっているか?
  - 傾向はないか? : 傾向があると言うことは正規性の仮定が崩れていること
  - 外れ値? 80Kg より重い 4名程度が吟味対象?
[演習1] : 「男性のみ」で、かつ「体重の大きい数名を除外」して実行してみよ。
- プログラム : les0503.sas、出力結果 : les0503_Results.txt , les0503-output.pdf
  where sex='M' and taijyuu<80;
- 当てはまりは良くなったか? : 異常値と外れ値の意味するもの
- 残差の正規性はどのように変化したか?

[回帰分析:要点] 解析する上での注意点
- 説明する変量の係数を求める手法
- 残差解析を行なう: 残差が正規分布をしていることを仮定している手法
- 残差に傾向はないか? : バナナカーブ、ラッパ状、...
- 対象集団を単質のものに絞る : 男性だけ、18歳だけ、...
- 外れ値を持つ個体を除外する : 標準体型の者だけ、...
  - ただし、恣意的な除外は謹むべきである
- 「外れ値」と「異常値」は区別して取り扱う・識別する : 「外れて」いても「異常」かどうかは別
  - 外れ値 : 大多数のデータから大きく離れた値
  - 異常値 : 通常ではあり得ない入力ミス等の異常な値
- 上記以外に「多重共線性(multi-colinearity)」の問題もある。逆行列が取れなくなる(ランク落ち)。不安定な解。
4つの尺度と回帰分析
- 名義尺度(nominal scale) : 性別、職業、支持政党。分類の項目。
- 順序尺度(ordinal scale) : 成績の優・良・可・不可。順序関係あり。「優と良の差」が「可と不可の差」に等しいわけではない。
- 間隔尺度(interval scale) : セ氏やカ氏で測定された温度。原点は任意。
- 比率尺度(ratio scale) : 長さ、重さ。間隔尺度に加えて絶対原点が存在。
- 質的データ : 名義尺度、順序尺度。離散変量とも。
- 量的データ : 間隔尺度、比率尺度。連続変量とも。
- 離散変量を回帰分析で用いるには注意が必要である。線形関係と言えるのか? 割り当てた値を変更したら影響力が変化する。
身近な統計の話題から
1. 誤用?!
  - 構造がシンプルで理解し易い分、"変な"利用も散見される
  - 内挿(観測点の内側)はまだしも、外挿(観測点の外側)を予測するのは難しい or 無謀。
  - 予測範囲でデータの構造が一定という"条件"が必要。
  - 線形で表現できる関係なのか? 線形の関係式が有効なのか? 非線型?
  - 自分で判断できる能力を持つこと。疑うこと。
  　[例1] 人間の成長曲線
  　[例2] 将来のプログラマ必要数予測 : 21世紀(?)には国民全員がプログラマ ('80s)
  　[例3] オリンピック 100m 走の男女記録 : 2156年には女性の方が速い (2004.09.30) :
  　　　　　　　 Japan Journal LTD の記事 , 朝日新聞の記事
  　　　　　[究極の命題!] 100m に 0.00秒要する(!?)ようになるのは何時?
2. 教育関係: 教育効果と所得の関係
  - 全国学力テスト　初めて英語も実施　小中２１２万人参加 _ : 産経ニュース 19年4月18日
  - 2019年度全国学力テスト: 公立学校の都道府県別・政令市別の平均正答率 _ : 毎日新聞 19年7月
  - 【全国学力テスト】理解力はあるのに…活用力や表現力課題 _ : 産経ニュース 19年7月31日
  - 全国学力テスト小中、応用力なお課題　正答率「基礎」と大差 _ : 毎日新聞2018年8月1日
  - 小４、中２全４教科で過去最高点…国際学力調査 2016年11月29日読売新聞
  - 読書量減少が原因？日本の「読解力」８位に低下 _ _ 2016年12月7日読売新聞
  - 国際数学・理科教育動向調査（TIMSS）の調査結果: 国際教育到達度評価学会（ＩＥＡ）
  - 国際学力調査（PISA、TIMSS）
  - 教育課程研究センター「全国学力・学習状況調査」 : 国立教育政策研究所
  - 就学援助率と成績の関係: 一見、負の相関のように見えるが...

回帰分析における変数選択 : 　回帰分析では回帰係数や重相関係数を知ることだけでなく残差解析も重要であることを強調したつもりである。次に説明変数の取捨選択(変数選択)について説明する。

今までのモデル : 指定した説明変数全てを取り込んで計算 : フルモデル

影響力の大きい説明変数だけをモデルに取り込みたい。

たくさんの説明変数の中から影響力の大きい、いくつかの変数を選別する。

幾つかの選択方法がある :

逐次増減法(stepwise)
前進選択法(forward)
後退選択法(backward)、...

総当たり法(rsquare)

全部の組合わせを概観してみたい : 総当たり
モデルの探索用 : どの組合わせを使うと説明力が大きいか

例題 : アメリカの大気汚染データ : usair2.txt
- 41都市、7変数。
- 7変量 : 都市名, SO2(二酸化硫黄)濃度, 気温, 製造業数, 人口, 風速, 降雨量, 降雨日数
- SO2 濃度を他の 6変数で説明してみよう。
- 6変数の中には影響力が大きいものがあるはずだ。<==> 影響力が小さいもの
- 参考 : 大気汚染物質広域監視システム「そらまめ君」とその測定項目

プログラム : les0504.sas

/* Lesson 05-4 */ /* File Name = les0504.sas 11/05/19 */ options nocenter linesize=78 pagesize=30; proc printto log = 'Kougi/les0504_log.txt' print = 'Kougi/les0504_Results.txt' new; ods listing gpath='Kougi/SAS_ODS99'; data air; infile '/folders/myfolders/Kougi/usair2.prn'; input id $ y x1 x2 x3 x4 x5 x6; /* label y='SO2 of air in micrograms per cubic metre' x1='Average annual temperature in F' x2='Number of manufacturing enterprises employing 20 or more workers' x3='Population size (1970 census); in thousands' x4='Average annual wind speed in miles per hour' x5='Average annual precipitation in inches' x6='Average number of days with precipitation per year' ; */ proc print data=air(obs=10); run; proc corr data=air; run; proc reg data=air; : model y=x1 x2 x3 x4 x5 x6; : フルモデル output out=outreg1 predicted=pred1 residual=resid1; : run; : proc plot data=outreg1; plot resid1*pred1 /vref=0; : plot resid1*x1 /vref=0; : ズラズラと列記 plot resid1*x2 /vref=0; : plot resid1*x3 /vref=0; : plot resid1*x4 /vref=0; : plot resid1*x5 /vref=0; : plot resid1*x6 /vref=0; : plot resid1*y /vref=0; : run; proc reg data=air; : model y=x1-x6 / selection=stepwise; : 逐次増減法 output out=outreg1 predicted=pred1 residual=resid1; : 連続変数の指定方法 run; : proc print data=outreg1(obs=15); run; proc plot data=outreg1; plot resid1*pred1 /vref=0; : plot resid1*(x1 x2 x3 x4 x5 x6) /vref=0; : 簡略形(上と比較せよ) plot resid1*(x1-x6) /vref=0; : 簡略形(これも同じ意味) plot resid1*y /vref=0; : run; proc reg data=air; : model y=x1-x6 / selection=rsquare; : 総当たり法 run; :

出力結果 : les0504_Results.txt , les0504-output.pdf
2019年11月 4日月曜日 15時26分46秒 90 Obs id y x1 x2 x3 x4 x5 x6 1 Phoenix 10 70.3 213 582 6.0 7.05 36 2 Little_R 13 61.0 91 132 8.2 48.52 100 3 San_Fran 12 56.7 453 716 8.7 20.66 67 4 Denver 17 51.9 454 515 9.0 12.95 86 5 Hartford 56 49.1 412 158 9.0 43.37 127 6 Wilmingt 36 54.0 80 80 9.0 40.25 114 7 Washingt 29 57.3 434 757 9.3 38.89 111 8 Jacksonv 14 68.4 136 529 8.8 54.47 116 9 Miami 10 75.5 207 335 9.0 59.80 128 10 Atlanta 24 61.5 368 497 9.1 48.34 115 2019年11月 4日月曜日 15時26分46秒 91 CORR プロシジャ 7 変数 : y x1 x2 x3 x4 x5 x6 単純統計量変数 N 平均標準偏差合計 y 41 30.04878 23.47227 1232 x1 41 55.76341 7.22772 2286 x2 41 463.09756 563.47395 18987 x3 41 608.60976 579.11302 24953 x4 41 9.44390 1.42864 387.20000 x5 41 36.76902 11.77155 1508 x6 41 113.90244 26.50642 4670 単純統計量変数最小値最大値 y 8.00000 110.00000 x1 43.50000 75.50000 x2 35.00000 3344 x3 71.00000 3369 x4 6.00000 12.70000 x5 7.05000 59.80000 x6 36.00000 166.00000 2019年11月 4日月曜日 15時26分46秒 93 CORR プロシジャ Pearson の相関係数, N = 41 H0: Rho=0 に対する Prob > |r| y x1 x2 x3 x4 x5 x6 y 1.00000 -0.43360 0.64477 0.49378 0.09469 0.05429 0.36956 0.0046 <.0001 0.0010 0.5559 0.7360 0.0174 x1 -0.43360 1.00000 -0.19004 -0.06268 -0.34974 0.38625 -0.43024 0.0046 0.2340 0.6970 0.0250 0.0126 0.0050 x2 0.64477 -0.19004 1.00000 0.95527 0.23795 -0.03242 0.13183 <.0001 0.2340 <.0001 0.1341 0.8405 0.4113 x3 0.49378 -0.06268 0.95527 1.00000 0.21264 -0.02612 0.04208 0.0010 0.6970 <.0001 0.1819 0.8712 0.7939 x4 0.09469 -0.34974 0.23795 0.21264 1.00000 -0.01299 0.16411 0.5559 0.0250 0.1341 0.1819 0.9357 0.3052 x5 0.05429 0.38625 -0.03242 -0.02612 -0.01299 1.00000 0.49610 0.7360 0.0126 0.8405 0.8712 0.9357 0.0010 x6 0.36956 -0.43024 0.13183 0.04208 0.16411 0.49610 1.00000 0.0174 0.0050 0.4113 0.7939 0.3052 0.0010 2019年11月 4日月曜日 15時26分47秒 94 REG プロシジャモデル : MODEL1 従属変数 : y 読み込んだオブザベーション数 41 使用されたオブザベーション数 41 分散分析要因自由度平方和平均平方 F 値 Pr > F Model 6 14755 2459.10601 11.48 <.0001 Error 34 7283.26641 214.21372 Corrected Total 40 22038 Root MSE 14.63604 R2 乗 0.6695 従属変数の平均 30.04878 調整済み R2 乗 0.6112 変動係数 48.70761 2019年11月 4日月曜日 15時26分47秒 95 REG プロシジャモデル : MODEL1 従属変数 : y パラメータの推定パラメータ変数自由度推定値標準誤差 t 値 Pr > |t| Intercept 1 111.72848 47.31810 2.36 0.0241 x1 1 -1.26794 0.62118 -2.04 0.0491 x2 1 0.06492 0.01575 4.12 0.0002 x3 1 -0.03928 0.01513 -2.60 0.0138 x4 1 -3.18137 1.81502 -1.75 0.0887 x5 1 0.51236 0.36276 1.41 0.1669 x6 1 -0.05205 0.16201 -0.32 0.7500 2019年11月 4日月曜日 15時26分47秒 96 Plot of resid1*pred1. Legend: A = 1 obs, B = 2 obs, etc. 50 + A | | | | | A 残 25 + 差 | | A A | A A | A A A A | B AAA A A 0 +---------------AA------A---A---------A-----------------------A------ | A ABA | A AA | AA B | C A A | -25 + A ---+--------+--------+--------+--------+--------+--------+--------+-- -20 0 20 40 60 80 100 120 予測値 y 2019年11月 4日月曜日 15時26分47秒 97 Plot of resid1*x1. Legend: A = 1 obs, B = 2 obs, etc. 50 + A | | | | | A 残 25 + 差 | | A A | A A | A A A A | AA A A A A A 0 +---------------AAAA------A--------------------------A-------- | B A A A | A A A | A A A A | A A A A A | -25 + A ---+-------------+-------------+-------------+-------------+-- 40 50 60 70 80 x1 2019年11月 4日月曜日 15時26分47秒 98 Plot of resid1*x2. Legend: A = 1 obs, B = 2 obs, etc. 50 + A | | | | | A 残 25 + 差 | | A A | A A | A A A A | BA AA A A 0 +----AAA-A--------A-------------------------------A----------- | A B B | A A A | AAAA | A B A A | -25 + A ---+-------------+-------------+-------------+-------------+-- 0 1000 2000 3000 4000 x2 2019年11月 4日月曜日 15時26分47秒 99 Plot of resid1*x3. Legend: A = 1 obs, B = 2 obs, etc. 50 + A | | | | | A 残 25 + 差 | | A A | A A | A A A A | AAA A AA A 0 +----A-AA--A------------A-------------------------A----------- | A AA B | A A A | AAA A | A AA B | -25 + A ---+-------------+-------------+-------------+-------------+-- 0 1000 2000 3000 4000 x3 2019年11月 4日月曜日 15時26分47秒 100 Plot of resid1*x4. Legend: A = 1 obs, B = 2 obs, etc. 50 + A | | | | | A 残 25 + 差 | | A A | A A | A A B | A AAA A B 0 +--------------------A--B-------A-A---------------A----------- | A A A A A | A A A | A A AA | A A A A A | -25 + A ---+-------------+-------------+-------------+-------------+-- 6 8 10 12 14 x4 2019年11月 4日月曜日 15時26分47秒 101 Plot of resid1*x5. Legend: A = 1 obs, B = 2 obs, etc. 50 + A | | | | | A 残 25 + 差 | | A A | A A | A AA A | A A A AA A A 0 +-----------------------------------B---A--A------A-----------------A- | A A A A A | A A A | A A A A | A A AA A | -25 + A --+----------+----------+----------+----------+----------+----------+- 0 10 20 30 40 50 60 x5 2019年11月 4日月曜日 15時26分47秒 102 Plot of resid1*x6. Legend: A = 1 obs, B = 2 obs, etc. 50 + A | | | | | A 残 25 + 差 | | A A | A A | A A B | A A A A AA A 0 +--------------------------A-----------------A-BA-A------------------- | A B A A | A A A | A A A A | B A A A | -25 + A --+----------+----------+----------+----------+----------+----------+- 25 50 75 100 125 150 175 x6 2019年11月 4日月曜日 15時26分47秒 103 Plot of resid1*y. Legend: A = 1 obs, B = 2 obs, etc. 50 + A | | | | | A 残 25 + 差 | | A A | A A | A A A A | AA BA A A 0 +------AA-------A----A---------A--------------------------A------- | A AA AA | AA A | B B | ABA A | -25 + A ---+---------+---------+---------+---------+---------+---------+-- 0 20 40 60 80 100 120 y 2019年11月 4日月曜日 15時26分47秒 104 REG プロシジャモデル : MODEL1 従属変数 : y 読み込んだオブザベーション数 41 使用されたオブザベーション数 41 ステップワイズ法: ステップ 1 変数 x2 の追加 : R2 乗 = 0.4157 C(p) = 23.1089 分散分析要因自由度平方和平均平方 F 値 Pr > F Model 1 9161.74469 9161.74469 27.75 <.0001 Error 39 12876 330.15789 Corrected Total 40 22038 2019年11月 4日月曜日 15時26分47秒 105 REG プロシジャモデル : MODEL1 従属変数 : y ステップワイズ法: ステップ 1 パラメータ Type II 変数推定値標準誤差平方和 F 値 Pr > F Intercept 17.61057 3.69159 7513.50474 22.76 <.0001 x2 0.02686 0.00510 9161.74469 27.75 <.0001 条件数における境界 : 1, 1 ------------------------------------------------------------------------------ ステップワイズ法: ステップ 2 変数 x3 の追加 : R2 乗 = 0.5863 C(p) = 7.5586 2019年11月 4日月曜日 15時26分47秒 106 REG プロシジャモデル : MODEL1 従属変数 : y ステップワイズ法: ステップ 2 分散分析要因自由度平方和平均平方 F 値 Pr > F Model 2 12921 6460.63359 26.93 <.0001 Error 38 9116.63526 239.91145 Corrected Total 40 22038 パラメータ Type II 変数推定値標準誤差平方和 F 値 Pr > F Intercept 26.32508 3.84044 11273 46.99 <.0001 x2 0.08243 0.01470 7548.02378 31.46 <.0001 x3 -0.05661 0.01430 3759.52248 15.67 0.0003 条件数における境界 : 11.434, 45.735 ------------------------------------------------------------------------------ ステップワイズ法: ステップ 3 2019年11月 4日月曜日 15時26分47秒 107 REG プロシジャモデル : MODEL1 従属変数 : y ステップワイズ法: ステップ 3 変数 x6 の追加 : R2 乗 = 0.6174 C(p) = 6.3610 分散分析要因自由度平方和平均平方 F 値 Pr > F Model 3 13606 4535.41173 19.90 <.0001 Error 37 8431.66725 227.88290 Corrected Total 40 22038 パラメータ Type II 変数推定値標準誤差平方和 F 値 Pr > F Intercept 6.96585 11.77691 79.72552 0.35 0.5578 x2 0.07433 0.01507 5547.32154 24.34 <.0001 x3 -0.04939 0.01454 2628.36952 11.53 0.0016 x6 0.16436 0.09480 684.96801 3.01 0.0913 2019年11月 4日月曜日 15時26分47秒 108 REG プロシジャモデル : MODEL1 従属変数 : y ステップワイズ法: ステップ 3 条件数における境界 : 12.65, 78.633 ------------------------------------------------------------------------------ モデル内のすべての変数は水準 0.1500 で有意です。モデルへの変数追加で、他の変数は有意水準 0.1500 で満たされていません。 2019年11月 4日月曜日 15時26分47秒 109 REG プロシジャモデル : MODEL1 従属変数 : y ステップワイズ法の要約変数の変数の取り込んだステップ追加削除変数の数偏 R2 乗モデル R2 乗 1 x2 1 0.4157 0.4157 2 x3 2 0.1706 0.5863 3 x6 3 0.0311 0.6174 ステップワイズ法の要約ステップ C(p) F 値 Pr > F 1 23.1089 27.75 <.0001 2 7.5586 15.67 0.0003 3 6.3610 3.01 0.0913 2019年11月 4日月曜日 15時26分48秒 110 Obs id y x1 x2 x3 x4 x5 x6 pred1 resid1 1 Phoenix 10 70.3 213 582 6.0 7.05 36 -0.032 10.0316 2 Little_R 13 61.0 91 132 8.2 48.52 100 23.646 -10.6461 3 San_Fran 12 56.7 453 716 8.7 20.66 67 16.285 -4.2849 4 Denver 17 51.9 454 515 9.0 12.95 86 29.410 -12.4103 5 Hartford 56 49.1 412 158 9.0 43.37 127 50.661 5.3392 6 Wilmingt 36 54.0 80 80 9.0 40.25 114 27.698 8.3020 7 Washingt 29 57.3 434 757 9.3 38.89 111 20.079 8.9208 8 Jacksonv 14 68.4 136 529 8.8 54.47 116 10.011 3.9887 9 Miami 10 75.5 207 335 9.0 59.80 128 26.844 -16.8439 10 Atlanta 24 61.5 368 497 9.1 48.34 115 28.673 -4.6731 11 Chicago 110 50.6 3344 3369 10.4 34.44 122 109.181 0.8191 12 Indianap 28 52.3 361 746 9.7 38.74 121 16.840 11.1603 13 Des_Moin 17 49.0 104 201 11.2 30.85 103 21.697 -4.6973 14 Wichita 8 56.6 125 277 12.7 30.58 82 16.053 -8.0528 15 Louisvil 30 55.6 291 593 8.3 43.11 123 19.522 10.4776 2019年11月 4日月曜日 15時26分48秒 111 Plot of resid1*pred1. Legend: A = 1 obs, B = 2 obs, etc. | 50 + A | | | A 25 + 残 | B 差 | A AAA A | A A AAA A B 0 +---------------A--AA-----------------------------------------A------- | CAAA A A A A | B AA | AA A A A -25 + | A | | -50 + | --+----------+----------+----------+----------+----------+----------+- 0 20 40 60 80 100 120 予測値 y 2019年11月 4日月曜日 15時26分48秒 112 Plot of resid1*x1. Legend: A = 1 obs, B = 2 obs, etc. | 50 + A | | | A 25 + 残 | A A 差 | A A A A A | A AA A A B A 0 +----------------AA-----A---A--------------------------------- | AA B B A B A | A A A A | A A A A A -25 + | A | | -50 + | ---+-------------+-------------+-------------+-------------+-- 40 50 60 70 80 x1 2019年11月 4日月曜日 15時26分48秒 113 Plot of resid1*x2. Legend: A = 1 obs, B = 2 obs, etc. | 50 + A | | | A 25 + 残 | A A 差 | A AAA A | BB A B A 0 +--A--A-----------A-------------------------------A----------- | AAAACA A A | A A AA | A A AAA -25 + | A | | -50 + | ---+-------------+-------------+-------------+-------------+-- 0 1000 2000 3000 4000 x2 2019年11月 4日月曜日 15時26分48秒 114 Plot of resid1*x3. Legend: A = 1 obs, B = 2 obs, etc. | 50 + A | | | A 25 + 残 | A A 差 | A B A A | ABA AAA A 0 +---A--A----------------A-------------------------A----------- | AAAAB AAA A | A A B | AA B A -25 + | A | | -50 + | ---+-------------+-------------+-------------+-------------+-- 0 1000 2000 3000 4000 x3 2019年11月 4日月曜日 15時26分48秒 115 Plot of resid1*x4. Legend: A = 1 obs, B = 2 obs, etc. | 50 + A | | | A 25 + 残 | A A 差 | A A AA A | AABB A A 0 +------A------A-----------------A-A--------------------------- | A A B A C A A | A A A A | A A A A A -25 + | A | | -50 + | ---+-------------+-------------+-------------+-------------+-- 6 8 10 12 14 x4 2019年11月 4日月曜日 15時26分48秒 116 Plot of resid1*x5. Legend: A = 1 obs, B = 2 obs, etc. | 50 + A | | | A 25 + 残 | A A 差 | A AA A A | A A A A AA A A 0 +-----------------------------------A---A------A-A-------------------- | A AB A A A BA | A A A A | A A A A A -25 + | A | | -50 + | --+----------+----------+----------+----------+----------+----------+- 0 10 20 30 40 50 60 x5 2019年11月 4日月曜日 15時26分48秒 117 Plot of resid1*x6. Legend: A = 1 obs, B = 2 obs, etc. | 50 + A | | | A 25 + 残 | A A 差 | A B AA | A A A AAA A A 0 +-----------------------------------------A--A--A-------A------------- | A A A BA AA A A | A B A | A A A A A -25 + | A | | -50 + | --+----------+----------+----------+----------+----------+----------+- 25 50 75 100 125 150 175 x6 2019年11月 4日月曜日 15時26分48秒 118 Plot of resid1*x1. Legend: A = 1 obs, B = 2 obs, etc. | 50 + A | | | A 25 + 残 | A A 差 | A A A A A | A AA A A B A 0 +----------------AA-----A---A--------------------------------- | AA B B A B A | A A A A | A A A A A -25 + | A | | -50 + | ---+-------------+-------------+-------------+-------------+-- 40 50 60 70 80 x1 2019年11月 4日月曜日 15時26分48秒 119 Plot of resid1*x2. Legend: A = 1 obs, B = 2 obs, etc. | 50 + A | | | A 25 + 残 | A A 差 | A AAA A | BB A B A 0 +--A--A-----------A-------------------------------A----------- | AAAACA A A | A A AA | A A AAA -25 + | A | | -50 + | ---+-------------+-------------+-------------+-------------+-- 0 1000 2000 3000 4000 x2 2019年11月 4日月曜日 15時26分48秒 120 Plot of resid1*x3. Legend: A = 1 obs, B = 2 obs, etc. | 50 + A | | | A 25 + 残 | A A 差 | A B A A | ABA AAA A 0 +---A--A----------------A-------------------------A----------- | AAAAB AAA A | A A B | AA B A -25 + | A | | -50 + | ---+-------------+-------------+-------------+-------------+-- 0 1000 2000 3000 4000 x3 2019年11月 4日月曜日 15時26分48秒 121 Plot of resid1*x4. Legend: A = 1 obs, B = 2 obs, etc. | 50 + A | | | A 25 + 残 | A A 差 | A A AA A | AABB A A 0 +------A------A-----------------A-A--------------------------- | A A B A C A A | A A A A | A A A A A -25 + | A | | -50 + | ---+-------------+-------------+-------------+-------------+-- 6 8 10 12 14 x4 2019年11月 4日月曜日 15時26分48秒 122 Plot of resid1*x5. Legend: A = 1 obs, B = 2 obs, etc. | 50 + A | | | A 25 + 残 | A A 差 | A AA A A | A A A A AA A A 0 +-----------------------------------A---A------A-A-------------------- | A AB A A A BA | A A A A | A A A A A -25 + | A | | -50 + | --+----------+----------+----------+----------+----------+----------+- 0 10 20 30 40 50 60 x5 2019年11月 4日月曜日 15時26分48秒 123 Plot of resid1*x6. Legend: A = 1 obs, B = 2 obs, etc. | 50 + A | | | A 25 + 残 | A A 差 | A B AA | A A A AAA A A 0 +-----------------------------------------A--A--A-------A------------- | A A A BA AA A A | A B A | A A A A A -25 + | A | | -50 + | --+----------+----------+----------+----------+----------+----------+- 25 50 75 100 125 150 175 x6 2019年11月 4日月曜日 15時26分48秒 124 Plot of resid1*y. Legend: A = 1 obs, B = 2 obs, etc. | 50 + A | | | A 25 + 残 | AA 差 | A AAA A | AA AAA A B 0 +---------------A--A-A------------------------------------A------- | ABAA B A A A | A B A | B A A A -25 + | A | | -50 + | ---+---------+---------+---------+---------+---------+---------+-- 0 20 40 60 80 100 120 y 2019年11月 4日月曜日 15時26分48秒 125 REG プロシジャモデル : MODEL1 従属変数 : y R2 乗選択法読み込んだオブザベーション数 41 使用されたオブザベーション数 41 取り込んだ変数の数 R2 乗モデルの独立変数 1 0.4157 x2 1 0.2438 x3 1 0.1880 x1 1 0.1366 x6 1 0.0090 x4 1 0.0029 x5 ----------------------------------------------------- 2 0.5863 x2 x3 2 0.5161 x1 x2 2 0.4981 x2 x6 2 0.4214 x2 x5 2 0.4194 x2 x4 2 0.4066 x1 x3 2 0.3657 x3 x6 2 0.2483 x3 x5 2 0.2458 x1 x5 2 0.2439 x3 x4 2 0.2291 x1 x6 2 0.1917 x1 x4 2 0.1587 x5 x6 2 0.1378 x4 x6 2 0.0120 x4 x5 ----------------------------------------------------- 3 0.6174 x2 x3 x6 3 0.6125 x1 x2 x3 3 0.5930 x2 x3 x5 3 0.5930 x2 x3 x4 3 0.5622 x1 x2 x5 3 0.5452 x1 x2 x6 3 0.5452 x1 x2 x4 3 0.5083 x2 x4 x6 3 0.5047 x2 x5 x6 3 0.4649 x1 x3 x5 3 0.4446 x1 x3 x6 3 0.4320 x1 x3 x4 3 0.4250 x2 x4 x5 3 0.3808 x3 x5 x6 3 0.3702 x3 x4 x6 3 0.2550 x1 x4 x5 3 0.2484 x3 x4 x5 3 0.2462 x1 x5 x6 3 0.2332 x1 x4 x6 3 0.1590 x4 x5 x6 ----------------------------------------------------- 4 0.6396 x1 x2 x3 x5 4 0.6329 x1 x2 x3 x4 4 0.6291 x1 x2 x3 x6 4 0.6285 x2 x3 x4 x6 4 0.6176 x2 x3 x5 x6 4 0.6028 x1 x2 x4 x5 4 0.5997 x2 x3 x4 x5 4 0.5747 x1 x2 x4 x6 4 0.5622 x1 x2 x5 x6 4 0.5164 x2 x4 x5 x6 4 0.5035 x1 x3 x4 x5 4 0.4708 x1 x3 x4 x6 4 0.4649 x1 x3 x5 x6 4 0.3871 x3 x4 x5 x6 4 0.2550 x1 x4 x5 x6 ----------------------------------------------------- 5 0.6685 x1 x2 x3 x4 x5 5 0.6501 x1 x2 x3 x4 x6 5 0.6396 x1 x2 x3 x5 x6 5 0.6290 x2 x3 x4 x5 x6 5 0.6040 x1 x2 x4 x5 x6 5 0.5043 x1 x3 x4 x5 x6 ----------------------------------------------------- 6 0.6695 x1 x2 x3 x4 x5 x6

結果の見方
- フルモデル
- 逐次選択法(stepwise)
  - 変量増減法。
  - 一度取り込まれても、組合わせによっては削除される。
- 総当たり法(rsquare)
  - 説明変数の組合わせ毎の決定係数(R^2)が表示される : 大きい順に
  - モデルの探索用。
  - 決定係数 : R-Square : 1 に近いほど当てはまりが良いと言える
  - 説明変数が増えると大きくなるのが一般的
  - 興味のある組合わせを見つけ出して、このあと計算させる。残差解析も行うこと。
- 他に、前進選択法(forward)、後退選択法(backward)、...
- 「数値計算上の最適モデル」と「その分野の知識からの最適モデル」には違いがあることを知っておくこと。
- 残差解析はいつの場合でも必要
  - 残差の性質 ===> 正規性 : 残差プロット、残差解析
  - 均等に散らばっているか?
  - 傾向はないか? : 傾向があると言うことは正規性の仮定が崩れていること
  - ...
- ...

次回は、... : 11月12日 16:20-17:50
- 主成分分析
- 因子分析
- ...

講義のホームページへ戻ります