Using -margins- to get SE of dydx with terms involving X and sqrt(X) in a nonlinear model

[Mike Lacy]

I want to use -margins- to get the derivative with respect to a variable X and the standard error of that derivative, where X appears in the model as both a linear term and a square root term. I know how to get the derivative, but I don't know how to get the standard error.

I would like to do this in the context of a -heckprobit- model, where the derivative of interest is for the probability in the outcome model with respect to X. To my understanding, the problem in using -margins- would be essentially the same if I had a plain -logit- model and a function other than sqrt(). On that count, I'll hope that people who don't use -heckprobit- will nevertheless take a look at this. I'll stick with -heckprobit- and sqrt() to be precise.

-margins- will give the derivative of the probability w.r.t. to each term involving X, which can then be combined into the overall derivative with respect to X according to conventional rules, and I have code below to show that, using the example data from -help heckprobit-, and using "years" as the "X" variable.

Can someone suggest how -margins- or other built-in means could be used to estimate the standard error of this overall derivative? Worst case, I could use -bootstrap-, but that seems like a hard way to go. As I understand it, the dydx option on -margins- only accepts varlists, not expressions.

Thanks, of course, for whatever thoughts you might have.

use http://www.stata-press.com/data/r15/school, clear
gen years_sqrt = sqrt(years)
heckprobit private years years_sqrt, select(vote=years years_sqrt loginc logptax)
//
// -margins- gives dydx w.r.t. linear and sqrt terms involving -year-
local yrval = 9  // Suppose years = 9 is of interest
local yrvalrt = sqrt(`yrval')
margins, dydx(years) atmeans at(years == `yrval' years_sqrt = `yrvalrt')
local d1 = el(r(b), 1, 1)     // linear
//
margins, dydx(years_sqrt) atmeans at(years == `yrval' years_sqrt = `yrvalrt')
local d2 = el(r(b), 1, 1)     // sqrt term
//
// Combine derivative terms
local dydx_years = `d1' + `d2'/(2 * `yrvalrt')
di "dydx w.r.t. years = `dydx_years'"

[Andrew Musau]

Can you not reparametrize this so that you have a quadratic term? In this way, using factor variables, margins understands that years and its square root are related.

use http://www.stata-press.com/data/r15/school, clear
gen years_sqrt = sqrt(years)
heckprobit private years years_sqrt, select(vote=years years_sqrt loginc logptax) nolog
*REPARAMETERIZATION
heckprobit private c.years_sqrt##c.years_sqrt, select(vote=c.years_sqrt##c.years_sqrt loginc logptax) nolog
margins, dydx(years_sqrt) atmeans at(years_sqrt=3)

Res.:

. heckprobit private years years_sqrt, select(vote=years years_sqrt loginc logptax) nolog

Probit model with sample selection              Number of obs     =         95
                                                      Selected    =         59
                                                      Nonselected =         36

                                                Wald chi2(2)      =       3.22
Log likelihood = -63.94751                      Prob > chi2       =     0.2003

------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
private      |
       years |  -3.302895   2.063555    -1.60   0.109     -7.34739    .7415988
  years_sqrt |   12.38636   7.396465     1.67   0.094    -2.110444    26.88317
       _cons |  -11.84635   6.623069    -1.79   0.074    -24.82733    1.134623
-------------+----------------------------------------------------------------
vote         |
       years |   .2901484   .0849602     3.42   0.001     .1236294    .4566674
  years_sqrt |  -2.316325    .609562    -3.80   0.000    -3.511044   -1.121605
      loginc |   1.039763   .5801362     1.79   0.073    -.0972827    2.176809
     logptax |  -1.604211   .6688682    -2.40   0.016    -2.915169   -.2932534
       _cons |   4.730536    4.74375     1.00   0.319    -4.567044    14.02812
-------------+----------------------------------------------------------------
     /athrho |  -2.013566   3.309914    -0.61   0.543    -8.500879    4.473747
-------------+----------------------------------------------------------------
         rho |  -.9649736   .2278082                     -.9999999    .9997399
------------------------------------------------------------------------------
LR test of indep. eqns. (rho = 0):   chi2(1) =     0.72   Prob > chi2 = 0.3966

