clear all // prog mata to compute the interval each duration belongs to cap program drop Inter mata: void Inter(string scalar t) { // import variables and matrices into mata st_view(vt,.,t) Tcut = st_matrix("Tcut") N = rows(vt) // interval to which t belongs (consider all intervals strictly included in [0,t] and add 1) Tcut_sup = Tcut[.,2..cols(Tcut)] nint = rowsum((vt :> (J(N,1,1) * Tcut_sup))) + J(N,1,1) // compute t minus the lower bound of the interval t belongs to Tcut_inf = Tcut[.,1..cols(Tcut) - 1] Ikt = (vt :> (J(N,1,1) * Tcut_inf)) :* ((vt :<= (J(N,1,1) * Tcut_sup))) Ik_length = Tcut_sup-Tcut_inf,0 lint = rowsum((vt * J(1,cols(Tcut) - 1,1) - J(N,1,1) * Tcut_inf) :* Ikt) // send results to stata st_addvar(("double", "double"), ("nint", "lint")) st_store(.,"nint",nint) st_store(.,"lint",lint) } end // prog that declares model and maximizes likelihood cap program drop maxlik program define maxlik syntax [, NV(real 1) MINIT] // declare model - parameters are ordered as follows: ba ea bp bb ga gp gb bz da dp db by I_a I_p I_b I_z I_y vp vz vy pzy local K_ALL = 0 foreach vv of varlist A1 P B1 Z Y{ local K_Tcut = colsof(Tcut_`vv')-2 foreach kk of numlist 1(1)`K_Tcut'{ local dkk = 12+`K_ALL'+`kk' local deltaIk_`vv' = "`deltaIk_`vv'' (delta`dkk':)" } local K_ALL = `K_ALL'+`K_Tcut' } scalar nv = `nv' local nv_all = (4*`nv'-1)*(`nv'>1) + 3*(`nv'==1) foreach hh of numlist 1(1)`nv_all'{ local dhh = 12+`K_ALL'+`hh' local deltahet = "`deltahet' (delta`dhh':)" } ml model lf lik (delta1: X*, nocons) (delta2:) (delta3: X*, nocons) (delta4: X*, nocons) (delta5:) (delta6:) (delta7:) /// (delta8: X*, nocons) (delta9:) (delta10:) (delta11:) (delta12: X*, nocons) /// `deltaIk_A1' `deltaIk_P' `deltaIk_B1' `deltaIk_Z' `deltaIk_Y' `deltahet', technique(bhhh 25 bfgs) // starting values if "`minit'"~=""{ ml init M, copy } else{ ml search } // maximize likelihood ml max end // prog to compute hazard rates and add the different bits of the likelihood //clear mata mata: function LL(string bvars, string dvars, string inta1, string intp, string intb1, string intz, string inty) { me = st_matrix("effects") st_view(NINT = ., ., tokens(inta1)) PCUT = st_matrix("Pcut_A1") SPCUT = st_matrix("SPcut_A1") LINT = NINT[.,2] NINT = NINT[.,1] H = PCUT[NINT,1] HINT = SPCUT[NINT:-1+(NINT:==1),1]:*(NINT:>1) + PCUT[NINT,1]:*LINT st_view(NINT = ., ., tokens(intp)) PCUT = st_matrix("Pcut_P") SPCUT = st_matrix("SPcut_P") m = cols(NINT)/2 LINT = NINT[.,m+1..2*m] NINT = NINT[.,1..m] H = H, PCUT[NINT[.,m],1] NINT = colshape(NINT',1) LINT = colshape(LINT',1) T = exp((0, me[1,1])) * ( I(m) + ( J(1,m,0)\(-I(m-1),J(m-1,1,0)) ) ) HINT = HINT, ( rowshape(SPCUT[NINT:-1+(NINT:==1),1]:*(NINT:>1) + PCUT[NINT,1]:*LINT, m)' * T' ) st_view(NINT = ., ., tokens(intb1)) PCUT = st_matrix("Pcut_B1") SPCUT = st_matrix("SPcut_B1") LINT = NINT[.,2] NINT = NINT[.,1] H = H, PCUT[NINT,1] HINT = HINT, SPCUT[NINT:-1+(NINT:==1),1]:*(NINT:>1) + PCUT[NINT,1]:*LINT st_view(NINT = ., ., tokens(intz)) PCUT = st_matrix("Pcut_Z") SPCUT = st_matrix("SPcut_Z") m = cols(NINT)/2 LINT = NINT[.,m+1..2*m] NINT = NINT[.,1..m] H = H, PCUT[NINT[.,m],1] NINT = colshape(NINT',1) LINT = colshape(LINT',1) T = exp((0, me[1,(2, 3, 4)])) * ( I(m) + ( J(1,m,0)\(-I(m-1),J(m-1,1,0)) ) ) HINT = HINT, ( rowshape(SPCUT[NINT:-1+(NINT:==1),1]:*(NINT:>1) + PCUT[NINT,1]:*LINT, m)' * T' ) st_view(NINT = ., ., tokens(inty)) PCUT = st_matrix("Pcut_Y") SPCUT = st_matrix("SPcut_Y") m = cols(NINT)/2 LINT = NINT[.,m+1..2*m] NINT = NINT[.,1..m] H = H, PCUT[NINT[.,m],1] NINT = colshape(NINT',1) LINT = colshape(LINT',1) T = exp((0, me[1,(5, 6, 7)])) * ( I(m) + ( J(1,m,0)\(-I(m-1),J(m-1,1,0)) ) ) HINT = HINT, ( rowshape(SPCUT[NINT:-1+(NINT:==1),1]:*(NINT:>1) + PCUT[NINT,1]:*LINT, m)' * T' ) st_view(BETA = ., ., tokens(bvars)) st_view(D = ., ., tokens(dvars)) BETA = BETA D = D mp = st_matrix("probas") mv = st_matrix("vhet") R = cols(mp) N = rows(D) L = rowshape( rowsum( (mv#J(N,1,1)) :* (J(R,1,1)#D) - (J(R,1,1)#HINT) :* exp(J(R,1,1)#BETA + mv#J(N,1,1)) ) , R)' L = log(rowsum(exp(mp:+L))) + rowsum(D:*(log(H)+BETA)) :- log(rowsum(exp(mp))) L = L + D[.,2] :* me[1,1]:*D[.,1] /// + D[.,4] :* (me[1,2]:*D[.,1]:*(1:-D[.,2]) + me[1,3]:*D[.,2]:*(1:-D[.,3]) + me[1,4]:*D[.,3]) /// + D[.,5] :* (me[1,5]:*D[.,1]:*(1:-D[.,2]) + me[1,6]:*D[.,2]:*(1:-D[.,3]) + me[1,7]:*D[.,3]) st_store(.,"lnvrais",L) } end // likelihood cap program drop lik program define lik local nv = nv foreach vv in "A1" "P" "B1" "Z" "Y"{ local K_`vv' = colsof(Tcut_`vv') - 1 local I_`vv' = "I2_`vv'" foreach kk of numlist 3(1)`K_`vv''{ local I_`vv' = "`I_`vv'' I`kk'_`vv'" } } foreach bb of numlist 1(1)`nv'{ local vphet = "`vphet' vp`bb'" local vzhet = "`vzhet' vz`bb'" local vyhet = "`vyhet' vy`bb'" if `bb' ~= `nv'{ local pzyhet = "`pzyhet' pzy`bb'" } } local pzy`nv' = 0 args lnf ba ea bp bb ga gp gb bz da dp db by `I_A1' `I_P' `I_B1' `I_Z' `I_Y' `vphet' `vzhet' `vyhet' `pzyhet' foreach vv of varlist A1 P B Z Y{ matrix Pcut_`vv' = .0001 matrix SPcut_`vv' = (.0001 * LTcut_`vv'[1,1]) \ J(`K_`vv''-1,1,0) foreach kk of numlist 2(1)`K_`vv''{ matrix Pcut_`vv' = Pcut_`vv' \ (exp(-`I`kk'_`vv''[1]) / (1+exp(-`I`kk'_`vv''[1]))) matrix SPcut_`vv'[`kk',1] = SPcut_`vv'[`kk'-1,1] + Pcut_`vv'[`kk',1]*LTcut_`vv'[1,`kk'] } } local int_P = "nint_PtA1 nint_P lint_PtA1 lint_P" local int_Z = "nint_ZtA1 nint_ZtP nint_ZtB1 nint_Z lint_ZtA1 lint_ZtP lint_ZtB1 lint_Z" local int_Y = "nint_YtA1 nint_YtP nint_YtB1 nint_Y lint_YtA1 lint_YtP lint_YtB1 lint_Y" matrix effects = (`ea'[1],`ga'[1],`gp'[1],`gb'[1],`da'[1],`dp'[1],`db'[1]) matrix probas = `pzy`nv'' matrix vhet = (J(1,3,`vp`nv''[1]), `vz`nv''[1], `vy`nv''[1]) if `nv'>1{ local nv1 = `nv'-1 foreach rr of numlist `nv1'(1)1{ matrix probas = `pzy`rr''[1], probas matrix vhet = (J(1,3,`vp`rr''[1]), `vz`rr''[1], `vy`rr''[1]) \ vhet } } gen double lnvrais = 0 mata: LL("`ba' `bp' `bb' `bz' `by'", "A1 P B1 Z Y", "nint_A1 lint_A1", "`int_P'", "nint_B1 lint_B1", "`int_Z'", "`int_Y'") qui replace `lnf' = lnvrais drop lnvrais end use table_75 drop *A2 *A3 *A4 *A5 *B2 *B3 *B4 *B5 // create intervals for piecewise constant hazards foreach vv of varlist P A1 B1 Z Y{ su t`vv' if `vv'==1, det _pctile t`vv' if `vv'==1, nq(11) matrix Tcut_`vv' = (0, r(r1), r(r2), r(r3), r(r4), r(r5), r(r6), r(r7), r(r8), r(r9), r(r10),.) } // length of each interval (used later to compute integrated hazard rates) qui foreach vv of varlist A1 P B1 Z Y{ local K_`vv' = colsof(Tcut_`vv')-2 matrix LTcut_`vv' = Tcut_`vv'[1,2] foreach cc of numlist 2(1)`K_`vv''{ matrix LTcut_`vv' = LTcut_`vv', Tcut_`vv'[1,`cc'+1]-Tcut_`vv'[1,`cc'] } noi matrix list Tcut_`vv' noi matrix list LTcut_`vv' } // tZ (resp. tY) is censored by tY (resp. tZ) replace Y = 0 if tZ < tY replace tY = tZ if tZ < tY // find which interval each duration belongs to foreach vv of varlist A1 P B1 Z Y{ matrix Tcut = Tcut_`vv' mata: Inter("t`vv'") rename nint nint_`vv' rename lint lint_`vv' } foreach vv of varlist P Z Y{ matrix Tcut = Tcut_`vv' mata: Inter("tA1") rename nint nint_`vv'tA1 rename lint lint_`vv'tA1 } foreach vv of varlist Z Y{ matrix Tcut = Tcut_`vv' mata: Inter("tP") rename nint nint_`vv'tP rename lint lint_`vv'tP } foreach vv of varlist Z Y{ matrix Tcut = Tcut_`vv' gen double tPB = tP+tB1*P mata: Inter("tPB") rename nint nint_`vv'tB1 rename lint lint_`vv'tB1 drop tPB } matrix M = ( .0601013, -.6818966, .5763798, -.0219434, .0352129, .0466404, -.0153628, .1542297, -.0490291, -.0709146, -.0589766, -.1894159, .0097903, /// -.0007275, -.0139814, -.0464572, -.051971, -.0087275, 2.124328, -.1050517, /// -.1284812, /// -.4705276, .8425404, -.8578912, -.2810217, .1874907, -.7291221, -.1494358, -.110895, -.285446, -.185316, -.3541172, -.2885243, -.2593461, /// -.1783272, .2036339, .0892493, .1059123, .1599462, -15.60995, -1.427539, /// -.132381, .1893963, -.2045503, -.1263874, .0782094, -.4638715, -.0152003, -.1352692, -.3380791, -.3934296, -.3768191, -.427929, -.2219454, /// -.2094316, .0260245, -.083015, -.0781566, .0543719, -8.583396, -1.285607, /// .8959077, 2.644384, 2.831976, /// -.1355234, 1.501476, -1.545338, -.1777239, .2031262, .2989467, -.0550441, -.0103289, -.2737245, -.266387, -.1852578, -.1424156, .1174956, /// -.4710776, .0685089, -.1538684, .0199283, .1429505, -.2121854, .7531427, /// .4091654, -.0459475, .5142501, /// .0313308, -1.106682, .9860802, -.2523569, .1860936, -.1885734, -.0098875, .0071584, .0462545, .0501849, -.0378257, -.1699541, -.0561349, /// .1793531, .1111822, .1173911, .0203909, .0034032, 2.807206, -.6820767, /// 9.168471, 9.248176, 9.085564, 9.367116, 9.109801, 8.947102, 9.552071, 9.978703, 10.8058, 11.97151, /// 10.25289, 10.35687, 10.56689, 11.10702, 11.31853, 11.36468, 11.02592, 11.41242, 11.51648, 12.14401, /// 9.250672, 9.359858, 9.367563, 9.428748, 9.590561, 9.618017, 9.438458, 8.672611, 8.634726, 10.62414, /// 9.332054, 9.202862, 9.245157, 9.221721, 9.226328, 9.394481, 9.490944, 9.704647, 9.800531, 10.76038, /// 8.566902, 9.227665, 9.416571, 9.568667, 9.508635, 9.488841, 9.668784, 9.808923, 9.768417, 9.603993, /// 5.033718, .2963769, 3.359013 ) maxlik, nv(1) minit matrix eb1 = e(b) matrix M = eb1[1,1..colsof(eb1)-3], 3, 5.0331126, 0, .34459065, 3, 3.3575577, 0 maxlik, nv(2) minit matrix list e(b) matrix eb2 = e(b) matrix eV2 = e(V) // gP-gA matrix bb = eb2[1,63]-eb2[1,62] matrix VV = (eV2[62,62], eV2[62,63]) \ (eV2[63,62], eV2[63,63]) matrix sd2 = (-1,1)*VV*(-1,1)' matrix list bb matrix list sd2 // dP-dA matrix bb = eb2[1,86]-eb2[1,85] matrix VV = (eV2[85,85], eV2[85,86]) \ (eV2[86,85], eV2[86,86]) matrix sd2 = (-1,1)*VV*(-1,1)' matrix list bb matrix list sd2 // gB-gP - this effect are called gamma_B in the paper matrix bb = eb2[1,64]-eb2[1,63] matrix VV = (eV2[63,63], eV2[63,64]) \ (eV2[64,63], eV2[64,64]) matrix sd2 = (-1,1)*VV*(-1,1)' matrix list bb matrix list sd2 // dB-dP - this effect are called delta_B in the paper matrix bb = eb2[1,87]-eb2[1,86] matrix VV = (eV2[86,86], eV2[86,87]) \ (eV2[87,86], eV2[87,87]) matrix sd2 = (-1,1)*VV*(-1,1)' matrix list bb matrix list sd2