私は、現在のプログラムのどれも管理できない問題を最終的に解決できるようにすることを目的として、Cox 比例ハザード モデルを実行する小さなプログラムを作成しようとしています。の基本的な機能を再現しようとしましたが、最初のハードルで立ち往生しましたstcox。を使用stcoxして適切に機能していることを確認すると、 を使用して最大値が同じ場所にあるため、可能性が正しいことを示すことができますml plot。しかし、この時点から、それは崩壊します。変数drugは値 1 または 0 を取り、これは Cox モデルであるため適合する定数がないため、これは方程式で指定されます。ただし、Xbを観察することによって対数尤度が計算されるたびに、対照群と治療群の両方の値が変化していることがわかります。実際、それらの差は固定されており、定数だけが変化しています。
clear all
webuse drugtr
*Show st settings
stset
*Fit Cox proportional hazards model
stcox drug, breslow
gen timem = -_t
sort timem
capture program drop mlcox
program mlcox
version 10
args lnf Xb
tempvar risk sumrisk sumrisk2
qui gen double `risk'=exp(`Xb')
qui gen double `sumrisk'=0
local N=_N
forval j=1/`N' {
qui replace `sumrisk'=`sumrisk' + exp(`Xb'[`j']) if _t[`j']>=_t
}
tab `Xb'
qui replace `lnf' = 0 if $ML_y1==0
qui replace `lnf' = `Xb' - log(`sumrisk') if $ML_y1==1
end
ml model lf mlcox ( died = drug, noconstant)
ml plot -4 1 30
ml maximize
編集: 最尤のコーディングを修正しましたが、まだ収束に問題があります。