.
. *REPARAMETERIZATION

.
. heckprobit private c.years_sqrt##c.years_sqrt, select(vote=c.years_sqrt##c.years_sqrt loginc logptax) nolog

Probit model with sample selection              Number of obs     =         95
                                                      Selected    =         59
                                                      Nonselected =         36

                                                Wald chi2(2)      =       3.22
Log likelihood = -63.94751                      Prob > chi2       =     0.2003

-------------------------------------------------------------------------------------------
                          |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
--------------------------+----------------------------------------------------------------
private                   |
               years_sqrt |   12.38635   7.396458     1.67   0.094    -2.110439    26.88314
                          |
c.years_sqrt#c.years_sqrt |  -3.302893   2.063553    -1.60   0.109    -7.347383    .7415974
                          |
                    _cons |  -11.84635   6.623063    -1.79   0.074    -24.82731    1.134619
--------------------------+----------------------------------------------------------------
vote                      |
               years_sqrt |  -2.316325   .6095619    -3.80   0.000    -3.511044   -1.121605
                          |
c.years_sqrt#c.years_sqrt |   .2901484   .0849602     3.42   0.001     .1236294    .4566674
                          |
                   loginc |   1.039763   .5801364     1.79   0.073    -.0972833     2.17681
                  logptax |  -1.604211   .6688683    -2.40   0.016    -2.915169    -.293253
                    _cons |   4.730536   4.743751     1.00   0.319    -4.567046    14.02812
--------------------------+----------------------------------------------------------------
                  /athrho |  -2.013566   3.309918    -0.61   0.543    -8.500887    4.473754
--------------------------+----------------------------------------------------------------
                      rho |  -.9649736   .2278081                     -.9999999    .9997399
-------------------------------------------------------------------------------------------
LR test of indep. eqns. (rho = 0):   chi2(1) =     0.72   Prob > chi2 = 0.3966

.
. margins, dydx(years_sqrt) atmeans at(years_sqrt=3)

Conditional marginal effects                    Number of obs     =         95
Model VCE    : OIM

Expression   : Pr(private=1), predict()
dy/dx w.r.t. : years_sqrt
at           : years_sqrt      =           3
               loginc          =    9.971017 (mean)
               logptax         =    6.939496 (mean)

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
  years_sqrt |  -.0001748   .0023449    -0.07   0.941    -.0047707    .0044212
------------------------------------------------------------------------------

.

You are interested in the marginal effect of years, but that is simply a matter of scaling.

[Mike Lacy]

Thanks Andrew Musau I get your basic idea here -- simply clever -- but I'm not sure I'm thinking right about what you describe as rescaling, which perhaps I might be overcomplicating. I'm thinking like this:

We use -margins- to get us dy/dw, where w = x^1/2. We want dy/dx, so we use the chain rule:
dy/dx = dy/dw * dw/dx = dy/dw * 1/[2sqrt(x)]

Is this what you mean? And, presuming that's right, what this would mean for translating the standard error from the sqrt() scale to the linear scale.

[Andrew Musau]

As the marginal effect represents a unit change, the marginal effect of the square root term and square term should be exactly the same as the square of one is one. It should be possible to show this, so I will post the math once I am able to write it down. But let me know if you are finding a fault in my logic.

[Andrew Musau]

This replicates the calculation of the AME by hand, so doing so for the square root term in the second parametrization and the square term in the first results in about the same estimated AME.

use http://www.stata-press.com/data/r15/school, clear
gen years_sqrt = sqrt(years)
heckprobit private c.years_sqrt##c.years_sqrt, select(vote=c.years_sqrt##c.years_sqrt loginc logptax) nolog
margins, dydx(years_sqrt)
qui sum years_sqrt if e(sample)
scalar h= (abs(r(mean))*.001)*.001
predict double private_0 if e(sample)
rename years_sqrt years_s
gen double years_sqrt= years_s+scalar(h)
predict double private_1 if e(sample)
gen double dydx= (private_1-private_0)/scalar(h)
sum dydx


