StataでのTrajectory解析
こちらのリンクから。全文訳ではないので注意。
要約
グループごとの軌道解析は対象とする集団の発達過程や行動過程の差を見るために用いられる（経済格差による小児の問題行動など）。この記事ではStataのコマンド”traj”についての解説になる。基本的に用いるのはlongデータであり、用いることのできるモデルとしてはCesored nomal, Poisson, zero-inflated Poisson, そしてベルヌーイ分布になる。
導入
発達の軌跡は年齢や時間経過に伴うアウトカムの経過を測定するものである。発達の軌跡の研究は発達心理学や異常心理学、精神医学、社会学や犯罪学におけるライフコース研究、医学や老年学における身体的・生物学的成果の中心的テーマである。これらの現象を研究するために様々な統計的手法が用いられている。ここではグループベースの軌跡モデルを推定するためのStataプラグインを導入する。
適切に定義された確率分布の有限な混合を用いて、発達の軌跡をモデル化するグループベースのアプローチは集団内の個人の軌跡の特徴的なクラスタを識別するための柔軟で容易に適用できる方法を提供する。
興味の主要な概念は年齢あるいはタイムポイントによるアウトカムの分布である。そこで、アウトカムの軌道分布を\(P(Y_i\mid age_i)\)で表す。ランダムベクトル\(Y_i\)は各個人\(i\)における縦方向の行動アウトカムを表し、ベクトル\(age_i\)は各個人\(i\)のアウトカムが測定された時点での年齢を表す。
グループベースの軌道モデルは、軌道の母集団分布が未知の次数\(J\)の有限の混合から生じると仮定している。つまり、個人\(i\)の尤度はグループ数\(J\)によって条件付けられており、以下のように$$P(Y_i\mid age_i) = \sum_{j = 1}^{J}{{\pi}^j \cdot P(Y_i\mid age_i, j:{\beta}^j)} \tag{1}\label{1}$$と表される。ここで\({\pi}^j\)はグループ\(j\)に属する確率であり、あるグループ\(J\)に属する個人\(i\)の条件付き分布は特にグループ固有の軌道の形状を決定する未知のパラメータベクトル\({\beta}^j\)によって指標化される。軌道は年齢（あるいはタイムポイント）の最大\(5\)次の多項式によってモデル化される。与えられた\(j\)について、条件付き独立は、\(T\)ポイントの期間それぞれに置ける測定の要素\(Y_i, y_{it}\)を逐次合成することによって仮定される。これを以下のように記載する。$$P(Y_i\mid age_i, j; {\beta}^j) = \prod_{i=1}^{T}{p(y_{it}\mid age_{it}, j;{\beta}^j)} \tag{2}\label{2}$$ここで\(p(\cdot )\)はグループ\(j\)に属する個人\(i\)のタイムポイント\(t\)における変数\(y_{it}\)の条件付き確率分布である。Stataの”traj”ではこの\(p(\cdot)\)のモデルについて、\(3\)の確率分布が用意されている。ひとつは打ち切り正規分布で、繰り返し測定され、おおよそ連続的な分布に対して使用される。スケールが最小か最大になったときに打ち切りとされる。zero-inflatedポアソン分布は、縦のカウントデータ（年齢による逮捕など）を解析するときに用いられる。アウトカムが二値を取る場合に、ロジスティック分布を用いる。
共変量をモデルに含めることも可能で、\(x_i\)を時間によらない共変量として、この効果は以下のソフトマックス関数で与えられる。$$\pi_j(x_i) = \frac{e^{x_i \theta_j}}{\sum_{j}{e^{x_i \theta_j}}}$$時間依存の共変量が軌道自体に及ぼす影響は、他の共変量を含まない基本モデルにおいて軌道の形状を定義する年齢（またはタイムポイント）の多項式関数の仕様を一般化し、時変の有無にかかわらずそのような共変量を含めることによってモデル化される。すべてのパラメータ効果推定値は、軌道群に特化したものである。これにより、年齢や時間に関するパラメータ推定値が軌道グループ間で自由に変化するだけでなく、軌道の仕様に含まれる他の共変量に関するパラメータ推定値も自由に変化させることができる。
インストール
前置きはここまで。ここからは実際のコマンドを触る。インストール方法は以下。
net from https://www.andrew.cmu.edu/user/bjones/traj
net install traj, replaceこのままインストールしようとしたが、
package does not contain traj.plugin for MACARM64 platform
could not load traj.pkg from http://www.andrew.cmu.edu/user/bjones/traj/
r(601);とエラーが出る。解決方法は以下のリンクを参照。
要するに、intel mac環境で使えるようにStataを右クリックして「Rosettaを使用して開く」のチェックをオンにすれば良い。
基本的な文法
traj [if exp], var (varlist) indep(varlist) model (string) order (numlist)
[min (rela) max (rela)] iorder (numlist) risk (varlist) tcov (varlist) plottov (matrix)
start (matrix) weight (varname) exposure (varlist) refgroup (integer) dropout (numlist)
dcov (varlist) obsmar (varname) outcome (varname) omodel (string) detail]となる。ひとつひとつ解説していく。ちなみに、
[if exp]とはもしifを用いるならここで挿入しろという意味である(if age >= 30、など)。
Trajectory Variables
var (varlist)こではアウトカムで、異なる年齢（あるいはタイムポイント）で測定されている。
indep (varlist)暴露因子、説明変数である。
Model
model (string)アウトカムに対する分布を指定する。指定できるのは
cnorm, zip, logitの\(3\)種類である。
order
order (numlist)それぞれのグループの軌道に対する多項式のタイプ（\(0\)は切片モデル、\(1\)は線形モデル、\(2\)は二次式、\(3\)は三次式）である。ここで
order (2 2 3)などとすると、一つ目のグループは\(2\)次、二つ目のグループも\(2\)次、三つ目のグループは\(3\)次式でモデルを推定する。
その他のオプション
min (real)
max (real)打ち切り正規分布モデルを用いるときの、最小値あるいは最大値。”cnorm”を指定するときは必要。
iorder (numlist)これは各グループについてゼロ過剰モデルを適合するときに用いるオプションで、orderと同じく（\(0\)は切片モデル、\(1\)は線形モデル、\(2\)は二次式、\(3\)は三次式）となる。
exposure (varlist)ゼロ過剰ポアソンモデルを用いるときのオプション。
weight (varname)確率重み付けの変数。
Time-Stable Covariates for Group Membership
risk (varlist)グループメンバーの共変量。
refgroup (integer)リファレンスグループ（デフォルトで\(1\)になっている）。リスクオプションが用いられたときに指定する。
Time-varying Covariates for Group Membership
tcov (varlist)グループ軌道の時間によって変化する共変量。
plotcov (matrix)軌道を時間によって変化する共変量とともにプロットするときのオプション。
Dropout model
dropout (numlist)waveごとのドロップアウト確率のロジスティックモデルである。\(0\)は一定、\(1\)は前期に依存、\(2\)は\(2\)期前に依存、というように指定する。
dcov (varlist)ドロップアウトモデルを採用するとき、時間依存性共変量のラグを指定する。
obsmar (varname)どの観測がドロップアウトモデルに含まれ、missing at randomとして扱われるかを指定する二値変数。\(1\)を指定すると観測はMARとして扱われ、\(0\)を指定するとドロップアウトされる観測として扱われる。Distal Outcome Modelは省略。
プロットについて
基本的な構文はこちら。
trajplot, [xtitle (string) ytitle (string) model (integer) ci]xtitleは\(x\)軸の名称。
ytitleは\(y\)軸の名称。
modelははジョイント軌道モデルのうち、どちらを指定するか（\(1\)か\(2\)で、デフォルトは\(1\)）
ciは\(95%\)信頼区間。
実例
データセットは横データwide形式にしておく。やり方は簡単で、アウトカムが計測されたタイミングの変数を予め生成しておく。（以下リンクも参照）。
foreach x in 3 6 12 18 {
gen time_`x' = `x'
}
こうして、
traj, traj, var(outcome3 outcome6 outcome12 outcome18) indep(t_3 t_6 t_12 t_18) model(logit) order(3 3 2)などというようにモデルと作ればよい。結果ではBICとAICが出力される。グループ数やorderの次数などを色々試してみて、BIC/AICの低いモデルを採用すればよい。モデルを回転させると、”_traj_Group”, “_traj_ProbG1”, “_traj_ProbG2”, \(\cdots\)というグループ変数が生成される。これはある個人がGroup_xに属する確率を表している。”_traj_Group”がどのグループに属するかの変数なので、これを用いてグループごとの解析を行えば良い。
コメント