前回の残り、回帰分析(前編)

データサイエンス : 第3回 (10/22/18)

　今回は、前回残した部分を説明した後、多変量解析の代表的な手法である回帰分析について解説する。
　工学系や農学系の実験等を行う領域では頻繁に使用される手法であるが、日常的な話題の中でも概念は広く利用されているので、取っ付き易い手法ではないだろうか。過去のデータからその構造を把握し、新規に測定されたデータに対する予測を行ないたいと言うときに、回帰分析は有用である。構造のシンプルな単回帰分析でこの手法の原理を理解し、複数の説明変量を用いた重回帰分析に拡張する。残差の取り方や、その二乗和を最少にするという考えは同じである。
　ここでは SAS プログラムとその出力を例示しながら説明するが、対象としたデータも掲載しておくので、ご自身の手慣れたプログラム言語を使って各自実行してみてほしい。

● 前回やり残したところ

[利用データ] 講義中の例示として解析対象として取り扱うデータ : 体格 & 小遣いデータ : all07au.txt

今までの講義で収集してきた 15回分をマージ。8変量、375ケース。
性別(M/F)、身長、体重、胸囲、仕送り/小遣いの別(G/J)、仕送り/小遣い額、携帯電話会社名、月平均通話料
1行目は、各変量の説明が入っているので、計算時には読み飛ばす必要がある。

【敢えての注意喚起】データ収集について [再掲&最終]
　第1回(10月1日)でも、第2回(10月15日)でもお伝えしましたが、各自で収集してきたデータに対して多変量解析を適用して分析し、それを報告してもらうことを考えています。よって、各自でデータを見つけてきてもらわないとレポートが書けません。
　これも既にお伝えしているように、この講義は前半と後半に分かれており、前半部で1回評価が入ります。課題を課す段階で提出期限を設けますが、指定した期日を過ぎたものは評価対象としませんのでご注意下さい。初回のアンケートへの回答状況を観ていて若干不安を感じますので予めお伝えしておきます。
　参考までに「多変量解析に適用する」ということは、余りサンプル数が少ない(少数例)や変量数の少ないデータだと適用できなかったり、適用しても興味深い知見が得られないことが多いので、収集する際に注意してください。加えて、集約されたデータではなく、個票データでないと細かな分析できません。
　引用元に取り立てての制限はないので、身の回りにある新聞、雑誌、書籍、ご自身の調査等いろいろなところから持ち込んでくれればと思うが、例えばWeb上にもデータを収録・掲載しているサイトがあり、以下のようなものもその取っ掛かりになるのではないか。
- 政府統計ポータルサイト
- ユーザー提供サンプルデータ: jmp限定
【蛇足のコメント】
1. メールにはタイトル(Subject)を付けてください。
2. 連絡等の送信先は stat.nitech@gmail.com とします。

散布図にもっともらしい直線を当てはめたい

何に使う?
何に使える?

散布図: 縦軸=体重、横軸=身長

    : taijyuu*shintyou. A=1, B=2, ...
taijyuu |
    100 +                                                  A
        |                                                  A
        |                                       A                A
        |                                            A
     80 +                                  A     A      AAAA          A
        |                                     B  A    A B    A  A
        |                                       AAB EDB C AAA F  A A B
        |                                B B BC AGC BBB DDCB CEA A A
     60 +                      A   AA  BB  AACCABIFBHIFACCDC AAA AA
        |                  A    B   EC AA  ECHGACDF EBB DB     A
        |               A  A  CDAC DE  DDB CFBAA AA A   A
        |           A AAA BABA BAB AB ABA   A
     40 +        A AA    A B   B
        |                 A
        |
        |
     20 +
        |
        --+------------+------------+------------+------------+------------+--
         140          150          160          170          180          190
                                       shintyou

どの様に直線の位置決めをすれば良い?
一意に決めるにはどの様なアイディアがある?
そもそも、どういう考えから導き出された?
回帰分析とは?

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

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

プログラム : DSles0301.sas
/* Lesson 3-1 */ /* File Name = DSles0301.sas 10/22/18 */ options linesize=72 pagesize=20; options nocenter linesize=78 pagesize=30; proc printto log = 'Kougi/DSles0301_log.txt' print = 'Kougi/DSles0301_Results.txt' new; ods listing gpath='Kougi/SAS_ODS99'; data gakusei; infile 'Kougi/all07au.txt' 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; : 変量を指定 output out=outreg1 predicted=pred1 residual=resid1; : 結果項目の保存 run; : : proc print data=outreg1(obs=15); : 表示してみる run; : : proc plot data=outreg1; : 散布図を描く where shintyou^=. and taijyuu^=.; : 欠損値ではないサンプルに対して plot taijyuu*shintyou/vaxis=20 to 100 by 20; : 体重と身長(縦軸指定) plot pred1*taijyuu; : 予測値と観測値 plot resid1*pred1 /vref=0; : 残差と予測値(残差解析)(水平軸指定) plot resid1*shintyou/vref=0; : 残差と説明変数(残差解析) plot resid1*taijyuu /vref=0; : 残差と目的変数(残差解析) run; : : proc univariate data=outreg1 plot normal; : 残差を正規プロットして確かめる var resid1; : run; :
[備考] 上記のコロン以降は説明のためのものであり、 SAS のプログラムではありません。