use http://www.stata-press.com/data/r15/school, clear
gen years_sqrt = sqrt(years)
heckprobit private years years_sqrt, select(vote=years years_sqrt loginc logptax) nolog
qui sum years if e(sample)
scalar h= (abs(r(mean))*.001)*.001
predict double private_0 if e(sample)
rename (years years_sqrt) (year years_s)
gen double years= year+scalar(h)
gen double years_sqrt = sqrt(years)
predict double private_1 if e(sample)
gen double dydx= (private_1-private_0)/scalar(h)
sum dydx

Res.:

.
. margins, dydx(years_sqrt)

Average marginal effects                        Number of obs     =         95
Model VCE    : OIM

Expression   : Pr(private=1), predict()
dy/dx w.r.t. : years_sqrt

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
  years_sqrt |   .0305817   .0744436     0.41   0.681     -.115325    .1764884
------------------------------------------------------------------------------

.
. sum dydx

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
        dydx |         95    .0305809    .4751426  -.7630771   .7843245

.
.
. sum dydx

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
        dydx |         95    .0292057    .1433681  -.1859099   .2863702

[FernandoRios]

Hi Mike
not on my computer right now, but I have the command f_able
check it out as well as the stata journal article
It can do exactly what you want as long as x is strictly positive
there is also my presentation of couple of years ago at the stata conference
Fernando

[FernandoRios]

Here is the example

ssc install f_able
use http://www.stata-press.com/data/r15/school, clear

* This creates the variable and keeps the info on how it was created 
. fgen years_sqrt = sqrt(years)

. heckprobit private years years_sqrt, select(vote=years years_sqrt loginc logptax)

Probit model with sample selection              Number of obs     =         95
                                                      Selected    =         59
                                                      Nonselected =         36

                                                Wald chi2(2)      =       3.22
Log likelihood = -63.94751                      Prob > chi2       =     0.2003

------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
private      |
       years |  -3.302893   2.063553    -1.60   0.109    -7.347384    .7415974
  years_sqrt |   12.38635   7.396459     1.67   0.094    -2.110439    26.88315
       _cons |  -11.84635   6.623064    -1.79   0.074    -24.82731    1.134619
-------------+----------------------------------------------------------------
vote         |
       years |   .2901484   .0849602     3.42   0.001     .1236294    .4566674
  years_sqrt |  -2.316325   .6095619    -3.80   0.000    -3.511044   -1.121605
      loginc |   1.039763   .5801364     1.79   0.073    -.0972833     2.17681
     logptax |  -1.604211   .6688683    -2.40   0.016    -2.915169   -.2932529
       _cons |   4.730536   4.743751     1.00   0.319    -4.567045    14.02812
-------------+----------------------------------------------------------------
     /athrho |  -2.013566   3.309918    -0.61   0.543    -8.500886    4.473754
-------------+----------------------------------------------------------------
         rho |  -.9649736   .2278081                     -.9999999    .9997399
------------------------------------------------------------------------------
LR test of indep. eqns. (rho = 0): chi2(1) = 0.72         Prob > chi2 = 0.3966

** This adds this info back to e() so it can estimate the derivatives
. f_able years_sqrt, auto
This is an experimental feature in f_able
Using this option, you do not need to add nochain or numerical options in margins
All variables in -nlvar- have been declared

. margins, dydx(years)

Average marginal effects                                    Number of obs = 95
Model VCE: OIM

Expression: Pr(private=1), predict()
dy/dx wrt:  years

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       years |   .0269328   .0241411     1.12   0.265    -.0203829    .0742486
------------------------------------------------------------------------------

Hope this helps
F

[Mike Lacy]

Thanks to both of you.

Andrew Musau : I'm going to have to think about what you wrote and see if I understand or at least ask an intelligent question <grin>. There's something I'm not quite getting here, but I've read enough of your postings to know that I'd implicitly trust your mathematical reasoning over mine.

FernandoRios : This looks great, very easy and general. Look forward to checking out your article to see what's inside your routine.

I'll follow up here with questions/comments when I've digested what's here.

[FernandoRios]

Absolutely! Always happy to share and answer this kind of questions!

[Mike Lacy]

Thanks again to both of you for your interest and efforts.

I have some summary comments here, with questions, and some example code:

