Stata 11 whitepaper mwp 4 mwp mwp-028 / 41 mwp mwp mwp-079 functions 72 mwp-076 insheet 89 mwp-030 recode 94 mwp-033 reshape wide

【基盤機能編】

目  次

機能区分 コマンド whitepaperタイトル ページ mwp番号 一般 − 推定機能 4 mwp-027 − 因子変数 23 mwp-028 − 日付/時間情報の入力 41 mwp-001 − プログラミング機能 51 mwp-078 − サーベイデータの扱い 62 mwp-079 データ管理 functions 関数群 72 mwp-076 insheet スプレッドシートの読込み 89 mwp-030 recode コードの変換 94 mwp-033 reshape wide/long形式変換 100 mwp-036 推定機能 ivregress 操作変数法による回帰 110 mwp-082 logistic/logit ロジスティック回帰 127 mwp-039 logistic/logit 推定後機能 137 mwp-040 margins マージン分析機能 151 mwp-029 mlogit 多項ロジスティック回帰 167 mwp-090 ologit 順序ロジスティック回帰 184 mwp-088 poisson ポアソン回帰 199 mwp-087 regress 線形回帰 217 mwp-037 regress 推定後機能 235 mwp-038 検定機能 anova/oneway 分散分析 247 mwp-042 sdtest χ2_{検定、F 検定 272 mwp-043} ttest t検定 279 mwp-041 テーブル作成 table 要約統計量のテーブル化 286 mwp-070 tabstat 要約統計量のテーブル化 294 mwp-071 tabulate 一元度数分布表 299 mwp-072 tabulate 二元度数分布表 305 mwp-073

c

⃝ 2011 Math工房 一部 ⃝ 2011 StataCorp LPc

Math工房 web: www.math-koubou.jp

mwp-076 functions - 関数群 Stataには数多くの関数が用意されており、generateコマンドや数式記述の中で利用されます。詳細は [D] functionsを参照いただくとして、ここでは代表的な関数に限定した形でその用法を紹介します。なお、 [D] egen (mwp-077 )についても併せてご参照ください。 1. 乱数の発生 1.1一様乱数 1.2正規乱数 2. 数学関数 2.1整数値への変換 2.2 Running sum 3. 統計分布関数 3.1 t検定の実行 3.2 t分布のグラフ化 3.3面積の算出 3.4限界値の算出 4. プログラミング関数 4.1コード化関数 5. 文字列関数 6. 日付/時間関数 7. 行列操作関数 補足１ c

1.

乱数の発生

統計解析を行う上で擬似乱数を発生させたいといったニーズが生ずる場合があります。Stataには種々の分布

に従って乱数を発生させる関数が用意されていますが、本セクションでは一様分布と正規分布のケースを例に とってその用例を紹介します。詳細については[D] functions のセクション“Random-number functions” をご参照ください。

1.1

一様乱数

関数runiform()は[0, 1)上で一様に分布する擬似乱数を発生させます。例えば観測データ (observations) 数を10に設定した上で、変数x₁中に一様乱数を発生させてみます。なお、repeatabilityを確保する意味で 最初に擬似乱数発生用のseedを設定しておきます。 . set obs 10 obs was 0, now 10 . set seed 2 . generate x1 = runiform()*1 . list*2 10. .409848 9. .3874338 8. .9104601 7. .1747266 6. .3360431 5. .518447 4. .6991696 3. .6303533 2. .0515642 1. .850512 x1 区間[a, a + b)上での一様乱数を発生させたい場合には a + b*runiform() という数式を用いれば良いわけです。例えば区間[10, 20)上での一様乱数を変数x₂中に発生させるには次の ように操作します。

*1メニュー操作：Data ◃ Create or change data ◃ Create new variable

. generate x2 = 10 + 10*runiform() . list x2 10. 12.13104 9. 10.89719 8. 18.96779 7. 14.25147 6. 17.27702 5. 14.2031 4. 14.25097 3. 16.53234 2. 11.52971 1. 19.59647 x2 整数値が欲しい場合には数学関数round()（四捨五入）、またはfloor()（切捨て）を使用します。 . generate x3 = round(x2) . list x2 x3 10. 12.13104 12 9. 10.89719 11 8. 18.96779 19 7. 14.25147 14 6. 17.27702 17 5. 14.2031 14 4. 14.25097 14 3. 16.53234 17 2. 11.52971 12 1. 19.59647 20 x2 x3

このround()やfloor()をruniform()と組み合わせて使用することもできます。 . generate x4 = round(10 + 10*runiform())

. list x4 10. 13 9. 18 8. 16 7. 15 6. 18 5. 15 4. 20 3. 20 2. 14 1. 16 x4

1.2

正規乱数

評価版では割愛しています。

2.

数学関数

Stataには数多くの数学関数が用意されています。全体像については[D] functions のセクション

“Mathe-matical functions”を参照いただくとして、ここでは注意を要する事項についてのみ記載しておきます。

