Stataでのspline補完
https://www.lightstone.co.jp/stata/files/201607spline.pdf
https://www.lightstone.co.jp/stata/files/201608spline.pdf
https://www.lightstone.co.jp/stata/files/201702spline.pdf
スプライン曲線というのは、区間ごとに多項式で表現された曲線のこと。
単純な直線で近似できるのが一番だが、実際のデータ分析ではそううまくいかないことも多い。Stataでは”mkspline”コマンドで、線形あるいは制限3次スプライン(restricted cubic spline, RCS)を行うことができる。
線形スプライン
まずはデータセットの読み込み。
use http://www.stata-press.com/data/r14/mksp1, clearIDと年齢、教育歴（年数）と収入、log（収入）が横データで与えられている。
以下のコードで、区間を7個の分割する。
mkspline age1 24 age2 30 age3 40 age4 48 age5 59 age6 65 age7 = age区間を\(n\)個のノット（結び目）で分割すると、\(n+1\)個の変数が生じる。構文としては、
mkspline 新変数1 ノット1 新変数2 ノット2 … 新変数n+1 = x
とする。新変数の値は、以下のようになる。
$$\begin{eqnarray}\text{age}1 & = & \min(\text{age}, k_1) \\ \text{age}2 & = & \max(\min(\text{age}, k2), k_1)-k_1 \\ \text{age}3 & = & \max(\min(\text{age}, k_3), k_2)-k_2\\ \cdots \\ \text{age}6 & = & \max(\min(\text{age}, k_6), k_5)-k_5 \\ \text{age}7 & = & \max(\text{age}, k_6)- k_6\end{eqnarray}$$
これは結構複雑で、詳しくは上記のpdfをグラフを見るとよく分かる。あるいは、コマンドを実行してみて、どのような変数が作られるのかを確認すると良い。
この変数をもとに、回帰を行う。
reg lninc age*
* あるいは
reg lninc age1-age7ワイルドカードを用いても良いし、\(\text{age}1-\text{age}7\)としても良い。そして、予測を行う。上記のpdfには数式を用いて、何を行っているのか、あるいはどのように変数の連続性が担保されているのかが詳細に記載されている。
ノット間で有意に傾きが変化したかどうかは、以下のように”lincom”を使って推定を行う。
lincom age4-age3predict yhat, xbここで、最悪なことに、日本語マニュアルのpdfの指示に従って図を描くと、以下のような図が出てくる。
これが欲しいものでないことは明らかなので、こちらを参考に、正しいコードは以下のようになる。
twoway scatter lninc age || mspline lninc age
このグラフは制限3次スプラインになっている。直線で補完する方法を色々調べたが、見つけられなかった。
ノットの数の選択について
上の2番目のpdfに記載されている。基本的に臨床的に妥当なものを選択すれば良いが、AICやBICなどの数値で決定することもできる。
重回帰やロジスティック回帰を行う場合
この場合、”xblc”というadoファイルをインストールする必要がある。
search xblcで”st0215_1″をクリック、さらに”click here to install”でOK。
復習として、以下のように、制約3次スプラインを実行する。
mkspline ages = age, cubic nknot(4)この場合、ノット数は\(4\)になる。つまり、作成される変数は\(3\)個である。これを用いて回帰、予測まで一気に行う。
reg lninc ages*
predict lninchat, xb
sort age
twoway line lninchat age || scatter lninc age xline(16 35 48 68), xlabel(, grid)
重回帰
ここまでは前と同様。まずは重回帰分析を行う。”educ”を加え、センタリングを行った上で回帰を行う。ここで、センタリングを行うのは、\(= 0\)で平均値を示すようにするためである。スプライン曲線では\(2\)軸に表さない共変量は、連続変数では値が\(0\)のとき、あるいはカテゴリー変数では値が最小値のときのグラフがそれぞれ描かれる。すなわち、連続変数についてはセンタリング後の変数を用いることで、センタリング前の変数の平均値でのグラフを描くことができるようになる。カテゴリー変数でも、推定コマンド（ib#.）などを用いてベースを指定することで、その値でのグラフを描写できる。
sum educ
gen educ_c = educ-r(mean)こうしておいて、回帰分析をこのセンタリング後の変数で行う。sumでeudcを確認しておかないと、r(mean)が欠損値になってしまうので注意。
reg lninc ages* educ_c後は、通常と同様。
levelsof age
xblc ages*, covname(age) at(`r(levels)') line
levelsofコマンドで指定する変数がdouble型の小数の場合、xblcを用いるとエラーが起こる。この場合、以下のようにもとの変数をfloat型に変換すると良い（double型の方は別名の変数を作り、退避させておく）。
gen age_double = age
recast float age, force
levelsof ageロジスティック回帰
今までと違うデータを用いる。
use http://nicolaorsini.altervista.org/data/pa_luts, clearはサーバーエラーで使えないので、以下を用いる。
webuse lbw, clearここでは”age”をノットで区切ることにする。
mkspline ages = age, nknots(4) cubic以下は一気にやってしまう。
logit low ages* i.smoke
levelsof age
xblc ages1-ages3, covname(age) at(`r(levels)') reference(25) eform lineオッズ比は基準となるオッズが必要になるので、reference()オプションで基準となる観測点をデータセット内にある数値から指定する。
コメント