# Median of Predictions With Standard Error using Margins command in Stata

I understand that by default, the "margins" command in Stata calculates the predicted value of the dependent variable for each observation, then reports the mean value of the predictions. Is it possible for "margins" to calculate the median value of the predictions instead? I do know how to calculate this manually, but I need the standard errors using the delta method that is provided by the "margins" command. Alternatively, if anyone knows how to manually compute the same standard errors as "margins" for the median of the predictions, I could implement that also.

Thanks for your help.

• Why not use margins after quantile regression for this? Commented Mar 17, 2020 at 5:22
• Thanks for the suggestion. I need to estimate the predictions from the user-written "twopm" command where the first part is a logistic regression and the second is a GMM with gamma distribution and log-link. If I understand correctly, this differs from what I would get from a quantile regression. Commented Mar 17, 2020 at 15:46
• I think I misunderstood what you had in mind, so please disregard my QR suggestion. Commented Mar 18, 2020 at 6:41

## 1 Answer

I wonder if you can just bootstrap the median prediction. I don't have my econometrics books with me and can't look this up because of pangolins, but the basic idea is below. If this is a stupid one, perhaps others can chime in on why that is.

Here we will fit the TPM, use margins to get the average prediction, and calculate the median of the predictions by hand. Then we will both bootstrap the mean and the median to see if the mean's SE agrees with the delta method approach. This might give us more confidence in that this works with the median.

. webuse womenwk, clear

. replace wage = 0 if wage==.
(657 real changes made)

. twopm wage educ age married children, firstpart(logit) secondpart(glm) nolog

Fitting logit regression for first part:

Iteration 0:   log likelihood = -1266.2225
Iteration 1:   log likelihood = -1040.6658
Iteration 2:   log likelihood = -1027.9567
Iteration 3:   log likelihood = -1027.9145
Iteration 4:   log likelihood = -1027.9144

Fitting glm regression for second part:

Iteration 0:   log likelihood = -4159.7605

Two-part model
------------------------------------------------------------------------------
Log pseudolikelihood =  -5187.675                 Number of obs   =       2000

Part 1: logit
------------------------------------------------------------------------------
Number of obs   =       2000
LR chi2(4)      =     476.62
Prob > chi2     =     0.0000
Log likelihood = -1027.9144                       Pseudo R2       =     0.1882

Part 2: glm
------------------------------------------------------------------------------
Number of obs   =      1343
Deviance         =  38542.35905                    (1/df) Deviance =  28.80595
Pearson          =  38542.35905                    (1/df) Pearson  =  28.80595

Variance function: V(u) = 1                        [Gaussian]
Link function    : g(u) = u                        [Identity]

AIC             =  6.202175
Log likelihood   = -4159.760537                    BIC             =   28905.2
------------------------------------------------------------------------------
wage |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
logit        |
education |   .0982513   .0186522     5.27   0.000     .0616936     .134809
age |   .0579303    .007221     8.02   0.000     .0437773    .0720833
married |   .7417775   .1264705     5.87   0.000     .4938998    .9896552
children |   .7644882   .0515289    14.84   0.000     .6634935     .865483
_cons |  -4.159247   .3320401   -12.53   0.000    -4.810034   -3.508461
-------------+----------------------------------------------------------------
glm          |
education |   .8750694    .050243    17.42   0.000     .7765949    .9735438
age |   .1514818   .0192717     7.86   0.000     .1137099    .1892537
married |  -.5395024   .3574519    -1.51   0.131    -1.240095    .1610904
children |  -.6862982   .1032256    -6.65   0.000    -.8886166   -.4839798
_cons |   7.934369   .9264515     8.56   0.000     6.118558    9.750181
------------------------------------------------------------------------------

. margins

Predictive margins                              Number of obs     =      2,000

Expression   : twopm combined expected values, predict()

------------------------------------------------------------------------------
|            Delta-method
|     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
_cons |   15.90929   .2388088    66.62   0.000     15.44123    16.37735
------------------------------------------------------------------------------

. predict yhat

. qui sum yhat, detail

. di "Median yhat is " r(p50)
Median yhat is 16.573483

.
. capture program drop my_margins

.
. program define my_margins, eclass
1.         twopm wage educ age married children, firstpart(logit) secondpart(glm)
2.         tempvar yhat
3.         predict yhat'
4.         sum yhat', detail
5.         ereturn scalar mean_yhat   = r(mean)
6.         ereturn scalar median_yhat = r(p50)
7. end

.
. bootstrap median_yhat=e(median_yhat) mean_yhat=e(mean_yhat), reps(100) seed(100179) nodots: my_margins

Bootstrap results                               Number of obs     =      2,000
Replications      =        100

command:  my_margins
median_yhat:  e(median_yhat)
mean_yhat:  e(mean_yhat)

------------------------------------------------------------------------------
|   Observed   Bootstrap                         Normal-based
|      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
median_yhat |   16.57348   .3778301    43.86   0.000     15.83295    17.31402
mean_yhat |   15.90929   .2695621    59.02   0.000     15.38096    16.43762
------------------------------------------------------------------------------


The mean and median match the output above. The mean SE is close, but does not quite match the delta method from margins, though to be fair, I did not expect it to match exactly.

Here's the code for those playing at home:

cls
webuse womenwk, clear
replace wage = 0 if wage==.
twopm wage educ age married children, firstpart(logit) secondpart(glm) nolog
margins
predict yhat
qui sum yhat, detail
di "Median yhat is " r(p50)

capture program drop my_margins

program define my_margins, eclass
twopm wage educ age married children, firstpart(logit) secondpart(glm)
tempvar yhat
predict yhat'
sum yhat', detail
ereturn scalar mean_yhat   = r(mean)
ereturn scalar median_yhat = r(p50)
end

bootstrap median_yhat=e(median_yhat) mean_yhat=e(mean_yhat), reps(100) seed(100179) nodots: my_margins

• This is extremely helpful, thanks Dimitriy! Commented Mar 18, 2020 at 20:19