Delta Method to obtain robust standard errors for discrete increases in each variables holding - the other constants - LOGIT model

[Linda Luciani]

Hello,

I have a logit regression that contains both categorical variables and continuous variables as regressors. I am interested in DISCRETE EFFECTS/INCREASES at the MEDIANS and therefore I have computed those. Now I want to compute robust standard errors and I was suggested to use the DELTA METHOD, although I am finding to understand how to apply it in my case. Could you please provide me with a code that can make me obtain the correct robust standard errors?

logit repair sav_rate inc_a age HighEd1 fam price_ratio c.year [pweight=weighta]

For example, this is the way I have computed the Discrete effect at the median of HIghEd1, for a discrete increase of delta_increase=+4

*median of HighEd1
gen educ=4
gen double xb0 = _b[sav_rate]*savr + _b[inc_a]* ina + _b[age]* ag + _b[fam]*siz+ _b[HighEd1]*4 + _b[price_ratio]* pri + _b[year]*2009 + _b[_cons]

generate double pe = logistic(xb0 -_b[HighEd1]*educ + _b[HighEd1]*(educ+4)) - logistic(xb0)
di pe
*.01164981

Thank you in advance!

[Mike Lacy]

You're on the edge of where I consider myself entirely competent, but let me suggest a few things that might help.

1) You are calculating the discrete effect *for each observation in your data set using its observed vector of predictor variables,* except that you have specified that the value of education to be 4 in your first equation, and 8 in the second one. Is this what you want, i.e., different estimates of the discrete effect for each different observation?

Note that when you issue the command "di pe," Stata will display the value of the variable pe for the first observation. The -display- command is not designed to be applied to whole variables, like pe, just single (scalar) values, so Stata defaults to the value of the variable for the first observation. Therefore, with what you show, you will obtain the discrete effect of education, comparing educ = 4 vs. educ = 8, for an observation that has the vector of predictors of your first observation, whatever those might be. I would guess that's not really what you want.

2) The -margins- command in Stata would likely do you want. If I'm correct, you're really doing things the hard way. The -margins- command allows you to specify the exact vector of predictor values you want. And, it allows you to obtain predictions while varying one of the predictors. It also gives standard errors of the predictions. Presuming you agree with what I said in 1), you'll want to stop and first learn about using -margins-. A useful source for learning about -margins- is https://www.stata-journal.com/articl...article=st0260
Presuming some familiarity with -margins-, you will find that you can use it to produce predictions for discrete changes in a predictor. For example:

sysuse auto, clear
logit foreign price weight, robust
// Compare effect on predicted Prob(foreign==1) for an observation
// with price = 6000 and weight = 3000 vs. 3500
margins, at(price = (6000) weight = (3000 3500))
// These quantities are available for further programming
mat list r(b) // predicted values
mat list r(V) // variance matrix for predictions

You could then pick values out of those matrices to compute the variance and standard error of the difference of predicted values for the discrete change associated with the two values of weight.

3) I'm understanding you here to want a Delta-Method formula to calculate robust standard errors for a *predicted* value. The -margins- command does give Delta-Method standard errors for predicted values, but I don't know whether they are based on robust methods. (In fact, I don't even know -- my ignorance -- whether robust standard errors is even a concept that applies to functions of predicted values, such as you are calculating. You can estimate your model with robust standard errors as I did above, and from a brief experiment I did, I found that it does change the Delta-Method estimates of the standard error of the predicted values.)

I hope I understood correctly what you were looking for, as you may well know more about the theory and of course your intentions than I do.

[Linda Luciani]

Dear Professor Mike Lacy ,

thank you for your kind reply.

1) Just to be clear I will explain what I think I have obtained though my code. First of all savr, ina, age , siz , pri, 2009 for the year and educ=4 are all the medians values of my sample and have been previously computed.
"pe" returns the change in the predicted probability when all the variables are fixed at their level but education increases from HighEd1=4 (ie the median value) to HighEd1=8. This said I want to see if this increase in predicted probability (ie pe= .01164981) that I have obtained is significant and therefore I need the standard errors I suppose.

2) As far as I have understood, your code is alternative to mine to the extent that .044892 is the predicted probability at price=6000 and weight=3000, while .0024804 is the predicted probability at price=6000 and weight=3500. Therefore (.0024804- .044892) returns the decrease in probability from an increase in weight=3500-3000=+ 500. Is the this point correct? Although I am not sure I have understood how you would compute the robust standard errors from results in mat list r(V)?
Could you actually show it to me through numbers and mention what the method is called in this case?

. margins, at(price = (6000) weight = (3000 3500))

Adjusted predictions Number of obs = 74
Model VCE : Robust

Expression : Pr(foreign), predict()

1._at : price = 6000
weight = 3000

2._at : price = 6000
weight = 3500

------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_at |
1 | .044892 .0386933 1.16 0.246 -.0309454 .1207294
2 | .0024804 .0041676 0.60 0.552 -.005688 .0106489
------------------------------------------------------------------------------


