二変量の関係、単回帰分析

統計解析 03 クラス : 第09回 (12/04/03)

　ここまで、SAS の使い方と，単変量(一変量)を取り扱う統計手法を主に紹介してきた。今後は、二変量以上、つまり、多変量解析を紹介していくことにする。まずは二変量の関係を説明する方法について紹介する。

複数変量の関係
- 身長と体重の関係 : どのように表現すれば良いのだろうか?
- 数値 : 相関係数、回帰分析 <=== 単変量の場合は「基礎統計量」であった。
- 図 : 散布図 <=== 単変量の場合は「ヒストグラム」であった。
  [例] 「「インターネット人口」をどのように読んだらいいか(Impressの記事, 2000/07/28) 」内の図
  - 情報通信白書 :2003年7月発表(平成15年度版)の統計, 国内のインターネット人口
    - インターネットの着実な普及 : ネット利用者は6942万人
    - 概要 (PDF)

散布図と相関係数

どの様に分布しているのだろう?
分布の仕方を数値的に表現すると?

プログラム : les0901.sas

/* Lesson 09-1 */ /* File Name = les0901.sas 12/04/03 */ data gakusei; infile 'all03b.prn' firstobs=2; input sex $ height weight chest jitaku $ kodukai carrier $ tsuuwa; proc print data=gakusei(obs=10); run; : proc plot data=gakusei; : 散布図を描く plot height*weight; : 散布図の変量を指定(縦軸、横軸の順) plot weight*height; : run; : : proc corr data=gakusei; : 相関係数(相関行列)を計算 run; :

出力結果 : les0901.lst
SAS システム 2 11:25 Thursday, November 20, 2003 プロット : HEIGHT*WEIGHT. 凡例: A = 1 OBS, B = 2 OBS, ... (NOTE: 36 オブザベーションが欠損値です．) HEIGHT | 200 + | | B A 180 + A ADCDDCBCA A B A A | AACDJHRMGFCDDCB BA | ADAGHEEDCBBDAA A A 160 + ADBDBHCCAABB | A D DCCA A A | A AAA 140 + ---+-----------+-----------+-----------+-----------+-- 20 40 60 80 100 WEIGHT SAS システム 3 11:25 Thursday, November 20, 2003 プロット : WEIGHT*HEIGHT. 凡例: A = 1 OBS, B = 2 OBS, ... (NOTE: 36 オブザベーションが欠損値です．) 100 + A | A A WEIGHT | A A B A A | B CBDDC CBEAC CBD B AA | A AA C B CABBF HBMHKBIFFCC BADBB A 50 + AAA CABCA CBF F EBBFF DAAAB A | A B B A BA | | | 0 + --+-----------+-----------+-----------+-----------+-----------+-- 140 150 160 170 180 190 HEIGHT SAS システム 4 11:25 Thursday, November 20, 2003 Correlation Analysis 5 'VAR' Variables: HEIGHT WEIGHT CHEST KODUKAI TSUUWA Simple Statistics Variable N Mean Std Dev Sum Minimum Maximum HEIGHT 255 168.0 8.0840 42846.9 145.0 186.0 WEIGHT 229 58.7852 9.0852 13461.8 35.0000 100.0 CHEST 87 86.9080 7.9864 7561.0 56.0000 112.0 KODUKAI 241 51566.4 52037.2 12427500 0 300000 TSUUWA 48 7961.5 5011.7 382150 200.0 30000.0 SAS システム 5 11:25 Thursday, November 20, 2003 Correlation Analysis Pearson Correlation Coefficients / Prob > |R| under Ho: Rho=0 / Number of Observations HEIGHT WEIGHT CHEST KODUKAI TSUUWA HEIGHT 1.00000 0.70774 0.36532 0.02347 0.05957 0.0 0.0001 0.0005 0.7204 0.6875 255 229 87 235 48 WEIGHT 0.70774 1.00000 0.66125 -0.03849 0.11510 0.0001 0.0 0.0001 0.5773 0.4515 229 229 87 212 45 CHEST 0.36532 0.66125 1.00000 -0.11120 0.53420 0.0005 0.0001 0.0 0.3169 0.0905 87 87 87 83 11 KODUKAI 0.02347 -0.03849 -0.11120 1.00000 -0.11297 0.7204 0.5773 0.3169 0.0 0.4600 235 212 83 241 45 TSUUWA 0.05957 0.11510 0.53420 -0.11297 1.00000 0.6875 0.4515 0.0905 0.4600 0.0 48 45 11 45 48