2.1

整数値への変換

実数値から整数値への変換という操作はしばしば必要となります。具体的には ◦ round() –四捨五入 ◦ floor() –切捨て ◦ ceil() –切上げ といった関数が用いられるわけですが、これ以外に ◦ int() – 0方向への切捨て という関数もあるので説明を補足しておきます。 (1) xが正の実数値の場合

round(x) floor(x) ceil(x) int(x)

x = 4.8 5 4 5 4

x = 5.2 5 5 6 5

この場合、int()の挙動はfloor()と等しくなります。

(2) xが負の実数値の場合

round(x) floor(x) ceil(x) int(x)

x = −4.8 −5 −5 −4 −4

x = −5.2 −5 −6 −5 −5

2.2 Running sum

数学関数の中にsum()という関数がありますが、これはrunning sumを算出する機能を提供します。例えば

1, 2, . . . , 10という値を持つ変数xを生成してみましょう。 . clear . set obs 10 . generate x = n*3 この変数xに対してrunning sumを計算すると次のようになります。 . generate y1 = sum(x) . list 10. 10 55 9. 9 45 8. 8 36 7. 7 28 6. 6 21 5. 5 15 4. 4 10 3. 3 6 2. 2 3 1. 1 1 x y1 変数y₁にはPi_j=1x_j (i = 1, 2, . . . , 10)という形の和がセットされていることがわかります。これに対して egenのtotal()関数を使用すると総和の値P10_j=1x_jが算出されます。 . egen y2 = total(x)*4 . list 10. 10 55 55 9. 9 45 55 8. 8 36 55 7. 7 28 55 6. 6 21 55 5. 5 15 55 4. 4 10 55 3. 3 6 55 2. 2 3 55 1. 1 1 55 x y1 y2

*3 _{n というシステム変数については [U] 13.4 System variables をご参照ください。} *4メニュー操作：Data ◃ Create or change data ◃ Create new variable (extended)

3.

統計分布関数

Stataには各種統計分布に関連した関数が一式用意されています。具体的には正規分布とか二項分布といった 統計分布ごとに ◦ 確率密度関数 ◦ 累積分布関数 ◦ 逆関数 等の値を算出する関数が提供されています。詳細については[D] functionsのセクション“Probability

distributions and density functions”をご参照ください。ここではStudentのt分布を例にとってその用例を

紹介します。

3.1 t

検定の実行

最初にExampleデータセットfuel.dtaを用いてt検定を実行してみます。 . use http://www.stata-press.com/data/r11/fuel, clear*5

このデータセット中にはある燃料添加剤を加えたときと加えなかったときでの1ガロン当りの走行距離が、12 台の車に関しwide形式で記録されています。変数mpg1は燃料添加剤を加えなかったときのデータに、mpg2 は燃料添加剤を加えたときのデータに対応しています。 . list, separator(0) 12. 19 23 11. 23 21 10. 24 27 9. 20 24 8. 24 28 7. 18 17 6. 17 18 5. 18 23 4. 25 22 3. 21 21 2. 23 25 1. 20 24 mpg1 mpg2 ここでmpg1 = mpg2と言えるかどうかについてt検定を実行してみましょう。ただし有意水準はデフォルト の5%を用います。

*5_{メニュー操作：File ◃ Example Datasets ◃ Stata 11 manual datasets と操作、Base Reference Manual [R] の ttest の項よ}

. ttest mpg1 == mpg2*6

Pr(T < t) = 0.0232 Pr(|T| > |t|) = 0.0463 Pr(T > t) = 0.9768 Ha: mean(diff) < 0 Ha: mean(diff) != 0 Ha: mean(diff) > 0 Ho: mean(diff) = 0 degrees of freedom = 11 mean(diff) = mean(mpg1 - mpg2) t = -2.2444 diff 12 -1.75 .7797144 2.70101 -3.46614 -.0338602 mpg2 12 22.75 .9384465 3.250874 20.68449 24.81551 mpg1 12 21 .7881701 2.730301 19.26525 22.73475 Variable Obs Mean Std. Err. Std. Dev. [95% Conf. Interval] Paired t test . ttest mpg1 == mpg2 この場合、mpg1とmpg2の間には1.75の差があり、検定統計量であるt値が−2.2444、自由度11、両側検定 時のp値が0.0463とレポートされています。なお、ttestコマンドについては[R] ttest (mwp-041 ) をご 参照ください。

3.2 t

分布のグラフ化

関数tden(n, t)を用いるとt分布の確率密度関数の値を算出できますが、ここではこの関数を用いてグラフを 作成してみます。なお、t分布の関数式は自由度nに依存するわけですが、今の場合、nの値としては11を 指定します。

. twoway (function tden(11, x), range(-4 4)), ytitle("") xtitle(t)

> title(t(11) distribution)*7

*6メニュー操作：Statistics ◃ Summaries, tables and tests ◃ Classical tests of hypotheses ◃ Mean-comparison test, paired data

