重回帰分析、外れ値の処理

統計解析 06 クラス : 第11回 (12/16/04)

　前回は、説明変量が一つである単回帰分析を紹介した。単回帰分析における説明変量が複数になった手法が重回帰分析であり、残差(予測誤差)の二乗和を最小にするという考え方は同じなので、その原理は容易に理解できると期待している。

重回帰分析 : 2変量以上の説明する変量(説明変量)で 1変量(目的変量)を説明

説明変量が複数になる : 単 ===> 重
体重を身長と胸囲で説明したい。予測したい。
[体重]=a[身長]+b[胸囲]+c : 回帰係数を求めたい。
単回帰とアイディアは同じ : 残差(予測誤差)の二乗和を最小にする(最小二乗法)
説明される変量(目的変量)と平行に残差を取る。

プログラム : les1101.sas

/* Lesson 11-1 */ /* File Name = les1101.sas 12/16/04 */ data gakusei; infile 'all04b.prn' firstobs=2; input sex $ shintyou taijyuu kyoui jitaku $ kodukai carryer $ tsuuwa; if sex^='M' & sex^='F' 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; :

出力結果 : les1101.lst
SAS システム 2 19:47 Tuesday, December 7, 2004 Model: MODEL1 Dependent Variable: TAIJYUU Analysis of Variance Sum of Mean Source DF Squares Square F Value Prob>F Model 2 8432.71135 4216.35567 119.383 0.0001 Error 99 3496.47532 35.31793 C Total 101 11929.18667 Root MSE 5.94289 R-square 0.7069 Dep Mean 58.86667 Adj R-sq 0.7010 C.V. 10.09551 SAS システム 3 19:47 Tuesday, December 7, 2004 Parameter Estimates Parameter Standard T for H0: Variable DF Estimate Error Parameter=0 Prob > |T| INTERCEP 1 -112.109751 11.52738966 -9.726 0.0001 SHINTYOU 1 0.694385 0.07197938 9.647 0.0001 KYOUI 1 0.632583 0.08287933 7.633 0.0001 SAS システム 4 19:47 Tuesday, December 7, 2004 S H T K C I A J O A T R N I K I D R S P E T J Y T U R U R S O S Y Y O A K Y U E I B E O U U K A E W D D S X U U I U I R A 1 1 1 F 145.0 38.0 . J 10000 . . . 2 F 146.7 41.0 85 J 10000 Vodafone 6000 43.5262 -2.52616 3 F 148.0 42.0 . J 50000 . . . 4 F 148.0 43.0 80 J 50000 DoCoMo 4000 41.2659 1.73405 5 F 148.9 . . J 60000 . . . 6 F 149.0 45.0 . G 60000 . . . 7 F 150.0 46.0 86 40000 . 46.4502 -0.45022 8 F 151.0 50.0 . G 60000 J-PHONE . . . 9 F 151.7 41.5 80 J 35000 . 43.8352 -2.33517 10 F 152.0 35.0 77 J 60000 DoCoMo 2000 42.1457 -7.14574 11 F 152.0 43.0 . J 20000 au 3500 . . 12 F 152.0 44.0 . 45000 DoCoMo 4000 . . 13 F 153.0 41.0 . J 125000 No . . . 14 F 153.0 42.0 . G 0 Vodafone 1000 . . 15 F 153.0 46.5 87 G 10000 . 49.1660 -2.66596 SAS システム 6 19:47 Tuesday, December 7, 2004 プロット : TAIJYUU*SHINTYOU. 凡例: A = 1 OBS, B = 2 OBS, ... 100 + A | A A TAIJYUU | A A A | B BABAB AACAA A B A AA | A A B A B ABA BAFBC ABA AABBA 50 + A A ADA BBE C BBACB A | A A B A A | | | 0 + --+-----------+-----------+-----------+-----------+-----------+- 140 150 160 170 180 190 SHINTYOU SAS システム 7 19:47 Tuesday, December 7, 2004 プロット : TAIJYUU*KYOUI. 凡例: A = 1 OBS, B = 2 OBS, ... 100 + A | A A TAIJYUU | AA A | A A C BBF BABA A A | A A C C AAE FBI AAA A 50 + A A AA B ICHBBA | A A B B | | | 0 + ---+-------+-------+-------+-------+-------+-------+-------+-- 50 60 70 80 90 100 110 120 KYOUI SAS システム 8 19:47 Tuesday, December 7, 2004 プロット : TAIJYUU*PRED1. 凡例: A = 1 OBS, B = 2 OBS, ... 100 + A | A A TAIJYUU | A A A | A B ACBB BA ABAABA A | A B BA BAA ABBAACFBABAA AB 50 + AA B CBBADCACDB B | ABAA A | | | 0 + --+-----------+-----------+-----------+-----------+-----------+- 40 50 60 70 80 90 Predicted Value of TAIJYUU SAS システム 9 19:47 Tuesday, December 7, 2004 プロット : RESID1*PRED1. 凡例: A = 1 OBS, B = 2 OBS, ... | R 40 + e | s | A i 20 + A d | A A A A u | A B AA B A B BBA A A A a 0 +----AAAA--B-CBB-CCAAABA-ABBAACEB-AA--B-ABA---------A------------- l | A AAAA BCA B A AA BAAAAC A | A -20 + ---+-----------+-----------+-----------+-----------+-----------+-- 40 50 60 70 80 90 Predicted Value of TAIJYUU SAS システム 10 19:47 Tuesday, December 7, 2004 プロット : RESID1*SHINTYOU. 凡例: A = 1 OBS, B = 2 OBS, ... | R 40 + e | s | A i 20 + A d | A A A A u | B B A BABAB ABAA a 0 +----------A-A-A-A-AAACA-BAD-B-BABBC-AACAB-BAC-A-C-A-A--A--------- l | A A AA B AA CB A BAA A ACAA A | A -20 + ---+-----------+-----------+-----------+-----------+-----------+-- 140 150 160 170 180 190 SHINTYOU SAS システム 11 19:47 Tuesday, December 7, 2004 プロット : RESID1*KYOUI. 凡例: A = 1 OBS, B = 2 OBS, ... | R 40 + e | s | A i 20 + A d | A A A A u | B A B A C ABD B a 0 +-----------------------C-A-D-CCDKCBAG-A-BC---B--------A---------- l | AA B BA FACBC A B A | A -20 + -+--------+--------+--------+--------+--------+--------+--------+- 50 60 70 80 90 100 110 120 KYOUI SAS システム 12 19:47 Tuesday, December 7, 2004 プロット : RESID1*TAIJYUU. 凡例: A = 1 OBS, B = 2 OBS, ... | R 40 + e | s | A i 20 + A d | AA A A u | A AAAAB B BBBA AA a 0 +----------------BAACCBFCAC-CFDBCAA-E----A---------------- l | A A BDABB A ADAAD A | A -20 + ---+------------+------------+------------+------------+-- 20 40 60 80 100 TAIJYUU SAS システム 13 19:47 Tuesday, December 7, 2004 Univariate Procedure Variable=RESID1 Residual Moments N 102 Sum Wgts 102 Mean 0 Sum 0 Std Dev 5.883755 Variance 34.61857 Skewness 1.217662 Kurtosis 2.157941 USS 3496.475 CSS 3496.475 CV . Std Mean 0.582579 T:Mean=0 0 Pr>|T| 1.0000 Num ^= 0 102 Num > 0 39 M(Sign) -12 Pr>=|M| 0.0223 Sgn Rank -341.5 Pr>=|S| 0.2563 W:Normal 0.922178 Pr SAS システム 17 19:47 Tuesday, December 7, 2004 Univariate Procedure Variable=RESID1 Residual Stem Leaf # Boxplot 2 4 1 * 1 8 1 0 1 01224 5 0 0 5567777778888 13 | 0 0000111112223334444 19 +--+--+ -0 44444333333333332222222222222111111100000 41 *-----* -0 998777776666665555555 21 | -1 0 1 | ----+----+----+----+----+----+----+----+- Multiply Stem.Leaf by 10**+1 SAS システム 18 19:47 Tuesday, December 7, 2004 Univariate Procedure Variable=RESID1 Residual Normal Probability Plot 22.5+ * | * | ****++++++ | *******+ | +++******* | ************* | * *********** -12.5+*++++++++ +----+----+----+----+----+----+----+----+----+----+ -2 -1 0 +1 +2

