注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

凯恩斯的博客

做一个无畏的行者

 
 
 

日志

 
 
关于我

Permanent Head Damage

网易考拉推荐

STATA-probit模型中如何做marginal effect  

2011-12-18 11:17:05|  分类: 默认分类 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

Note: This FAQ is for Stata 10 and older versions of Stata. In Stata 11, the margins command replaced mfx.

I am using mfx after an estimation that has an offset. How does mfx take that into account?

Title Marginal effects after estimations with offsets
AuthorMay Boggess, StataCorp
DateApril 2004; minor revisions October 2005

The command mfx evaluates at the mean value of the offset. Let’s see how this works in an example:

 . sysuse auto, clear  (1978 Automobile Data)     . probit foreign weight mpg, offset(turn) nolog    Probit estimates                                  Number of obs   =         74                                                    Wald chi2(2)    =     440.82  Log likelihood = -130.53661                       Prob > chi2     =     0.0000    ------------------------------------------------------------------------------       foreign |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]  -------------+----------------------------------------------------------------        weight |   -.008765   .0005981   -14.65   0.000    -.0099373   -.0075928           mpg |  -.1921153   .0527769    -3.64   0.000    -.2955561   -.0886746         _cons |  -10.73888   2.721338    -3.95   0.000     -16.0726   -5.405153          turn |   (offset)  ------------------------------------------------------------------------------  Note: 8 failures and 0 successes completely determined.    . mfx, predict(p)     Marginal effects after probit        y  = Pr(foreign) (predict, p)           =  .04973232  ------------------------------------------------------------------------------  variable |      dy/dx    Std. Err.     z    P>|z|  [    95% C.I.   ]      X  ---------+--------------------------------------------------------------------    weight |  -.0009001      .00031   -2.89   0.004  -.001511  -.00029   3019.46       mpg |  -.0197293      .00802   -2.46   0.014   -.03544 -.004019   21.2973      turn |  (offset)                                                   39.6486  ------------------------------------------------------------------------------    . summarize turn        Variable |       Obs        Mean    Std. Dev.       Min        Max  -------------+--------------------------------------------------------          turn |        74    39.64865    4.399354         31         51    . replace turn=r(mean)  turn was int now float  (74 real changes made)    . quietly summarize weight    . quietly replace weight=r(mean)    . quietly summarize mpg    . quietly replace mpg=r(mean)    . predict y, p   . list y in 1         +----------+       |        y |       |----------|    1. | .0497323 |       +----------+ 

In the above example, we demonstrated that the value y (above the table in the mfx output), is the value you obtain when you predict at the average values of the covariates and the average value of the offset. In the next example, we will go further and calculate the marginal effects of two dichotomous variables by hand:

 . clear   . set obs 50  obs was 0, now 50    . set seed 85642    . generate y=uniform()>0.5    . generate x1=uniform()>0.5    . generate x2=uniform()>0.5    . generate off=uniform()*30+50    . probit y x1 x2, offset(off) nolog    Probit estimates                                  Number of obs   =         50                                                    Wald chi2(2)    =      58.13  Log likelihood = -286.69439                       Prob > chi2     =     0.0000    ------------------------------------------------------------------------------             y |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]  -------------+----------------------------------------------------------------            x1 |   3.081539   .4533974     6.80   0.000     2.192896    3.970181            x2 |  -1.103905   .4571051    -2.41   0.016    -1.999815   -.2079954         _cons |  -68.40088   .4041975  -169.23   0.000    -69.19309   -67.60867           off |   (offset)  ------------------------------------------------------------------------------  Note: 10 failures and 3 successes completely determined.   . mfx, predict(p)     Marginal effects after probit        y  = Pr(y) (predict, p)           =  .08459144  ------------------------------------------------------------------------------  variable |      dy/dx    Std. Err.     z    P>|z|  [    95% C.I.   ]      X  ---------+--------------------------------------------------------------------        x1*|   .6777208      .11014    6.15   0.000   .461848  .893594        .4        x2*|  -.1832137      .09128   -2.01   0.045  -.362116 -.004311       .52       off |  (offset)                                                   66.3675  ------------------------------------------------------------------------------  (*) dy/dx is for discrete change of dummy variable from 0 to 1    . summarize off        Variable |       Obs        Mean    Std. Dev.       Min        Max  -------------+--------------------------------------------------------           off |        50    66.36746    8.437674   50.19214   79.75008    . replace off=r(mean)  (50 real changes made)    . quietly {  . summarize x1  . local mean1=r(mean)  . summarize x2  . local mean2=r(mean)  . replace x1=`mean1'  . replace x2=`mean2'  . predict p, p  . noisily display _n "y = " p _n      y = .08459145    . replace x1=0  . predict p0,p  . replace x1=1  . predict p1, p  . noisily display  _n "marginal effect x1 = " p1-p0  _n     marginal effect x1 = .67772082    . drop p1 p0  . replace x1=`mean1'  . replace x2=0  . predict p0,p  . replace x2=1  . predict p1, p  . noisily display  _n "marginal effect x2 = " p1-p0      marginal effect x2 = -.18321373  . } 