ttest出力の主張するところは、|t| > 2.2444の部分における曲線下の面積が双方足して0.0463であるとい うことです。

3.3

面積の算出

評価版では割愛しています。

3.4

限界値の算出

評価版では割愛しています。

4.

プログラミング関数

Stataにはプログラミングを支援するための関数が多数用意されています。詳細は [D] functions のセク ション“Programming functions”を参照いただくとして、ここでは有用性の高いコード化に関連した関数につ いてその用例を紹介します。

4.1

コード化関数

ある連続変数（もしくはそれに準ずる離散変数）が与えられたときに、その値に応じてデータを区分化/ コード化したいというニーズはしばしば起ります。ここでは年齢に関するデータを例にとって関連する関数 recode(), autocode(), irecode()の用法を説明します。使用するデータセットはage01.dtaです。

. use http://www.math-koubou.jp/stata/data11/age01.dta データセットの内容は次のとおりです。

. list, separator(0) 8. 82 7. 73 6. 65 5. 58 4. 44 3. 36 2. 21 1. 18 age (1) recode()の用例 関数recode()は(x, x₁, x₂, . . . , x_n)という引数を取ります。x₁, x₂, . . . , x_nが変数xに関する境界を意味す ることになります。例えば

. generate code1a = recode(age, 20, 30, 40, 50, 60, 70, 80, 90) と入力すると次のような結果が得られます。

. list age code1a, separator(0)

8. 82 90 7. 73 80 6. 65 70 5. 58 60 4. 44 50 3. 36 40 2. 21 30 1. 18 20 age code1a 境界値x₁, x₂, . . . , x_nは等間隔である必要はなく、またx > x_nとなるようなデータがあっても構いません。

. generate code1b = recode(age, 25, 40, 60) . list age code1b, separator(0)

8. 82 60 7. 73 60 6. 65 60 5. 58 60 4. 44 60 3. 36 40 2. 21 25 1. 18 25 age code1b x > xnとなるようなデータに対しては値x_nがセットされている点に注意してください。

(2) autocode()の用例

評価版では割愛しています。

(3) irecode()の用例

評価版では割愛しています。

5.

文字列関数

Stataには文字列を操作する関数が多数用意されています。全容は [D] functions のセクション“String

functions”を参照いただくとして、ここでは和暦を西暦に変換する操作を一例として紹介します。

評価版では割愛しています。

6.

日付

/

時間関数

Data Editor上で10/15/2011や01jan2011ような日付データを入力した場合、これらは文字列データとし て認識されてしまうため、このままでは時系列分析等で利用することができません。このため最初に数値デー

タ、具体的には01jan1960を起点とした整数値に変換してやる必要があります。そのために用いられるのが

daily()とかmonthly()といった関数です。詳細については [D] functionsのセクション“Date and time

functions”を、用例についてはmwp-001 をご参照ください。

7.

行列操作関数

Stataには行列演算に関連した関数も多数用意されています。詳細は [D] functionsのセクション“Matrix

functions returning a matrix”、及び“Matrix functions returning a scalar”をご参照いただくとして、ここ では一部の用例のみを紹介します。最初に次のような行列 A =  1 1 1_{1 2 2} 1 2 3   を定義します。 . matrix input A = (1,1,1\1,2,2\1,2,3) . matrix list A r3 1 2 3 r2 1 2 r1 1 c1 c2 c3 symmetric A[3,3] この場合、Aは対称行列のため、下三角要素しか表示されていない点に注意してください。Aの逆行列を求め るには関数inv()を使用します。

. matrix B = inv(A) . matrix list B c3 0 -1 1 c2 -1 2 c1 2 r1 r2 r3 symmetric B[3,3] すなわち A−1 ₌  ₋₁2 −1_{2 −1}0 0 −1 1   と求まったわけです。 次にAをCholesky分解してみます。使用するのは関数cholesky()です。 . matrix C = cholesky(A) . matrix list C r3 1 1 1 r2 1 1 0 r1 1 0 0 c1 c2 c3 C[3,3] 実際、CCT を計算してみると . matrix D = C*C’ . matrix list D r3 1 2 3 r2 1 2 r1 1 r1 r2 r3 symmetric D[3,3] より A = CCT となっていることが確認できます。

補足１

–

グラフ作成コマンド操作

評価版では割愛しています。 ¥

mwp-082 ivregress - 操作変数法による回帰 説明変数が誤差項と相関を持つ場合、OLS推定量には不偏性も一致性も期待できなくなります。そのような 場合には操作変数を用いた推定が必要となります。 1. OLSの前提条件 2. OLS推定量の性質 3. 確率的説明変数 4. 操作変数法 5. 2段階最小2乗法 6. ivregressの用例 6.1データセット 6.2 2SLS推定法 6.3 LIML推定法 6.4 GMM推定法 7. ivregress postestimationの用例 7.1 estat endogenous 7.2 estat firststage 7.3 estat overid