結果の見方
- 対象になったのは 102名。
- 説明変量群が予測に役立っているか?
  - 回帰に役立っているか : Prob>F : 小さいと有意
  - 「役立っている」と言える : 0.01% だから 1% で有意
- 決定係数 : R-Square ( 相関係数 : R )
  - 目的変量が説明変量でどの程度説明しているかの割合。
  - 1 に近いほど当てはまりが良いと言える。: 70.7%
  - 説明変量数が増えると大きくなるのが一般的。
- 回帰係数 : Parameter Estimate
  - 回帰式: a=0.694, b=0.633, c=-112
- ある特定の説明変量が予測に役立っているか?
  - 回帰係数の検定(帰無仮説:係数=0 か?) : Prob>|T| : 小さいと有意
  - 両方とも(身長も胸囲も)有意
  - 「各係数は 0ではない」と言える : 0.01% だから 1% で有意
- 残差の性質 ===> 正規性 : 残差プロット、残差解析
  - 残差(予測誤差)は正規分布をしていると仮定してモデルが構築されている。
  - この仮定が覆ると、回帰分析として成立していないことになる。
  - 残差が正規分布をしているか確認する必要がある。
  - 均等に散らばっているか?
  - 傾向はないか? : もし傾向があると言うことになれば正規性の仮定が崩れている
  - 体重の大きい 2例程度が外れ値と考えられるか要確認 ===> [演習](第3節)
  - ...