The above example gives you the idea of what to do if you want to evaluate marginal effects at a value of the offset that is not the mean. You can replace the offset by the value you want before running mfx, as follows:

 . sysuse auto, clear  (1978 Automobile Data)    . replace weight=weight/1000  weight was int now float  (74 real changes made)    . probit foreign weight length, offset(turn) nolog    Probit estimates                                  Number of obs   =         74                                                    Wald chi2(2)    =     511.38  Log likelihood = -129.03886                       Prob > chi2     =     0.0000    ------------------------------------------------------------------------------       foreign |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]  -------------+----------------------------------------------------------------        weight |  -3.277497   .8475066    -3.87   0.000    -4.938579   -1.616414        length |  -.1330356   .0294001    -4.52   0.000    -.1906588   -.0754124         _cons |  -5.978285   3.276507    -1.82   0.068    -12.40012    .4435497          turn |   (offset)  ------------------------------------------------------------------------------  Note: 1 failure and 0 successes completely determined.    . mfx, predict(p)     Marginal effects after probit        y  = Pr(foreign) (predict, p)           =  .10979816  ------------------------------------------------------------------------------  variable |      dy/dx    Std. Err.     z    P>|z|  [    95% C.I.   ]      X  ---------+--------------------------------------------------------------------    weight |  -.6154742      .19077   -3.23   0.001  -.989375 -.241573   3.01946    length |  -.0249825      .00869   -2.87   0.004  -.042023 -.007942   187.932      turn |  (offset)                                                   39.6486  ------------------------------------------------------------------------------    . replace turn=41  (70 real changes made)    . mfx, predict(p)     Marginal effects after probit        y  = Pr(foreign) (predict, p)           =  .54924357  ------------------------------------------------------------------------------  variable |      dy/dx    Std. Err.     z    P>|z|  [    95% C.I.   ]      X  ---------+--------------------------------------------------------------------    weight |  -1.297558      .34523   -3.76   0.000   -1.9742 -.620918   3.01946    length |  -.0526687      .01159   -4.54   0.000  -.075384 -.029954   187.932      turn |  (offset)                                                        41  ------------------------------------------------------------------------------ 

If you have used exposure() rather than offset(), then mfx will evaluate at the exp(average of ln(exposure)). For example,

 . webuse airline, clear     . poisson injur XYZ, exposure(n) nolog    Poisson regression                                Number of obs   =          9                                                    LR chi2(1)      =       1.77                                                    Prob > chi2     =     0.1836  Log likelihood = -23.027177                       Pseudo R2       =     0.0370    ------------------------------------------------------------------------------      injuries |      Coef.   Std. Err.      z    p>|z|     [95% Conf. Interval]  -------------+----------------------------------------------------------------      XYZowned |   .3808084   .2780192     1.37   0.171    -.1640993    .9257161         _cons |   4.061204    .147442    27.54   0.000     3.772223    4.350185             n | (exposure)  ------------------------------------------------------------------------------    . mfx, predict(n)     Marginal effects after poisson        y  = predicted number of events (predict, n)           =  6.4864947  ------------------------------------------------------------------------------  variable |      dy/dx    Std. Err.     z    p>|z|  [    95% C.I.   ]      X ---------+--------------------------------------------------------------------  XYZowned*|   2.647899     2.14322    1.24   0.217  -1.55274  6.84854   .333333     ln(n) |  (offset)                                                  -2.31842  ------------------------------------------------------------------------------  (*) dy/dx is for discrete change of dummy variable from 0 to 1    . generate offset=ln(n)    . summarize offset        Variable |       Obs        Mean    Std. Dev.       Min        Max  -------------+--------------------------------------------------------        offset |         9   -2.318418    .5339422   -2.98975  -1.571179    . replace n=exp(r(mean))  (9 real changes made)    . replace XYZ=0  (3 real changes made)    . predict n0, n    . replace XYZ=1  (9 real changes made)    . predict n1, n    . display _n "marginal effect XYZ = " n1-n0    marginal effect XYZ = 2.6478992 

