Mplusでロジスティック回帰分析をする方法である。炭鉱夫の年齢と炭鉱で働くと生じる健康被害の息切れの関係を推定した例を使用する。
解説するのはこちらの例。
https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.367.3127&rep=rep1&type=pdf
データはこちらからダウンロードできる。 https://www.statmodel.com/mplusbook/chapter5.shtml
title: logit regression for coal miners data: file =coalminer.dat; variable: names = x u w; categorical = u; freqweight = w; define: x = x/10; model: u on x; analysis: estimator = ml; output: tech1 sampstat standardized;
ロジステック回帰分析は従属変数がカテゴリカル変数(この例では2値)であるため、CATEGORICALオプションで指定する。このコードの場合にはcategorical = u
がそれに相当する。
また、いくつかオプションが書き加えられている。freqweightオプションは、頻度(ケース)重み付けの変数指定でする。頻度のウェイトの値は整数でなければならない。(参照)
コードの解説だけではよくわからないので、元データをみてみよう。
データをみると炭鉱夫22歳のところはYesは1952,Noが16となっている。
分析のデータであるcoalminer.dat
をみると、スタック形式のデータセットになっている。
22 0 1936 22 1 16 27 0 1759 27 1 32 32 0 2040
つまり、ケース数を重みづけとして掛けているわけだ。ロジスティック回帰分析を習う時にはクロス表の拡張ですよ、と習うはずなので基本的な勉強をした人にとっては、それほど不思議なデータではないはずだ。スタック形式のデータ表現は三重クロス表も列を増やせばよいので、1ケース1行のデータ形式でなくても、クロス表があれば分析ができる。もともとロジスティック回帰分析とはそういうものなので、再分析をする時などに知っておくとよいのかもしれない。
普段よく使用する1ケース1行のデータセットの場合は頻度の重み付けが不要なのでfreqweight
を指定しなければよいだけである。
DEFINEコマンドは新しい変数を作成するコマンドだ。独立変数の値を10分の1にするという定義がされているが、これをすると推定結果が10倍になる。推定結果が異なるのでオッズ比も違いが出てくる。1年刻みであったものを10分の1にするということは10年でどのくらい違ってくるか、というデータの解釈になる。
結果は以下のように出る。
MODEL RESULTS Two-Tailed Estimate S.E. Est./S.E. P-Value U ON X 1.025 0.025 41.758 0.000 Thresholds U$1 6.564 0.124 52.873 0.000 LOGISTIC REGRESSION ODDS RATIO RESULTS U ON X 2.787
オッズ比はロジスティック回帰なので、対数でと計算できる。
xが1単位(10年)増えると、炭鉱夫の息切れ症状が現れる確率が2.79倍になるという解釈だ。1年、つまりdefine: x = x/10;
を指定しない場合は1年のオッズ比が出て、結果は1.108倍になる。10年単位が知りたければ、1.10810の10乗を計算すれば、同じ計算をしていることになる。ただ、小数点の四捨五入の関係で、実際に計算すると、2.788673となり、出力と0.001の差が出る。通常使う分には問題ないようには思うが、DEFINEを知っておくといいこともありそうだ。
通常のロジスティック回帰分析は最小二乗法で推定するが、Mplusの例は最尤法なので、推定値は変わらないが、標準誤差が異なり、場合によってはP値も異なる時がある。どちらが良い推定方法なのかは、標準誤差が小さい方がよいとされている。要するにデータによって異なるということである。なお、Mplusでは切片ではなく閾値(Thresholds)が出力されるが、これは符号が違うだけである。