1. OLS

の前提条件

ここでは煩雑さを避ける意味で、単回帰モデル yi= β0+ β1xi+ ui, i = 1, 2, . . . , n (M1)

を対象に議論を進めて行くことにします。通常、β₀ とβ₁ の推定には最小2 乗法 (OLS: ordinary least

squares)が使用されますが、その推定が適切に行われるためには次の4条件が前提となる点に注意する必要

があります。

c

前提１： 説明変数xは確率変数ではなく固定変数(fixed variable)である。 前提２： 誤差項uは確率変数で、期待値は0である。 E(ui) = 0, i = 1, 2, . . . , n (M2) 前提３： 誤差項の分散は均一である。 V (ui) = E(u2i) = σ2, i = 1, 2, . . . , n (M3) 前提４： 異なった誤差項は無相関である。 Cov(ui, uj) = E(uiuj) = 0, i ̸= j, i, j = 1, 2, . . . , n (M4)

これらの前提のもとでOLS推定法は残差2乗和(RSS: residual sum of squares)、すなわち

RSS =Xn i=1 (yi− ˆyi)2= n X i=1 (yi− ˆβ0− ˆβ1xi)2 (M5) を最小化する推定量 ˆβ₀, ˆβ₁を求めます。具体的には ∂RSS ∂ ˆβ0 = −2 n X i=1 (yi− ˆβ0− ˆβ1xi) = 0 (M6a) ∂RSS ∂ ˆβ1 = −2 n X i=1 (yi− ˆβ0− ˆβ1xi)xi= 0 (M6b) を連立させて解けばよく、ˆβ₀, ˆβ₁は次式で与えられることになります。 ˆβ1= P_n

i=1_P(xi− ¯x)(yi− ¯y) n i=1(xi− ¯x)2 µ = P_n i=1(xi− ¯x)yi P_n i=1(xi− ¯x)2 ¶ (M7a) ˆβ0= ¯y − ˆβ1¯x (M7b)

2. OLS

推定量の性質

(1)不偏性(unbiasedness) ある母集団パラメータをθ、その推定量をˆθとします。何度も標本抽出と推定を繰り返したとき、ˆθの期待値 であるE(ˆθ)が真の値θに一致するとき、すなわち E(ˆθ) = θ (M8) が成り立つとき、ˆθは不偏推定量 (unbiased estimator)であると言います。 前提条件1から4のもとで得られるOLS推定量 ˆβ₀, ˆβ₁は不偏推定量です。またそれらは線形不偏推定量の中

で最も分散の小さな推定量、すなわち最良線形不偏推定量(BLUE: best linear unbiased estimator)でもあ ります[ Gauss-Markovの定理 ]。

(2)一致性(consistency) 標本数nを∞にしたとき、推定量の値ˆθがθに収束する、すなわち任意のϵ > 0について lim n→∞P (|ˆθ − θ| ≥ ϵ) = 0 が成立するとき、ˆθはθの一致推定量(consistent estimator)であると言います。またˆθがθの一致推定量で あることをˆθがθに確率収束すると言い、 plim ˆθ = θ (M9) のように書きます。 OLS推定量 ˆβ₀, ˆβ₁は一致推定量でもあります。

3.

確率的説明変数

説明変数が非確率変数であることを主張する前提1が崩れた場合、推定量の重要な資質である不偏性と一致性 は共に失われることになります。 (1)不偏性

今、説明変数xが確率変数であることに注意して(M7a)式を変形します。y_i− ¯y = β₁(x_i− ¯x) + (u_i− ¯u)と 書けることから

ˆβ1=

P_n

i=1_P(xi− ¯x)(yi− ¯y) n

i=1(xi− ¯x)2

= β1+

P_n

i=1_P(xi− ¯x)(ui− ¯u) n

i=1(xi− ¯x)2 (M10)

従ってその期待値は

E[ ˆβ1] = β1+ E

·Pn

i=1_P(xi− ¯x)(ui− ¯u) n i=1(xi− ¯x)2 ¸ (M11) で与えられることになります。この第2項は変数xが誤差項uと統計的に独立であれば0となりますが、一 般的にはxとuの間の相関によって0とはならず、ˆβ₁はβ₁の不偏推定量とは言えなくなります。 (2)一致性 今、 plim(V (xi)) = plim " _n X i=1 (xi− ¯x)2/n # = σ2 x plim(Cov(xi, ui)) = plim " _n X i=1 (xi− ¯x)(ui− ¯u)/n # = σxu

であるとします。このとき(M10)より plim( ˆβ1) = plim

·

β1+

P_n

i=1_P(xi− ¯x)(ui− ¯u) n

i=1(xi− ¯x)2

¸

= β1+ plim

·P_n

i=1_P(xi− ¯x)(ui− ¯u)/n n

i=1(xi− ¯x)2/n

¸

