mwp-037 regress -機能概要と用例 regressは線形回帰モデルのフィットを行います。 1. 結果の解釈 1.1シミュレーションモデル 1.2良好なフィット 1.3小標本時のフィット 2. ベータ係数 3. 定数項の抑止 4. 重み付き回帰 5. ロバスト推定 補足１
1.
結果の解釈
regressは回帰分析系のコマンドの中で最も頻繁に使用されるコマンドです。このコマンドを実行すると例 えば次のような推定結果が得られます。 _cons 41.6797 2.165547 19.25 0.000 37.36172 45.99768 foreign -1.650029 1.075994 -1.53 0.130 -3.7955 .4954422 weight -.0065879 .0006371 -10.34 0.000 -.0078583 -.0053175 mpg Coef. Std. Err. t P>|t| [95% Conf. Interval] Total 2443.45946 73 33.4720474 Root MSE = 3.4071 Adj R-squared = 0.6532 Residual 824.171761 71 11.608053 R-squared = 0.6627 Model 1619.2877 2 809.643849 Prob > F = 0.0000 F( 2, 71) = 69.75 Source SS df MS Number of obs = 74 . regress mpg weight foreignc
ここに示した結果はStataのExampleデータセットauto.dtaに対し mpg = β0+ β1· weight + β2· foreign + ϵ という線形モデルをフィットさせたときの出力を示しています。Stataのマニュアルを見ればコマンドの使い 方は細かくわかるわけですが、残念ながら得られた結果をどう解釈・評価したら良いかに関する記述は十分で はありません。結果の評価は使用するデータやモデルに深く依存するため、一般的に言えることは表面的な側 面に限られてしまうという事情があるからです。 ここではやや変則的なアプローチですが、人為的に用意されたデータを用いてモデルフィットを行い、十分な サンプルがある場合とそうでない場合とでフィット結果がどう変化するかを見てみることにします。regress のアウトプットを評価する上でどこに注意しなくてはならないかは、このプロセスによって自ずと明白になる ものと期待されます。
1.1
シミュレーションモデル
データセットregress1.dtaには1,000件のデータが含まれています。 use http://www.math-koubou.jp/stata/data11/regress1.dta 先頭の5個のデータをリスト出力しておくと次のようになります。 . list in 1/5*1 5. 16.8 45 8 3 4. 11.6 26 2 3 3. 36.2 73 9 6 2. 35.8 54 7 10 1. 25.4 49 3 3 y x1 x2 x3 この変数yの値は . generate y = 0.5*x1 + 2*x2 - 10 + rnormal(0, 10) のような形で変数x1とx2の値をベースに生成したものです。rnormal(m, s)というのは正規分布に従う乱 数発生関数であり、mは平均値を、sは標準偏差を意味します（[D] functions参照）。すなわち y = 0.5x1+ 2.0x2− 10 (1) という線形の関係を基調に据えて乱数を発生させているわけです。なお、データセット中にはx1, x2とは全く 無関係の変数x3も含めてあります。yに対する線形フィットを実行するに際し、説明変数にx3も加えたとき にどのような現象が生ずるかを確認するためです。1.2
良好なフィット
まずデータが1,000件ある状態でregressコマンドを実行してみます。
• Statistics ◃ Linear models and related ◃ Linear regressionと操作 • Modelタブ: Dependent variable: y
Independent variables: x1 x2 x3 図1 regressダイアログ- Modelタブ _cons -8.906289 1.697041 -5.25 0.000 -12.23647 -5.576103 x3 .0072188 .1083894 0.07 0.947 -.205479 .2199165 x2 1.974713 .1090362 18.11 0.000 1.760746 2.18868 x1 .476147 .0303301 15.70 0.000 .4166288 .5356652 y Coef. Std. Err. t P>|t| [95% Conf. Interval] Total 160820.539 999 160.98152 Root MSE = 9.8954 Adj R-squared = 0.3917 Residual 97526.2879 996 97.9179598 R-squared = 0.3936 Model 63294.2507 3 21098.0836 Prob > F = 0.0000 F( 3, 996) = 215.47 Source SS df MS Number of obs = 1000 . regress y x1 x2 x3 (1)診断プロット regressの出力結果に関する考察を始める前に、フィットの状況を視覚化しておきましょう。これが単回帰で あれば散布図とフィット直線の組合せで視覚化は容易に行えるのですが、重回帰の場合にはそのような直接的 なアプローチは取れません。n(n ≥ 3)次元空間中における平面のフィット状況を何らかの形で2次元のグラ フに焼き直す必要があるわけですが、ここではrvfplotコマンドを用いてグラフ化を行ってみます。
. rvfplot, yline(0)*2
モデルをフィットさせると、それぞれの観測点(x1i, x2i, . . .)ごとにモデルによる推定値byiが算出されます。 その場合、実際の観測値yiとの差yi− byiが残差となるわけですが、rvfplot (residual-versus-fitted plot)
はこの残差の値をフィット値byiに対向させる形でプロットしたものです。標準偏差10の正規分布を前提に乱
数を発生させているために、それなりにバラツキを持った分布となっていることがわかります。
(2)決定係数R2
それではregressのアウトプットについて見て行くことにしましょう。regress出力の上半分にはANOVA
表が掲載されています。これは変数yの変動（バラツキ）に関する情報を集約したもので、SSは平方和(sum of squares)を意味します。ここでyi(i = 1, . . . , n)の平均値を¯y、最小2乗法によるyiの推定値をbyiと書く ことにしたとき、yの変動は X (yi− ¯y)2= X (yi− byi)2+ X (byi− ¯y)2 (2) のように分解することができます。今、 P
(yi− ¯y)2をTSS (total sum of squares)
P
(byi− ¯y)2をMSS (model sum of squares)
P
(yi− byi)2をRSS (residual sum of squares)
と書くことにすれば、上の例の場合、TSSは160820.539、MSSは63294.2507、RSSは97526.2879という ことになります。
モデルフィットの良し悪しを判断する上で良く用いられる決定係数(coefficient of determination)は R2= MSS TSS = 1 − RSS TSS (3) として定義されます。上の例では63294.2507/160820.539 = 0.3936となるわけですが、この値はANOVA 表の右側にR-squaredとしてレポートされています。すなわち全変動のうち39%がフィットされたモデルに よって表現されたことをこの数値は示しているわけです。 (3)モデル全体に対するp値
R-squaredという出力の1つ上にp値がProb > Fという形で示されています。これは分散分析 (ANOVA:
analysis of variance)の結果に基づくF検定の結果を示しているわけですが、その場合の帰無仮説H0は H0: β1= β2= · · · = 0
とするものです（β0は含まれません）。全く無意味な説明変数を設定しない限り、この帰無仮説は棄却される
のが普通です。
(4)係数推定値
ANOVA表の下には係数表が示されています。Coef. (coefficients)と書かれた列に係数の推定値が示されて
いるわけですが、ここでの例の場合、 b β1= 0.48 βb2= 1.97 βb3= 0.01 βb0= −8.91 という結果が得られています。背景にあるモデル式(1)と比べて十分良好な点推定値が得られていることがわ かります。 (5)係数に対するp値 それぞれの係数推定値に対してもp値が示されています。その場合の帰無仮説はβj = 0とするもので、そ れに対するt検定の結果がp値として示されているわけです。ここで特筆すべきは変数x3に対するp値が 0.947となっている点です。すなわちβ3= 0という帰無仮説が全く棄却できないことをこれは示しています。
実際、β3に関する95%信頼区間（CI: confidence interval） は[−0.21, 0.22]とレポートされており、0とい
う値がその中に含まれている点に注意してください。変数x3は有意ではないと判断せざるを得ません。
これに対し変数x1, x2、及び定数項の場合、p値はそれらの係数値が0とは言えないと主張しているだけなの
で、regressが出力してきたp値は余り参考にはならないでしょう。もし仮に理論モデル等によりx1の係数
値としては0.5が予想されるとしたとき、postestimation機能の1つであるtestコマンドを使用すればその
仮説を次のようにして検定することができます。
• Statistics ◃ Postestimation ◃ Tests ◃ Test linear hypothesesと操作 • Mainタブ: Test type for specification 1: Linear expressions are equal
図2 testダイアログ- Mainタブ Prob > F = 0.4318 F( 1, 996) = 0.62 ( 1) x1 = .5 . test (x1 = 0.5) (6)標準誤差
Coef.の隣には標準誤差(standard error)の値がStd. Err.として示されています。標準誤差とは標本平均の
標準偏差(standard deviation)を意味するもので、信頼区間を規定する上での重要な統計量であると言えま す。実際、変数x1の場合について与えられた標準誤差の値から95% CIを算出してみましょう。残差の自由 度(degrees of freedom)が996とレポートされていることに注意してt分布の逆関数 invttail(n, p)を使用すると（[D] functions参照）、 . display invttail(996, 0.025)*3 1.9623486 という値が求められます。これが95% の面積に対応した棄却限界値(critical value)であるわけですから、信 頼区間は点推定値と標準誤差の値から . display 0.476147 - 1.9623486*0.0303301 .41662877 . display 0.476147 + 1.9623486*0.0303301 .53566523 のように算出することができるわけです。係数表においてx1に対する95% CIは[0.4166288, 0.5356652]と レポートされているわけですが、その内容がこれによって検証できたことになります。
いずれにせよ、得られた係数推定値の信頼度はこれら標準誤差と信頼区間の値から判断することになります。 ここでの例の場合、x1, x2に対する信頼区間は十分狭いものとなっていますが、定数項についてはやや幅の 広いものとなっています。データ生成のプロセスを勘案するなら、これ以上は期待できないことは明らかで しょう。
1.3
小標本時のフィット
セクション1.2では1,000件のデータに対してフィットを行いましたが、今度は標本サイズを1/10に減らし て同じ操作を行ってみます。このresamplingにはランダムなプロセスが伴うため、ここでは再現性を期して 最初にseedの設定を行います。繰り返すごとに異なる結果が出ても構わないということであればseedの設定 は不要です。 . set seed 111*4 次にランダムに100件のデータを抽出します。 . sample 100, count*5 (900 observations deleted) これでメモリ上のデータは100件のみとなったわけです。この状態でregressを実行してみます。 . regress y x1 x2 x3 _cons -12.28143 6.422559 -1.91 0.059 -25.03011 .4672456 x3 -.1063396 .3448411 -0.31 0.758 -.7908438 .5781646 x2 1.565079 .3726188 4.20 0.000 .8254362 2.304721 x1 .583658 .1117937 5.22 0.000 .3617492 .8055668 y Coef. Std. Err. t P>|t| [95% Conf. Interval] Total 15850.8101 99 160.109193 Root MSE = 10.554 Adj R-squared = 0.3043 Residual 10693.3387 96 111.388945 R-squared = 0.3254 Model 5157.47143 3 1719.15714 Prob > F = 0.0000 F( 3, 96) = 15.43 Source SS df MS Number of obs = 100 . regress y x1 x2 x3 R2値にはそれほど大きな変化はありません。また今回も変数x 3は有意ではないと判定されています。一方、 係数の推定値は微妙に変化していますが、その要因は標準誤差の値がセクション1.2の場合に比べて3 ∼ 4倍 に増えている点にあります。結果的に信頼区間の幅がそれだけ広くなってしまっています。例えばx1の係数 値としては0.37とか0.8といった値も有意水準5%では否定できないことになります。実際、testコマンド を用いて検定を行ってみましょう。 *4[R] set seed 参照. test (x1 = 0.37) Prob > F = 0.0590 F( 1, 96) = 3.65 ( 1) x1 = .37 . test (x1 = 0.37) . test (x1 = 0.8) Prob > F = 0.0559 F( 1, 96) = 3.74 ( 1) x1 = .8 . test (x1 = 0.8) 共にp値が0.05よりも大きいため、仮説を棄却することはできないわけです。 有意水準(significance level) としては5%という値がデフォルトとして仮定されますが、それは
Reportingタブ上のConfidence levelの設定によって調整できます（図3参照）。
2.
ベータ係数
今度はStataのExampleデータセットauto.dtaを用いて線形回帰モデルのフィットを行ってみます。
. sysuse auto.dta*6
(1978 Automobile Data)
このデータセット中には1978年に米国で販売された74車種のデータが収納されていますが、ここでは走行
燃費性能を表すmpg (miles per gallon)を車重weightとその2乗weight2とによって説明する線形回帰モ
デルを考えます。2乗項については交互作用演算子を用いてc.weight#c.weightと記述する方法もあります
が（mwp-028参照）、ここでは従来通り、新たな変数を生成するアプローチで対応することにします。
. generate weightsq = weight^2*7 まずデータの内容を確認しておきます。
*6メニュー操作：File ◃ Example datasets ◃ Example datasets installed with Stata と操作しダウンロードする。
. format weightsq %10.0g*8
. list mpg weight weightsq in 1/5*9
5. 15 4,080 16646400 4. 20 3,250 10562500 3. 22 2,640 6969600 2. 17 3,350 11222500 1. 22 2,930 8584900 mpg weight weightsq このデータの場合、説明変数の桁数に大きな違いがある点に注意する必要があります。あらかじめweightの 単位をKポンドにでもしておけば、その2乗項も2桁程度の数値でおさまり扱いやすくなるわけですが、こ こでは敢えて現状のままフィットを行ってみます。
. regress mpg weight weightsq*10
_cons 51.18308 5.767884 8.87 0.000 39.68225 62.68392 weightsq 1.32e-06 6.26e-07 2.12 0.038 7.67e-08 2.57e-06 weight -.0141581 .0038835 -3.65 0.001 -.0219016 -.0064145 mpg Coef. Std. Err. t P>|t| [95% Conf. Interval] Total 2443.45946 73 33.4720474 Root MSE = 3.3587 Adj R-squared = 0.6630 Residual 800.937487 71 11.2808097 R-squared = 0.6722 Model 1642.52197 2 821.260986 Prob > F = 0.0000 F( 2, 71) = 72.80 Source SS df MS Number of obs = 74 . regress mpg weight weightsq
この場合、
d
mpg = −1.42e-2 · weight + 1.32e-6 · weight2+ 51.18
というモデル式が得られたわけです。係数値だけを見れば、mpgdに対する影響の大きさはweightの方が weight2より104倍も大きいように見えるわけですが、これは正しい解釈とは言えません。weight2の変数 値自体が106とか107といった大きな値を取るからです。 説明変数が及ぼす効果を相対的に比較・評価したいといった場合に良く利用されるのがベータ係数です。これ はすべての変数値を規格化(normalize)、すなわち平均が0、標準偏差が1となるように標準化した上で求め た係数値であり、被説明変数に対する効果を相対的に比較する上で有用な尺度であると言えます。このベータ 係数はbetaオプションを指定することで算出することができます。
*8メニュー操作：Variables ウィンドウ内で weightsq を右クリック、Format ’weightsq’ を選択
*9メニュー操作：Data ◃ Describe data ◃ List data
• regressダイアログ: Reportingタブ: Standardized beta coefficients: X
図3 regressダイアログ- Reportingタブ
_cons 51.18308 5.767884 8.87 0.000 . weightsq 1.32e-06 6.26e-07 2.12 0.038 1.104148 weight -.0141581 .0038835 -3.65 0.001 -1.901918 mpg Coef. Std. Err. t P>|t| Beta Total 2443.45946 73 33.4720474 Root MSE = 3.3587 Adj R-squared = 0.6630 Residual 800.937487 71 11.2808097 R-squared = 0.6722 Model 1642.52197 2 821.260986 Prob > F = 0.0000 F( 2, 71) = 72.80 Source SS df MS Number of obs = 74 . regress mpg weight weightsq, beta
betaオプションが指定された場合には信頼区間に関する情報に代ってベータ係数の推定値が出力されます。
ここでの例で言えば、weightの効果はweight2の効果の約1.7倍であることがわかります。
3.
定数項の抑止
引続き同じauto.dtaを使用しますが、今度は車重weightを車長lengthの線形関数として表現したいとし
ます。この場合、特別な指定を行わなければ定数項β0（切片）を含んだ形のモデル式がフィットされます。す
なわちlength = 0のときにβ0という重量の残るモデルとなるわけです。これに対しβ0 = 0、すなわち原点
を通る直線をフィットさせたいという場合にはnoconstantオプションを指定します。
• Statistics ◃ Linear models and related ◃ Linear regressionと操作 • Modelタブ: Dependent variable: weight
Independent variables: length Suppress constant term: X
length 16.29829 .2774752 58.74 0.000 15.74528 16.8513 weight Coef. Std. Err. t P>|t| [95% Conf. Interval] Total 718762200 74 9713002.7 Root MSE = 451.68 Adj R-squared = 0.9790 Residual 14892897.8 73 204012.299 R-squared = 0.9793 Model 703869302 1 703869302 Prob > F = 0.0000 F( 1, 73) = 3450.13 Source SS df MS Number of obs = 74 . regress weight length, noconstant
車長が1インチ増加するごとに車重は16.3ポンド増加するというモデル式が誘導されたことになります。
4.
重み付き回帰
本セクションではExampleデータセットcensus9.dtaを使用することにします。
. use http://www.stata-press.com/data/r11/census9, clear*11
(1980 Census data by state)
このデータセット中には米国における1980年国勢調査の結果が州別に記録されています。次に示すのは先頭
5州のデータです。
*11メニュー操作：File ◃ Example Datasets ◃ Stata 11 manual datasets と操作、Base Reference Manual [R] の regress の項
. list state drate pop medage region in 1/5, nolabel*12 5. California 79 23,667,902 29.90 4 4. Arkansas 99 2,286,435 30.60 3 3. Arizona 78 2,718,215 29.20 4 2. Alaska 40 401,851 26.10 4 1. Alabama 91 3,893,888 29.30 3
state drate pop medage region
ここでは死亡率 (drate) のデータを州別の年齢中央値 (medage)と地域区分 (region)とによって説明し
ようとする線形回帰モデルをフィットさせます。ただしregionはカテゴリ変数であり、Northeast, North
Central, South, Westという米国内地域区分がそれぞれ1, 2, 3, 4というコードによって表現されています。
このフィットに際して注意すべきは、drateもmedageも共に平均化された値であるという点です。そのため
州の人口を表す変数popを重みとして用いた形のフィットがより適切であると判断されます。その際、重み
の選択には注意を払う必要があります（mwp-027参照）。
図4 regressダイアログ - Weightsタブ
今の場合、Analytic weightsが正しい選択です。間違えてFrequency weightsを選択した場合、例えばAlabama
州であれば州人口の3, 893, 888人すべてが29.3才という扱いとなってしまうので、極めて奇妙な結果が導か
れてしまう結果となります。
• Statistics ◃ Linear models and related ◃ Linear regressionと操作 • Modelタブ: Dependent variable: drate
Independent variables: medage i.region*13
• Weightsタブ: Analytic weights: pop
_cons -39.14727 17.23613 -2.27 0.028 -73.86262 -4.431915 4 -10.90629 2.681349 -4.07 0.000 -16.30681 -5.505777 3 -1.438452 2.320244 -0.62 0.538 -6.111663 3.234758 2 .3138738 2.456431 0.13 0.899 -4.633632 5.26138 region medage 4.283183 .5393329 7.94 0.000 3.196911 5.369455 drate Coef. Std. Err. t P>|t| [95% Conf. Interval] Total 5335.01916 49 108.877942 Root MSE = 5.246 Adj R-squared = 0.7472 Residual 1238.40987 45 27.5202192 R-squared = 0.7679 Model 4096.6093 4 1024.15232 Prob > F = 0.0000 F( 4, 45) = 37.21 Source SS df MS Number of obs = 50 (sum of wgt is 2.2591e+08)
. regress drate medage i.region [aweight = pop]
コマンド構文上は[aweight = pop]という指定が追加されている点に注意してください。
5.
ロバスト推定
特に何も指定しなかった場合、regressは最小2乗法(OLS: ordinary least squares)による回帰を実行しま
す。しかしOLSの場合には分散の均一性(homoskedasticity) が一つの前提となっている点に注意する必要
があります。例えばExampleデータセットauto.dtaに対して次のようなフィットを行ってみます。
. sysuse auto, clear (1978 Automobile Data)
. replace weight = weight/1000*14
*13i. 演算子については mwp-028 を参照。
*14メニュー操作：Data ◃ Create or change data ◃ Change contents of variable 適正な係数値を得るため、車重の単位をポンドからK ポンドに変更します。
. regress mpg weight
_cons 39.44028 1.614003 24.44 0.000 36.22283 42.65774 weight -6.008687 .5178782 -11.60 0.000 -7.041058 -4.976316 mpg Coef. Std. Err. t P>|t| [95% Conf. Interval] Total 2443.45946 73 33.4720474 Root MSE = 3.4389 Adj R-squared = 0.6467 Residual 851.469221 72 11.8259614 R-squared = 0.6515 Model 1591.99024 1 1591.99024 Prob > F = 0.0000 F( 1, 72) = 134.62 Source SS df MS Number of obs = 74 . regress mpg weight
このモデルは単回帰であるため、フィットの様子を直接視覚化することが可能です。
. twoway (scatter mpg weight) (lfit mpg weight), ytitle(mpg)*15
フィットされた直線周囲のデータ点のバラツキは余り均等とは言えないように見えます。この点は診断プロッ トの1つであるrvpplot (residual-versus-predictor plot) を使うとより鮮明なものにできます。
. rvpplot weight, yline(0)*16
*15メニュー操作：Graphics ◃ Twoway graph (scatter, line, etc.) 詳細については補足１を参照。
rvpplotは残差を特定の説明変数（今の場合、weight）と対向させた形でプロットする機能を提供します。
この分散不均一性(heteroskedasticity)はフォーマルな検定によっても確かめることができます。
• Statistics ◃ Postestimation ◃ Reports and statisticsと操作
• estatダイアログ: Reports and statistics: Tests for heteroskedasticity (hettest)
Prob > chi2 = 0.0009 chi2(1) = 11.05
Variables: fitted values of mpg Ho: Constant variance
Breusch-Pagan / Cook-Weisberg test for heteroskedasticity . estat hettest
estat hettestの出力に示されているように、均一な分散を仮定した帰無仮説H0に対するp値が0.0009で
すから、H0は棄却せざるを得ないわけです。
このように今回設定したモデルについては分散の均一性が主張できないため、OLSの結果を疑ってかかる必
要がありそうです。OLSに代る推定法を選択するにはregressダイアログのSE/Robustタブを利用します
（SEは標準誤差(standard error)の略です）。
• Statistics ◃ Linear models and related ◃ Linear regressionと操作 • Modelタブ: Dependent variable: mpg
Independent variables: weight • SE/Robustタブ: Robust
_cons 39.44028 1.98832 19.84 0.000 35.47664 43.40393 weight -6.008687 .5840839 -10.29 0.000 -7.173037 -4.844337 mpg Coef. Std. Err. t P>|t| [95% Conf. Interval]
Robust
Root MSE = 3.4389 R-squared = 0.6515 Prob > F = 0.0000 F( 1, 72) = 105.83 Linear regression Number of obs = 74 . regress mpg weight, vce(robust)
vce(robust)オプションを指定した場合には分散不均一性を許容する推定法が用いられます。OLSの場合と 比べ点推定値に変化は生じませんが、標準誤差や信頼区間の推定値が変わってきます。たまたま今回の例では 次の表に見られるように95% CIに大きな違いは見られませんでした。 OLS Robust weight [−7.04, −4.98] [−7.17, −4.84] cons [36.22, 42.66] [35.48, 43.40]
SE/Robustタブ上でClustered robustを選択した場合には、グループ（クラスタ）内での相関を許 容した推定法が実行されます。
補足１
–
グラフ作成コマンド操作
セクション5では散布図にフィット直線を重ね合わせる形でグラフを作成しましたが、その操作の詳細を記し
ておきます。
• Graphics ◃ Twoway graph (scatter, line, etc.) と操作
• Plotsタブ上でCreateボタンをクリック、表示されるPlot 1ダイアログ上で次の設定を行う。 ◦ Choose a plot category and type: Basic plots（デフォルト）
◦ Basic plots: Scatter ◦ Y variable: mpg ◦ X variable: weight
• Plotsタブ上で再度Createボタンをクリック、表示されるPlot 2ダイアログ上で次の設定を行う。 ◦ Choose a plot category and type: Fit plots
◦ Fit plots: Linear prediction（デフォルト） ◦ Y variable: mpg
◦ X variable: weight • Y axisタブ: Title: mpg
. twoway (scatter mpg weight) (lfit mpg weight), ytitle(mpg)
