After running a probit, how can I generate the margins for the whole distribution?

Question

I'm using Stata.

I ran a probit of the form $$ \text{outcome}_i = \beta \ f(\text{income}_i) + \gamma\text{ Controls}_i $$

Where $f(\text{income}_i)$ is a fractional polynomial. I'm interested on the marginal probability $\frac{\delta \ P \ (\text{outcome}_i | \text{mean (Controls}_i))}{\delta \ \text{income}_i}$ for all values of $i$. When I run margins, dydx(varname) I only get the average marginal effect, not for the whole distribution. Moreover, I have to calcualte bootstrap standard errors, so I need some way to get the marginal probability that doesn't rely on specifying the functional form beforehand (i.e. _b[var1] + 2*_b[var2]*var2 for a 2nd order polynomial), as $f(\cdot)$ will change for each of the 500 bootstrap estimations.

A bit of context in case this is an instance of an xy problem.

I'm implementing an estimator of the form:

$$ \beta(\mathbf{x}) = \triangledown \Psi(\mathbf{x}) + \frac{\Psi(\mathbf{x})}{\mathrm{G_M}(\mathbf{x})}\triangledown \mathrm{G_M}(\mathbf{x}) $$

Where $\mathbf{x}$ is a vector of size $N$ (number of observations), $\Psi(\mathbf{x})$ is the prediction of the conditional mean, $\mathrm{G_M}(\mathbf{x})$ is the prediction of the probability that the outcome is greater than zero, and $\triangledown \Psi(\mathbf{x})$, $\triangledown \mathrm{G_M}(\mathbf{x})$ refer to their respective marginal effects.

The paper I'm following specifies:

But it doesn't include more information on this topic.

I hope I have been clear enough. Thanks in advance.

Answer 1

First of all, thanks to Dimitriy V. Masterov for pointing me towards the fact that I can save the margins for the whole distribution using margins, dydx(var1 var2), gen(newvar), which will generate newvar1 and newvar2 for var1 and var2 respectively.

Concerning the automation of marginal probability from the fp command, I wrote some code that can handle it, although only for two fp terms.

A fractional polynomial $x^{f(p)}+x^{f(q)}$ of $x$ is simply $x^p + x^q$, with the excepton that $x^{f(0)}=\ln(x)$. Also, if any term is repeated, say $x^{f(p)}+x^{f(p)}$ then the second term is $ln(x^{f(p)})$.

Our outcome is $y= x^{f(p)}+x^{f(q)}$. Then, $$ \frac{\delta y}{\delta x} = \frac{\delta y}{\delta x^{f(p)}}(\frac{\delta x^{f(p)}}{\delta x})+ \frac{\delta y}{\delta x^{f(q)}}(\frac{\delta x^{f(q)}}{\delta x}) $$

Where $\frac{\delta y}{\delta x^{f(p)}}$ and $\frac{\delta y}{\delta x^{f(q)}}$ are the results from margins, dydx(var1 var2) on the transformed variables,

and $\frac{\delta x^{f(p)}}{\delta x}$, $\frac{\delta x^{f(q)}}{\delta x}$ are the derivatives of the transformed variables with respect to $x$. Thus the code is as follows:

fp <ingpre_noL>, scale(1 1) replace: probit anytransfer <ingpre_noL> $controls i.id_entidad [pw=factor], robust cluster(upm) iterate(30)

margins, dydx(ingpre_noL_1 ingpre_noL_2) gen(me)

tempvar beta1p beta2p

scalar r = e(fp_fp)[1,1]

if r==0 {
    gen `beta1p' = 1/(ingpre_noL)
    }
else if r!=0 {
    gen `beta1p' = r*(ingpre_noL)^(r-1)
}

scalar t = e(fp_fp)[1,2]

if r!=t & t==0{
    gen `beta2p' = 1/(ingpre_noL+1)
}
else if r!=t & t!=0 {
    gen `beta2p' = t*(ingpre_noL+1)^(t-1)
}
else if r==t & t==0 {
    gen `beta2p' = 1/(ingpre_noL+1)*ln(ingpre_noL)
}
else if r==t & t!=0{
    gen `beta2p' = t/(ingpre_noL+1)
}

gen beta_p = me1*`beta1p' + me2*`beta2p'