1
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ Why not use margins after quantile regression for this? $\endgroup$
    – dimitriy
    Commented Mar 17, 2020 at 5:22
  • $\begingroup$ 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. $\endgroup$ Commented Mar 17, 2020 at 15:46
  • $\begingroup$ I think I misunderstood what you had in mind, so please disregard my QR suggestion. $\endgroup$
    – dimitriy
    Commented Mar 18, 2020 at 6:41

1 Answer 1

0
$\begingroup$

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
$\endgroup$
1
  • $\begingroup$ This is extremely helpful, thanks Dimitriy! $\endgroup$ Commented Mar 18, 2020 at 20:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.