05
Jul
2018
Example Stata Code for Hanmer and Kalkan (2013, AJPS)
by Shawna K. Metzger
shawnakmetzger.com/wp/stata_hk2013/
Accessed April 27, 2024. Last modified April 2, 2021
Quick code to demonstrate how to implement Hanmer and Kalkan’s observed value approach when generating predicted probabilities in Stata. If you find errors in the code, let me know (email button’s at bottom of page). I might generalize the code some day, but until then:
Click the code block to expand it.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 |
* Dichotomous IV global nSims = 1000 mata mata clear webuse lbw, clear estsimp logit low age lwt smoke ptl ht ui, sims($nSims) // ^ Could also generate betas using drawnorm // Simulate using H&K's observed values approach. // NOTE: requires moremata package. (ssc install moremata) // Written entirely in Mata because, for super huge datasets, non-Mata code // will take a while. Mata's 7.8 times faster, for this specific example. // Let's say moving smoke from 0 to 1 is the effect we want to see local keyIV smoke qui{ // Make the matrix of simulated betas putmata b=(b1-b7), replace matrix b = e(b) // Make the matrix filled with the data // ** make sure covariates appear in same order as simulated betas putmata X=(age lwt smoke ptl ht ui), replace omit mata: X = X, J(rows(X), 1, 1) // adding col for constant local col = colnumb(b, "`keyIV'") // column of matrix containing key IV // loop over covar values forvalues sm=0/1{ mata: X[,`col'] = J(rows(X), 1, `sm') // Insert the covariate value // Go over each sim draw, calculate pred prob forvalues s=1/$nSims{ mata: temp = invlogit(rowsum(b[`s',] :* X)) // pred prob (vector of 1 x n) if(`s'==1) mata: meanIntm`sm' = mean(temp) // running list of each sim draw's average else mata: meanIntm`sm' = meanIntm`sm' \ mean(temp) } } // All simulations done. Need to now average across all the sim draws + pull quantiles. forvalues sm=0/1{ // can comment out, if we want to keep mata mem tidy mata: mean`sm' = mean(meanIntm`sm') mata: lb`sm' = mm_quantile(meanIntm`sm', 1, 0.025) mata: ub`sm' = mm_quantile(meanIntm`sm', 1, 0.975) // put the means and CIs into local macros mata: st_local("mean`sm'", strofreal(mean(meanIntm`sm'))) mata: st_local("lb`sm'", strofreal(mm_quantile(meanIntm`sm', 1, 0.025))) mata: st_local("ub`sm'", strofreal(mm_quantile(meanIntm`sm', 1, 0.975))) } // Display results forvalues sm=0/1{ noi di as gr _n "`=strproper("`keyIV'")' = " as ye "`sm'" noi di as gr " Mean: " as ye %5.4f "`mean`sm''" noi di as gr " LB95: " as ye %5.4f "`lb`sm''" noi di as gr " UB95: " as ye %5.4f "`ub`sm''" } } // end of qui // mata mata clear // if you want mata memory emptied. |
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 |
// requires moremata package (ssc install moremata) * Continuous IV, with different SDs away from mean global nSims = 1000 webuse fullauto, clear global dv rep77 // DV local ivs foreign length weight // All IVs // Estimate oprobit $dv `ivs' , robust keep if(e(sample)==1) // for simplicity's sake // Doing the betas manually. { // convenience collapse bracket mat b = e(b) mat V = e(V) local params = colsof(b) local pString = "" // for drawnorm local aux = "" // for cutpoint betas local msn = "" // for covar betas local rhs = "" // for the actual rhs varnames local cNames: colnames b local colEqNames: coleq b forvalues b = 1/`params'{ local pString = "`pString' b`b'" // get eqname for cut (Stata 15 does this differently) if(`c(version)'<15) local stStr = "cut" else local stStr = "/" // get the colname, and if it's _cut, add to the aux list local cN: word `b' of `cNames' local cEq: word `b' of `colEqNames' if(substr("`cEq'", 1, 3)=="`stStr'") local aux = "`aux' b`b'" else{ local msn = "`msn' b`b'" local rhs = "`rhs' `cN'" } } preserve drawnorm `pString', n($nSims) mean(b) cov(V) d clear // can force to be (semi)pos definite if needed with forcespd option // go through and label the vars, for sanity's sake forvalues b = 1/`params'{ local cN: word `b' of `cNames' local cEq: word `b' of `colEqNames' if(substr("`cEq'", 1, 3)=="cut'") label variable b`b' "Simulated _`cEq' parameter" else label variable b`b' "Simulated `cN' parameter" } tempfile sims save `sims', replace restore // Bring back in to main dataset merge 1:1 _n using `sims', nogen order `pString', first } // end of beta convenience collapse bracket *********** // We're now to the point where we can do the simulations local keyIV length local sdRange "-2 -1 0 1 2" // This will take longer to run (compared to Clarify's usual blink-and-done). Prepare yourself. qui{ sum `keyIV' local mean = `r(mean)' local sd = `r(sd)' levelsof $dv, local(dvLvs) // Make the matrix of simulated betas putmata b=(`msn'), replace omit putmata cut=(`aux'), replace omit // Make the matrix filled with the data * make sure covariates appear in same order as simulated betas putmata X=(`rhs'), replace omit ** NOTE: no constant for oprobit // col of matrix containing kevIV local col = colnumb(b, "`keyIV'") local cnt = 1 // since SD range can have neg values, bypass that problem altogether with a counter foreach sm of numlist `sdRange'{ noi di as gr _n "SD value: `sm'" // get the keyIV's value local val = `mean' + `sm'*`sd' mata: X[,`col'] = J(rows(X), 1, `val') // key IV local cutptCnt = 1 // first column to start in cut matrix foreach o of local dvLvs{ noi di as gr " DV = `o'" // Go over each sim draw, calculate pred prob forvalues s=1/$nSims{ // for percentage complete updates if(mod(`s',$nSims*0.05)==0){ if(`s'==$nSims*0.05) noi di _c as ye " " `s'/$nSims*100 "%" else noi di _c as ye "..." `s'/$nSims*100 "%" } // first outcome in list if("`o' `ferest()'"=="`dvLvs'") mata: temp = normal(cut[`s',`cutptCnt'] :- rowsum(b[`s',] :* X)) // cutlo - xb // not first else{ // not last outcome in list if("`ferest()'"!="") mata: temp = normal(cut[`s',`cutptCnt'] :- rowsum(b[`s',] :* X)) /// // cuthi - cutlo - normal(cut[`s',`cutptCnt'-1] :- rowsum(b[`s',] :* X)) // last outcome in list else mata: temp = normal(rowsum(b[`s',] :* X) :- cut[`s',`=`cutptCnt'-1']) // xb - cutlo } if(`s'==1) mata: meanIntm`cnt'_`o' = mean(temp) // running list of each sim draw's average else mata: meanIntm`cnt'_`o' = meanIntm`cnt'_`o' \ mean(temp) } noi di "" local `cutptCnt++' } local `cnt++' } // All simulations done. Need to now average across all the sim draws + pull quantiles. local cnt = 1 foreach sm of numlist `sdRange'{ foreach o of local dvLvs{ // can comment out, if we want to keep mata mem tidy mata: mean`cnt'_`o' = mean(meanIntm`cnt'_`o') mata: lb`cnt'_`o' = mm_quantile(meanIntm`cnt'_`o', 1, 0.025) mata: ub`cnt'_`o' = mm_quantile(meanIntm`cnt'_`o', 1, 0.975) // put the means and CIs into local macros mata: st_local("mean`cnt'_`o'", strofreal(mean(meanIntm`cnt'_`o'))) mata: st_local("lb`cnt'_`o'", strofreal(mm_quantile(meanIntm`cnt'_`o', 1, 0.025))) mata: st_local("ub`cnt'_`o'", strofreal(mm_quantile(meanIntm`cnt'_`o', 1, 0.975))) } local `cnt++' } // Display results, save to variables gen double xVal= . foreach o of local dvLvs{ gen double pr`o' = . gen double pr`o'_lb = . gen double pr`o'_ub = . noi di as gr _n _n ">> Pr(`e(depvar)') = " as ye "`o'" as gr " <<" local cnt = 1 foreach sm of numlist `sdRange'{ // get the x value local val = `mean' + `sm'*`sd' noi di as gr _n "`keyIV' = " as ye "`val'" noi di as gr " Mean: " as ye %5.4f "`mean`cnt'_`o''" noi di as gr " LB95: " as ye %5.4f "`lb`cnt'_`o''" noi di as gr " UB95: " as ye %5.4f "`ub`cnt'_`o''" replace xVal = `val' in `cnt' replace pr`o' = `mean`cnt'_`o'' in `cnt' replace pr`o'_lb = `lb`cnt'_`o'' in `cnt' replace pr`o'_ub = `ub`cnt'_`o'' in `cnt' local `cnt++' } } } // end of qui // mata mata clear // if you want mata memory emptied. |