= β1+plim [

P_n

i=1(xi− ¯x)(ui− ¯u)/n]

plim [Pn_i=1(xi− ¯x)2/n] = β1+σ_σxu₂ x (M12) となるため、σ_xu̸= 0の場合、すなわちxとuの間に相関がある場合には、ˆβ₁はβ₁の一致推定量とは言えな くなります。 (3)観測誤差を伴うモデル 今、2つの物理量ξとηの間の関係が ηi= β0+ β1ξi+ ϵi, i = 1, 2, . . . , n という単回帰モデルで表現できるものとします。これらの物理量を観測するに際し、ηについては誤差のない 値が観測できるが、ξについては観測誤差を伴うものとします。この場合、観測値x_iとy_iの間の関係式を yi= β0+ β1xi+ ui, i = 1, 2, . . . , n と書いた場合、x_iは実際にはξ_i+ v_i（v_iは観測誤差）と表現されるため、u_i= ϵ_i− β₁v_iであることになりま す。従って Cov(xi, ui) = Cov(ξi+ vi, ϵi− β1vi) = β1σv2̸= 0 すなわち、説明変数xと誤差項uの間には相関が生じ、OLS推定量は不偏性も一致性も持たないことになり ます。 (4)同時方程式モデル Y は所得、Cは消費、Iは投資を表す変数としたとき、次のようなマクロ経済モデルを考えることにします。 ( Yi = Ci+ Ii Ci = β0+ β1Yi+ ui ただしIは外生変数(exogenous variable)であるとします。このとき Yi= _{1 − β}β0 1+ 1 1 − β1Ii+ 1 1 − β1ui Ci= _{1 − β}β0 1+ β1 1 − β1Ii+ 1 1 − β1ui と表現できることから、所得Y は確率変数であることがわかります。

回帰式においてY とuの相関をチェックすると次のようになります。 Cov(Yi, ui) = E · u2 i 1 − β1 ¸ =_{1 − β}σ2 1 ̸= 0

このことから、内生変数(endogenous variable)が説明変数となる場合にも、OLS推定量は不偏性も一致性

も持たないことがわかります。

4.

操作変数法

評価版では割愛しています。

5. 2

段階最小

2

乗法

評価版では割愛しています。

6. ivregress

の用例

評価版では割愛しています。

7. ivregress postestimation

の用例

評価版では割愛しています。 ¥

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 foreign

c

ここに示した結果は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) のような形で変数x₁とx₂の値をベースに生成したものです。rnormal(m, s)というのは正規分布に従う乱 数発生関数であり、mは平均値を、sは標準偏差を意味します（[D] functions参照）。すなわち y = 0.5x1+ 2.0x2− 10 (1) という線形の関係を基調に据えて乱数を発生させているわけです。なお、データセット中にはx₁, x₂とは全く 無関係の変数x₃も含めてあります。yに対する線形フィットを実行するに際し、説明変数にx₃も加えたとき にどのような現象が生ずるかを確認するためです。

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

モデルをフィットさせると、それぞれの観測点(x_1i, x_2i, . . .)ごとにモデルによる推定値by_iが算出されます。 その場合、実際の観測値y_iとの差y_i− by_iが残差となるわけですが、rvfplot (residual-versus-fitted plot)

はこの残差の値をフィット値by_iに対向させる形でプロットしたものです。標準偏差10の正規分布を前提に乱

数を発生させているために、それなりにバラツキを持った分布となっていることがわかります。

(2)決定係数R2

それではregressのアウトプットについて見て行くことにしましょう。regress出力の上半分にはANOVA

表が掲載されています。これは変数yの変動（バラツキ）に関する情報を集約したもので、SSは平方和(sum of squares)を意味します。ここでy_i(i = 1, . . . , n)の平均値を¯y、最小2乗法によるy_iの推定値をby_iと書く ことにしたとき、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検定の結果を示しているわけですが、その場合の帰無仮説H₀は

H0: β1= β2= · · · = 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値として示されているわけです。ここで特筆すべきは変数x₃に対するp値が 0.947となっている点です。すなわちβ₃= 0という帰無仮説が全く棄却できないことをこれは示しています。

実際、β₃に関する95%信頼区間（CI: confidence interval） は[−0.21, 0.22]とレポートされており、0とい

う値がその中に含まれている点に注意してください。変数x₃は有意ではないと判断せざるを得ません。

これに対し変数x₁, x₂、及び定数項の場合、p値はそれらの係数値が0とは言えないと主張しているだけなの

で、regressが出力してきたp値は余り参考にはならないでしょう。もし仮に理論モデル等によりx₁の係数

値としては0.5が予想されるとしたとき、postestimation機能の1つであるtestコマンドを使用すればその

仮説を次のようにして検定することができます。

• Statistics ◃ Postestimation ◃ Tests ◃ Test linear hypothesesと操作

• Mainタブ: Test type for specification 1: Linear expressions are equal Specification 1, linear expression: x1 = 0.5

