st: mata competing risks with discrete distribution of frailty

["Amparo Nagore Garcia" <amparo.nagore@uv.es>]

Dear statalist,

I am beginner programming in mata,and I am programming the evaluator for a likelyhood function of a competing risks model with three destination states, the baseline hazard follows exponential and piecewise distribution and the frailty follows a correlated discrete distribution with two points of support.
My data include multiple spells per individual and data are split in order to include time varying covariates. 
I am using chapter 11 of Maximum likelihood estimation with Stata (fouth edition from W Gould, J Pitblado, W Sribney) but as there is no examples of competing risks I have no clear how to program for each panel (corresponding to each individual) the hazard and survival function for each competing risks (from "for(i=1; i<= npanels;++i) in my code). Should I create a panel submatrix for each competing risks? At the bottom of this e-mail is my code.
Any help is really appreciated.
Thank you very much for your help.
Amparo

mata
//evaluator function
void CR_mata(transmorphic scalar M, real scalar todo,real rowvector b,real scalar lnf, real rowvector g, real matrix H)
{
//dependent variables for three competing risks
dt0_j=moptimize_util_depvar(M,1)
dt1_j=moptimize_util_depvar(M,2)
dj=moptimize_util_depvar(M,3)

dt0_u=moptimize_util_depvar(M,4)
dt1_u=moptimize_util_depvar(M,5)
du=moptimize_util_depvar(M,6)

dt0_n=moptimize_util_depvar(M,7)
dt1_n=moptimize_util_depvar(M,8)
dn=moptimize_util_depvar(M,9)

//parameters of three competing risks
xbj=moptimize_util_xb(M,b,1)
xbu=moptimize_util_xb(M,b,2)
xbn=moptimize_util_xb(M,b,3)

//scalar parameters (unobserved heterogeneity following a discrete distribution with two mass points)
a1=moptimize_util_eq_indices(M,1)
Vj1=moptimize_util_eq_indices(M,2)
Vu1=moptimize_util_eq_indices(M,3)
Vn1=moptimize_util_eq_indices(M,4)

//Terms of mass points:
// p1, p2 y Vj2 Vu2 Vn2
real scalar p1
p1=exp(a1)/(1+exp(a1))
real scalar p2
p2=1/(1+exp(a1))

real scalar Vj2
Vj2=-(p1*Vj1)/(1-p1)

real scalar Vu2
Vu2=-(p1*Vu1)/(1-p1)

real scalar Vn2
Vn2=-(p1*Vn1)/(1-p1)

//panvar is a matrix including id and the value of all explanatory variables
st_view (panvar=.,.,tokens("`panel'"))
paninfo=panelsetup(panvar,1)
npanels=panelstats(paninfo)[1]

//X is a matrix including all the explanatory variables
st_view (X=.,.,tokens("`xlist'"))
cols(X)
rows(X)

//lnfj is a matrix with the number of rows equal to the number of ids and one column 
lnfj=J(npanels,1,0)

// Since it is a multi-spell panel and the panel is split for including time-varying covariates and the term of frailty is common for
// all the spells of the same individual, we need to sum for each individual and for each competing risks xb
//So that I have created zi, a matrix with all the explanatory variables for each i.

for (i=1; i<=npanels; ++i)
zi=panelsubmatrix(X,i,paninfo)

sumbj1=sum((xbj:+Vj1):*dj)
sumbj2=sum((xbj:+Vj2):*dj) 
sumbu1=sum((xbu:+Vu1):*du)
sumbu2=sum((xbu:+Vu2):*du) 
sumbn1=sum((xbn:+Vn1):*dn)
sumbn2=sum((xbn:+Vn2):*dn) 

//To construct the survival function of each id
inthazj1=sum(exp(xbj:+Vj1')*(dt1_j:-dt0_j))
inthazj2=sum(exp(xbj:+Vj2')*(dt1_j:-dt0_j))
inthazu1=sum(exp(xbu:+Vu1')*(dt1_u:-dt0_u))
inthazu2=sum(exp(xbu:+Vu2')*(dt1_u:-dt0_u))
inthazn1=sum(exp(xbn:+Vn1')*(dt1_n:-dt0_n))
inthazn2=sum(exp(xbn:+Vn2')*(dt1_n:-dt0_n))

//likelihood for each id
lnfj[i]=p1*(exp(sumbj1)*exp(sumbu1)*exp(sumbn1)*exp(-inthazj1)*exp(-inthazu1)*exp(-inthazn1)) +
p2*(exp(sumbj2)*exp(sumbu2)*exp(sumbn2)*exp(-inthazj1)*exp(-inthazu1)*exp(-inthazn1))

lnf=moptimize_util_sum(M,ln(lnfj))
}

st_view (d_j=.,.,"`d_j'")
cols(d_j)
rows(d_j)
st_view(t0_j=.,.,"`t0_j'")

st_view(t_j=.,.,"`t_j'")
cols(t_j)
rows(t_j)
st_view(d_u=.,.,"`d_u'")
st_view(t0_u=.,.,"`t0_u'")
st_view(t_u=.,.,"`t_u'")

st_view(d_n=.,.,"`d_n'")
st_view(t0_n=.,.,"`t0_n'")
st_view(t_n=.,.,"`t_n'")

end

ml model d0 CR_mata() (cr_05_j: t0_j t_j d_j t0_u t_u d_u t0_n t_n d_n = tp* $personal $labour_bis if id_job05==1, nocons) ///
(cr_05_u:tp* $personal $labour_bis if id_job05==1, nocons) ///
(cr_05_n:tp* $personal $labour_bis if id_job05==1, nocons) ///

ml maximize, trace difficult

     
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/