Some multiple-equation estimators allow more than one offset. Let's look at an example using hetprob. First, we will create some data and run the model:

 . clear     . set obs 1000  obs was 0, now 1000    . set seed 1234567    . generate x = 1-2*uniform()    . generate xhet=uniform()    . generate sigma=exp(1.5*xhet)    . generate p=normal((0.3+2*x)/sigma)    . generate y = cond(uniform()<=p,1,0)    . generate off1=uniform()+0.5    . generate off2=uniform()+0.6    . hetprob y x, het(xhet, offset(off2))  offset(off1) nolog    Heteroskedastic probit model                    Number of obs     =       1000                                                  Zero outcomes     =        452                                                  Nonzero outcomes  =        548                                                    Wald chi2(1)      =      79.18  Log likelihood = -576.5979                      Prob > chi2       =     0.0000    ------------------------------------------------------------------------------             y |      Coef.   Std. Err.      z    p>|z|     [95% Conf. Interval]  -------------+----------------------------------------------------------------  y            |             x |   6.658445   .7482918     8.90   0.000      5.19182     8.12507         _cons |  -.3720305   .2473793    -1.50   0.133     -.856885    .1128241          off1 |   (offset)  -------------+----------------------------------------------------------------  lnsigma2     |          xhet |   1.724798   .2680221     6.44   0.000     1.199484    2.250112          off2 |   (offset)  ------------------------------------------------------------------------------  Likelihood-ratio test of lnsigma2=0: chi2(1) =   101.70   Prob > chi2 = 0.0000 

Now we will run mfx with the predict(xb) option, which is the linear predictor from the first equation:

 . mfx, predict(xb)     Marginal effects after hetprob        y  = Linear prediction (predict, xb)           =  .76705778  ------------------------------------------------------------------------------  variable |      dy/dx    Std. Err.     z    p>|z|  [    95% C.I.   ]      X  ---------+--------------------------------------------------------------------         x |   6.658445      .74829    8.90   0.000   5.19182  8.12507   .020114      xhet |          0           0       .       .         0        0   .502716      off1 |  (offset1)                                                  1.00516      off2 |  (offset2)                                                  1.09709  ------------------------------------------------------------------------------    . summarize x if e(sample)        Variable |       Obs        Mean    Std. Dev.       Min        Max  -------------+--------------------------------------------------------             x |      1000    .0201138    .5846345  -.9999619   .9977288    . scalar meanx=r(mean)    . summarize off1  if e(sample)        Variable |       Obs        Mean    Std. Dev.       Min        Max  -------------+--------------------------------------------------------          off1 |      1000    1.005162    .2875233   .5014016   1.497421    . scalar meanoff1=r(mean)    . display _n "pred xb = " meanx*_b[x]+_b[_cons] + meanoff1    pred xb = .76705774 

Let’s finish with a more complicated prediction option. By differentiating the formula for the probability of success, we can verify that mfx is correctly calculating the marginal effect. The formula for the probability of success is in [R] hetprob:

                  meanx*_b[x]+_b[_cons] + meanoff1   p = normal( ------------------------------------------  )               exp(meanxhet*_b[lnsigma2:xhet] + meanoff2) 

where normal() is the cumulative distribution function for the standard normal distribution. The derivative is

 dp                  _b[x]   -- = ------------------------------------------     dx   exp(meanxhet*_b[lnsigma2:xhet] + meanoff2)                             meanx*_b[x]+_b[_cons] + meanoff1         *   normalden( ------------------------------------------  )                        exp(meanxhet*_b[lnsigma2:xhet] + meanoff2) 

where normalden() is the probability density function for the standard normal distribution. Let’s use these formulas to check mfx:

 . mfx, predict(p)      Marginal effects after hetprob        y  = Pr(y) (predict, p)           =  .54284206  ------------------------------------------------------------------------------  variable |      dy/dx    Std. Err.     z    p>|z|  [    95% C.I.   ]      X  ---------+--------------------------------------------------------------------         x |   .3704576      .03237   11.44   0.000   .307015    .4339   .020114      xhet |  -.0736092      .02423   -3.04   0.002  -.121095 -.026124   .502716      off1 |  (offset1)                                                  1.00516      off2 |  (offset2)                                                  1.09709  ------------------------------------------------------------------------------    . summarize xhet if e(sample)         Variable |       Obs        Mean    Std. Dev.       Min        Max  -------------+--------------------------------------------------------          xhet |      1000    .5027165    .2938006   .0004071   .9983656    . scalar meanxhet=r(mean)    . summarize off2  if e(sample)         Variable |       Obs        Mean    Std. Dev.       Min        Max  -------------+--------------------------------------------------------          off2 |      1000    1.097091    .2822951   .6011194   1.597888    . scalar meanoff2=r(mean)    . scalar predy=normal((meanx*_b[x]+_b[_cons] +   > meanoff1)/exp(_b[lnsigma2:xhet]*meanxhet +meanoff2))    . display _n "pred y = " predy     pred y = .54284206    . scalar dpdx=_b[x]*normalden((meanx*_b[x]+_b[_cons] +   > meanoff1)/exp(_b[lnsigma2:xhet]*meanx  > het  +meanoff2))/exp(_b[lnsigma2:xhet]*meanxhet +meanoff2)    . display _n "dpdx = " dpdx     dpdx = .37045757
  评论这张
 
阅读(2189)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017