図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)を意味するもので、信頼区間を規定する上での重要な統計量であると言えま す。実際、変数x₁の場合について与えられた標準誤差の値から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 のように算出することができるわけです。係数表においてx₁に対する95% CIは[0.4166288, 0.5356652]と レポートされているわけですが、その内容がこれによって検証できたことになります。

いずれにせよ、得られた係数推定値の信頼度はこれら標準誤差と信頼区間の値から判断することになります。 ここでの例の場合、x₁, x₂に対する信頼区間は十分狭いものとなっていますが、定数項についてはやや幅の 広いものとなっています。データ生成のプロセスを勘案するなら、これ以上は期待できないことは明らかで しょう。

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倍 に増えている点にあります。結果的に信頼区間の幅がそれだけ広くなってしまっています。例えばx₁の係数 値としては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_{倍も大きいように見えるわけですが、これは正しい解釈とは言えません。weight}2_の変数 値自体が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の線形関数として表現したいとし

ます。この場合、特別な指定を行わなければ定数項β₀（切片）を含んだ形のモデル式がフィットされます。す

なわちlength = 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

*13_{i. 演算子については mwp-028 を参照。}

*14メニュー操作：Data ◃ Create or change data ◃ Change contents of variable

. 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の出力に示されているように、均一な分散を仮定した帰無仮説H₀に対するp値が0.0009で

すから、H₀は棄却せざるを得ないわけです。

このように今回設定したモデルについては分散の均一性が主張できないため、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)

mwp-042

anova/oneway - 機能概要と用例

anovaは分散分析の機能を提供する汎用的なコマンドです。これに対しonewayは一元配置ANOVAに特化

し、多重比較機能も包含するなど、より使いやすさを追求したコマンドです。 1. 分散分析と多重比較 2. 分散分析の前提条件 3. 一元配置ANOVA – oneway 4. 一元配置ANOVA – anova 5. 線形回帰モデル 6. 二元配置ANOVA 7. 二元配置ANOVA後の多重比較 8. 反復測定ANOVA 補足１ 補足２ 補足３

1.

分散分析と多重比較

平均値が等しいと言えるかどうかを検定する場合、対象とする標本の数が2つ以下の場合には通常t検定が用 いられます（mwp-041参照）。それでは3つの標本A, B, Cが与えられたとき、t検定を繰返し使用したら何 が悪いのでしょうか？今、検定の有意水準αを5%とすると (1) A-B間の比較で有意差が検出されない確率= 0.95 (2) A-C間の比較で有意差が検出されない確率= 0.95 (3) B-C間の比較で有意差が検出されない確率= 0.95 ですから、3回の検定を通じて有意差が検出されない確率は0.953= 0.86となります。逆に言えば(1), (2), (3)のいずれかで有意差が検出される確率は1 − 0.86 = 0.14となり、正しい帰無仮説を棄却してしまう過誤 （第1種過誤）の確率が本来の5%よりも大きく増大するといった問題が生じます。 c

このため標本数が3以上の場合には分散分析 (ANOVA: analysis of variance) という手法が用いられます。 これはすべての平均値間に差がないことをF 検定によって確認しようとするものです。しかしこの仮説が棄 却された場合に、どの平均とどの平均の間に有意差が認められるかについて、分散分析自体は何ら情報をもた らしません。このため、有意確率の補正を伴う多重比較(multiple comparison)検定を併用することが通常行 われます。

2.

分散分析の前提条件

分散分析の実行に際しては次の要件が満たされていることが前提となります。 (1) 従属変数（応答変数）は量的（区間尺度）データであること (2) 従属変数（応答変数）は正規分布に従うこと (3) 観測データは互いに独立であること (4) 各グループの分散は均一であること 反復測定(repeated-measures) ANOVAの場合には独立性に関する前提条件が成り立たなくなり ます。

3.

一元配置

ANOVA – oneway

要因として想定する因子(factor) の数が1つの場合を一元配置ANOVA (one-way ANOVA)と言います。一

元配置ANOVAに対してはanova, onewayコマンド双方を使用することができますが、本セクションでは

onewayを用いて分析を行ってみます。使用するデータセットはanova1.dtaです。 . use http://www.math-koubou.jp/stata/data11/anova1.dta, clear

全部で24個の血圧(blood pressure)データが記録されていますが、ここではその一部のデータをリスト表示 しておきます。 . list if n <= 3 | n >= 22, separator(3)*1 24. 123 4 23. 139 4 22. 137 4 3. 115 1 2. 121 1 1. 126 1 bp drug 変数drugは薬剤の種別を表すカテゴリ変数で1, 2, 3, 4という4つの値を取ります。このdrugの値ごとに データを整理し、平均値を算出すると次のようになります*2。

*1メニュー操作：Data ◃ Describe data ◃ List data