結果の見方
- 縦軸と横軸の該当部分が交差したところにマークを付置
- データが1つなら「Aマーク」、2つなら「Bマーク」、...
- データ全体がどこに分布しているかが判る
- 縦軸と横軸を交換するだけで印象が異なる
- 各変量の平均値との比較
- 外れ値(Outlier)を見つける <===> 異常値
- サンプルサイズ、平均、標準偏差、最大値、最小値 <=== proc means だけでなく proc corr でも得られる。
- 相関係数(R) / 仮説「相関係数(R)=0」の起る確率 / サンプルサイズ
- -1 ≦ 相関係数(R)≦ 1
- R＝0 : 無相関。R＞0 : 正の相関、右肩上がり。R＜0 : 負の相関、右肩下がり。
- 相関係数(R)が 0 かの検定 : 値が小さいと有意(相関係数が 0 ではない)
  この例 : 身長と体重、身長と胸囲、体重と胸囲の間には有意な関係があると言える(5%, 1%)。
[注意] 相関行列は細切れに表示されるので、不要部分を削除することによって整形しレポート等に使うこと。

単回帰分析 : 予測等に使う、連続変量の関係

体重を身長で説明(回帰)したい : [体重]=a[身長]+b : 回帰係数
説明される変量と説明する変量 : 説明する変数が一つ＝「単」 <===>「重」
説明される変量 : 目的変数、従属変数、dependent variable
説明する変量 : 説明変数、独立変数、independent variable
どうやって直線を決める? : 予測誤差の2乗和を最小にする

プログラム : les0902.sas
/* Lesson 09-2 */ /* File Name = les0902.sas 12/04/03 */ data gakusei; infile 'all03b.prn' firstobs=2; input sex $ height weight chest jitaku $ kodukai carrier $ tsuuwa; proc print data=gakusei(obs=10); run; : proc reg data=gakusei; : 回帰分析 model weight=height; : 変量を指定 output out=outreg1 predicted=pred1 residual=resid1; : 結果項目の保存 run; : : proc print data=outreg1(obs=15); : まずは表示 run; : : proc plot data=outreg1; : 散布図を描く plot weight*height/vaxis=20 to 100 by 20; : 体重と身長(縦軸指定) plot pred1*weight; : 予測値と観測値 plot resid1*pred1/vref=0; : 残差と予測値(残差解析)(水平軸指定) plot resid1*height/vref=0; : 残差と説明変数(残差解析) plot resid1*weight/vref=0; : 残差と目的変数(残差解析) run;
[補足] proc plot の下に以下の行を追加した方がより正確ではある。欠損値を含むデータを解析対象から除外する事を指示する命令文である。「欠損値です」の表示が無くなるだけで、得られる図は同じ(欠損値は描画できないから)。試しに追加する/しないの両方で実行してみよ。
where weight^=. and height^=.;