- ...

特定グループでの解析

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

プログラム : les1102.sas

/* Lesson 11-2 */ /* File Name = les1102.sas 12/16/04 */ data gakusei; infile 'all04b.prn' firstobs=2; 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;

出力結果 : les1102.lst
SAS システム 2 19:47 Tuesday, December 7, 2004 Correlation Analysis 5 'VAR' Variables: SHINTYOU TAIJYUU KYOUI KODUKAI TSUUWA Simple Statistics Variable N Mean Std Dev Sum Minimum Maximum SHINTYOU 64 172.5 6.1180 11041.1 156.0 185.0 TAIJYUU 64 64.6828 9.0923 4139.7 46.0000 100.0 KYOUI 64 88.5000 8.6189 5664.0 56.0000 112.0 KODUKAI 60 55100.0 57720.5 3306000 0 300000 TSUUWA 8 9250.0 3453.8 74000.0 5000.0 15000.0 SAS システム 3 19:47 Tuesday, December 7, 2004 Correlation Analysis Pearson Correlation Coefficients / Prob > |R| under Ho: Rho=0 / Number of Observations SHINTYOU TAIJYUU KYOUI KODUKAI TSUUWA SHINTYOU 1.00000 0.41073 0.19351 0.10468 -0.18380 0.0 0.0007 0.1255 0.4261 0.6631 64 64 64 60 8 TAIJYUU 0.41073 1.00000 0.66124 -0.05998 0.15571 0.0007 0.0 0.0001 0.6489 0.7127 64 64 64 60 8 KYOUI 0.19351 0.66124 1.00000 -0.12114 -0.16606 0.1255 0.0001 0.0 0.3565 0.6943 64 64 64 60 8 KODUKAI 0.10468 -0.05998 -0.12114 1.00000 0.39786 0.4261 0.6489 0.3565 0.0 0.3290 60 60 60 60 8 TSUUWA -0.18380 0.15571 -0.16606 0.39786 1.00000 0.6631 0.7127 0.6943 0.3290 0.0 8 8 8 8 8 SAS システム 6 19:47 Tuesday, December 7, 2004 Model: MODEL1 Dependent Variable: TAIJYUU Analysis of Variance Sum of Mean Source DF Squares Square F Value Prob>F Model 2 2709.84404 1354.92202 33.082 0.0001 Error 61 2498.32705 40.95618 C Total 63 5208.17109 Root MSE 6.39970 R-square 0.5203 Dep Mean 64.68281 Adj R-sq 0.5046 C.V. 9.89398 SAS システム 7 19:47 Tuesday, December 7, 2004 Parameter Estimates Parameter Standard T for H0: Variable DF Estimate Error Parameter=0 Prob > |T| INTERCEP 1 -67.062537 23.09102487 -2.904 0.0051 SHINTYOU 1 0.436590 0.13432815 3.250 0.0019 KYOUI 1 0.637583 0.09535081 6.687 0.0001 SAS システム 10 19:47 Tuesday, December 7, 2004 プロット : TAIJYUU*SHINTYOU. 凡例: A = 1 OBS, B = 2 OBS, ... TAIJYUU | 100 + A | A A | 75 + A A A A A AA | B B A C A A A C A A D A A A | A A A B A B A D B C A AAA A A AA A 50 + A B A | | 25 + --+---------+---------+---------+---------+---------+---------+- 155 160 165 170 175 180 185 SHINTYOU SAS システム 11 19:47 Tuesday, December 7, 2004 プロット : TAIJYUU*KYOUI. 凡例: A = 1 OBS, B = 2 OBS, ... TAIJYUU | 100 + A | A A | 75 + AA BA A A | A A C BAH BAAB A | A A B C AAC EBE AA A 50 + A A A A | | 25 + ---+-------+-------+-------+-------+-------+-------+-------+-- 50 60 70 80 90 100 110 120 KYOUI SAS システム 12 19:47 Tuesday, December 7, 2004 プロット : TAIJYUU*PRED1. 凡例: A = 1 OBS, B = 2 OBS, ... TAIJYUU | 100 + A | A A | 75 + AA AAAA A | AAAA DABBBABB AAA | A B AA CAAACDDAABA 50 + A A A A | | 25 + --+-----------+-----------+-----------+-----------+-----------+- 40 50 60 70 80 90 Predicted Value of TAIJYUU SAS システム 13 19:47 Tuesday, December 7, 2004 プロット : RESID1*PRED1. 凡例: A = 1 OBS, B = 2 OBS, ... | R 40 + e | s | i 20 + A A d | A A u | A A AAAA CBA AA a 0 +---------------AA--A----A-BAA-BBBABB--B-BA---------A------------- l | A A A AADBAABBB AA | -20 + ---+-----------+-----------+-----------+-----------+-----------+-- 40 50 60 70 80 90 Predicted Value of TAIJYUU SAS システム 14 19:47 Tuesday, December 7, 2004 プロット : RESID1*SHINTYOU. 凡例: A = 1 OBS, B = 2 OBS, ... | R 40 + e | s | i 20 + A A d | A A u | A B A B A B A B AA a 0 +----A-------A-----------A-B---A-C-A-A-A--AC---A-BA--B-A-A---A---- l | A B A A B C A B A A A BA A A | -20 + ---+---------+---------+---------+---------+---------+---------+-- 155 160 165 170 175 180 185 SHINTYOU SAS システム 15 19:47 Tuesday, December 7, 2004 プロット : RESID1*KYOUI. 凡例: A = 1 OBS, B = 2 OBS, ... | R 40 + e | s | i 20 + A A d | A A u | A B A B BD B a 0 +------------------A----B-A-B---ABABAE-A-AC---A--------A---------- l | A B BBBAE AAB A A | -20 + -+--------+--------+--------+--------+--------+--------+--------+- 50 60 70 80 90 100 110 120 KYOUI SAS システム 16 19:47 Tuesday, December 7, 2004 プロット : RESID1*TAIJYUU. 凡例: A = 1 OBS, B = 2 OBS, ... | R 40 + e | s | i 20 + A A d | A A u | A A B AAAC A A AA a 0 +----------A------AA--ADACA-DA-A-CB------A------------------------ l | A A BA B FABABA A | -20 + ---+---------+---------+---------+---------+---------+---------+-- 40 50 60 70 80 90 100 TAIJYUU SAS システム 17 19:47 Tuesday, December 7, 2004 Univariate Procedure Variable=RESID1 Residual Moments N 64 Sum Wgts 64 Mean 0 Sum 0 Std Dev 6.2973 Variance 39.65598 Skewness 1.20336 Kurtosis 1.801302 USS 2498.327 CSS 2498.327 CV . Std Mean 0.787162 T:Mean=0 0 Pr>|T| 1.0000 Num ^= 0 64 Num > 0 26 M(Sign) -6 Pr>=|M| 0.1686 Sgn Rank -130 Pr>=|S| 0.3889 W:Normal 0.913681 Pr SAS システム 20 19:47 Tuesday, December 7, 2004 Univariate Procedure Variable=RESID1 Residual Stem Leaf # Boxplot 2 2 1 0 1 8 1 0 1 24 2 | 0 55556677789 11 | 0 00001133344 11 +--+--+ -0 44433333322221111100 20 *-----* -0 998887665555555555 18 +-----+ ----+----+----+----+ Multiply Stem.Leaf by 10**+1 SAS システム 21 19:47 Tuesday, December 7, 2004 Univariate Procedure Variable=RESID1 Residual Normal Probability Plot 22.5+ * | * ++ | **++++++ 7.5+ ++****** | +++******* | ********** -7.5+ * * ********** +----+----+----+----+----+----+----+----+----+----+ -2 -1 0 +1 +2

