Capture in iteration loop

[giacomomason]

Hi,

I'm using Stata 14.1, and in particular the "optimise" mata routine.
I'm writing some stata and mata code to perform permutation-based inference with inverse probability weights to correct for attrition. I'm generating these weights by performing simple maximum likelihood estimation of a logit model on R samples with randomly permuted treatment assignment vectors. This code is embedded in a mata function which is then called by an external stata routine.

My code is currently similar to the following:

for (r=1; r<=R; r++) { /* permutation loop */

/* computing OLS starting values (X are covariates, A is outcome) */
    data = A,X
    XpXi = quadcross(X, X)
    XpXi = invsym(XpXi)
    binit  = XpXi*quadcross(X, A) 

/* Optimisation ------------------------------- */
S=optimize_init()
optimize_init_evaluator(S,&logitf())
optimize_init_evaluatortype(S,"v0")
optimize_init_argument(S,1,data)
optimize_init_params(S,binit')
optimize(S)

/* get coefficients */
beta=optimize_result_params(S)'
/* compute linear index */
linind = X * beta

/* ... more stuff with the estimated coefficients */

} /* permutation loop */

My problem is that, in my empirical application, the optimiser fails in certain instances, for two main errors:
"convergence not achieved"
"Hessian is not negative semidefinite"
This is expected given the small sample size (around 250 observations, with lower numbers due to attrition).

I'm fine with the optimisation failing for some permutations, but I have two questions:
1) How can I manage these errors and avoid the whole stata routine that wraps my function to stop?
2) How can I keep track of how many of the permutations failed to optimise?

Thanks so much for your help!

Giacomo

[Christophe Kolodziejczyk]

You should use _optimize() which returns an error code instead of the parameter vector. The return code is 0 if the optimization converged. Based on these return codes you can take the appropriate actions (as recommended in the manual). And you can count the number of times the return codes is different from 0.
See the help optimize for further details.

nFailures = 0
for (r=1; r<=R; r++) 
{

...


rc = _optimize(S)

if (rc!=0) nFailures = nFailures + 1

...
}

[giacomomason]

Dear Cristophe,

thanks so much for your help, that is exactly what I was looking for! I ended up combining your help with the following to achieve the result I wanted and return both the optimisation result and the number of failures
http://www.stata.com/statalist/archi.../msg00633.html

Best,

Giacomo

[giacomomason]

Hi Cristophe,

Sorry for having to ask for your advice again. I thought everything was running smoothly, but I still can't get past the Hessian error.

I've switched to a while loop, and am now using _optimize(). My code now looks as follows:

r=1
fails = 0
while (r<=R) {

/* computing OLS starting values (X are covariates, A is outcome) */    
data = A,X    
XpXi = quadcross(X, X)    
XpXi = invsym(XpXi)    
binit  = XpXi*quadcross(X, A)
 
/* Optimisation ------------------------------- */
S=optimize_init()
optimize_init_evaluator(S,&logitf())
optimize_init_evaluatortype(S,"v0")
optimize_init_argument(S,1,data)
optimize_init_params(S,binit')
opt = _optimize(S)
conv = optimize_result_converged(S)  

if (opt!=0 | conv==0) { /* if failed */        
     fails = fails+1    
} else {     /* if successful */        
     beta=optimize_result_params(S)'        
     /* other operations ... */      
}
 
r=r+1  
} /* permutation loop */

I'm finding that this successfully counts the times the code fails to converge (increasing the fails count). However at the first occurrence of a "Hessian is not negative semidefinite" error the loop still stops!

Any help is greatly appreciated.

Best,

Giacomo

[Christophe Kolodziejczyk]

I wonder if you shouldn't rather code this way, i.e. _optimize is defined as void and not returning a return code and then retrieve the return code with optimize_result_code().

The manual pp.26-27 gives the following example. The following code stops the optimization in case of an error code, so you will obviously have to change this in order to fit into your problem.

(void) _optimize(S)
. . .
if (ec = optimize_result_code(S)) {
errprintf("{p}\n")
errprintf("%s\n", optimize_result_errortext(S))
errprintf("{p_end}\n")
exit(optimize_result_returncode(S))
/*NOTREACHED*/
}

[giacomomason]

Hi Cristophe,

thanks for all your help. I managed to solve the problem by switching optimisation algorithm. Apparently if numerical errors occur _optimize() with Newton-Raphson stops regardless of what I do, while other stuff like BHHH seems to manage to move to the next iteration. This is puzzling, but at least I'm happy it works now!


Best,

Giacomo