drug bp 平均値 1 126 121 115 123 125 113 120.5 2 112 123 115 129 106 108 115.5 3 123 112 133 124 130 121 123.8 4 122 132 125 137 139 123 129.7 薬剤の種類によって平均値は微妙に異なるわけですが、これらの差が統計的に有意と言えるかどうかをoneway を使って検定してみます。なお、有意水準αの値としてはデフォルトの5%を用いることにします。また要 約統計量を示す表の作成と多重比較検定の実行に関するオプションも指定することにします。

• Statistics ◃ Linear models and related ◃ ANOVA/MANOVA ◃ One-way ANOVAと操作

• Mainタブ: Response variable: bp Factor variable: drug

Multiple-comparison tests: Bonferroni Output: Produce summary table: X

Bartlett's test for equal variances: chi2(3) = 1.1493 Prob>chi2 = 0.765 Total 1719.625 23 74.7663043 Within groups 1083.16667 20 54.1583333 Between groups 636.458333 3 212.152778 3.92 0.0237 Source SS df MS F Prob > F Analysis of Variance Total 122.375 8.6467511 24 4 129.66667 7.3665913 6 3 123.83333 7.3598007 6 2 115.5 8.9162773 6 1 120.5 5.3572381 6

drug Mean Std. Dev. Freq. Summary of bp . oneway bp drug, bonferroni tabulate 0.260 0.020 1.000 4 9.16667 14.1667 5.83333 1.000 0.383 3 3.33333 8.33333 1.000 2 -5 Col Mean 1 2 3 Row (Bonferroni) Comparison of bp by drug (1) ANOVA表の解釈 ここでの操作ではtabulateオプションを指定しているので、最初にまず因子水準ごとに要約統計量を整理し た表が出力されています。この例では度数(frequency)がいずれも6ということでバランスの取れたデータと

なっていますが、onewayにしろanovaにしろ、バランスの取れていないデータ (unbalanced data)も扱う

ことができます。

要約統計量の表に続く形で出力されているのがANOVA表です。

Bartlett's test for equal variances: chi2(3) = 1.1493 Prob>chi2 = 0.765 Total 1719.625 23 74.7663043

Within groups 1083.16667 20 54.1583333

Between groups 636.458333 3 212.152778 3.92 0.0237 Source SS df MS F Prob > F

SSは平方和 (sum of squares)を意味します。regressの項（mwp-037参照）ではyの変動を X (yi− ¯y)2= X (yi− byi)2+ X (byi− ¯y)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)

と表現しましたが、今の場合、MSSに相当するのが群間変動(between groups) 636.46、RSSに相当するの が群内変動 (within groups) 1083.17、TSSに相当するのが全変動(total) 1719.63です。このとき、それぞれ の変動を自由度(df: degrees of freedom)で割ることによってMS (mean square)すなわち平均平方が各々 212.15, 54.16と算出されます。これらの比を取ったものがF 値で212.15/54.16 = 3.92という値になります。 このF 値を使ったF検定の結果がp値0.0237と示されているわけですが、この場合の帰無仮説はすべての 平均値が等しいとするものです。今の場合、p値は< 0.05ですから帰無仮説は棄却されることになります。 なお、ANOVA表の末尾にBartlett検定の結果が表示されていますが、これは等分散性に関するものです。 等分散性はANOVAの前提条件の一つであるわけですが、今の場合p値は≫ 0.05ですから、この前提条件 に関する問題はないと言えます。 (2)多重比較 上記ANOVA表からはµ₁= µ₂= µ₃= µ₄が主張できないということはわかったわけですが、どのペア間に 有意差があるかを見るためには多重比較検定のプロセスが必要となります。上記操作ではbonferroniとい うオプションを指定しているので、Bonferroni補正を施した形の多重比較検定の結果がANOVA表に引き続 く形で出力されています。 0.260 0.020 1.000 4 9.16667 14.1667 5.83333 1.000 0.383 3 3.33333 8.33333 1.000 2 -5 Col Mean 1 2 3 Row (Bonferroni) Comparison of bp by drug 表示されている行列要素をM_ijと表記することにすれば、M_ij にはµ_i− µ_jの値とそれが0と言えるか否か に関するBonferroni補正後のp値が示されています。この結果からするとµ₄–µ₂間の差のみ有意と判定され ています。多重比較補正の手法としてはBonferroniの他にScheffe, ˇSid´akが選択できます。

4.

一元配置

ANOVA – anova

評価版では割愛しています。

5.

線形回帰モデル

評価版では割愛しています。

6.

二元配置

ANOVA

評価版では割愛しています。

7.

二元配置

ANOVA

後の多重比較

評価版では割愛しています。

8.

反復測定

ANOVA

評価版では割愛しています。

補足１

–

表の作成

評価版では割愛しています。

補足２

–

繰返しのない二元配置

ANOVA

評価版では割愛しています。

補足３

–

二元配置データのグラフ化

評価版では割愛しています。 ¥