. mat list r(V)

symmetric r(V)[2,2]
1. 2.
_at _at
1._at .00149717
2._at .00015704 .00001737

3) I know that the robust standard errors through - margins- function can be obtained through the specification- vce(unconditional)-
in the following way:

. margins, at(price = (6000) weight = (3000 3500)) vce(unconditional)

Adjusted predictions Number of obs = 74

Expression : Pr(foreign), predict()

1._at : price = 6000
weight = 3000

2._at : price = 6000
weight = 3500

------------------------------------------------------------------------------
| Unconditional
| Margin Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_at |
1 | .044892 .0386933 1.16 0.246 -.0309454 .1207294
2 | .0024804 .0041676 0.60 0.552 -.005688 .0106489
------------------------------------------------------------------------------

.
Although the results are identical to the your code without specification, and therefore I suppose -margins- calculates robust standard errors by default

. margins, at(price = (6000) weight = (3000 3500))

Adjusted predictions Number of obs = 74
Model VCE : Robust

Expression : Pr(foreign), predict()

1._at : price = 6000
weight = 3000

2._at : price = 6000
weight = 3500

------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_at |
1 | .044892 .0386933 1.16 0.246 -.0309454 .1207294
2 | .0024804 .0041676 0.60 0.552 -.005688 .0106489
------------------------------------------------------------------------------



Thank you so much for your help!


Best,
Linda

[Mike Lacy]

You had said:

1) Just to be clear I will explain what I think I have obtained though my code. First of all savr, ina, age , siz , pri, 2009 for the year and educ=4 are all the medians values of my sample and have been previously computed.

We had no way to know that. Apparently savr, ina, etc. are Stata scalar values? What you did, then, carries out exactly the same calculation, with exactly the same scalars, for every single observation in your file, because you used -generate-, which creates values for all variables. You can do this, but no one would normally do that, since why would you want the same calculation done +_N times?
For that reason, I assumed that those values you had were variables, not scalar values. What you were doing was reasonable, but doesn't fit with how Stata and similar languages work.
To do what you want to do, you'd just do that computation and assign the result to a Stata scalar or local. I'd suggest you get someone to explain this kind of thing to you a bit, or read more extensively about Stata scalars and locals. Here's how you could use the "do it yourself" calculation that you were trying to do:

scalar xbo = _b[sav_rate]*savr + _b[inc_a]* ina + _b[age]* ag + _b[fam]*siz+ _b[HighEd1]*4 + _b[price_ratio]* pri + _b[year]*2009 + _b[_cons]
scalar pe = logistic(xb0 -_b[HighEd1]*educ + _b[HighEd1]*(educ+4)) - logistic(xb0)
display "Discrete change in predicted value = " pe - xb0

2) As far as I have understood, your code is alternative to mine to the extent that .044892 is the predicted probability at price=6000 and weight=3000, while .0024804 is the predicted probability at price=6000 and weight=3500

Yes.

(.0024804- .044892) returns the decrease in probability from an increase in
weight=3500-3000=+ 500. Is the this point correct?

Yes. What I showed used some built-in features of -margins- to do what is done in a hard and primitive way with the scalar computation. When you check out the documentation for -margins-, you'll also find that it has ways for you to ask for the median (or other percentile values) so you don't need to pre-compute them.

Although I am not sure I have understood how you would compute the robust standard errors from results in mat list r(V)?

As I said, I don't know for sure whether those standard errors are robust or not, but I do know how to calculate the standard error of the difference of those two predicted values. I didn't describe that because I had assumed you had more knowledge about math stat because you were talking about pretty "heavy-duty" stuff, Delta-Method estimators and the like. It sounds like what you might be missing here is that there is a basic theorem in statistics stating that the variance of the difference of two random variables, say X and Y, is:

var(X - Y) = var(X) + var(Y) - 2*cov(X,Y)

After -margins-, r(V) holds the variance/covariance matrix of your predicted values, which -margins- leaves behind. The diagonal elements of r(V) give the variance of each estimate, and the off-diagonal elements give the covariance. The notation is easier if you assign r(V) to an actual Stata matrix, so you could do this:

logistic....
margins ... ... (with two different values of educ)
mat V = r(V)
scalar vardiff = V[1,1]^2 + V[2,2]^2  - 2*V[1,2]
di "SE of diff of predicted values = " sqrt(vardiff)

It sounds like you're in over your head somewhat here, both with regard to Stata and to math/stats (no shame in that, by the way), so I'd encourage you to pick up some of the background you need to go further.

[Linda Luciani]

Dear Professor Mike Lacy ,

thank you again for your kind reply!

For that reason, I assumed that those values you had were variables, not scalar values.

Yes, you are right about that. I have simply called the scalar values of the medians after those names: ina, siz etc. . Sorry about that and thanks for guidance about using scalar rather than generate function.

