Re: st: ml methods d1 and d2 and robust / clustered standard errors

[jpitblado@stata.com (Jeff Pitblado, StataCorp LP)]

Timothee Carayol <timothee.carayol@gmail.com> has an -ml- method -d1-
evaluator program designed prior to Stata 11, and asks why -ml method d1- no
longer works with the -vce(robust)- option:

> I am left rather confused with Stata's response to my attempts to
> obtain robust standard errors using the -ml- command with method d1 or
> d2 evaluator in Stata 11.
> I tried to reproduce (word by word) the Stata 11 base reference
> manual's example on this (page 1081); with no success.
> 
> My input:
> 
> ************
> cap program drop weib1
> program  weib1
> args  todo  b  lnf  g  negH  g1  g2 //  negH,  g1,  and  g2  are  new
> tempvar  t1  t2
> mleval  `t1'  =  `b',  eq(1)
> mleval  `t2'  =  `b',  eq(2)
> local  t  "$ML_y1"
> local  d  "$ML_y2"
> tempvar  p  M  R
> quietly  gen  double  `p'  =  exp(`t2')
> quietly  gen  double  `M'  =  (`t'*exp(-`t1'))^`p'
> quietly  gen  double  `R'  =  ln(`t')-`t1'
> mlsum  `lnf'  =  -`M'  +  `d'*(`t2'-`t1'  +  (`p'-1)*`R')
> if  (`todo'==0  |  `lnf'>=.)  exit
> tempname  d1  d2
> quietly  replace  `g1'  =  `p'*(`M'-`d') /*  <--  new         */
> quietly  replace  `g2'  =  `d'  -  `R'*`p'*(`M'-`d')       /*  <--
> new         */
> mlvecsum  `lnf'  `d1'    =  `g1',  eq(1) /*  <--  changed  */
> mlvecsum  `lnf'  `d2'    =  `g2',  eq(2) /*  <--  changed  */
> matrix  `g'  =  (`d1',`d2')
> end
> 
> use  http://www.stata-press.com/data/r11/cancer,  clear
> 
> gen drug2 = drug==2
> gen drug3 = drug==3
> 
> 
> ml  model  d1  weib1  (studytime  died  =  drug2  drug3  age)  /s,  vce(robust)
> ml  maximize
> *************
> 
> 
> Stata's output:
> 
> option vce(robust) is not allowed with evaltype d1
> 
> 
> In my understanding, the g1 and g2 variables should allow Stata to
> calculate the robust standard errors. It does work in Stata 10 and
> 10.1; but fails with Stata 11 and 11.1. So I can fix it by using
> -version 10-, but out of curiosity mainly I'd like to get it to work
> in versions 11 and 11.1.

With the release of Stata 11, -ml-'s optimization engine is implemented in
Mata through the -optimize()- routines.  With this change, the -d1- and -d2-
evaluator methods no longer work with -vce(robust)-, -pweight-s, or
-vce(cluster ...)-.  The old behavior continues to work only under version
control since the original -ml- engine was perserved and is callable only
under version control.

Although short-lived in the official on-line documentation, the -d1- and -d2-
methods supporting equation-level scores were renamed to -e1- and -e2- (see
-help ml_11-).  Note also that the modern -ml- with -e2- will need the -negh-
option because -ml- now assumes that the programmer computed second derivative
matrix is the Hessian matrix instead of the negative Hessian matrix (-negh-
implies that the evaluator returns the negative Hessian matrix).

Timothee can use the -e1- evaluator type in Stata 11 with -weib1- above.

Since the release of the updates in Stata 11.1, the official evaluator types
that work with equation level scores are named

	short name		long name
	---------------------------------------------------------------------
	lf			linearform

	lf0			linearform0
	lf1			linearform1
	lf2			linearform2

The -lf- evaluator type is as it was before.

Syntactically, the -lf#- evaluators are similar to the -d#- evaluators.

-lf0- is like -d0- except the evaluator is expected to return the
observation-level log-likelihood values instead of the overall log-likelihood
value.

-lf1- is like -d1- except the evaluator is expected to return the
observation-level log-likelihood values instead of the overall log-likelihood
value, and equation-level score (first derivatives) variables instead of the
gradient vector.

-lf2- is like -d2- except the evaluator is expected to return the
observation-level log-likelihood values instead of the overall log-likelihood
value, equation-level score (first derivatives) variables instead of the
gradient vector, and the Hessian matrix (second derivatives)

Here is a version of Timothee's -weib1- program modified to be an -lf1-
evaluator:

***** BEGIN:
program weib_lf1
	version 11.1
	args todo b lnf g1 g2
	tempvar t1 t2
	mleval `t1' = `b', eq(1)
	mleval `t2' = `b', eq(2)
	local t "$ML_y1"
	local d "$ML_y2"
	tempvar p M R
	quietly gen double `p' = exp(`t2')
	quietly gen double `M' = (`t'*exp(-`t1'))^`p'
	quietly gen double `R' = ln(`t')-`t1'
	quietly replace `lnf' = -`M' + `d'*(`t2'-`t1' + (`p'-1)*`R')
	if (`todo'==0) exit
	quietly replace `g1' = `p'*(`M'-`d')
	quietly replace `g2' = `d' - `R'*`p'*(`M'-`d') 
end
***** END:

Here is the Stata 11 version of the example Timothee posted using this new
evaluator and the 'i.' operator on the -drug- factor variable:

***** BEGIN:
. webuse cancer
(Patient Survival in Drug Trial)

. ml model lf1 weib_lf1 (studytime died = i.drug age) /s, vce(robust) maximize

initial:       log pseudolikelihood =       -744
alternative:   log pseudolikelihood = -356.14276
rescale:       log pseudolikelihood = -200.80201
rescale eq:    log pseudolikelihood = -136.69232
Iteration 0:   log pseudolikelihood = -136.69232  (not concave)
Iteration 1:   log pseudolikelihood = -124.13044  
Iteration 2:   log pseudolikelihood = -113.99048  
Iteration 3:   log pseudolikelihood = -110.30732  
Iteration 4:   log pseudolikelihood = -110.26748  
Iteration 5:   log pseudolikelihood = -110.26736  
Iteration 6:   log pseudolikelihood = -110.26736  

. ml display

                                                  Number of obs   =         48
                                                  Wald chi2(3)    =      30.69
Log pseudolikelihood = -110.26736                 Prob > chi2     =     0.0000

------------------------------------------------------------------------------
             |               Robust
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
eq1          |
        drug |
          2  |   1.012966   .2801321     3.62   0.000     .4639171    1.562015
          3  |    1.45917   .2878814     5.07   0.000      .894933    2.023407
             |
         age |  -.0671728   .0189346    -3.55   0.000    -.1042839   -.0300616
       _cons |   6.060723   1.023302     5.92   0.000     4.055089    8.066358
-------------+----------------------------------------------------------------
s            |
       _cons |   .5573333   .1359451     4.10   0.000     .2908859    .8237808
------------------------------------------------------------------------------
***** END:

--Jeff
jpitblado@stata.com
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/