mata, maximum likelihood using ml maximize, and nonconvergence

[William Rhodes]

I have written a Stata maximum likelihood routine to mimic stpm2, a user contributed routine to estimate survival models with restricted cubic splines. My routine reproduces results from stpm2. I have rewritten my routine with Mata using ml model lf as a shell. When I use results from stpm2 as starting values, and set ml maximize, difficult iterate(0) search(off), I get the correct likelihood, suggesting to me that my likelihood function is written correctly. But if I use these starting values and allow the program to iterate, I cannot achieve convergence. Nor can I achieve convergence if I use different starting values. Any suggestions for a next step?

mata:
void my_po_lf(transmorphic scalar ML, real rowvector b, real colvector lnfj)
{
real vector event
real vector XB, first_sum, second_sum, survivor, hazard
real vector parmindex, parms
real matrix DER
event = moptimize_util_depvar(ML,1)
XB = moptimize_util_xb(ML,b,1)
first_sum = moptimize_util_xb(ML,b,2)
parmindex = moptimize_util_eq_indices(ML, 2)
parms = b[parmindex[1,2]..parmindex[2,2]]
parms = parms[1,1..4]'
DER = st_data(.,"drcs1 drcs2 drcs3 drcs4")
second_sum = DER * parms
survivor = (1 :+ exp(first_sum) :* exp(XB)) :^ (-1)
hazard = survivor :* (exp(first_sum) :* exp(XB)) :* second_sum
lnfj = ((event :== 1) :* ln(hazard)) :+ ln(survivor)
}
end
ml maximize, difficult

[Paul Lambert]

There are two things needed to get this working

1. ml needs some help getting reasonable starting values. In the code below I use how stpm2 would obtain initial values for this model.
2. I am not sure of the exact reason, but it seems you can't extract and use the beta's as you do in your program when using method lf. I have thus changed it to a gf0 evaluator in the code below.

clear all
webuse brcancer
stset rectime, f(censrec=1) scale(365.24)

mata:
void my_po_gf0(transmorphic scalar ML, real scalar todo, real rowvector b, real colvector lnfj, real matrix S, real matrix H)
{
real vector event
real vector XB, first_sum, second_sum, survivor, hazard
real vector parmindex, parms
real matrix DER
event = moptimize_util_depvar(ML,1)
XB = moptimize_util_xb(ML,b,1)
first_sum = moptimize_util_xb(ML,b,2)
parmindex = moptimize_util_eq_indices(ML, 2)
parms = b[parmindex[1,2]..parmindex[2,2]]
parms = parms[1,1..4]'
DER = st_data(.,"_d_rcs1 _d_rcs2 _d_rcs3 _d_rcs4")
second_sum = DER * parms
survivor = (1 :+ exp(first_sum) :* exp(XB)) :^ (-1)
hazard = survivor :* (exp(first_sum) :* exp(XB)) :* second_sum
lnfj = ((event :== 1) :* ln(hazard)) :+ ln(survivor)
}
end

// stpm2 mmodel
stpm2 hormon, scale(odds) df(4)
estimates store stpm2

// obtain initial values 
stcox hormon, 
predict S0, basesurv
predict xb, xb
gen S = S0^exp(xb)
gen Z = log((1-S)/S)
regress Z hormon _rcs1 _rcs2 _rcs3 _rcs4 if _d
matrix initb = e(b)

// define and fit model
ml model gf0 my_po_gf0() (xb: _d = hormon, nocons) (time: _rcs1 _rcs2 _rcs3 _rcs4)
ml init initb, copy
ml maximize 
estimates store ml

// compare parameters
est tab stpm2 ml, keep(xb: time:)

[William Rhodes]

Thank you, Paul. I did use STPM2 to derive starting values. I will try your recommendation to use a gf0 evaluator. I have been struggling with this problem and I greatly appreciate your recommendations.