First, speaking to Andrew: I like the simplicity of your reparameterization suggestion, but to the best of my understanding, it gives the derivative with respect to sqrt(years), while I want the derivative with respect to years itself. The only way I understand to get that is to apply the chain rule to the derivative your method gives for the sqrt(years) term. That result agrees with what I get from my original clumsy approach (illustrated in code below). I still don't see how either of our approaches can give the SE of that derivative, and I still don't understand what you had in mind about rescaling--no criticism intended; I'm just missing something. Presumably the SE could be obtained by a Delta Method treatment of the variance of the dydx of the sqrt(years) term, but that manipulation is beyond my ken.

Speaking to Fernando: Your approach, not surprisingly, conveniently yielded a result that agrees with what I got in my original approach, and also with what I obtained from Andrew's approach combined with the chain rule. (See code at the end below). And, using your -f_able- command *does* give a standard error.

Thinking now that Fernando's method gave me everything I needed, I decided to compare the standard error it gives what I could get from -bootstrap-. I wrote a wrapper for Fernando's method, and fed it to -bootstrap-, but it kept resulting in impossibly long iterations and no solution. I also tried the same bootstrap approach with my own method and had the same result. I have diagnosed such problems before, but couldn't figure this one out. There seems to be something about the combination of -heckprobit- and -margins- and -bootstrap- that causes problems. This is a bit of a side note here, but for whatever interest that might offer, here's some code to illustrate that problem:

cap prog drop dydxyears
prog dydxyears, rclass
// wrapper for FernandoRios approach
syntax, atspec(string)
cap drop years_sqrt
fgen years_sqrt = sqrt(years)
heckprobit private years years_sqrt, ///
   select(vote=years years_sqrt loginc logptax)