出力結果 : les0902.lst
SAS システム 2 11:25 Thursday, November 20, 2003 Model: MODEL1 Dependent Variable: WEIGHT Analysis of Variance Sum of Mean Source DF Squares Square F Value Prob>F Model 1 9426.48780 9426.48780 227.815 0.0001 Error 227 9392.78172 41.37789 C Total 228 18819.26952 Root MSE 6.43257 R-square 0.5009 Dep Mean 58.78515 Adj R-sq 0.4987 C.V. 10.94250 SAS システム 3 11:25 Thursday, November 20, 2003 Parameter Estimates Parameter Standard T for H0: Variable DF Estimate Error Parameter=0 Prob > |T| INTERCEP 1 -77.661171 9.05004324 -8.581 0.0001 HEIGHT 1 0.808135 0.05354181 15.094 0.0001 SAS システム 4 11:25 Thursday, November 20, 2003 K C H W J O A T R E E C I D R S P E I I H T U R U R S O S G G E A K I U E I B E H H S K A E W D D S X T T T U I R A 1 1 1 F 145.0 38.0 . J 10000 . 39.5184 -1.5184 2 F 148.0 42.0 . J 50000 . 41.9428 0.0572 3 F 148.0 43.0 80 J 50000 DoCoMo 4000 41.9428 1.0572 4 F 148.9 . . J 60000 . 42.6701 . 5 F 149.0 45.0 . G 60000 . 42.7509 2.2491 6 F 150.0 46.0 86 40000 . 43.5590 2.4410 7 F 151.0 50.0 . G 60000 J-PHONE . 44.3672 5.6328 8 F 151.7 41.5 80 J 35000 . 44.9329 -3.4329 9 F 152.0 35.0 77 J 60000 DoCoMo 2000 45.1753 -10.1753 10 F 153.0 41.0 . J 125000 No . 45.9835 -4.9835 11 F 153.0 46.5 87 G 10000 . 45.9835 0.5165 12 F 153.0 50.0 . G 70000 DoCoMo 10000 45.9835 4.0165 13 F 153.0 55.0 78 J 30000 . 45.9835 9.0165 14 F 153.0 . . G 120000 DoCoMo 200 45.9835 . 15 F 153.5 46.0 . J 30000 J-PHONE 8000 46.3875 -0.3875 SAS システム 6 11:25 Thursday, November 20, 2003 プロット : WEIGHT*HEIGHT. 凡例: A = 1 OBS, B = 2 OBS, ... (NOTE: 36 オブザベーションが欠損値です．) WEIGHT | 100 + A | A A 80 + A A B A A | B CBDDC CBEAC CBD B AA 60 + A AA C B CABBF HBMHKBIFFCC BADBB A | AAA CABCA CBF F EBBFF DAAAB A 40 + A B B A BA | 20 + | --+-----------+-----------+-----------+-----------+-----------+-- 140 150 160 170 180 190 HEIGHT SAS システム 7 11:25 Thursday, November 20, 2003 プロット : PRED1*WEIGHT. 凡例: A = 1 OBS, B = 2 OBS, ... (NOTE: 36 オブザベーションが欠損値です．) 80 + | PRED1 | A B A | A ADAAFAB D A A A A | A BBBKFDD GAA A BB 60 + BEBHHGGJBGBAADAB A | AE EHBF BCAAC | BBBCCDAB AAA | BA BBABA A A | A B AB B A 40 + A AA ---+------------+------------+------------+------------+-- 20 40 60 80 100 WEIGHT SAS システム 8 11:25 Thursday, November 20, 2003 プロット : RESID1*PRED1. 凡例: A = 1 OBS, B = 2 OBS, ... (NOTE: 36 オブザベーションが欠損値です．) | R 50 + e | s | A A i 25 + d | A A A A AB A u | A A A A BBAA BBBCDCCAB AA A A a 0 +-------------A--BAA-ABBBCABBF-DEBBCGIBKHHHFIDBBF-A-AA------------ l | AA BAA B BA AFDCACCEB EBBAABBA | A -25 + ---+-----------+-----------+-----------+-----------+-----------+-- 30 40 50 60 70 80 Predicted Value of WEIGHT SAS システム 9 11:25 Thursday, November 20, 2003 プロット : RESID1*HEIGHT. 凡例: A = 1 OBS, B = 2 OBS, ... (NOTE: 36 オブザベーションが欠損値です．) | R 50 + e | s | A A i 25 + d | A A A A A B A u | A A A A B BAA B BBCDC BBB B A A a 0 +--------A---BAA-A-CABCA-BBF-D-EBBCG-IBKFIAHFGBD-BBF-A--AA-------- l | A A BA A B B A AEE CACCDAB CBB BABAB A | A -25 + ---+-----------+-----------+-----------+-----------+-----------+-- 140 150 160 170 180 190 HEIGHT SAS システム 10 11:25 Thursday, November 20, 2003 プロット : RESID1*WEIGHT. 凡例: A = 1 OBS, B = 2 OBS, ... (NOTE: 36 オブザベーションが欠損値です．) | R 50 + e | s | A A i 25 + d | A B BB A u | A AAABBAHBDBFAB A A a 0 +--------------A-BAADCDJHDJGISOFKAH-C--------------------- l | A CABCF CKAHCCECAA | A -25 + ---+------------+------------+------------+------------+-- 20 40 60 80 100 WEIGHT