結果の見方
- 単変量毎の相関が有意なのは、身長と体重、体重と胸囲の間。
- 対象になったのは 64名。
- 回帰に役立っているか : 役立っている : 0.01% だから 1% で有意
- 決定係数(R-square)は 0.520
- 個々の説明変量が予測に役立っているか?
  - 係数がゼロか? : 定数項も身長も胸囲も有意(1% で有意)
- 残差の性質 ===> 正規性 : 残差プロット、残差解析
  - 均等に散らばっているか?
  - 傾向はないか? : 傾向があると言うことは正規性の仮定が崩れていること
  - 外れ値? 85Kg より重い 2名程度が吟味対象?

[演習1] : 「男性のみ」で、かつ「体重の大きい 3名を除外」して実行してみよ。
- プログラム : les1103.sas、出力結果 : les1103.lst
  where sex='M' and taijyuu<85;
- 当てはまりは良くなったか? : 異常値と外れ値の意味するもの
- 残差の正規性はどのように変化したか?
[演習2] : 「男性のみ」で、かつ「体重の大きい 3名を除外」の集団に対して、「胸囲」の代わりに「胸囲の2乗」を使って実行してみよ。 ===> [おまけ1]参照のこと
- プログラム : les1104.sas、出力結果 : les1104.lst
  kyo2=kyoui**2;