scalar xbo = _b[sav_rate]*savr + _b[inc_a]* ina + _b[age]* ag + _b[fam]*siz+ _b[HighEd1]*4 + _b[price_ratio]* pri + _b[year]*2009 + _b[_cons] scalar pe = logistic(xb0 -_b[HighEd1]*educ + _b[HighEd1]*(educ+4)) - logistic(xb0) display "Discrete change in predicted value = " pe - xb0

This seems correct to me as well but I think that the discrete change is simply: "Discrete change in predicted value = " pe , therefore the "-xb0" should not be there. Indeed doing so I obtain the same result as the difference of the two beta obtained after the following command you suggested (here below):

margins, at(sav_rate= 3.522739 price_ratio= 1.865185 age= 49 inc_a = 27616.16 fam= 2 year=2009 HighEd1=(4 8)) vce(unconditional)

Your Code:

 logistic.... margins ... ... (with two different values of educ) mat V = r(V) scalar vardiff = V[1,1]^2 + V[2,2]^2  - 2*V[1,2] di "SE of diff of predicted values = " sqrt(vardiff)

I think there is a typo in the code above since in line #4 since it already extracts variances from the var-cov matrix and therefore that should be rather be without the squares i.e. scalar vardiff = V[1,1] + V[2,2] - 2*V[1,2]

Thank you for this and for the relating explanation, the reason I have I did ask is the following: When I run the logit regression: logit repair sav_rate inc_a age HighEd1 fam price_ratio c.year [pweight=weighta] the coefficient for sav_rate is highly insignificant while when I run its -margins-, SE (the way you suggested) and then the p-value, I get p-value=0.0000 and therefore high significance. This seems very weird to me given non-significance in the logit and makes me wonder if this way to compute the SE is correct.

Here I report the steps I took for SE and p-value for significance:

*SAV_RATE
*for an increase in sav_rate equal to +1% ; 3.522739*1.01=3.5579664
margins, at( sav_rate=( 3.522739 3.5579664 ) price_ratio= 1.865185 age= 49 inc_a = 27616.16 fam= 2 year=2009 HighEd1=4) vce(unconditional)

mat list r(b) // predicted values
*difference=.06423018 - .06423016= 2.000e-08
mat list r(V) // variance matrix for predictions
mat V = r(V)
scalar vardiff = V[1,1] + V[2,2]  - 2*V[1,2]
di vardiff
di sqrt(vardiff)
*sqrt(vardiff)=4.543e-08

*t-value  =  (x-μ) / (s/√n)
di (2.000e-08 - 0)/(4.543e-08/sqrt(9081))
*t-value=41.952138 ; cross-chacking with t-values table I get a p-value=0.0000

I am also a bit skeptical since a professor suggested me to go with the Delta Method and I have found this source (comment #4) in which they compute the SE for the Average Partial Effects (APEs) of a logit model. The way they do it is not only hard for me to understand but I also wonder if the further steps they take in there are due to the fact that SE of APE is much more complex to compute than the SE of a Discrete Increase at the Medians (like it is my case), since in the former derivates have to be computed for every observed point of the distribution.

Again, thank you very much for helping me!

[Mike Lacy]

Whoops, I got so wrapped up in various things that I did make some mistakes here, I'm sorry. First, here's a correct example.

sysuse auto, clear
logit foreign price weight
// predicted probability at price = 6000 and weight = 3000
scalar xb3000 = _b[_cons] + _b[price]*6000 + _b[weight]*3000  // predicted logit
scalar p3000 = invlogit(xb3000) // probability
di "p3000 = " p3000
//
// predicted probability at price = 6000 and weight = 3500
scalar xb3500 = _b[_cons] + _b[price]*6000 + _b[weight]*3500  // predicted logit
scalar p3500 = invlogit(xb3500)  // probability
di "p3500 = " p3500
// Same thing from margins.
margins, at(price = (6000) weight = (3000 3500))

.... I think there is a typo in the code above since in line #4 since it already extracts variances from the var-cov matrix

Yes, of course; my apologies.

...the coefficient for sav_rate is highly insignificant while when I run its -margins-. etc.
t-value = (x-µ) / (s/vn)
di (2.000e-08 - 0)/(4.543e-08/sqrt(9081))

This is distinctly the wrong formula: You treated the *sampling variance* of the difference in predicted values as though it was a variance of a variable for which a mean had been calculated, and therefore divided by sqrt(N) as you would in calculating the standard error of the mean. Simple means and related formulae are not at all relevant here. The t-statistic should be calculated as: 2.000e-08/4.54e-08 = 0.44

Cautions: 1) I would not be certain that this quantity nicely follows a t-distribution for the purpose of producing confidence intervals or p-values. 2) Numbers like 2e-8 are at the edges of the precision of a calculation in Stata, and might give you wrong results when you do further calculations with them.

I can't comment on the other posting you cite, as digging into that is beyond the time I want to put in here, and may also be beyond my competence level.

One last comment: Rather than putting so much attention onto "significance" and the like, I'd pay more attention to "Is the difference big enough to care about," and only secondarily worry about "how precisely have I measured that difference."