結果の見方
- 説明変量が予測に役立っているか?
  - 回帰に役立っているか : Prob>F : 小さいと有意(役立っている)
    [この例] 1% 未満(0.01%) なので役に立っていると言える。
- 決定係数 : R-Square ( 相関係数 : R )
  - 目的変量が説明変量でどの程度説明しているかの割合。
  - 1 に近いほど当てはまりが良いと言える。
    [この例] 50% 程(半分, 50.09)を説明できている。
- 回帰係数 : Parameter Estimate
  [この例] a=0.808, b=-77.7
  説明変数が予測に役立っているか?
  回帰係数の検定(係数=0 か?) : Prob>|T| : 小さいと有意(ゼロではないと言える)
  [この例] 両者とも 1% 未満(0.01%) なので回帰係数はゼロではない(何らかの意味がある数字と言える)。
  残差の性質 ===> 正規性 : 残差プロット、残差解析
  
  残差(予測誤差)は正規分布をしていると仮定してモデルが構築されている。
  この仮定が覆ると、回帰分析として成立していないことになる。
  残差が正規分布をしているか確認する必要がある。
  均等に散らばっているか?
  傾向はないか?
  ...
  [この例] 残差には概ね傾向は見られない。ただし体重の大きい 3例程度は要確認。場合によっては外れ値として除外も。
[注意] 誤差は「説明変量」の軸と垂直に取ることに注意せよ。誤差は測定時に混入していると考えてモデルが構築されているから。

[先取り説明] : 12月18日に説明します
　今回は単回帰分析を紹介した。SAS プログラムでは「model weight=height;」と書けば良いことが解ってもらえたと思う。一つの変量を別のいくつか(複数)の変量で説明しようと考えると、これが重回帰分析となり、SAS プログラム上は「model weight=height chest;」のように、model 式の右辺に変量名をズラズラと書くだけで実現できる。興味があれば実行してみよ。
[発展的演習] : 12月11日にでもやってみてください
　各自が収集したデータ、もしくは、「 The Data and Story Library 」(http://lib.stat.cmu.edu/DASL/)で公開されているデータの中から、回帰分析(Regression)に適用することが適当と思われるデータを一つ選んで、解析してみよ。
　後者について具体的に手順を書いておく。
1. 「 The Data and Story Library 」を開き「 List all methods 」のタイルを選択すると、統計手法の一覧が表示される。
2. その中の「 Regression 」を選ぶと、「回帰分析に適用することが適当と思われるデータ」の一覧が表示される(50ほど表示されるであろう)。
3. データに付けられた表題や説明文から、自分の興味に合致したデータを見つけよ。
4. 特定のデータを決めたら、そのページの「Datafile Name」欄をクリックし、データの説明の最後に記述されている「The Data」の部分を「カット & ペースト」で切り出して Windows のファイルとして保存する。ファイル名は自分の覚えやすいものを付ける。
5. 保存したファイルを stat システムに転送して回帰分析等、各種の解析を SAS で実行する。
次回は、... : 12月18日 14:45 (今年最後)
- 多変量解析 : 重回帰分析、...

おまけ : 単変量、二変量を視覚的に捉えると? by Mathematica

$1 dim. Normal Distribution$ [式(a)] 1次元正規分布 N(0,1)
$2 dim. Normal Distribution$ [式(b)] 2次元正規分布 N({0,0},{1,1}, ρ=0.0)
$2 dim. Normal Distribution$ [式(c)] 2次元正規分布 N({0,0},{1,1}, ρ=0.7)
$2 dim. Normal Distribution$ [式(d)] 2次元正規分布 N({0,0},{1,1}, ρ=0.7)、y=1 で切り出し
$2 dim. Normal Distribution$ [式(e)] 2次元正規分布 N({0,0},{1,1}, ρ=0.7)、x+y=2 で切り出し

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