mwp-070 table -要約統計情報のテーブル化 tableコマンドはカテゴリ変数の組合せごとに要約統計情報を算出し、結果をテーブル形式で出力します。要 約統計情報の種類としては数多くのものが用意されています。Superrows, supercolumnsを用いたテーブル 構成にも対応しています。 1. tableコマンド 2. 一元表 3. 二元表 4. 三元表

1. table

コマンド

tableコマンドを使用する場合にキーとなるのはカテゴリ変数です。その数がn個の場合にはn元の表 (n-way table)が作成されます（ただしn ≤ 7）。一方、カテゴリ変数の値の組合せに対して個々のテーブルセ ルが対応することになるわけですが、そこに表示できる要約統計情報としては次のようなものが用意されてい ます。 ◦ 度数(frequency) ◦ 平均値 (mean) ◦ 標準偏差 (standard deviation) ◦ 最大値 (maximum) ◦ 最小値 (minimum) ◦ 中央値 (median) ◦ 四分位範囲(interquantile range) ◦ パーセンタイル(percentile) それぞれのセルに複数の統計情報を出力させることもできます。 c

2.

一元表

最初に一元の場合についてその用例を示します。使用するExampleデータセットはauto.dtaです。 . sysuse auto.dta*1 (1978 Automobile Data) こ の デ ー タ セ ッ ト 中 に は 1978 年 に 米 国 で 販 売 さ れ た 74 車 種 の デ ー タ が 収 納 さ れ て い ま す 。変 数 displacementは排気量を表す連続変数ですが単位が立方インチでわかりにくいので、最初にそれをリッター に変換します。

. generate disp1 = 16.4*displacement/1000*2

変数disp1の最大値と最小値をチェックしておくと

. summarize disp1*3

disp1 74 3.235676 1.50613 1.2956 6.97 Variable Obs Mean Std. Dev. Min Max . summarize disp1 より、排気量の値は1.3lから7.0lまでの範囲に及ぶことがわかります。この変数disp1から次のようなカテ ゴリ変数classを生成してみます。 class 排気量区分 0 disp1 ≤ 2.0 1 2.0 < disp1 ≤ 3.0 2 3.0 < disp1 ≤ 4.0 3 4.0 < disp1

. generate class = irecode(disp1, 2.0, 3.0, 4.0)*4

参考までに先頭から5台の車種についてdisp1とclassの値をリスト出力しておきます。

. list make disp1 class in 1/5

5. Buick Electra 5.74 3 4. Buick Century 3.2144 2 3. AMC Spirit 1.9844 0 2. AMC Pacer 4.2312 3 1. AMC Concord 1.9844 0 make disp1 class

*1メニュー操作：File ◃ Example datasets ◃ Example datasets installed with Stata と操作しダウンロードする。

*2メニュー操作：Data ◃ Create or change data ◃ Create new variable

*3メニュー操作：Statistics ◃ Summaries, tables and tests ◃ Summary and descriptive statistics ◃ Summary statistics

このカテゴリ変数classを用いて一元の表を作成してみましょう。その際、変数mpg (miles per gallon)に 関する要約統計情報をいくつか出力させることにします。

• Statistics ◃ Summaries, tables, and tests ◃ Tables ◃ Table of summary statistics (table)と操作

• Mainタブ: Row variable: class

Statistics 1: Count nonmissing mpg Statistics 2: Mean mpg

Statistics 3: Standard deviation mpg

図1 tableダイアログ- Mainタブ 3 19 16.4211 3.48514 2 19 19.4211 2.734873 1 13 20.6154 3.819652 0 23 27.2609 5.100957 class N(mpg) mean(mpg) sd(mpg) . table class, contents(count mpg mean mpg sd mpg )

定義された排気量クラスごとに変数mpgの度数、平均値、標準偏差が出力されていることがわかります。な

お、変数classに対し値ラベルを設定してやればよりわかりやすい出力とすることができます*5。

. label define class 0 "<2000cc" 1 "2000-3000cc" 2 "3000-4000cc" 3 ">4000cc" . label values class class

. table class, contents(count mpg mean mpg sd mpg )

>4000cc 19 16.4211 3.48514

3000-4000cc 19 19.4211 2.734873

2000-3000cc 13 20.6154 3.819652

<2000cc 23 27.2609 5.100957 class N(mpg) mean(mpg) sd(mpg) . table class, contents(count mpg mean mpg sd mpg )

さらにtableダイアログのOptionsタブ上で表示書式を明示してやれば数値の表示様式も調整できます。

• Optionsタブ: Override display format for numbers in cells: %9.2f

図2 tableダイアログ- Optionsタブ >4000cc 19 16.42 3.49 3000-4000cc 19 19.42 2.73 2000-3000cc 13 20.62 3.82 <2000cc 23 27.26 5.10 class N(mpg) mean(mpg) sd(mpg)

3.

二元表

評価版では割愛しています。

4.

三元表

評価版では割愛しています。