Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.
From | "Amparo Nagore Garcia" <amparo.nagore@uv.es> |
To | statalist@hsphsun2.harvard.edu |
Subject | st: mata competing risks with discrete distribution of frailty |
Date | Wed, 19 Mar 2014 14:45:57 +0100 |
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/