Stataの「使える」adoファイル Stataの「使える」adoファイル(ver1.03 2006.4.8)
統計をするためのソフトウェアです．SPSSよりも値段が安く，機能も豊富です．アカデミックで600ドルくらいだっただろうか（Intercooled Stataの場合）．SPSSで同じ機能をそろえるよりもずいぶん安上がりになります．欠点も少しあります．日本語が使えません．一番いいところは，最新バージョンでMac版があるところです．その他，コマンドライン入力を受け付けているだけあって，データの加工がやりやすい．というのが利点だと思います．
Stataは，ユーザがボランタリーに提供してくれるプログラムをインターネット経由でインストールすることで，多機能化します．ただ追加分のプログラムに関する情報は日本語ではあまり見かけませんので，英語が苦手な人向けにここで少しだけ紹介します．
インターネットに接続した状態で，以下のように打ち込み，あとは指示に従えば，本体プログラムとヘルプファイルが自動的にインストールされます．
. findit _____
「 ____ 」には下で紹介するプログラム名を入れてください．
変数名を一括して変換するためのコマンドです．以下のようにすれば二つ同時に変数名を変更できます．
. sysuse auto,clear (1978 Automobile Data) . renames make price \ maker carprice
また，複数の変数に一括して同じprefix（接頭辞?）あるいはsuffix（接尾辞?）をつけることもできます．
. sysuse auto, clear (1978 Automobile Data) . renames make price mpg, suffix(0)
とやると，上の三つの変数の接尾に「0」をつけてくれます．変数のリコードする前準備などに使えるコマンドだと思います．
とやれば，欠損値のカウント，比率を返します．変数名を入れなければ，メモリ中の全ての変数の欠損値のリストを返します．この場合，質問デザイン上欠損値としているものもカウントしますので，注意しましょう．
変数三つのクロス表を作ります．
クロス表で，各セルの期待値を観察値の下に返します．aオプションでadjusted residual（調整残差？）も付けてくれますので，どのセルが非独立に貢献しているのかが一目で分かります．
連続変数を自由に区切ってカテゴリー化するためのegenコマンドの関数+オプションです．実際の使い方は（説明するとややこしいので）ヘルプファイルをご覧になってください．
変数を中心化します．たとえば
. sysuse auto, clear (1978 Automobile Data) . center price, prefix(c)
とすれば，cpriceという，priceを中心化（値から平均値を引いた）変数を作ります．prefixオプションをつけない場合，自動的にc_が変数名の頭につけられます．standardizeオプションをつければ，標準化された変数（zスコア）を作ります．同じ機能でzscoreというのもあります．prefixはこの場合stub(z_)というふうに付けることができます．
職業分類のユーティリティ．ILOによるISCO（International Classification of Occupations）分類（四桁コード）を他のインデックス（ISEIなど）や分類に変換してくれます．Erikson，GoldthorpeらによるEGP分類への変換などが便利．
iscoというプログラムもありますが，こちらはISCO88には対応していません．iskoを使いましょう．
使用法については，ヘルプをじっと見てください．
回帰分析の後にモデルの期待値を出力するには最初からStataに備わっている「predict」コマンドを使いますが，「prvalue」コマンドでは，説明変数の個々の値に対応するモデル期待値と95%信頼区間を返します．
. u http://www.ats.ucla.edu/stat/stata/notes/hsb2, clear
. quietly reg read write
. predict pread
(option xb assumed; fitted values)
. scatter pread write
. prvalue, x(write=50)
regress: Predictions for read
Predicted value of y: 50.43866 95% ci: ( 49.24738, 51.62994)
write
x= 50
つまりライティングスコア50点の場合のリーディングスコアの期待値は50.44で，95%の確率で49.25点から51.63点である．
回帰分析の後に，推定結果をより見やすい形にアレンジして返してくれます．特にダミーを使った回帰分析の場合，デフォルトでは結果に「I_variable」のかたちで返されて，どのダミーがどのカテゴリー値を示しているか分かりませんが，reformatを使うとカテゴリー値のラベルを表示してくれます．
. sysuse auto, clear
(1978 Automobile Data)
. xi:reg price i.foreign
i.foreign _Iforeign_0-1 (naturally coded; _Iforeign_0 omitted)
Source | SS df MS Number of obs = 74
-------------+------------------------------ F( 1, 72) = 0.17
Model | 1507382.66 1 1507382.66 Prob > F = 0.6802
Residual | 633558013 72 8799416.85 R-squared = 0.0024
-------------+------------------------------ Adj R-squared = -0.0115
Total | 635065396 73 8699525.97 Root MSE = 2966.4
------------------------------------------------------------------------------
price | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_Iforeign_1 | 312.2587 754.4488 0.41 0.680 -1191.708 1816.225
_cons | 6072.423 411.363 14.76 0.000 5252.386 6892.46
------------------------------------------------------------------------------
. reformat
REGRESS formatted output
Outcome variable: price (Price), n=74
-----------------------------------------------------------------------------
Covariate Coef. Std. Err. P>|t| 95% Conf. Interval
-----------------------------------------------------------------------------
Car type
Domestic* 0
Foreign 312.259 754.449 0.680 (-1.2e+03 to 1816.225)
Constant 6072.423 411.363 <0.001 (5252.386 to 6892.460)
-----------------------------------------------------------------------------
* Baseline category
他のオプションについてはヘルプファイルを参照してください．
因子分析を行うときに便利なコマンドを集めたパッケージです．回転後に（デフォルトで）0.3以上の因子負荷量のみ表示してくれるfadisp，斜交回転後に因子間の相関を返すfacorrなど．サンプルログを載せておきます．
. u http://www.ats.ucla.edu/stat/stata/examples/sws5/planets, clear
(Solar system data)
. factor rings logdsun - logdense, ml nolog
(obs=9)
number of factors adjusted to 3
(maximum likelihood factors; 3 factors retained)
Factor Variance Difference Proportion Cumulative
------------------------------------------------------------------
1 3.81295 1.92230 0.6596 0.6596
2 1.89065 1.81394 0.3271 0.9867
3 0.07671 . 0.0133 1.0000
Test: 3 vs. no factors. Chi2( 18) = 141.77, Prob > chi2 = 0.0000
Test: 3 vs. more factors. Chi2( 0) = 2.26, Prob > chi2 = .
note: Above tests may not apply, see manual.
Factor Loadings
Variable | 1 2 3 Uniqueness
-------------+-------------------------------------------
rings | 0.88934 0.41910 0.18282 0.00000
logdsun | 0.26048 0.88924 0.04017 0.13979
lograd | 0.99284 0.11224 -0.04094 0.00000
logmoons | 0.84310 0.46130 0.01433 0.07619
logmass | 0.99645 -0.08279 0.01506 0.00000
logdense | -0.51450 -0.83187 0.19891 0.00372
. quietly rotate, promax factors(3)
. facorr
Promax Factor Correlations
1 2 3
1 1.0000
2 0.4695 1.0000
3 0.2719 0.5875 1.0000
. fadisp rings logdsun - logdense
Rotated Factor Loadings
1 2 3
rings 0.71595 . 0.30884
logdsun . 0.86710 .
lograd 0.95975 . .
logmoons 0.64038 0.40886 .
logmass 1.06033 . .
logdense . -0.99707 .
fadispの表示下限を変更したいときは，
. fadisp varnames, blank(.4)
のようにオプション指定します．
3つ以上のダミー変数のインタラクションを許可するコマンドで，ログリニア分析をする際などに便利なコマンドです．
xi3: logit y i.cate1*i.cate2 ...
のように，いくらでも複雑にログリニアできます．普通は，lrtestや次に紹介するfitstatなどの診断コマンドと組み合わせて最適なモデルを探します．
判別分析determinant analysisを行うためのファイルです．
Stataでパス解析を行うためのコマンドです．例えば「var1→var2→var3」というモデルで推定したければ，
. pathreg (var2 var1) (var3 var2)
のように打ち込んでください．ただし，なぜか「if」オプションが使えないので，条件に合わないcaseを最初に削除しておく必要があるようです．試しにEQSを使って同じモデルを作ってみたところ，同じパラメータ推定を返しました．
ブール代数分析を行うためのプログラム．
χ二乗分布を利用した適合度検定のためのコマンド．ある母数を仮定し，そのモデルからの乖離の度合いをχ二乗検定する．
. sysuse auto
(1978 Automobile Data)
. ta foreign
Car type | Freq. Percent Cum.
------------+-----------------------------------
Domestic | 52 70.27 70.27
Foreign | 22 29.73 100.00
------------+-----------------------------------
Total | 74 100.00
. csgof foreign, exp(50 50)
+----------------------------------------+
| foreign expperc expfreq obsfreq |
|----------------------------------------|
| Domestic 50 37 52 |
| Foreign 50 37 22 |
+----------------------------------------+
chisq(1) is 12.16, p = .0005
つまり「国産車と外車の割合は半々である」というモデル（仮説）の適合度は5%水準で，有意ではない．
. csgof foreign, exp(70 30) +----------------------------------------+ | foreign expperc expfreq obsfreq | |----------------------------------------| | Domestic 70 51.8 52 | | Foreign 30 22.2 22 | +----------------------------------------+ chisq(1) is 0, p = .9595
これに対して「国産車と外車の割合は7:3である」という仮説の適合度は，当然高くなります．
このような場合だとふつうにttest foreign=0.5 (or 0.3)でもいいのでしょうが，値が三つ以上ある場合だとcsgofは便利かもしれません．ただ，たぶん一つの変数の分布しか検定できないので，多変数の場合は（たいてい）線形モデルを作った後でlrtest等をすることになると思います．
LS重回帰モデルの多重共線性multicollinearityの程度を計算します．デフォルトではパラメータ推定後にvifコマンドで見ることができますが，collinの方は固有値や相関係数も返してくれます．モデル全体の決定係数が大きいのに個々の偏回帰係数が有意にならないと思ったら，とりあえずtoleranceが最も小さい変数をモデルから除いてみましょう．
モデルの診断コマンドです．デフォルトでlrtestが使えますが，fitstatの方がいろんな情報を返してくれるみたいです．回帰分析の後にオプション抜きで打ち込むと，主効果（と交互作用項）の回帰係数が全部ゼロのモデルと比べたときの値を返します（尤度比検定）．このときp値は小さい方が良いです．p値が小さいと，「モデルの説明変数の少なくとも一つは有意である」ということになります．
他方，あるモデル（モデル0とする）の推定のあとに
fitstat, saving(m0)
と打ち込み，さらに別のモデル(モデル1とする）の推定を行い，その後で
fitstat, using(m1)
とすると，モデル0とモデル1のdeviance（尤離度?）の比較検定をします．
例えば，モデル0にシンプルなモデルを持ってきて，モデル1で説明変数を追加してからfitstatをし，尤度比検定のp値が大きければ二つのモデルの説明力にたいした差がないことになり，「余計なことしないでシンプルなモデル（m0）を採用した方がいい」ということになります．ログリニアの場合，飽和モデルsaturated modelとの差を考えなくてはなりませんから，説明変数が多い場合，複雑なモデルからはじめて尤度比検定のp値が大きいモデルを選ぶ，というふうにすると思います．この場合のp値は「飽和モデルと当該モデルの適合度が離れている危険度」ということになります．離れていては拙いので，この場合はp値は大きい方がよい，ということです．
ちなみにオプション抜きで回帰モデルの推定後にlrtestをすると，飽和モデルとの比較（適合度検定）を行います．
これもログリニア・モデルで使います．
ipf, fit(cate1*cate2+cate2*cate3+cate3*cate1)
のように「+」でつなげて使います．この場合，たとえばカテゴリー変数がx1，x2，x3の三つである場合，x1+x2+x3は「三つの変数が独立」なモデルです．x1とx2のインタラクション項を加えたい場合は，
ipf, fit(x1*x2+x3)
のようにします．ちなみにipfはモデルのフィットの具合を見るためのコマンドなので，何もしないでも適合度検定の数値を返してくれますが，save(ファイル名)というオプションを追加することで，期待値（モデルセルカウント）を加えたデータファイルを別に作ってくれます．ベストのモデルが作れたと思ったら，最後にsaveオプションをやってみて，観察値との一致度を目で確かめてみましょう．
これは別に適合度を計算するものではないのですが，ログリニア・モデルをStataのグラフィック機能で視覚的に描画してくれるものです．詳しくはヘルプをご覧になってください．
J検定のためのプログラム．非常に簡単な文法だが、ただしifオプションなどは使えない．（下の例は適当．）
. sysuse auto . nnest price mpg (weight) M1 : Y = a + Xb with X = [mpg] M2 : Y = a + Zg with Z = [weight] J test for non-nested models H0 : M1 t(71) 2.72324 H1 : M2 p-val 0.00813 H0 : M2 t(71) 0.57468 H1 : M1 p-val 0.56732 Cox-Pesaran test for non-nested models H0 : M1 N(0,1) -3.51752 H1 : M2 p-val 0.00022 H0 : M2 N(0,1) -0.61988 H1 : M1 p-val 0.26767
回帰分析の後に，結果をASCIIテクストファイル形式の別ファイルに出力する機能．
. reg y x . outreg x using model1, se nolabel bdec(3) rdec(3) sigsymb(*,**,+) 10pct coefastr ctitle(I); . reg y x z using model1, se nolabel bdec(3) rdec(3) sigsymb(*,**,+) 10pct coefastr ctitle(I) append;
のように書けば，二つの回帰モデルの結果を並列して____.outというファイルに書き出してくれます（ファイルの在処が不安な場合，やる前にpwdでカレント・ディレクトリを確認しましょう）．最後のappendオプションは，すでにあるファイル（ここでは「model1.out」）に追加して書き出しますよ，ということ．replaceオプションをつけると，追加せずに上書きされます．
「outreg」と同じく，いろんな回帰分析の後に結果をまとめてテクスト形式に出力します．オプションが多すぎて使い方が難しそうなので，興味があればインストールした後にヘルプファイルを参照してください．mkcorrという仲間もいるみたいです．説明は省きます．
「summary」のLaTeX出力．単純ですので普通にやってみてください．
ほぼtabulateコマンドと同じように機能しますが，クロス表をLaTex形式に変換してから返します．（resultウィンドウ内にテクストデータの形で返しますので，そのままエディタにコピペです．）latabstatは，同じことのtabstat版です．例えば次のように打ち込めば，論文でも使いやすいでしょう．
. latabstat price, by(foreign) stat(n mean sd) f(%9.1fc)
なにげにこれが一番おいしい機能かなあと思います．回帰分析（OLS回帰でもロジットでもポワッソンでもよし）を行った後にこれを打ち込むと，結果をLaTeXのテーブルにしてくれます．levelおよびdetailオプションで，論文に載せるときのようなテーブルが出来上がります．
. sysuse auto (1978 Automobile Data) . quietly regress price mpg . outtex, level detail
のように打ち込み，出力されたものをLaTeXでレンダリングすると，以下のような見栄えになります．（t値がほしいところですが，どうにかしたら出力してくれるのだろうか？）
以上．
私は数式を見ただけでウッと来てしまうよくいる社会学者のひとりで，統計に関してはほとんど素人のようなものですが，周囲にStataユーザが増えるといいこと（使い方を教えてもらえたり）があると思い，このような一覧を作ることにしました．間違い等あれば，お知らせください．