*** Last modified 11th October 2005 program define intcens_ll version 8.2 local t0 "$ML_y1" local t1 "$ML_y2" local f "$ML_y3" quietly{ if "$intcens_maxm"=="lf"{ if "$intcens_dist"=="logn"{ args lnf xb logsigma tempvar sigma gen double `sigma'=exp(`logsigma') if $ML_samp replace `lnf'=log(normden(log(`t0'), `xb', `sigma')) if `f'==1&$ML_samp replace `lnf'=log(norm(-(log(`t0')-`xb')/`sigma')) if `f'==2&$ML_samp replace `lnf'=log(norm((log(`t1')-`xb')/`sigma')) if `f'==3&$ML_samp replace `lnf'=log(norm((log(`t1')-`xb')/`sigma')-norm((log(`t0')-`xb')/`sigma')) if `f'==4&$ML_samp } else if "$intcens_dist"=="logl"{ args lnf xb loggamma tempvar gamma lambda gen double `gamma'=exp(`loggamma') if $ML_samp gen double `lambda'=exp(-`xb'/`gamma') if $ML_samp replace `lnf'=-`xb'/`gamma'+(1/`gamma')*log(`t0')-`loggamma'-2*log(1+`lambda'*(`t0')^(1/`gamma')) if `f'==1&$ML_samp replace `lnf'=-log(1+`lambda'*(`t0')^(1/`gamma')) if `f'==2&$ML_samp replace `lnf'=-log(1+(`lambda')^(-1)*(`t1')^(-1/`gamma')) if `f'==3&$ML_samp replace `lnf'=log(1/(1+`lambda'*(`t0')^(1/`gamma'))-1/(1+`lambda'*(`t1')^(1/`gamma'))) if `f'==4&$ML_samp } else if "$intcens_dist"=="exp"{ args lnf xb tempvar lambda gen double `lambda'=exp(-`xb') if $ML_samp replace `lnf'=-`xb'-`lambda'*`t0'+log(`t0') if `f'==1&$ML_samp replace `lnf'=-`lambda'*`t0' if `f'==2&$ML_samp replace `lnf'=log(1-exp(-`lambda'*`t1')) if `f'==3&$ML_samp replace `lnf'=log(exp(-`lambda'*`t0')-exp(-`lambda'*`t1')) if `f'==4&$ML_samp } else if "$intcens_dist"=="weib"{ args lnf xb logp tempvar lambda p gen double `p'=exp(`logp') if $ML_samp gen double `lambda'=exp(-`xb'*`p') if $ML_samp replace `lnf'=-`xb'*`p'+`logp'+`p'*log(`t0')-`lambda'*(`t0')^`p' if `f'==1&$ML_samp replace `lnf'=-`lambda'*(`t0')^`p' if `f'==2&$ML_samp replace `lnf'=log(1-exp(-`lambda'*(`t1')^`p')) if `f'==3&$ML_samp replace `lnf'=log(exp(-`lambda'*(`t0')^`p')-exp(-`lambda'*(`t1')^`p')) if `f'==4&$ML_samp } else if "$intcens_dist"=="gam"{ args lnf xb logalpha tempvar lambda alpha gen double `alpha'=exp(`logalpha') if $ML_samp gen double `lambda'=exp(-`xb') if $ML_samp replace `lnf'=-`xb'*`alpha'+`alpha'*log(`t0')-`lambda'*`t0'-lngamma(`alpha') if `f'==1&$ML_samp replace `lnf'=log(1-gammap(`alpha', `lambda'*`t0')) if `f'==2&$ML_samp replace `lnf'=log(gammap(`alpha', `lambda'*`t1')) if `f'==3&$ML_samp replace `lnf'=log(gammap(`alpha', `lambda'*`t1')-gammap(`alpha', `lambda'*`t0')) if `f'==4&$ML_samp } else if "$intcens_dist"=="gen"{ args lnf xb logsigma kappa1 tempvar sigma gamma firstml tempname kappa gen byte `firstml'=0 replace `firstml'=sum($ML_samp) if `firstml'[_n-1]==0|_n==1 su `kappa1' if `firstml'==1, meanonly scalar `kappa'=r(min) gen double `sigma'=exp(`logsigma') if $ML_samp gen double `gamma'=1/(`kappa')^2 if $ML_samp if abs(`kappa')<=1E-6{ local sign=1 } else{ local sign=sign(`kappa') } forval i=0/1{ tempvar z`i' u`i' gen double `z`i''=`sign'*(log(`t`i'')-`xb')/`sigma' if $ML_samp gen double `u`i''=`gamma'*exp(abs(`kappa')*`z`i'') if $ML_samp } if abs(`kappa')>1E-6{ replace `lnf'=-2*(`gamma'-0.5)*log(abs(`kappa'))+`z0'/abs(`kappa')-`u0'-`logsigma'-lngamma(`gamma') if `f'==1&$ML_samp if `kappa'>1E-6{ replace `lnf'=log(1-gammap(`gamma', `u0')) if `f'==2&$ML_samp replace `lnf'=log(gammap(`gamma', `u1')) if `f'==3&$ML_samp replace `lnf'=log(gammap(`gamma',`u1')-gammap(`gamma',`u0')) if `f'==4&$ML_samp } else if `kappa'<-1E-6{ replace `lnf'=log(gammap(`gamma', `u0')) if `f'==2&$ML_samp replace `lnf'=log(1-gammap(`gamma', `u1')) if `f'==3&$ML_samp replace `lnf'=log(gammap(`gamma',`u0')-gammap(`gamma',`u1')) if `f'==4&$ML_samp } } else if abs(`kappa')<=1E-6{ replace `lnf'=log(normden(`z0'))-`logsigma' if `f'==1&$ML_samp replace `lnf'=log(norm(-`z0')) if `f'==2&$ML_samp replace `lnf'=log(norm(`z1')) if `f'==3&$ML_samp replace `lnf'=log(norm(`z1')-norm(`z0')) if `f'==4&$ML_samp } } else if "$intcens_dist"=="gomp"{ args lnf xb gamma1 tempvar lambda firstml gen double `lambda'=exp(`xb') if $ML_samp tempname gamma gen byte `firstml'=0 replace `firstml'=sum($ML_samp) if `firstml'[_n-1]==0|_n==1 su `gamma1' if `firstml'==1, meanonly scalar `gamma'=r(min) if abs(`gamma')<=1E-6{ replace `lnf'=`xb'-`lambda'*`t0'+log(`t0') if `f'==1&$ML_samp replace `lnf'=-`lambda'*`t0' if `f'==2&$ML_samp replace `lnf'=log(1-exp(-`lambda'*`t1')) if `f'==3&$ML_samp replace `lnf'=log(exp(-`lambda'*`t0')-exp(-`lambda'*`t1')) if `f'==4&$ML_samp } else{ replace `lnf'=`xb'+`gamma'*`t0'-(`lambda'/`gamma')*(exp(`gamma'*`t0')-1)+log(`t0') if `f'==1&$ML_samp replace `lnf'=-(`lambda'/`gamma')*(exp(`gamma'*`t0')-1) if `f'==2&$ML_samp replace `lnf'=log(1-exp(-(`lambda'/`gamma')*(exp(`gamma'*`t1')-1))) if `f'==3&$ML_samp replace `lnf'=log(exp(-(`lambda'/`gamma')*(exp(`gamma'*`t0')-1))-exp(-(`lambda'/`gamma')*(exp(`gamma'*`t1')-1))) if `f'==4&$ML_samp } } else if "$intcens_dist"=="invg"{ args lnf xb logphi tempvar eta phi gen double `eta'=exp(`xb') if $ML_samp gen double `phi'=exp(`logphi') if $ML_samp replace `lnf'=-(`logphi'+log(`t0'))/2-(`t0'-`eta')^2/(2*(`eta')^2*`t0'*`phi') if `f'==1&$ML_samp replace `lnf'=log(norm(sqrt(1/(`t0'*`phi'))*(1-`t0'/`eta')) /* */ -exp(2/(`eta'*`phi'))*norm(-sqrt(1/(`t0'*`phi'))*(`t0'/`eta'+1))) if `f'==2&$ML_samp replace `lnf'=log(1-norm(sqrt(1/(`t1'*`phi'))*(1-`t1'/`eta')) /* */ +exp(2/(`eta'*`phi'))*norm(-sqrt(1/(`t1'*`phi'))*(`t1'/`eta'+1))) if `f'==3&$ML_samp replace `lnf'=log(norm(sqrt(1/(`t0'*`phi'))*(1-`t0'/`eta'))- norm(sqrt(1/(`t1'*`phi'))*(1-`t1'/`eta')) /* */ +exp(2/(`eta'*`phi'))*(norm(-sqrt(1/(`t1'*`phi'))*(`t1'/`eta'+1))-norm(-sqrt(1/(`t0'*`phi'))*(`t0'/`eta'+1)))) if `f'==4&$ML_samp } else if "$intcens_dist"=="wien"{ args lnf xb logc tempvar mu c eta phi logphi gen double `mu'=exp(`xb') if $ML_samp gen double `c'=exp(`logc') if $ML_samp gen double `eta'=`c'/`mu' if $ML_samp gen double `phi'=1/(`c')^2 if $ML_samp gen double `logphi'=-2*`logc' if $ML_samp replace `lnf'=-(`logphi'+log(`t0'))/2-(`t0'-`eta')^2/(2*(`eta')^2*`t0'*`phi') if `f'==1&$ML_samp replace `lnf'=log(norm(sqrt(1/(`t0'*`phi'))*(1-`t0'/`eta')) /* */ -exp(2/(`eta'*`phi'))*norm(-sqrt(1/(`t0'*`phi'))*(`t0'/`eta'+1))) if `f'==2&$ML_samp replace `lnf'=log(1-norm(sqrt(1/(`t1'*`phi'))*(1-`t1'/`eta')) /* */ +exp(2/(`eta'*`phi'))*norm(-sqrt(1/(`t1'*`phi'))*(`t1'/`eta'+1))) if `f'==3&$ML_samp replace `lnf'=log(norm(sqrt(1/(`t0'*`phi'))*(1-`t0'/`eta'))- norm(sqrt(1/(`t1'*`phi'))*(1-`t1'/`eta'))+exp(2/(`eta'*`phi'))* /* */ (norm(-sqrt(1/(`t1'*`phi'))*(`t1'/`eta'+1))-norm(-sqrt(1/(`t0'*`phi'))*(`t0'/`eta'+1)))) if `f'==4&$ML_samp } else if "$intcens_dist"=="wienran"{ args lnf xb logc logtau tempvar mu c tau gen double `mu'=exp(`xb') if $ML_samp gen double `c'=exp(`logc') if $ML_samp gen double `tau'=exp(`logtau') if $ML_samp forval i=0/1{ tempvar u`i' gen double `u`i''=(`t`i''*`tau')^2+`t`i'' if $ML_samp } replace `lnf'=`logc'-0.5*log(`u0')-(`mu'*`t0'-`c')^2/(2*`u0') if `f'==1&$ML_samp replace `lnf'=log(norm((`c'-`mu'*`t0')/sqrt(`u0'))-exp(2*`c'*`mu'+2*(`c'*`tau')^2)* /* */ norm(-(`c'+2*`c'*`t0'*(`tau')^2+`mu'*`t0')/sqrt(`u0'))) if `f'==2&$ML_samp replace `lnf'=log(norm(-(`c'-`mu'*`t1')/sqrt(`u1'))+exp(2*`c'*`mu'+2*(`c'*`tau')^2)* /* */ norm(-(`c'+2*`c'*`t1'*(`tau')^2+`mu'*`t1')/sqrt(`u1'))) if `f'==3&$ML_samp replace `lnf'=log(norm((`c'-`mu'*`t0')/sqrt(`u0'))+exp(2*`c'*`mu'+2*(`c'*`tau')^2)*(norm(-(`c'+2*`c'*`t1'*(`tau')^2+`mu'*`t1')/sqrt(`u1')) /* */ -norm(-(`c'+2*`c'*`t0'*(`tau')^2+`mu'*`t0')/sqrt(`u0')))-norm((`c'-`mu'*`t1')/sqrt(`u1'))) if `f'==4&$ML_samp } exit } /* end of lf method */ else{ /* method d2, only written for 2-parameter gamma */ if "$intcens_dist"=="gam"{ args todo b lnf g negH g1 g2 tempvar ll theta1 lambda tempname theta2 alpha gen double `ll'=. mleval `theta1'=`b', eq(1) mleval `theta2'=`b', eq(2) scalar gen double `lambda'=exp(-`theta1') if $ML_samp scalar `alpha'=exp(`theta2') forval i=0/1{ tempvar P`i' lt`i' gen double `lt`i''=`lambda'*`t`i'' if $ML_samp gen double `P`i''=gammap(`alpha', `lt`i'') if $ML_samp } replace `ll'=-`theta1'*`alpha'+`alpha'*log(`t0')-`lt0'-lngamma(`alpha') if `f'==1&$ML_samp replace `ll'=log(1-`P0') if `f'==2&$ML_samp replace `ll'=log(`P1') if `f'==3&$ML_samp replace `ll'=log(`P1'-`P0') if `f'==4&$ML_samp mlsum `lnf'=`ll' if (`todo'==0|`lnf'>=.) exit /************************************************************************ ************************************************************************/ tempname d1 d2 tempvar dlambda dalpha gen double `dlambda'=. if $ML_samp gen double `dalpha'=. if $ML_samp replace `dlambda'=`alpha'/`lambda'-`t0' if `f'==1&$ML_samp replace `dalpha'=-`theta1'+log(`t0')-digamma(`alpha') if `f'==1&$ML_samp replace `dlambda'=-`t0'*dgammapdx(`alpha',`lt0')/(1-`P0') if `f'==2&$ML_samp replace `dalpha'=-dgammapda(`alpha',`lt0')/(1-`P0') if `f'==2&$ML_samp replace `dlambda'=`t1'*dgammapdx(`alpha',`lt1')/`P1' if `f'==3&$ML_samp replace `dalpha'=dgammapda(`alpha',`lt1')/`P1' if `f'==3&$ML_samp replace `dlambda'=(`t1'*dgammapdx(`alpha',`lt1')-`t0'*dgammapdx(`alpha',`lt0'))/(`P1'-`P0') if `f'==4&$ML_samp replace `dalpha'=(dgammapda(`alpha',`lt1')-dgammapda(`alpha',`lt0'))/(`P1'-`P0') if `f'==4&$ML_samp replace `g1'=-`lambda'*`dlambda' if $ML_samp replace `g2'=`alpha'*`dalpha' if $ML_samp mlvecsum `lnf' `d1'=`g1', eq(1) mlvecsum `lnf' `d2'=`g2', eq(2) matrix `g'=(`d1', `d2') if (`todo'==1|`lnf'>=.) exit /************************************************************************ ************************************************************************/ tempname d11 d12 d22 tempvar dtheta11 dtheta12 dtheta22 d2lambda dlambdaalpha d2alpha gen double `d2lambda'=. gen double `dlambdaalpha'=. gen double `d2alpha'=. replace `d2lambda'=-`alpha'/(`lambda')^2 if `f'==1&$ML_samp replace `dlambdaalpha'=1/`lambda' if `f'==1&$ML_samp replace `d2alpha'=-trigamma(`alpha') if `f'==1&$ML_samp replace `d2lambda'=-(`t0')^2*dgammapdxdx(`alpha',`lt0')/(1-`P0')-(`dlambda')^2 if `f'==2&$ML_samp replace `dlambdaalpha'=-`t0'*dgammapdadx(`alpha',`lt0')/(1-`P0')-`dlambda'*`dalpha' if `f'==2&$ML_samp replace `d2alpha'=-dgammapdada(`alpha',`lt0')/(1-`P0')-(`dalpha')^2 if `f'==2&$ML_samp replace `d2lambda'=(`t1')^2*dgammapdxdx(`alpha',`lt1')/`P1'-(`dlambda')^2 if `f'==3&$ML_samp replace `dlambdaalpha'=`t1'*dgammapdadx(`alpha',`lt1')/`P1'-`dlambda'*`dalpha' if `f'==3&$ML_samp replace `d2alpha'=dgammapdada(`alpha',`lt1')/`P1'-(`dalpha')^2 if `f'==3&$ML_samp replace `d2lambda'=((`t1')^2*dgammapdxdx(`alpha',`lt1')-(`t0')^2*dgammapdxdx(`alpha',`lt0'))/(`P1'-`P0')-(`dlambda')^2 if `f'==4&$ML_samp replace `dlambdaalpha'=(`t1'*dgammapdadx(`alpha',`lt1')-`t0'*dgammapdadx(`alpha',`lt0'))/(`P1'-`P0')-`dlambda'*`dalpha' if `f'==4&$ML_samp replace `d2alpha'=(dgammapdada(`alpha',`lt1')-dgammapdada(`alpha',`lt0'))/(`P1'-`P0')-(`dalpha')^2 if `f'==4&$ML_samp gen double `dtheta11'=-`g1'+(`lambda')^2*`d2lambda' if $ML_samp gen double `dtheta12'=-`lambda'*`alpha'*`dlambdaalpha' if $ML_samp gen double `dtheta22'=`g2'+(`alpha')^2*`d2alpha' if $ML_samp mlmatsum `lnf' `d11'=-`dtheta11', eq(1) mlmatsum `lnf' `d12'=-`dtheta12', eq(1,2) mlmatsum `lnf' `d22'=-`dtheta22', eq(2) matrix `negH'=(`d11', `d12' \ `d12'', `d22') } } } end