ループを使って繰り返しのコマンドをひとまとめにするコツを以前の記事でまとめていますが、今回は、ループをいくつも重ねてすっきり洗練されたプログラムを書くコツをまとめてみたいと思います．
1．多変量モデルをいくつも走らせて、メインの曝露因子の結果をまとめて表示させる
これまでの知識を総動員しながらプログラムを組んでいきましょう．生存時間分析を実施しますが、メインの曝露因子は4レベルの変数で、4つめを参照カテゴリーとします．aval_evはeventが生じるまでの時間、cnsr_evはイベントコードとします．後からでてきますが、複数のeventを検討するとき、イベントの種類をheart attack, stroke, deathの3つとします．Covariatesのセットはvar1 ~ var5までの組み合わせで、unadjusted, model 1, model 2, model 3を作るものとします．統計学の計算式はこちらを参照ください．
*** 多変量モデルをマクロに格納する ⇒ ここではglobal macroを使います．
global cov_0
global cov_1 var1 var2 var3
global cov_2 $cov_1 var4
global cov_3 $cov_2 var5
qui {
stset aval_ev, failure(cnsr_ev == 0) /* アウトカムをcnsr_evが0のときと定義 */
noi disp "Exposure", "Model", "category", "Coef", "Lowlim", "Uplim", "P"
forvalues n = 1/3 {
/* 曝露因子exposureが4レベルのカテゴリ変数とし、リファレンスカテゴリーを4とした時 */
stcox ib4.exposure ${cov_`n'} , strata(sitecd)
/* Cox比例ハザードモデルを、4番目のレベルをリファレンスとして実行 */
local coef_1 = _b[1.exposure]
local low_1 = _b[1.exposure] - invnormal(0.975)*(_se[1.exposure])
/* 信頼区間の計算は別ページを参照ください！ */
local high_1 = _b[1.exposure] + invnormal(0.975)*(_se[1.exposure])
local p_1 = 2*normprob(-abs(_b[1.exposure]/_se[1.exposure]))
noi disp "`v'", "`n'", 1, %3.2f exp(`coef_1'), %3.2f exp(`low_1'), %3.2f exp(`high_1'), %4.3f `p_1'
/* displayするときにコンマで区切るとスペースが空く */
/* 連続変数をdisplayするとき、頭に桁を指定すると小数点の位置を指定できる */
/* 以下は別のカテゴリ */
local coef_2 = _b[2.exposure]
local low_2 = _b[2.exposure] - invnormal(0.975)*(_se[2.exposure])
/* 信頼区間の計算は別ページを参照ください！ */
local high_2 = _b[2.exposure] + invnormal(0.975)*(_se[2.exposure])
local p_2 = 2*normprob(-abs(_b[2.exposure]/_se[2.exposure]))
noi disp "`v'", "`n'", 2, %3.2f exp(`coef_2'), %3.2f exp(`low_2'), %3.2f exp(`high_2'), %4.3f `p_2'
local coef_3 = _b[3.exposure]
local low_3 = _b[3.exposure] - invnormal(0.975)*(_se[3.exposure])
/* 信頼区間の計算は別ページを参照ください！ */
local high_3 = _b[3.exposure] + invnormal(0.975)*(_se[3.exposure])
local p_3 = 2*normprob(-abs(_b[3.exposure]/_se[3.exposure]))
noi disp "`v'", "`n'", 3, %3.2f exp(`coef_3'), %3.2f exp(`low_3'), %3.2f exp(`high_3'), %4.3f `p_3'
}
} ここで、複数のカテゴリーの結果を一つずつ書いているのがうっとうしいなあと気がついていただけたでしょうか？それが今回の主題の「ループを重ねる」につながるメンタリティです．もっと楽できないかな～と考えることこそが美しいプログラムを作る強力なドライブになると思っています（笑）．
コマンドを外に外に広げるのがコツです．今回はたまたま中身にループを組み込む形になりますが、繰り返しが見えたらその共通部分と違いを明確にするとプログラムに書き下ろしやすくなります．完成形の表を思い浮かべながら、どんなアウトプットにしたいのかを明確にすると間違えずに済むでしょう．
最初のステップはこのオレンジ色の部分でしたが、ここに青の繰り返しコマンドを組み込むと以下のようになります．
*** 繰り返しのコマンドを中に組み込みます
qui {
stset aval_ev, failure(cnsr_ev == 0) /* アウトカムをcnsr_evが0のときと定義 */
noi disp "Exposure", "Model", "category", "Coef", "Lowlim", "Uplim", "P"
forvalues n = 1/3 {
/* 曝露因子exposureが4レベルのカテゴリ変数とし、リファレンスカテゴリーを4とした時 */
forvalues m = 1/4 {
stcox ib4.exposure ${cov_`n'} , strata(sitecd)
/* Cox比例ハザードモデルを、4番目のレベルをリファレンスとして実行 */
local coef_`m' = _b[`m'.exposure]
local low_`m' = _b[`m'.exposure] - invnormal(0.975)*(_se[`m'.exposure])
local high_`m' = _b[`m'.exposure] + invnormal(0.975)*(_se[`m'.exposure])
local p_`m' = 2*normprob(-abs(_b[`m'.exposure]/_se[`m'.exposure]))
noi disp "`v'", "`n'", "`m'", %3.2f exp(`coef_`m''), %3.2f exp(`low_`m''), %3.2f exp(`high_`m''), %4.3f `p_`m''
/* displayするときにコンマで区切るとスペースが空く */
/* 連続変数をdisplayするとき、頭に桁を指定すると小数点の位置を指定できる */
}
}
} まとめるとすっきりしますね.
2．メインの曝露因子それぞれの結果を表示させる
さて、メインの曝露因子をそれぞれに表示させたい場合は以下の様にコマンドを一般化していくようにループコマンドを使いましょう．ここでは曝露因子は exposure1 と exposure2 の二つとします．シンプルに考えるため、カテゴリ数が同じとしてプログラムを書いていきます．それぞれの曝露因子につき別々に多変量解析の結果を表示したいので、この場合は曝露因子をより上位に置きます.
*** for many exposures
qui {
stset aval_ev, failure(cnsr_ev == 0) /* Set the outcome */
*reference=4
noi disp "Exposure", "Model", "category", "Coef", "Lowlim", "Uplim", "P"
foreach v of varlist exposure1 exposure2 {
forvalues n = 1/3 {
forvalues m = 1/4 {
stcox ib4.`v' ${cov_`n'} , strata(sitecd)
local coef_`m' = _b[`m'.`v']
local low_`m' = _b[`m'.`v'] - invnormal(0.975)*(_se[`m'.`v'])
local high_`m' = _b[`m'.`v'] + invnormal(0.975)*(_se[`m'.`v'])
local p_`m' = 2*normprob(-abs(_b[`m'.`v']/_se[`m'.`v']))
noi disp "`v'", "`n'", "`m'", %3.2f exp(`coef_`m''), %3.2f exp(`low_`m''), %3.2f exp(`high_`m''), %4.3f `p_`m''
}
}
}
} 3．別の種類のアウトカムそれぞれの結果を表示させる
ここで、複数のアウトカムについて一気に結果を出力したい、となったときの手順をご説明します．アウトカムを「ev+数字」という形のlocal macroにしてしまう、というのが味噌です．
それで、これまでのプログラムの外側にこれを追加する、という感じです．
*** 複数のアウトカムを設定
local ev1 "heartattack"
local ev2 "stroke"
local ev3 "death"
qui {
forvalues l = 1/3 {
noi disp "Outcome", "`ev`l''"
stset aval_`ev`l'' , failure(cnsr_`ev`l'' == 0) /* Set the outcome */
*reference=4
noi disp "Exposure", "Model", "category", "Coef", "Lowlim", "Uplim", "P"
foreach v of varlist exposure1 exposure2 {
forvalues n = 1/3 {
forvalues m = 1/4 {
stcox ib4.`v' ${cov_`n'} , strata(sitecd)
local coef_`m' = _b[`m'.`v']
local low_`m' = _b[`m'.`v'] - invnormal(0.975)*(_se[`m'.`v'])
local high_`m' = _b[`m'.`v'] + invnormal(0.975)*(_se[`m'.`v'])
local p_`m' = 2*normprob(-abs(_b[`m'.`v']/_se[`m'.`v']))
noi disp "`v'", "`n'", "`m'", %3.2f exp(`coef_`m''), %3.2f exp(`low_`m''), %3.2f exp(`high_`m''), %4.3f `p_`m''
}
}
}
}
} 4．層別解析結果を表示させる
モデルよりも上位に層別解析を置きたいので、以下のようなプログラムになります.
*** 層別解析を追加
local ev1 "heartattack"
local ev2 "stroke"
local ev3 "death"
qui {
forvalues l = 1/3 {
noi disp "Outcome", "`ev`l''"
stset aval_`ev`l'' , failure(cnsr_`ev`l'' == 0) /* Set the outcome */
*reference=4
noi disp "Strata", "Exposure", "Model", "category", "Coef", "Lowlim", "Uplim", "P"
foreach v of varlist exposure1 exposure2 {
foreach s of vrlist strata1 strata2 strata3 {
noi disp "`s'"
forvalues t = 0/1 {
sum `s' if `s' == `t'
noi disp `t', r(N)
forvalues n = 1/3 {
forvalues m = 1/4 {
stcox ib4.`v' ${cov_`n'} if `s' == `t' , strata(sitecd)
local coef_`m' = _b[`m'.`v']
local low_`m' = _b[`m'.`v'] - invnormal(0.975)*(_se[`m'.`v'])
local high_`m' = _b[`m'.`v'] + invnormal(0.975)*(_se[`m'.`v'])
local p_`m' = 2*normprob(-abs(_b[`m'.`v']/_se[`m'.`v']))
noi disp "`v'", "`n'", "`m'", %3.2f exp(`coef_`m''), %3.2f exp(`low_`m''), %3.2f exp(`high_`m''), %4.3f `p_`m''
stcox ib4.`v'##`s' ${cov_3} , strata(sitecd)
test 1.`v'#1.`s' 2.`v'#1.`s' 3.`v'#1.`s'
noi disp "P interaction" %4.3f r(p)
}
}
}
}
}
}
} 5．まとめ
複数のループを組み合わせたいときはより上位の概念が一番外にあり、順番に中に格納していきます.
注意すべきは、同じ記号のマクロを使わないようにすることです.
外側だけで有効な場合、一番内側まで有効な場合、などによってエラーがでてきますので、試行錯誤しながら微調整をしていきましょう.
Stataのプログラミングを深く学びたい方はコチラ．
コメント