出力結果 : DSles0301_Results.txt , DSles0301_out.pdf
2018年10月19日金曜日 14時27分07秒 50 Obs sex shintyou taijyuu kyoui jitaku kodukai carryer tsuuwa 1 F 145.0 38.0 . J 10000 . 2 F 146.7 41.0 85 J 10000 Vodafone 6000 3 F 148.0 42.0 . J 50000 . 4 F 148.0 43.0 80 J 50000 DoCoMo 4000 5 F 148.9 . . J 60000 . 6 F 149.0 45.0 . G 60000 . 7 F 150.0 46.0 86 40000 . 8 F 151.0 45.0 . J 20000 docomo 5000 9 F 151.0 50.0 . G 60000 J-PHONE . 10 F 151.7 41.5 80 J 35000 . 2018年10月19日金曜日 14時27分07秒 51 REG プロシジャモデル : MODEL1 従属変数 : taijyuu 読み込んだオブザベーション数 370 使用されたオブザベーション数 325 欠損値を含むオブザベーション数 45 分散分析要因自由度平方和平均平方 F 値 Pr > F Model 1 14055 14055 318.56 <.0001 Error 323 14251 44.12105 Corrected Total 324 28306 Root MSE 6.64237 R2 乗 0.4965 従属変数の平均 58.78092 調整済み R2 乗 0.4950 変動係数 11.30021 2018年10月19日金曜日 14時27分07秒 52 REG プロシジャモデル : MODEL1 従属変数 : taijyuu パラメータの推定パラメータ変数自由度推定値標準誤差 t 値 Pr > |t| Intercept 1 -79.35152 7.74804 -10.24 <.0001 shintyou 1 0.81883 0.04588 17.85 <.0001 2018年10月19日金曜日 14時27分09秒 53 Obs sex shintyou taijyuu kyoui jitaku kodukai carryer tsuuwa pred1 resid1 1 F 145.0 38.0 . J 10000 . 39.3789 -1.3789 2 F 146.7 41.0 85 J 10000 Vodafone 6000 40.7709 0.2291 3 F 148.0 42.0 . J 50000 . 41.8354 0.1646 4 F 148.0 43.0 80 J 50000 DoCoMo 4000 41.8354 1.1646 5 F 148.9 . . J 60000 . 42.5724 . 6 F 149.0 45.0 . G 60000 . 42.6542 2.3458 7 F 150.0 46.0 86 40000 . 43.4731 2.5269 8 F 151.0 45.0 . J 20000 docomo 5000 44.2919 0.7081 9 F 151.0 50.0 . G 60000 J-PHONE . 44.2919 5.7081 10 F 151.7 41.5 80 J 35000 . 44.8651 -3.3651 11 F 152.0 35.0 77 J 60000 DoCoMo 2000 45.1107 -10.1107 12 F 152.0 43.0 . J 20000 au 3500 45.1107 -2.1107 13 F 152.0 44.0 . 45000 DoCoMo 4000 45.1107 -1.1107 14 F 153.0 41.0 . J 125000 No . 45.9296 -4.9296 15 F 153.0 42.0 . G 0 Vodafone 1000 45.9296 -3.9296 2018年10月19日金曜日 14時27分09秒 54 : taijyuu*shintyou. A=1, B=2, ... taijyuu | 100 + A | A | A A | A 80 + A A AAAA A | B A A B A A | AAB EDB C AAA F A A B | B B BC AGC BBB DDCB CEA A A 60 + A AA BB AACCABIFBHIFACCDC AAA AA | A B EC AA ECHGACDF EBB DB A | A A CDAC DE DDB CFBAA AA A A | A AAA BABA BAB AB ABA A 40 + A AA A B B | A | | 20 + | --+------------+------------+------------+------------+------------+-- 140 150 160 170 180 190 shintyou 2018年10月19日金曜日 14時27分09秒 55 : pred1*taijyuu. A=1, B=2, ... 80 + 予 | 測 | 値 | A | A A B A t | A BABCB G A A a | BDAA EBB A A B A A i | A B BBADBCACBBA A AB j | A BBDBCIJBBB ECC A A y 60 + AAABHCDHFDBE BAA A A A u | A DCCGCDABA D AA u | BEACABDAAA ABA A | BB HC A BB | DAGDAAACAAA | BA B CB A A A | B ABABAA A | A AAAB A | AA A 40 + A A ---+-------------+-------------+-------------+-------------+-- 20 40 60 80 100 taijyuu 2018年10月19日金曜日 14時27分09秒 56 : resid1*pred1. A=1, B=2, ... 40 + | | A | A A | | A A 残 20 + A 差 | B A A | A A AA B | A A A AA A BDA A A A | A A B ABAA BBACAACAA A A | A A A A A CA A AACAGABABCAA AF 0 +-------------A-AB---AABABDACADBDA-C-BDEAEGFDCCBCCA---A--A------------ | AABA B BCC ADBB DEDBGEBFDCCB BE AA | BA B CA C AACAABCBAABBBA AAA A | A DA AC BC AA AAA | A A | -20 + --+------------+------------+------------+------------+------------+-- 30 40 50 60 70 80 予測値 taijyuu 2018年10月19日金曜日 14時27分09秒 57 : resid1*shintyou. A=1, B=2, ... 40 + | | A | A A | | A A 残 20 + A 差 | B A A | A A A AAA | A A A A A A B DA A A A | A A B ABA A BB ACA ACA A A A | AA A A A CA AAAC AFB BAB CAA A F 0 +--------A-AB---A-AAB-BDAC-ADB-DA--CBDCBAEFAFDBACBCC-A---A---A-------- | AAB A B B DB ADB BCFD BGDABFD CCB BE A A | BA B ABA CAABAAAB CBA ABBB AAA A A | A DA A C B CA A A AA | A A | -20 + --+------------+------------+------------+------------+------------+-- 140 150 160 170 180 190 shintyou 2018年10月19日金曜日 14時27分09秒 58 : resid1*taijyuu. A=1, B=2, ... 40 + | | A | A A | | A A 残 20 + A 差 | AA B | A A BB | A AA AA ACD A A A | A BABAAAH EB B | B AAAAC ABDEDCBDFA 0 +---------------A-BAACAGFDD-GDEIMDAGBAAA---------------------- | CAADAE ICEMCFLDCBC B | BA DCBADBBDDBDAAAA | A ADB DBB ABA | A A | -20 + ---+-------------+-------------+-------------+-------------+-- 20 40 60 80 100 taijyuu 2018年10月19日金曜日 14時27分09秒 59 UNIVARIATE プロシジャ変数 : resid1 (残差) モーメント N 325 重み変数の合計 325 平均 0 合計 0 標準偏差 6.63210957 分散 43.9848774 歪度 1.42132974 尖度 4.0008373 無修正平方和 14251.1003 修正済平方和 14251.1003 変動係数 . 平均の標準誤差 0.36788325 基本統計量位置ばらつき平均 0.00000 標準偏差 6.63211 中央値 -1.02372 分散 43.98488 最頻値 -2.30618 範囲 47.54351 四分位範囲 6.73181 Note: 3個の最頻値があります(度数: 5)。表では最頻値のなかで最も小さな値を表示します。 2018年10月19日金曜日 14時27分09秒 60 UNIVARIATE プロシジャ変数 : resid1 (残差) 位置の検定 H0: Mu0=0 検定 -統計量- ------p 値------- Student の t 検定 t 0 Pr > |t| 1.0000 符号検定 M -22.5 Pr >= |M| 0.0145 符号付順位検定 S -3248.5 Pr >= |S| 0.0552 正規性の検定検定 --統計量--- ------p 値------- Shapiro-Wilk W 0.915825 Pr < W <0.0001 Kolmogorov-Smirnov D 0.109794 Pr > D <0.0100 Cramer-von Mises W-Sq 0.914603 Pr > W-Sq <0.0050 Anderson-Darling A-Sq 5.307447 Pr > A-Sq <0.0050 2018年10月19日金曜日 14時27分09秒 61 UNIVARIATE プロシジャ変数 : resid1 (残差) 分位点 (定義 5) 水準分位点 100% 最大値 33.59967 99% 22.24447 95% 11.59967 90% 8.24447 75% Q3 2.69382 50% 中央値 -1.02372 25% Q1 -4.03799 10% -7.28364 5% -8.57436 1% -10.94384 0% 最小値 -13.94384 2018年10月19日金曜日 14時27分09秒 62 UNIVARIATE プロシジャ変数 : resid1 (残差) 極値 ----最小値---- ----最大値---- 値 Obs 値 Obs -13.9438 282 21.6938 267 -13.4474 346 22.2445 125 -11.4873 229 28.5997 326 -10.9438 284 29.2235 181 -10.9438 283 33.5997 327 欠損値 ---パーセント--- 欠損値カウント全体欠損値 . 45 12.16 100.00

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

幾つかのチェック項目
- 残差に一様性があるか? 残差プロット
- 残差が正規性をしているか? 診るには? ~~Q-Q~~ P-P plot
有効桁数に注意せよ : どこまでが「意味ある桁」か?
測定精度上回る計算結果は出せても、意味はない。プログラムの出力をそのまま書き写さないように。
[重要な注意] 統計ソフトは単なる道具。使いこなすのは各自。
　 [例1] 四捨五入の数値で考えてみれば : 精度(正確さ)が異なることに注意
　 [例2] 日本の観測史上の最高気温は、 2018(平成30)年7月23日に熊谷市で観測された41.1度であり、最低気温は、1902(明治35)年1月25日に北海道旭川市の-41度であった。===> -41.0度
　 [例3] 2001年のイチロー選手の打率は3割5分であった。 2006年は3割3分1厘であった。===> 3割5分0厘
次回は、... : 10月29日 13:00-14:30
- 今回の続き (回帰分析(後編))
- 「残差が大きい者を除外」して実行してみると?
- 重回帰分析
- ...
- 参考資料は、引き続き同じものを使うので持参ください

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