if (_rc == 0 ) { worked ok   
    f_able years_sqrt, auto
    margins, dydx(years) `atspec'
    return scalar dydx = el(r(b),1,1)
}
else {
   return scalar dydx = .
}
end
//
use http://www.stata-press.com/data/r15/school, clear
bootstrap dydx = r(dydx), noisily reps(1000) seed(234): dydxyears, atspec(atmeans at(years = 9))

The program executes fine on the original data, but the bootstrap iterations result in lots of "not concave" messages, wild LL values, and interminable iteration. While I have no reason to distrust the SE given by the -ftable- approach, being able to check it on a few examples would be nice.

Back to the main theme: For anyone wanting to compare our three approaches as regards the identity of the derivative estimates they yield, here's some example code:

// Comparing three approaches to dydx from -heckprobit- with a sqrt() term
use http://www.stata-press.com/data/r15/school, clear
local atval = 9
local atvalrt = sqrt(`atval')
// Per Andrew Musau's reparameterization approach
gen years_sqrt = sqrt(years)
heckprobit private c.years_sqrt##c.years_sqrt, select(vote=c.years_sqrt##c.years_sqrt loginc logptax) nolog
margins, dydx(years_sqrt) atmeans at(years_sqrt = `atvalrt' )
// Use his result and the chain rule to obtain dydx(years) from dydx(years_sqrt)
local dMusau = el(r(b),1,1)/(2*`atvalrt')
//
// Per FernandoRios, using his -f_able-
drop years_sqrt
fgen years_sqrt = sqrt(years)
heckprobit private years years_sqrt, select(vote=years years_sqrt loginc logptax)
f_able years_sqrt, auto
margins, dydx(years) atmeans at(years =`atval')
local dFernandoRios = el(r(b),1,1)
//
// Per Lacy's original approach
drop years_sqrt
gen years_sqrt = sqrt(years)
heckprobit private years years_sqrt, select(vote=years years_sqrt loginc logptax)
quiet margins, dydx(years) atmeans at(years = `atval' years_sqrt = `atvalrt')
local d1 = el(r(b), 1, 1)     // linear
quiet margins, dydx(years_sqrt) atmeans at(years = `atval' years_sqrt = `atvalrt')
local d2 = el(r(b), 1, 1)     // sqrt term
// Combine derivative terms
local dydx_years = `d1' + `d2'/(2 * `atvalrt')
local dLacy = `dydx_years'
foreach N in Musau FernandoRios Lacy {
   di "dydx_`N' = " %12.8f `d`N''
}

[FernandoRios]

Hi Mike
I think the bootstrap approach for the heckprobit may be asking a bit too much of the data because the rho parameter is hitting very close to 1. Which makes the identification difficult to estimate.
So, what i would suggest is something slightly different.
In my SJ paper, for example, I compare the results of using f_able, with those using factor notation. And compare to -nl-

sysuse auto
regress price c.mpg##c.mpg
margins, dydx(*)
fgen mpg2=mpg*mpg
regress price mpg mpg2
f_able, auto
f_able mpg2, auto
margins, dydx(*)

nl (price = {a0}+{a1}*mpg+{a2}* sqrt(mpg)), variables(mpg)
margins, dydx(mpg)
fgen mpgsqrt=sqrt(mpg)
reg price mpg mpgsqrt
f_able mpgsqrt, auto
margins, dydx(*)

In fact, the way f_able works is similar to -nl-. It tells margins to estimate numerical derivatives of variables that are transformations of others. Updating it internally everything it does.

F

[Mike Lacy]

Hi Fernando,
It would be interesting to try the bootstrap on some nicer data set. I can understand that this might be a small and brittle example data set, but at least I'd say that there is a decent univariate distribution on the response variable. It just puzzles me that that wrapper program behaved find on the original data set but immediately get lost while bootstrapping at N = _N. I tried a couple of different seeds, thinking that perhaps it was just hitting a funky sample, but that didn't cure the problem. The actual data set I want to use this on is probably not the best set either, so perhaps I can think of something else to try it on for a test.

[Andrew Musau]

First, speaking to Andrew: I like the simplicity of your reparameterization suggestion, but to the best of my understanding, it gives the derivative with respect to sqrt(years), while I want the derivative with respect to years itself. The only way I understand to get that is to apply the chain rule to the derivative your method gives for the sqrt(years) term. That result agrees with what I get from my original clumsy approach (illustrated in code below). I still don't see how either of our approaches can give the SE of that derivative, and I still don't understand what you had in mind about rescaling--no criticism intended; I'm just missing something. Presumably the SE could be obtained by a Delta Method treatment of the variance of the dydx of the sqrt(years) term, but that manipulation is beyond my ken.

Mike Lacy, first of all, thank you for the really insightful illustration. I went through the math and it agrees with your derivation in #3, so there is no need to reproduce it here. I think that we agree that when we are talking about the derivative of y with respect to either x or the square root of x, these are simply transformations of each other and not two separate derivatives (in the sense that you cannot change the square root of x without changing x). As a matter of fact, you have shown that if we define the derivative wrt the square root term as dy/dw, then dy/dx = dy/dw * 1/[2sqrt(x)] in #3, where the scaling parameter is 1/[2sqrt(x)]. Therefore, we can insert this scaling parameter within -expression()- in margins to obtain the desired derivative for years.

use http://www.stata-press.com/data/r15/school, clear
local atval = 9
local atvalrt = sqrt(`atval')
// Per Andrew Musau's reparameterization approach
gen years_sqrt = sqrt(years)
heckprobit private c.years_sqrt##c.years_sqrt, select(vote=c.years_sqrt##c.years_sqrt loginc logptax) nolog
margins, dydx(years_sqrt) atmeans at(years_sqrt = `atvalrt' ) expression(predict()*1/(2*`atvalrt'))

Res.:

