Re: st: Problems with ml max Survival analysis

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

Dear all again,

I have been working on the ml of the Likelihood for survival analysis and
I think that my problem is found with the gradient, because I get the following message:"Hessian is not negative semidefinite" when I run ml model d1, but running the same model with d0 "ml model d0" it works well. However, I need to use model d1 because in the following step I want to include cluster(id)
 
The model that I am trying to estimate:
h(t|x(t)=h0(t)*exp(beta'x(t) where the baseline hazard is piecewise (tp 28)
and beta'x include time-varying covariates and multi-spells (so far I am considering different spells of the individual as independent).
The Likelihood is: L= exp(Beta'X)^d*Exp(-sum(exp(beta'X)^(_t-_t0)
where d takes value 1 if failure is given and 0 otherwise.
I have defined the gradient for one individual as:
 d-sum(exp(beta'X)*(_t-_t0) taking only the last value for each spell.

Just below it is attached the code of the program, the function to maximize, the ml check and ml results. Please, any help is appreciated. Thank you very much for your help and time!

Amparo 

*****************************************
capture program drop CR_stable_3 
****Competing risk without uh
program CR_stable_3
version 11.2
args todo b lnf 

***Declare and define the arguments of the LL
tempvar beta1
mleval `beta1'=`b', eq(1)

***macros que tenemos que crear???
local dt0_1="$ML_y1"
local dt1_1="$ML_y2"
local d1="$ML_y3"

tempvar inthaz1 last  p1

*1) Likelihood calculation
tempvar sumb1
sort contadorbis quarter
by contadorbis: gen double `sumb1'=sum(`d1'*`beta1') if id_job05==1
tempvar haz1 
by contadorbis: gen double `haz1'=exp(`sumb1'[_N]) if  id_job05==1
by contadorbis: gen double `inthaz1'=sum(exp(`beta1')*(`dt1_1'-`dt0_1')) if id_job05==1
by contadorbis: gen double `p1'= `haz1'*exp(-`inthaz1'[_N]) if id_job05==1
by contadorbis: gen byte `last'=(_n==_N) if id_job05==1
mlsum `lnf'= ln(`p1') if `last'==1  & id_job05==1

*Calculate the gradient,I have to supply dlnl/dtheta and mlvecsum returns *dlnL/dbetai=sum dlnl/dtheta*Xij

if (`todo'==0 | `lnf'>=.) exit
mlvecsum `lnf' `g' = (`d1'-`inthaz1') if id_job05==1 & last==1, eq(1)

end

//--------------------------------------------------
// Construction of the duration variables
//-------------------------------------------------
gen t0_s=_t0
gen t_s=_t
ge d_s=_d

sort contadorbis quarter _t0

global ID_ "id_bis"

**************
**My model
ml model d1 CR_stable_3 (cr_05_stable:t0_s t_s d_s=tp* u_rate male h_skill m_skill non_manual municipio spanish_speakers no_spanish_speaker Aged05_16_19 /// 
Aged05_20_24 Aged05_25_29 Aged05_30_34 Aged05_35_39 Aged05_40_44 Aged05_45_51 Aged05_older61 responsabilities if id_job05==1, nocons)
*, cluster(`ID_')

ml check
ml init -5.70	-5.67	-5.66	-5.86	-5.93	-6.03	-6.01	-6.42	-6.41	-6.55	-6.68	-6.65	-6.64	-6.78	-6.86	-6.83	-6.96	-7.08	-7.06	-7.19	-7.43	-7.08	-7.05	-7.27	-6.29	-8.26	-7.86	-8.83	-3.98	0.04	0.12	-0.05	-0.04	-0.09	-0.08	-0.18	0.05	0.33	0.43	0.40	0.40	0.46	0.46	-0.62	0.13,copy

ml max, difficult
*************************
ml check

Test 1:  Calling CR_stable_3 to check if it computes log likelihood and
         does not alter coefficient vector...
         Passed.

Test 2:  Calling CR_stable_3 again to check if the same log likelihood value
         is returned...
         Passed.

------------------------------------------------------------------------------
The initial values are not feasible.  This may be because the initial values
have been chosen poorly or because there is an error in CR_stable_3 and it
always returns missing no matter what the parameter values.

Stata is going to search for a feasible set of initial values.
If CR_stable_3 is broken, this will not work and you will have to press Break
to stop the search.

Searching...
initial:       log likelihood =     -<inf>  (could not be evaluated)
searching for feasible values ....--Break--
r(1);


*********************************************
. ml max, difficult

initial:       log likelihood =  -247110.7
rescaling entire vector ..
rescale:       log likelihood =  -247110.7
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Iteration 0:
                                                   log likelihood =  -247110.7
Gradient length       =        .
Step length           =        .
Parameters + step -> new parameters
                                                   log likelihood =          0
                                                           (initial step good)
(1) Stepping forward, step length =        .
                                                   log likelihood =          0
                                                          (ignoring last step)
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Iteration 1:
                                                   log likelihood =          0
                                                             vtol =  .99999595
Gradient length       =        .
Step length           =        .
Parameters + step -> new parameters
                                                   log likelihood =          0
                                                           (initial step good)
(1) Stepping forward, step length =        .
                                                   log likelihood =          0
                                                          (ignoring last step)
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Iteration 2:
Hessian is not negative semidefinite
r(430);

end of do-file

r(430);





*
*   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/