capture program drop margeff program define margeff *! Obtain partial effects after estimation *! Version 2.2.0 (20 August 2009) (Revision of Stata Journal submission) *! Author: Tamas Bartus (Corvinus University, Budapest) version 9.2 ******************************************************************************* * * [1] DETERMINING WHAT TO DO * ******************************************************************************* if "`1'"=="?" { gettoken tmp 0 : 0 if `"`0'"'=="" SendVersion else `0' exit } if "`e(cmd)'"=="" exit 301 // Determining what to do // if no arguments specified: replay results, provided the -replace- option was used if `"`0'"'=="," | `"`0'"'=="" { if "`e(cmd)'"=="margeff" local doit = 0 else local doit = 1 } // if something specified: estimate, provided replay subcmd not found else { gettoken cmd 0 : 0 , parse(" ,") local l = length("`cmd'") if substr("replay",1,`l')==`"`cmd'"' { if "`e(cmd)'"=="margeff" & "`e(margeff_cmd)'"=="" { di as red "It is impossible to replay marginal effects" error 301 } local doit = 0 } else if substr("compute",1,`l')=="`cmd'" { local doit = 1 } else { local 0 `cmd' `0' local doit = 1 } } if `doit'==1 { preserve CheckSupport `0' Estimate `0' Replay `0' restore sret clear } else { Replay `0' } end program define CheckSupport , sclass version 8.2 syntax [if] [in] [fweight pweight iweight /] [ , Model(string) Link(string) eform * ] sret local model if "`e(cmd)'"=="margeff" exit 301 local cmd = e(cmd) // Taking care of xt commands if substr(`"`cmd'"',1,2)=="xt" { local family xt local cmd = substr("`cmd'",3,length("`cmd'")-2) } else local family if ("`cmd'"=="glm" | "`family'`cmd'"=="xtgee") & "`link'"=="" & "`model'"=="" { if "`family'"=="xt" local elink = e(link) else local elink = e(linkt) local elink = lower("`elink'") if "`elink'"=="logit" local model logit else if "`elink'"=="probit" local model probit else if "`elink'"=="complementary log-log" local model cloglog else if "`elink'"=="cloglog" local model cloglog else if "`elink'"=="log" local link log else if "`elink'"=="log-log" local link loglog else if "`elink'"=="log complement" local link logc else if "`elink'"=="neg. binomial" local link nbinomial else if "`elink'"=="nbinomial" local link nbinomial else if "`elink'"=="power" local link power `e(power)' else if "`elink'"=="odds power" local link opower `e(power)' else if "`elink'"=="odds p" local link opower `e(power)' } if "`eform'"!="" local link log if `"`link'"'!="" { local allow_link log loglog logc nbinomial power opower gettoken model power : link local temp : subinstr local allow_link "`model'" "" , count(local c) word if `c'==0 { di as error "margeff does not support the link -`model'-" error 198 } local subcmd link } else { if `"`model'"'=="" local model `cmd' local allow_normal probit logit logistic cloglog probit poisson nbreg local allow_ordered oprobit ologit local allow_multieq gologit2 mlogit heckprob biprobit zip zinb local allow_select heckman tobit truncreg cnreg foreach word in normal ordered multieq select { local temp : subinstr local allow_`word' "`model'" "" , count(local c) word local subcmd `word' if `c'==1 { continue , break } } if `c'==0 { di as error "margeff does not support `family'`model'; please use -mfx- instead" error 198 } } sret local model `model' sret local family `family' sret local link `link' sret local subcmd `subcmd' end ******************************************************************************* * * [2] ESTIMATION + REPLAYING * ******************************************************************************* program define Estimate , eclass version 8.2 syntax [if] [in] [fweight pweight iweight /] /// [ , at(string) Model(string) Link(string) /// Dummies(string) noOFFset noWght CONStant /// OUTcome(string) Percent Replace /// eform dx(string) count * ] // CHECKING AND PROCESSING OPTIONS foreach w in model family link subcmd { local `w' `s(`w')' } if `e(df_m)'==0 { di as error "There are no independent variables; running -margeff- makes no sense" exit 301 } if "`model'"=="cnreg" & `"`at'"'!="" { di as error "Only average partial effects can be calculated after cnreg" di as error "Option -at()- will be ignored" local at } if `"`dummies'"'!="" { DumList `dummies' if "`s(error)'"!="" { di as error "`s(error)'" error 198 } } * Weights if "`wght'"=="nowght" | "`e(wtype)'"=="" { local weight } else { local weight [iweight`e(wexp)'] local weight : subinstr local weight "= " "=" } if "`e(offset)'"!="" & "`offset'"=="" { tempvar offsetvar qui gen double `offsetvar' = `e(offset)' local offset offset(`offsetvar') } else local offset // Setting key S macros SetSMacros_`subcmd' `model' `outcome'`power' if "`outcome'"!="" { CheckOut `outcome' } else { numlist "1/`s(nout)'" } local outcome = r(numlist) // MARKING SAMPLE, GETTING E RESULTS, TEMPNAMES tempvar touse mark mark `mark' `if' `in' `weight' qui gen byte `touse' = e(sample) if `mark'==1 qui count if `touse'==1 local numobs = r(N) local depvar = e(depvar) local cmd = e(cmd) tempname bold vold // estimation results tempname bnew vnew // marginal effects tempname atmat // vector of evaluation points tempname dxmat // matrix of changes tempname at_dx // matrix of at values and changes tempname matcell // temporary matrix - checking whether a variable is dummy or not // Variables for which marginal effects should be computed mat `bold' = e(b) mat `vold' = e(V) local ncol = colsof(`bold') local colnames : colnames `bold' local varlist : list uniq colnames local varlist : subinstr local varlist "_cons" "" , all word count(local hascons) local varlist : subinstr local varlist "_se" "" , all word // this should eliminate anc parameters as well local nvar : word count `varlist' if `"`at'"'!="" { local method at AtList , at(`at') varlist(`varlist') local atvarlist = r(atvarlist) local atnumlist = r(atnumlist) local atstat = r(atstat) if "`atstat'"=="p50" { if `"`weight'"'!="" & "`wtype'"=="pweight" { di as error "Warning: pweights will be ignored during the calculation of medians" local weight } } if "`offset'"!="" { if "`atstat'"=="mean" { qui su `offsetvar' if `touse' `weight' local mean = r(mean) } else if "`atstat'"=="p50" { qui su `offsetvar' if `touse' `weight' , detail local mean = r(p50) } else local mean = 0 local offset offset(`mean') } } else { local method av local atstat mean // if "`offset'"!="" local offset offset(`offsetvar') } // Constructing the matrices of values at which marginal effects be computed // Making dx & at specification complete mat `atmat' = J(1,`nvar',0) mat `dxmat' = J(2,`nvar',.) mat colnames `dxmat' = `varlist' mat colnames `atmat' = `varlist' tempname dx1 dx2 forval x = 1/`nvar' { local var : word `x' of `varlist' local value = 0 // Atvalues if `"`at'"'!="" { local pos : list posof `"`var'"' in atvarlist if `pos'!=0 { local value : word `pos' of `atnumlist' } else if "`atstat'"!="zero" { if "`atstat'"=="p50" qui su `var' if `touse' , detail else qui su `var' if `touse' `weight' local value = r(`atstat') } } // end of if `"`at'"'!="" .... // From and To values of discrete change capture assert `var'==0 | `var'==1 if `touse' if _rc==0 { // dummy found scalar `dx1' = 1-`value' scalar `dx2' = 0-`value' } else { CalcUnit `var' `touse' scalar `dx1' = r(unit) scalar `dx2' = -`dx1' } * Filling in elements of atmat and dxmat matrices mat `atmat'[1,`x'] = `value' mat `dxmat'[1,`x'] = `dx1' mat `dxmat'[2,`x'] = `dx2' } // end of forvalue loop // Getting estimation results + removing anciliary parameters // Specifying xmean matrices with the help of atmat ProcessEstMat_`s(type)' `ncol' `bold' `vold' `method' mat `vold' = r(v) forval x = 1/`s(neq)' { tempname bold`x' mat `bold`x'' = r(bold`x') local boldlist `boldlist' `bold`x'' // At values tempname xmean`x' MakeXMeanMat_`s(type)' `x' `bold`x'' `atmat' mat `xmean`x'' = r(xmean) local xmeanlist `xmeanlist' `xmean`x'' if "`method'"=="av" { tempvar xbvar`x' qui gen double `xbvar`x'' = . local xbvarlist `xbvarlist' `xbvar`x'' } } // Defining some locals needed during calculations local neq = `s(neq)' local exe `s(model)' if "`s(model)'"=="gologit2" { local exe o`e(link)' } if "`method'"=="av" { local tempcmd tempvar local gencmd qui gen double local gencond if `touse'==1 local xlist xlist(`xbvarlist') local margeff_type Average partial effects } else { local tempcmd tempname local gencmd scalar local gencond local xlist xlist(`xmeanlist') local margeff_type Partial effects at fixed values } if "`family'"=="xt" { if e(predict)=="xtbin_p" { local genfxopt = sqrt(1-`e(rho)') } else local genfxopt = 1 } else local genfxopt = 1 tempname xmat bnew qui replace `touse'=1 // Calculation of partial effects and standard errors foreach out in `outcome' { forval i = 1/`nvar' { local var : word `i' of `varlist' // Setting mean / vars to zero if dummy found + Taking care of the dummies() option local deltamat `dxmat' local xmatlist `xmeanlist' // Prediction after variable under treatment increased then decreased by predefined value forval r = 1/2 { `tempcmd' hat`r' local dx`r' = `deltamat'[`r',`i'] GenXB_`s(type)' `var' `touse' `weight' , /// method(`method') change(`dx`r'') blist(`boldlist') `xlist' outcome(`out') `dlist' `offset' GenFx_`exe' `out' `genfxopt' `gencmd' `hat`r'' = `s(cdf)' `gencond' forval eq = 1/`neq' { `tempcmd' phi`r'`eq' `gencmd' `phi`r'`eq'' = `s(dens`eq')' `gencond' } } // Partial effect tempname me CalcME_B_`method' `hat1' `hat2' `dx1' `dx2' `touse' `weight' mat `me' = r(me) mat colnames `me' = `var' EqLabel `out' `touse' if `: word count `outcome''>1 { mat coleq `me' = "`s(lab)'" } else local category `s(cat)' // needed for labeling output mat `bnew' = nullmat(`bnew') , `me' // Derivatives of partial effect with respect to coefficients tempname pder add vec forval eq = 1/`neq' { local xmat : word `eq' of `xmatlist' CalcME_SE_`method' `var' `phi1`eq'' `phi2`eq'' `dx1' `dx2' `xmat' `touse' `weight' mat `vec' = r(vec) if "`s(type)'"=="ordered" { if `eq'==1 mat `pder' = `vec' else mat `pder' = `pder'+`vec' } else mat `pder' = nullmat(`pder') , `vec' } mat `vnew' = nullmat(`vnew') \ `pder' } // end of variable loop if `hascons'==1 & "`constant'"!="" { tempname consder GenXB_`s(type)' `touse' `weight' , cons method(`method') blist(`boldlist') `xlist' outcome(`out') `offset' GenFx_`exe' `out' `genfxopt' scalar `hat' = `s(cdf)' CalcME_B_`method' `hat' 0 1 0 `touse' `weight' mat `me' = r(me) mat colnames `me' = _cons if `: word count `outcome''>1 { EqLabel `out' `touse' mat coleq `me' = "`s(lab)'" } mat `bnew' = `bnew' , `me' forval eq = 1/`neq' { tempname phi`eq' scalar `phi`eq'' = `s(dens`eq')' local xmat : word `eq' of `xmatlist' mat `vec' = `xmat'*0 local pos = colsof(`xmat') mat `vec'[1,`pos'] = `phi`eq'' if "`s(type)'"=="ordered" { if `eq'==1 mat `consder' = `vec' else mat `consder' = `consder'+`vec' } else mat `consder' = nullmat(`consder') , `vec' } mat `vnew' = nullmat(`vnew') \ `consder' } } // end of outcome loop // matrix of standard errors * mat coleq `vold' = _ * mat roweq `vold' = _ mat `vnew' = (`vnew'*`vold')*`vnew'' local rownames : colnames `bnew' local roweq : coleq `bnew' , quoted mat rownames `vnew' = `rownames' mat colnames `vnew' = `rownames' mat roweq `vnew' = `roweq' mat coleq `vnew' = `roweq' mat rownames `dxmat' = plus minus // Returning results local title `s(title)' // if "`category'"!="" & `s(catdepv)'==1 local title `title' == `category' if "`percent'"!="" { mat `bnew' = 10^2 * `bnew' mat `vnew' = 10^4 * `vnew' } if "`replace'"!="" { tempname bpost vpost * Saves scalars local scalars : e(scalars) foreach el in `scalars' { tempname s`el' scalar `s`el'' = e(`el') } mat `bpost' = `bnew' mat `vpost' = `vnew' ereturn post `bpost' `vpost' , depname("variable") obs(`numobs') esample(`touse') * Restores scalars foreach el in `scalars' { eret scalar `el' = `s`el'' } ereturn local cmd margeff } ereturn mat margeff_dx `dxmat' ereturn mat margeff_at `atmat' ereturn mat margeff_V `vnew' ereturn mat margeff_b `bnew' * ereturn local margeff_def `title' ereturn local margeff_depname `s(title)' ereturn local margeff_cmd `e(cmd)' ereturn local margeff_title `margeff_type' end program define Replay , eclass version 8 syntax [if] [in] [fweight pweight iweight /] [, Level(integer `c(level)') * ] di di as text "`e(margeff_title)' after " as res "`e(margeff_cmd)'" di as text _col(7) as text "y = " as res "`e(margeff_depname)' " di tempname ehold b v mat `b' = e(margeff_b) mat `v' = e(margeff_V) nobreak { _est hold `ehold' ereturn post `b' `v' , depname("variable") ereturn disp , level(`level') } _est unhold `ehold' end *============================================================ * * [3] SETTING KEY S MACROS --- SetSMacros_ subroutines * *============================================================ program define SetSMacros_normal , sclass // probit logit logistic cloglog xtprobit poisson nbreg version 8 local depvar = e(depvar) local depvar : word 1 of `depvar' if "`1'"=="poisson" | "`1'"=="nbreg" { sret local title E(`depvar') (expected number of counts) sret local catdepv = 0 } else { sret local title probability of `depvar' sret local title Pr(`depvar') sret local catdepv = 1 } SMacReturn normal `depvar' `1' 1 1 end program define SetSMacros_link , sclass version 8 args link power local depvar = e(depvar) local depvar : word 1 of `depvar' if "`power'"!="" { capture confirm number `power' if _rc!=0 { di as error "A number must be specified in the link(power #) option" error 198 } sret local power `power' } sret local catdepv = 0 local fnc = e(linkf) if `"`fnc'"'!="" local title : subinstr local fnc "u" "`depvar'" , all local title "`1'(`depvar')" sret local title "`title'" SMacReturn normal `depvar' `1' 1 1 end program define SetSMacros_ordered , sclass // oprobit ologit version 8 local depvar = e(depvar) local nout = e(k_cat) local neq = `nout'-1 if "`2'"!="" local depvar `depvar'=`2' sret local title Pr(`depvar') sret local catdepv = 1 SMacReturn ordered `depvar' `1' `nout' `neq' end program define SetSMacros_multieq , sclass // mlogit gologit2 heckprob biprobit zip zinb version 8 local depvar = e(depvar) gettoken depvar depv2 : depvar local nout = 1 local neq = 2 local type multieq // Modification of default settings if "`1'"=="heckprob" { local neq = 1 sret local title probability of `depvar' sret local title Pr(`depvar') sret local catdepv = 1 } else if "`1'"=="biprobit" { sret local title joint probabilities of `depvar' and `depv2' if "`2'"=="" sret local title Pr(`depvar',`depv2') else { local prob : word `2' of p00 p01 p10 p11 local v1 = substr("`prob'",2,1) local v2 = substr("`prob'",3,1) sret local title Pr(`depvar'=`v1',`depv2'=`v2') } if substr(e(title),1,1)=="B" { local nout = 4 local type multieq sret local title Pr(`depvar'=`v1',`depv2'=`v2') } else { local nout = 4 local type multieq sret local title Pr(`depvar'!=0,`depv2'!=0) } } else if "`1'"=="zip" | "`1'"=="zinb" { sret local title E(`depvar') (expected number of counts) sret local catdepv = 0 } else { if "`e(k_cat)'"!="" local nout = e(k_cat) else local nout = e(k_out) local neq = `nout'-1 if "`2'"!="" local depvar `depvar'=`2' sret local title probability of `depvar' sret local title Pr(`depvar') sret local catdepv = 1 } SMacReturn `type' `depvar' `1' `nout' `neq' end program define SetSMacros_select , sclass // heckman tobit cnreg intreg truncreg version 8 local depvar = e(depvar) local depvar : word 1 of `depvar' sret local title E(`depvar'|`depvar' observed) if "`1'"=="heckman" | "`1'"=="heckprob" { SMacReturn multieq `depvar' `1' 1 2 } else { SMacReturn normal `depvar' `1' 1 1 } end program define SMacReturn , sclass version 8 sret local type `1' sret local depvar `2' sret local model `3' sret local nout = `4' sret local neq = `5' end *============================================================ * * [4] VECTORS OF COEFFICIENTS AND AT VALUES FOR EACH EQUATION * ProcessEstMat_ subroutines * MakeXMeanMat_ subroutines * *============================================================ program define ProcessEstMat_normal , rclass version 8.2 args ncol b v at tempname bold1 local eqlist : coleq `b' , quoted local eqlist : list uniq eqlist if `"`eqlist'"'!="_" { // take first eq local eqname : word 1 of `eqlist' mat `b' = `b'[1,"`eqname':"] mat `v' = `v'["`eqname':","`eqname':"] } // Taking care of constant or noconstant option in est model local pos = colnumb(`b',"_cons") if `pos'!=. { mat `b' = `b'[1,1..`pos'] mat `v' = `v'[1..`pos',1..`pos'] } return mat bold1 = `b' return mat v = `v' end program define ProcessEstMat_multieq , rclass version 8.2 args ncol b v at local eqlist : coleq `b' , quoted local eqlist : list uniq eqlist tempname vce local nvar = 0 // counter for n of cols containing var and constant forval eq = 1/`s(neq)' { local eqname : word `eq' of `eqlist' tempname bold`eq' mat `bold`eq'' = `b'[1,"`eqname':"] local ncol = colsof(`bold`eq'') // Taking care of constant or noconstant option in est model local nvar = `nvar'+`ncol' return mat bold`eq' = `bold`eq'' } mat `vce' = `v'[1..`nvar',1..`nvar'] return mat v = `vce' end program define ProcessEstMat_ordered , rclass version 8 args ncol b v method tempname coef bcut local ncol = colsof(`b') local ncons = e(k_cat)-1 local nvar = e(df_m) mat `coef' = `b'[1,1..`nvar'] mat `bcut' = J(1,`ncons',0) mat coleq `coef' = _ mat colnames `bcut' = _cons forval eq = 1/`s(neq)' { tempname bold`eq' mat `bold`eq'' = `bcut' local pos = `nvar'+`eq' local cut = `b'[1,`pos'] if "`method'"=="av" local cut = -1*`cut' mat `bold`eq''[1,`eq'] = `cut' mat `bold`eq'' = `coef', `bold`eq'' return mat bold`eq' = `bold`eq'' } return mat v = `v' end program define ProcessEstMat_difficult , rclass version 8 di as error "Subroutine ProcessEstMat_difficult not implemented yet" error 198 end program define MakeXMeanMat_normal , rclass version 8.2 args eq b at tempname xmean // Taking care of constant or noconstant option in est model local pos = colnumb(`b',"_cons") if `pos'!=. { tempname one mat `one' = J(1,1,1) mat colnames `one' = _cons mat `xmean' = `at' , `one' } else mat `xmean' = `at' return add // preserves matrices created by ProcessEstMat program return mat xmean = `xmean' end program define MakeXMeanMat_multieq , rclass version 8 args eq b at local ncol = colsof(`b') local colnames : colnames `b' tempname xmean mat `xmean' = `b'*. mat colnames `xmean' = `colnames' forval col = 1/`ncol' { local name : word `col' of `colnames' if "`name'"=="_cons" local value = 1 else { local pos = colnumb(`at',"`name'") local value = `at'[1,`pos'] } mat `xmean'[1,`col'] = `value' } return add // preserves matrices created by ProcessEstMat program return mat xmean = `xmean' end program define MakeXMeanMat_ordered , rclass args eq b at local ncol = colsof(`b') local colnames : colnames `b' tempname xmean mat `xmean' = `b'*. mat colnames `xmean' = `colnames' forval col = 1/`ncol' { local name : word `col' of `colnames' if "`name'"=="_cons" { local cons = `b'[1,`col'] local value = cond(`cons'==0,0,-1) } else { local value = `at'[1,`col'] } mat `xmean'[1,`col'] = `value' } return add // preserves matrices created by ProcessEstMat program return mat xmean = `xmean' end program define MakeXMeanMat_difficult , rclass version 8 di as error "Subroutine MakeXMeanMat_difficult not implemented yet" error 198 end *============================================================ * * [5] LINEAR PREDICTION --- GenXB_ subroutines * *============================================================ program define GenXB_normal , rclass version 8.2 syntax varlist [fweight pweight iweight /] [ , method(string) change(real 0) blist(string) xlist(string) offset(string) cons * ] if `"`weight'"'!="" local weight w(`exp') if `"`offset'"'!="" local offset offset(`offset') gettoken treat touse : varlist return clear tempname b x xb mat `b' = `blist' if "`method'"=="at" { GenXB_Mat `treat' `touse' , b(`b') x(`xlist') change(`change') `weight' `offset' `cons' scalar `xb' = r(xb) return scalar xb1 = `xb' } else { GenXB_Var `treat' `touse' `xlist' , b(`b') change(`change') `weight' `offset' `cons' `options' return local xb1 `xlist' } end program define GenXB_multieq , rclass version 8.2 syntax varlist [iweight /] [ , method(string) change(real 0) blist(string) xlist(string) offset(string) cons * ] if `"`weight'"'!="" local weight w(`exp') if `"`offset'"'!="" local offset offset(`offset') gettoken treat touse : varlist return clear tempname b x xb local neq : word count `blist' forval eq = 1/`neq' { local bvec : word `eq' of `blist' local xvec : word `eq' of `xlist' mat `b' = `bvec' if "`method'"=="at" { GenXB_Mat `treat' `touse' , b(`b') x(`xvec') change(`change') `weight' `offset' `cons' scalar `xb' = r(xb) return scalar xb`eq'= `xb' } else { GenXB_Var `treat' `touse' `xvec' , b(`b') change(`change') `weight' `offset' `cons' `options' eq(`eq') return local xb`eq' `xvec' } } end program define GenXB_ordered GenXB_multieq `0' end program define GenXB_select GenXB_multieq `0' end program define GenXB_Mat , rclass version 8 syntax varlist [ , change(real 0) b(string) x(string) offset(string) cons w(varname) * ] tempname xb mat if "`cons'"!="" { local ncol = colsof(`b') scalar `xb' = `b'[1,`ncol'] return scalar xb = `xb' exit } gettoken treat touse : varlist // Taking care of dummies option mat `mat'=`x'' if "`s(D_`treat')'"!="" { local dvarlist `s(D_`treat')' local dvarlist : subinstr local dvarlist "`treat'" "" , all word foreach dvar in `dvarlist' { local pos = rownumb(`mat', "`dvar'" ) if `pos'<=rowsof(`mat') & `pos'>0 mat `mat'[`pos',1]=0 } } mat `mat' = `b'*`mat' scalar `xb' = `mat'[1,1] if "`offset'"!="" { scalar `xb' = `xb' + `offset' } local beta = colnumb(`b',"`treat'") if `beta'!=. local beta = `b'[1,`beta'] else local beta = 0 scalar `xb' = `xb' + `change'*`beta' return scalar xb = `xb' end program define GenXB_Var , rclass version 8 syntax varlist [ , change(real 0) b(string) x(string) offset(varname) cons w(varname) eq(integer 0) * ] gettoken treat rest : varlist gettoken touse xbvar : rest capture assert `treat'==0 | `treat'==1 if `touse' local isdummy = (_rc==0) tempvar xb tempname mat coef // Baseline probability if "`cons'"!="" { local pos = colsof(`b') if "`s(type)'"=="ordered" & `eq'<`s(neq)' { local pos = `pos'-(`s(neq)'-`eq') di "Pos-(s(neq)-eq = `pos'-(`s(neq)'-`eq')" `pos'-(`s(neq)'-`eq') } local beta = `b'[1,`pos'] qui gen double `xb' = `beta' if `touse' return local xb `xb' exit } // Preparing the matrix for calculating the linear prediction mat `coef' = `b' // Taking care of dummies option if "`s(D_`treat')'"!="" { local dvarlist `s(D_`treat')' local dvarlist : subinstr local dvarlist "`treat'" "" , all word foreach dvar in `dvarlist' { local pos = colnumb(`coef', "`dvar'" ) if `pos'<=colsof(`coef') & `pos'>0 mat `coef'[1,`pos']=0 } } if "`offset'"!="" { tempname add mat `add' = J(1,1,1) mat colnames `add' = `offset' mat `coef' = `add' , `coef' } // Linear prediction mat score double `xb' = `coef' if `touse'==1 , forcezero // Manipulating the linear prediction local beta = colnumb(`b',"`treat'") if `beta'!=. local beta = `b'[1,`beta'] else local beta = 0 if `isdummy'==1 { if `change'==1 qui replace `xb' = `xb' + `beta' if `touse'==1 & `treat'==0 else qui replace `xb' = `xb' - `beta' if `touse'==1 & `treat'==1 } else { qui replace `xb' = `xb' + `change'*`beta' if `touse'==1 } qui replace `xbvar' = `xb' return add *return local xb `xb' end *============================================================ * * [6] CalcME_ SUBROUTINES CALCULATING PARTIAL EFFECTS AND STANDARD ERRORS * *============================================================ program define CalcME_B_at , rclass args hat1 hat0 x1 x0 tempname me scalar `me' = (`hat1' - `hat0')/(`x1'-`x0') return scalar me = `me' end program define CalcME_B_av , rclass args hat1 hat0 x1 x0 touse weight tempvar me tempname b mean qui gen double `me' = (`hat1' - `hat0')/(`x1'-`x0') if `touse'==1 qui su `me' if `touse'==1 `weight' scalar `mean' = r(mean) return scalar me = `mean' end program define CalcME_SE_at , rclass args var phi1 phi0 x1 x0 xmat tempname vec add mat `vec' = `xmat'*(`phi1'-`phi0')/(`x1'-`x0') local pos = colnumb(`xmat',"`var'") if `pos'!=. { scalar `add' = `vec'[1,`pos'] scalar `add' = `add'+(`x1'*`phi1'-`x0'*`phi0')/(`x1'-`x0') mat `vec'[1,`pos'] = `add' } return mat vec = `vec' end program define CalcME_SE_av , rclass args var phi1 phi0 x1 x0 xmat markvar weight tempname vec add cons tempvar phivar addvar local xvars : colnames `xmat' local xvars : subinstr local xvars "_cons" "" , word all local ncol = colsof(`xmat') local pos : word count `xvars' qui gen double `phivar' = (`phi1'-`phi0')/(`x1'-`x0') if `markvar' mat vecaccum `vec' = `phivar' `xvars' if `markvar' , nocons qui count if `markvar'==1 local numobs = r(N) mat `vec' = `vec' / `numobs' // Adding constants CalcMean `phivar' if `markvar' `weight' local pos = colnumb(`xmat',"_cons") if `pos'!=. { mat `cons' = `xmat'[1,`pos'..`ncol'] mat `cons' = `cons'*r(mean) mat colnames `cons' = _cons mat `vec' = `vec' , `cons' } // Taking care of diagonal elements local pos : list posof "`var'" in xvars if `pos'!=0 { scalar `add' = `vec'[1,`pos'] qui gen double `addvar' = (`x1'*`phi1'-`x0'*`phi0')/(`x1'-`x0') CalcMean `addvar' if `markvar' `weight' scalar `add' = `add'+ r(mean) mat `vec'[1,`pos'] = `add' } return mat vec = `vec' end program define CalcMean , rclass version 8 syntax varlist [if/] [fw pw iw aw /] [ , n(integer 0) ] tempvar varsum wgtsum tempname mean if "`exp'"=="" { qui gen double `varsum' = sum(`varlist'*`if') scalar `mean' = `varsum'[_N] if `n'==0 { qui count if `if'==1 scalar `mean' = `mean'/r(N) } else scalar `mean' = `mean'/`n' } else { qui gen double `varsum' = sum(`varlist'*`exp'*`if') qui gen double `wgtsum' = sum(`exp'*`if') scalar `mean' = `varsum'[_N]/`wgtsum'[_N] } return scalar mean = `mean' end *============================================================ * * [7] CDF AND DENSITY --- GenFx_ subroutines * *============================================================ program define GenFx_log , sclass sret local cdf exp(`2'*`r(xb1)') sret local dens1 `2'*exp(`2'*`r(xb1)') * sret local title "exp(xb) (xb: linear prediction)" end program define GenFx_loglog , sclass sret local cdf exp(-exp(-`r(xb1)')) sret local dens1 exp(-exp(-`r(xb1)'))*exp(-`r(xb1)') * sret local title "exp(-exp(-xb)) (xb: linear prediction)" end program define GenFx_logc , sclass sret local cdf 1-exp(`r(xb1)') sret local dens1 -exp(`r(xb1)') * sret local title "1-exp(xb) (xb: linear prediction)" end program define GenFx_nbinomial , sclass sret local cdf `cdf' sret local dens1 (`cdf')*(1+`cdf') * sret local title "exp(xb)/(1-exp(xb)) (xb: linear prediction)" end program define GenFx_power , sclass local power = 1/`s(power)' sret local cdf (`r(xb1)')^`power' sret local dens1 `power'*(`r(xb1)')^(`power'-1) * sret local title "(xb)^(1/`s(power)') (xb: linear prediction)" end program define GenFx_opower , sclass local power = 1/`s(power)' local xb (1+`power'*`r(xb1)')^(1/`power') sret local cdf (`xb')/(1+(`xb')) sret local dens1 (1/(1+(`xb'))^2)*(1+`power'*`r(xb1)')^((1/`power')-1) * sret local title "1/(1+(1+xb/`s(power)')^(-`s(power)') (xb: linear prediction)" end program define GenFx_probit , sclass sret local cdf normprob(`2'*`r(xb1)') sret local dens1 `2'*normd(`2'*`r(xb1)') end program define GenFx_heckprob , sclass local depvar = e(depvar) local depvar : word 1 of `depvar' //GenFx_probit `depvar' `e(rho)' GenFx_probit `depvar' 1 end program define GenFx_logit , sclass sret local cdf invlogit(`2'*`r(xb1)') sret local dens1 `2'*(invlogit(`2'*`r(xb1)'))*(1-invlogit(`2'*`r(xb1)')) end program define GenFx_logistic GenFx_logit `0' end program define GenFx_cloglog , sclass sret local cdf invcloglog(`2'*`r(xb1)') sret local dens1 -`2'*(1-invcloglog(`2'*`r(xb1)'))*exp(`2'*`r(xb1)') end program define GenFx_oprobit , sclass args out local ncat = e(k_cat) if `out'>`ncat' { di as error "Requested outcome larger than number of categories" error 198 } local neq = cond("`e(cmd)'"=="ologit2", `s(neq)', `s(neq)') forval i = 1/`neq' { sret local dens`i' 0 } local low = `out'-1 if `out'==1 { sret local cdf 1-normprob(`r(xb1)') sret local dens1 -normd(`r(xb1)') exit } else if `out'==`ncat' { sret local cdf normprob(`r(xb`low')') sret local dens`low' normd(`r(xb`low')') exit } else { sret local cdf normprob(`r(xb`low')')-normprob(`r(xb`out')') sret local dens`low' normd(`r(xb`low')') sret local dens`out' -normd(`r(xb`out')') } end program define GenFx_ologit , sclass args out depvar local ncat = e(k_cat) if `out'>`ncat' { di as error "Requested outcome larger than number of categories" error 198 } local neq = cond("`e(cmd)'"=="gologit2", `s(neq)', `s(neq)') forval i = 1/`neq' { sret local dens`i' 0 } local low = `out'-1 if `out'==1 { sret local cdf 1-invlogit(`r(xb1)') sret local dens1 -invlogit(`r(xb1)')*(1-invlogit(`r(xb1)')) exit } else if `out'==`ncat' { sret local cdf invlogit(`r(xb`low')') sret local dens`low' invlogit(`r(xb`low')')*(1-invlogit(`r(xb`low')')) exit } else { sret local cdf invlogit(`r(xb`low')')-invlogit(`r(xb`out')') sret local dens`low' invlogit(`r(xb`low')')*(1-invlogit(`r(xb`low')')) sret local dens`out' - invlogit(`r(xb`out')')*(1-invlogit(`r(xb`out')')) } end program define GenFx_ocloglog , sclass args out if `out'>`s(nout)' { di as error "Requested outcome larger than number of categories" error 198 } local neq = cond("`e(cmd)'"=="ologit2", `s(neq)', `s(neq)') forval i = 1/`neq' { sret local dens`i' 0 } local low = `out'-1 if `out'==1 { sret local cdf 1-invcloglog(`r(xb1)') sret local dens1 (1-invcloglog(`r(xb1)'))*exp(`r(xb1)') exit } else if `out'==`s(nout)' { sret local cdf invcloglog(`r(xb`low')') sret local dens`low' -(1-invcloglog(`r(xb`low')'))*exp(`r(xb`low')') exit } else { sret local cdf invcloglog(`r(xb`low')')-invcloglog(`r(xb`out')') sret local dens`low' -(1-invcloglog(`r(xb`low')'))*exp(`r(xb`low')') sret local dens`out' (1-invcloglog(`r(xb`out')'))*exp(`r(xb`out')') } end program define GenFx_ologlog , sclass args out if `out'>`s(nout)' { di as error "Requested outcome larger than number of categories" error 198 } local neq = cond("`e(cmd)'"=="ologit2", `s(neq)', `s(neq)') forval i = 1/`neq' { sret local dens`i' 0 } local low = `out'-1 if `out'==1 { sret local cdf invcloglog(-`r(xb1)') sret local dens1 -(1-invcloglog(-`r(xb1)'))*exp(-`r(xb1)') exit } else if `out'==`s(nout)' { sret local cdf 1-invcloglog(-`r(xb`low')') sret local dens`low' (1-invcloglog(-`r(xb`low')'))*exp(-`r(xb`low')') exit } else { sret local cdf invcloglog(-`r(xb`out')')-invcloglog(-`r(xb`low')') sret local dens`low' (1-invcloglog(-`r(xb`low')'))*exp(-`r(xb`low')') sret local dens`out' -(1-invcloglog(-`r(xb`out')'))*exp(-`r(xb`out')') } end program define GenFx_ocauchit , sclass args out if `out'>`s(nout)' { di as error "Requested outcome larger than number of categories" error 198 } local neq = cond("`e(cmd)'"=="ologit2", `s(neq)', `s(neq)') forval i = 1/`neq' { sret local dens`i' 0 } local low = `out'-1 if `out'==1 { sret local cdf .5+(1/_pi)*atan(-(`r(xb1)')) sret local dens1 -(1/_pi)/(1+(`r(xb1)')^2) exit } else if `out'==`s(nout)' { sret local cdf (1-.5-(1/_pi)*atan(-`r(xb`low')')) sret local dens`low' (1/_pi)/(1+(`r(xb`low')')^2) exit } else { sret local cdf (1/_pi)*atan(`r(xb`low')')-(1/_pi)*atan(`r(xb`out')') sret local dens`low' (1/_pi)/(1+(`r(xb`low')')^2) sret local dens`out' -(1/_pi)/(1+(`r(xb`out')')^2) } end program define GenFx_mlogit , sclass args out if "`e(ibasecat)'"!="" local base = e(ibasecat) else local base = e(ibaseout) local i = `out'-(`out'>`base') // denominator local denom 1 forval eq = 1/`s(neq)' { local denom `denom'+exp(`r(xb`eq')') } if `out'==`base' sret local cdf "(1/(`denom'))" else sret local cdf (exp(`r(xb`i')')/(`denom')) forval eq = 1/`s(neq)' { if `out'==`base' { sret local dens`eq' "(-exp(`r(xb`eq')')/(`denom')^2)" } else { if `i'==`eq' { sret local dens`eq' "(exp(`r(xb`i')')/(`denom')-exp(`r(xb`i')')^2/(`denom')^2)" } else { sret local dens`eq' "(-exp(`r(xb`i')')*exp(`r(xb`eq')')/(`denom')^2)" } } } end program define GenFx_poisson , sclass GenFx_log `0' end program define GenFx_nbreg , sclass GenFx_log `0' end program define GenFx_zip , sclass version 8 local mu "exp(`r(xb1)')" if "`e(inflate)'"=="logit" { local cuminf 1/(1+exp(-`r(xb2)')) local deninf (`cuminf')*(1-`cuminf') } else { local cuminf normprob(`r(xb2)') local deninf normd(`r(xb2)') } sret local cdf `mu'*(1-`cuminf') sret local dens1 `mu'*(1-`cuminf') sret local dens2 -`mu'*`deninf' end program define GenFx_zinb, sclass GenFx_zip `0' end program define GenFx_biprobit , sclass args out local rho = e(rho) local srh = sqrt(1-`rho'^2) local A local C local B if `out'<3 local A - if `out'==1 | `out'==3 local B - if `out'==2 | `out'==3 local C - sret local cdf binorm(`A'`r(xb1)',`B'`r(xb2)',`C'`rho') sret local dens1 (`A'normd(`A'`r(xb1)')*normprob((`B'`r(xb2)'-`C'`rho'*`A'`r(xb1)')/`srh')) sret local dens2 (`B'normd(`B'`r(xb2)')*normprob((`A'`r(xb1)'-`C'`rho'*`B'`r(xb2)')/`srh')) end ******* SELECTION MODELS ************** program define GenFx_truncreg , sclass version 8 args v1 local sig = `e(sigma)' local ul = e(ulopt) local ll = e(llopt) if `ll'==. local vll = -10000 else local vll (`ll'-`r(xb1)')/`sig' if `ul'==. local vul = 10000 else local vul (`ul'-`r(xb1)')/`sig' local mills "(normd(`vll')-normd(`vul'))/(normprob(`vul')-normprob(`vll')) " sret local cdf "`r(xb1)'+`sig'*`mills'" sret local dens1 "1+`mills'-(`mills')^2" end program define GenFx_heckman , sclass local sig = `e(sigma)'*`e(rho)' local mills normd(`r(xb2)')/normprob(`r(xb2)') sret local cdf (`r(xb1)'+`sig'*`mills') sret local dens1 1 sret local dens2 -`sig'*(`r(xb2)'*(`mills')+(`mills')^2 ) end program define GenFx_tobit , sclass // handling version 8 and 9 incompatibility capture di _b[sigma:_cons] if _rc==0 local sig = _b[sigma:_cons] else local sig = _b[_se] local ul = e(ulopt) local ll = e(llopt) if `ul'==. local ul = 0 if `ll'==. local ll = 0 if `ll'==0 local vll = -10000 else local vll (`ll'-`r(xb1)')/`sig' if `ul'==0 local vul = 10000 else local vul (`ul'-`r(xb1)')/`sig' GenFx_Censor , depvar(`depvar') ul(`ul') vul(`vul') ll(`ll') vll(`vll') sig(`sig') end program define GenFx_cnreg , sclass local depvar = e(depvar) local censor = e(censored) // handling version 8 and 9 incompatibility capture di _b[sigma:_cons] if _rc==0 local sig = _b[sigma:_cons] else local sig = _b[_se] local ll "((`censor'==-1)*`depvar')" local ul "((`censor'==1)*`depvar')" local vll "((`censor'==-1)*( `ll'-`r(xb1)')/`sig'+(1-`censor')*(-10000))" local vul "((`censor'==-1)*(`ull'-`r(xb1)')/`sig'+(1-`censor')*( 10000))" GenFx_Censor , depvar(`left' `right') ul(`ul') vul(`vul') ll(`ll') vll(`vll') sig(`sig') end program define GenFx_intreg , sclass args left right local sig = e(sigma) local censor "(`left'!=`right')" local ll "(`censor'*`left')" local ul "(`censor'*`right')" local vll "(`censor'*( `ll'-`r(xb1)')/`sig'+(1-`censor')*(-10000))" local vul "(`censor'*(`ull'-`r(xb1)')/`sig'+(1-`censor')*( 10000))" GenFx_Censor , depvar(`left' `right') ul(`ul') vul(`vul') ll(`ll') vll(`vll') sig(`sig') end program define GenFx_Censor , sclass // Source: Long Regression models.... p.213. version 8 syntax [, depvar(string) ll(string) ul(string) vll(string) vul(string) sig(real 0)] local cumu11 "`ll'*normprob(`vll')+`ul'*normprob(-`vul')" local cumu11 "`cumu11' + (normprob(`vul')-normprob(`vll'))*`r(xb1)'" local cumu11 "`cumu11' + `sig'*(normd(`vll')-normd(`vul'))" sret local cdf (`cumu11') sret local dens1 "normprob(`vul')-normprob(`vll')" end *============================================================ * * [8] ALL OTHER SUBROUTINES * *============================================================ program define CheckOut , rclass version 8 local nout : word count `0' if `nout'>1 { di as error "Only one number can be specified in the outcome() option" error 198 } capture confirm integer number `0' if _rc!=0 { di as red "the outcome() option must contain an integer number" error 198 } if `0'<1 { di as red "the number specified in the outcome() option must be positive" error 198 } local outcome = `0' if `s(nout)'==1 { di as text "outcome() option ignored because `s(model)' is a single-outcome model" di local outcome 1 } else { if `outcome'>`s(nout)' { di as error "Invalid outcome() option" error 125 } } return local numlist = `outcome' end program define AtList , rclass version 8 syntax , at(string asis) varlist(varlist) if "`s(dummies)'"!="" local dumlist = s(dummies) // Searching for mean median zero local at : subinstr local at "=" " " , all gettoken first : at if `"`first'"'=="mean" | `"`first'"'=="median" | `"`first'"'=="zero" { if `"`first'"'=="median" return local atstat p50 else return local atstat `first' gettoken first at : at } else return local atstat mean // Parsing var = value part while `"`at'"'!="" { gettoken varname at : at capture confirm var `varname' if _rc!=0 { di as error "Invalid at() option: "`varname' must be a numeric variable" error 198 } local test : subinstr local varlist "`varname'" "" , count(local c) all if `c'!=1 { di as error "Invalid at() option: "`varname' is not (or appears twice) among the explanatory variables" error 198 } gettoken value at : at if "`value'"=="=" | "`value'"=="==" { gettoken value at : at } capture confirm number `value' if _rc!=0 { di as error "Invalid at() option: `value' (which follows `varname') must be a number" error 198 } // Accumulating list local atvarlist `atvarlist' `varname' local atnumlist `atnumlist' `value' } if "`: word count `atvarlist' '"!="`: word count `atnumlist' '" { di as error "Invalid at() option: Number of variables not equal to number of values" error 198 } return local atvarlist `atvarlist' return local atnumlist `atnumlist' end program define DumList , sclass tempname b mat `b' = e(b) local xvars : colnames `b' local xvars : list uniq xvars // Step 1: Recognizing lists of dummies as tokens tokenize "`0'" , p(" \ ") local list `0' local n = 1 while `"`list'"'!="" { gettoken token list : list if "`token'"=="\" { if "`tokens`n''"!="" local ++n } else { local tokens`n' `tokens`n'' `token' } } // Step 2: Unabbreviating lists of dummies & eliminating dropped vars forval i = 1/`n' { capt tsunab list`i' : `tokens`i'' if _rc!=0 { sret local error "Invalid dummies() option: `tokens`i'' not among the regressors" exit } local list`i' : list list`i' & xvars } // Step 3: Adding the lists to s(D_) locals // Again a loop - reason: tsunab overwrites s macros! forval i = 1/`n' { foreach var in `list`i'' { sret local D_`var' `list`i'' local dummies `dummies' `var' } } local dummies : list uniq dummies sret local dummies `dummies' end program define CalcUnit , rclass version 8 args var touse // This part determining unit of measurement is taken from codebook.ado tempvar p scalar `p' = 1 capture assert float(`var') == float(round(`var',1)) if `touse' if _rc == 0 { while _rc == 0 { scalar `p' = `p'*10 capture assert float(`var') == float(round(`var',`p')) if `touse' } scalar `p' = `p'/10 } else { while _rc { scalar `p' = `p'/10 capture assert float(`var') == float(round(`var',`p')) if `touse' } } return scalar unit = `p' end program define EqLabel , sclass version 8 args out touse if "`s(model)'"=="biprobit" { local label : word `out' of p00 p01 p10 p11 sret local name "`label'" sret local cat "`label'" exit } if "`out'"=="" { di as error "Error when calling subroutine EqLabel: option out() not specified" error 198 } local depvar = e(depvar) if `: word count `depvar''>1 { local depvar : word 1 of `e(depvar)' } tempname mat qui tab `depvar' if `touse' , matrow(`mat') local cat = `mat'[`out',1] local label : label (`depvar') `cat' if `"`label'"'=="" { local label "p`cat'" } sret local cat "`cat'" sret local lab "`label'" end program define SendVersion , sclass version 8 sret local margeff_version = 218 end