. margins, dydx(years_sqrt) atmeans at(years_sqrt = `atvalrt' ) expression(predict()*1/(2*`atvalrt'))

Conditional marginal effects                    Number of obs     =         95
Model VCE    : OIM

Expression   : predict()*1/(2*3)
dy/dx w.r.t. : years_sqrt
at           : years_sqrt      =           3
               loginc          =    9.971017 (mean)
               logptax         =    6.939496 (mean)

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
  years_sqrt |  -.0000291   .0003908    -0.07   0.941    -.0007952    .0007369
------------------------------------------------------------------------------

I quite like your idea to verify this computation by bootstrapping in whatever dataset is possible.

[Enrique Pinzon (StataCorp)]

Andrew, Fernando, and Mike:

Needless to say, I found your discussion to be really fun. I want to add a quick suggestion that may be helpful for anyone thinking about adding -margins- in a bootstrap routine. -margins- has a -nose- option that tells -margins- not to compute standard errors. Depending on the number of parameters, given that -margins- uses the delta method, this will save a lot of computational time. So I would call -margins- with -nose- in your -bootstrap- routine.

[Mike Lacy]

The "TLDR;" here is that the new code from Andrew completely agrees with Fernando's, but both produce a standard error estimate about 40% *higher* than what I get from -bootstrap- using a data set modified to be friendly to -heckprobit-. That's a bigger difference than I would expect to arise from the asymptotic approximation of the Delta method. A summary of these results and code to produce them appears below.

A comment first:

Andrew: Your use of the predict(expression) to be "tacked onto" the dyd() option was just the kind of thing I was hoping for. Kudos to you for working out how to use what -margins- offers. I'm understanding you now re "rescaling". I was mentally focused more on the SE, and thinking of how in general var(f(x) != f(var(x)), so I wasn't thinking of the derivative as a rescaling.

Summary of results:
-margins- report from Andrew's method (identical results from Fernando's, which are omitted here)

Conditional marginal effects                    Number of obs     =        475
Model VCE    : OIM
 
Expression   : predict()*1/(2*3)
dy/dx w.r.t. : years_sqrt
at           : years_sqrt      =           3
               loginc          =    9.971017 (mean)
               logptax         =    6.939496 (mean)
 
------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
  years_sqrt |  -.0352885   .0246137    -1.43   0.152    -.0835304    .0129534
------------------------------------------------------------------------------

-bootstrap- report using Andrew's method (same to about 3 sig. digits as Fernando's, omitted here)

    command:  Mdydxyears, atspec(atmeans at(years = 9)) atvalrt(3)
         dydx:  r(dydx)

------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        dydx |  -.0352885   .0175207    -2.01   0.044    -.0696284   -.0009486
------------------------------------------------------------------------------
Note: One or more parameters could not be estimated in 1 bootstrap replicate;
      standard-error estimates include only complete replications.

Code to produce the preceding results:

cap prog drop FRdydxyears
cap prog FRdydxyears, rclass
// wrapper for FernandoRios approach
syntax, atspec(string) [wantse]
local nose = cond("`wantse'" != "", "", "nose")
cap drop years_sqrt
fgen years_sqrt = sqrt(years)
heckprobit private years years_sqrt, ///
   select(vote=years years_sqrt loginc logptax) nolog
if (_rc == 0 ) { worked ok  
    f_able years_sqrt, auto
    margins, dydx(years) `atspec' `nose'
    return scalar dydx = el(r(b),1,1)
}
else {
   return scalar dydx = .
}
end
//
cap prog drop Mdydxyears
prog Mdydxyears, rclass
// wrapper for Musau approach
syntax, atspec(string) atvalrt(real) [wantse]
local nose = cond("`wantse'" != "", "", "nose")
heckprobit private c.years_sqrt##c.years_sqrt, ///
   select(vote=c.years_sqrt##c.years_sqrt loginc logptax) nolog
if (_rc == 0 ) { worked ok  
    margins, dydx(years_sqrt) atmeans at(years_sqrt = `atvalrt' ) expression(predict()*1/(2*`atvalrt'))
    return scalar dydx = el(r(b),1,1)
}
else {
   return scalar dydx = .
}
end
//
//
use http://www.stata-press.com/data/r15/school, clear
set seed 46534
// Expand the data set, tweak the years variable to reduce rho and make
// the -heckprobit- estimation less brittle.  Keep iterations reasonable.
expand 5
replace years = years + ceil(runiform() * 5) if (years < 6)
set maxiter 100
//
local atval = 9
local atvalrt = sqrt(`atval')
gen years_sqrt = sqrt(years)
// With asymptotic SE
FRdydxyears, atspec(atmeans at(years = `atval')) wantse
Mdydxyears, atspec(atmeans at(years = `atval'))  atvalrt(`atvalrt') wantse
//
bootstrap dydx = r(dydx), reps(1000) seed(83956): ///
   FRdydxyears, atspec(atmeans at(years = `atval'))
//
bootstrap dydx = r(dydx), reps(1000) seed(83956): ///
   Mdydxyears, atspec(atmeans at(years = `atval'))  atvalrt(`atvalrt')