回帰分析における変数選択 :
- 今までのモデル : 指定した説明変数全てを取り込んで計算 : フルモデル
- 影響力の大きい説明変数だけをモデルに取り込みたい。
- たくさんの説明変数の中から影響力の大きい、いくつかの変数を選別する。
- 幾つかの選択方法がある : 逐次増減法(stepwise)、前進選択法(forward)、後退選択法(backward)、...
次回は、... : 01月13日 14:45
- 主成分分析
- レポート課題提示
- ...
[おまけ1] 変数変換 : 新しい変量の算出 : [注意] 以下はあくまでも計算できることの例です。
以下に示したような演算子や関数を使って、新しい変量を生成することができる。利用可能なものの一部を掲載しておく。
data gakusei; infile 'all04b.prn' firstobs=2; input sex $ shintyou taijyuu kyoui jitaku $ kodukai carryer $ tsuuwa; dekasa=shintyou+taijyuu+kyoui; : 変量間の加減乗除 kyo_2=kyoui**2; : 二乗 kyo_sr=sqrt(kyoui); : ルート
[算術演算子]
+ : 足し算を算出します。
- : 引き算を算出します。
* : 掛け算を算出します。
/ : 割り算を算出します。
** : 巾乗を算出します。
sqrt : 平方根(ルート)を算出します。
[数値関数]
arcos : 逆余弦(アークコサイン)を算出します。
arsin : 逆正弦(アークサイン)を算出します。
atan : 逆正接(アークタンジェント)を算出します。
cos : 三角関数の余弦(コサイン)を求めます。
cosh : 双曲線の余弦（コサイン）を求めます。
sin : 正弦（サイン）を算出します。
sinh : 双曲線正弦を算出します。
tan : 正接（タンジェント）を算出します。
tanh : 双曲線正接の値を算出します。

ceil : 引数より大きいかまたは等しい整数のうち最小の値を戻します。
floor : 引数より小さいかまたは等しい整数のうち最大のものを戻します。
fuzz : 引数と最も近い整数との差が10^-12以内のときに、その整数を戻します。
int : 小数部を切り捨てて、整数値を戻します。
round : 四捨五入します。
tranc : 指定された長さに切り捨てた数値を戻します。
[おまけ2] 1レコード(1行)のレコード長が長い場合 : 一行の長さの指定、例えば 230バイトだと
[おまけ3] 行末がそろってないデータの読み込み : 一行の長さを指定し、揃ってないことを明示
[おまけ4] 固定長データの読み込み : カラム位置の指定
[おまけ5] データファイルの先頭を読み飛ばす : コメント等が挿入されているような場合。先頭 3行を読み飛ばす場合は
[おまけ6] データの生成 : 乱数で(正規乱数) : 種(seed)の指定。種は任意の数値(普通は整数)。
[参考] : 乱数を生成する関数

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