Stataはマイクロデータの処理に優れているのですが,マイクロデータを使うとさまざ
まな推定を行うことが可能になります.サンプルが多い場合には自分で尤度関数を書いて 最尤法を使って推定することもあるかもしれません.もちろん,Stataでも自分で最尤法を 使って推定することができます(収束計算の方法がすくないとの指摘もありますが…).こ こでは,全体の対数尤度がそれぞれの観測値の対数尤度の和で表現できるような,もっと も簡単な場合の最尤法の,もっとも初歩的な使い方を紹介します.それゆえ,パネルデー
タのrandom effectのような,ちょっとややこしいケースは除きます.
Stataの最尤法プログラムとしては,lf, d0, d1, d2の4通りが用意されていますが,
使いやすさと計算の正確さ・速さから,lfの使用が推奨されています.もっとも,lfは全 体の対数尤度が観測値の対数尤度の和で表現できる場合にしか対応していません.いくつ かの観測値間に相関があるようなケースでは,それらをまとめてひとつの観測値としてし まい(つまり横に並べてしまい),対数尤度の定義の仕方を工夫したほうが楽かもしれま せん.
最尤法の使い方はさきほどの非線形最小2乗法とよく似ています.ただし,ここでは説 明変数の指定はあとで,関数の指定・被説明変数と説明変数の指定・推定,という順序を 採ります.ここで扱うような,観測値間の独立性が仮定されているケースではlfオプショ ンを使います.まず一般形を書いておくと次のようになります.
program define 式名
args lnf theta1 theta2 ...
tempvar 一時変数名 対数尤度関数の定義
quietly replace ‘lnf’ = ...
end
ml model lf 式名 (被説明変数 = 説明変数名)...
ml check ml maximize
これらの途中で初期値の設定などのオプション指定もできますがここでは省略します.
最尤法では式の名前に特に制限はありませんが,コマンド名になっているものは避けるべ きです.lnfは対数尤度です.programの定義文の最後でquietly replace ‘lnf’ = の 形で締めくくることになります.theta1, theta2,... は説明変数(と定数項)の線形和を 表します.推定されるのは,thetaたちの係数になります.標準偏差など,定数を推定す るばあいにもこれを使います.推定は収束計算をともなうので,推定される値の域が実数 全体になるように工夫します.被説明変数とは,ここでは係数が付かない変数を指します.
対数尤度関数を定義するときには$ML y1, $ML y2,... といった形で表現されます.一時変 数とは,対数尤度関数を定義するためだけに一時的に定義する変数です.対数尤度関数が 複雑な形をとるとき,一時変数を使って見通しをよくすることができます.対数尤度関数 を定義したあと,被説明変数と説明変数を指定します.=でつながれていますが,これに 特に意味はない(最も簡単なケースでは意味がありそうに見える)だけなので気にしない でください.
最尤法の使い方はやってみないとわかりにくいところがあるのですが,ここでひとつ例 をだします(Stataでの最尤法の入門書,Gould and Sribney 1999, 27ページ35).Weibull 分布を仮定したサバイバル分析を行うケースです.この場合,対数尤度関数は以下のよう にかけます.
lnfi= (tiexp(xiβ))exp(s)+di(s−xiβ+ (exp(s)−1)(lnti−xiβ))
ここで推定したいのは,説明変数の線形結合の係数β,Weibull分布のパラメタsです.ま た,ti, diは係数が付かないのでここでは被説明変数と呼べます.いま,説明変数の線形結 合xiβをtheta1,パラメタsをtheta2とし,tiを1つめの被説明変数,diを2つめの被 説明変数と考えてやれば,Stataで対数尤度関数を定義することができます.推定される係 数はいずれも(sも),すべての実数値を取りうることに注意してください.exp(s)に一時
35Gould, William and William Sribney. 1999. Maximum likelihood estimation with STATA. Stata Press, Texas.
変数p,右辺第1項に一時変数M,lnti−xiβに一時変数Rを割り当てて,式を読みやすく してやると,Stataによる尤度関数の定義はつぎのようになります.
program define myweib version 6.0
args lnf theta1 theta2 tempvar p M R
quietly generate double ‘p’ = exp(‘theta2’)
quietly generate double ‘M’ = ($ML_y1*exp(-‘theta1’))^‘p’
quietly generate double ‘R’ = ln($ML_y1)-‘theta1’
quietly replace ‘lnf’ = -‘M’+$ML_y2*(‘theta2’-‘theta1’+(‘p’-1)*‘R’) end
対数尤度関数を定義するときにはgenerate, replaceはよく使いますが,いずれもオプ ションquietlyとdoubleを使っていることに注意してください.これらがないと出力がひ どくなったり,収束計算の精度が非常に落ちたりします.また,パラメタや一時変数はシ ングルクォーテーションでくくっていることにも注意しましょう.ここではまだ説明変数 はなんでもいいことになっています.
次に説明変数と被説明変数を指定します.変数名はtiがtime,diがdied,xiとして定 数項とageであるとしましょう.すると,
ml model lf myweib (time died = age) ()
となります.2個目のカッコはtheta2に対応していますので省略はできません.コロン を使ってそれぞれに式の名前をつけることもでき,その場合には,
ml model lf myweib (timeeq: time died = age) (sigma: )
などとなります.被説明変数・説明変数の順序は,対数尤度関数の定義のときの$ML y1,
$ML y2,...,theta1, theta2,...,の順序に従います.被説明変数は=の左側,説明変数は 右側に来ますが,説明変数だけをカッコのなかにいれることもできます.被説明変数だけ をカッコのなかにいれることはできません.被説明変数をいれるカッコは順序をまちがわ なければどこでもいいことにはなっていますが,よくわからなければ,すべて最初のカッ コの等号の左側に書いておけばよいでしょう.もし,説明変数の線形結合が定数項を含ま ない場合には次のようになります.
ml model lf myweib (timeeq: time died = age, nocons) (sigma: )
ここまで書けば対数尤度の定義は終わりです.あとは最大化する命令を書けばいいこと になります.
ml check ml maximize