In Stata 16, the code below cannot run normally, but in Stata 15 and Stata 13, the code can run correctly, however, the results between them is different.
Which result should I believe?
Which result should I believe?
* Example generated by -dataex-. To install: ssc install dataex clear input double(dv id serial) float(iv med mod1 newid) 4 1 1 .12499996 .5486109 -.05694437 1 4 1 2 -.07500005 -.17361127 -.05694437 1 4 1 3 -.07500005 .5486109 -.05694437 1 3.8 1 4 -.07500005 -.06250016 -.05694437 1 4 1 5 -.07500005 .3263887 -.05694437 1 3 1 6 -.07500005 -.3958335 -.05694437 1 4 1 7 .325 -.3958335 -.05694437 1 4 1 8 -.07500005 -.3958335 -.05694437 1 3.2 2 1 1.725 -.4375001 .14305563 2 3.7 2 2 .725 -.3819445 .14305563 2 3.7 2 3 -.6750001 .006944391 .14305563 2 3.4 2 4 -.07500005 .2291666 .14305563 2 3.3 2 5 -1.075 -.54861116 .14305563 2 3 2 6 .325 .2291666 .14305563 2 4.7 2 7 -.27500004 .4513888 .14305563 2 5 2 8 -.6750001 .4513888 .14305563 2 4 3 1 -.125 -.013888836 -.3569444 3 3.5 3 2 .075 -.013888836 -.3569444 3 3.8 3 3 -.125 -.013888836 -.3569444 3 3.6 3 4 -.125 .09722228 -.3569444 3 3.8 3 5 .075 -.013888836 -.3569444 3 4 3 6 .075 -.013888836 -.3569444 3 3.8 3 7 .075 -.013888836 -.3569444 3 4 3 8 .075 -.013888836 -.3569444 3 3.2 4 1 -.17142864 -.007936531 -.4569444 4 3.7 4 2 -.7714286 .04761903 -.4569444 4 2.7 4 3 .02857137 -.06349209 -.4569444 4 2.2 4 4 .4285714 -.11904764 -.4569444 4 2.4 4 5 .6285714 .10317458 -.4569444 4 2.9 4 6 .22857137 .3253968 -.4569444 4 . 4 7 . . -.4569444 4 3.8 4 8 -.3714286 -.2857143 -.4569444 4 3.9 5 1 .025000047 .58333325 -.15694436 5 4.3 5 2 -.17499995 .02777767 -.15694436 5 4.1 5 3 -.17499995 .4166666 -.15694436 5 3.4 5 4 .025000047 -.25000012 -.15694436 5 3.5 5 5 .425 -.3611112 -.15694436 5 3.8 5 6 .22500005 -.19444455 -.15694436 5 3.1 5 7 -.17499995 -.25000012 -.15694436 5 4.3 5 8 -.17499995 .02777767 -.15694436 5 4 6 1 .7000001 -.08333354 .14305563 6 4.2 6 2 .3000001 .472222 .14305563 6 4 6 3 .1000001 .027777566 .14305563 6 3.6 6 4 -.2999999 -.08333354 .14305563 6 3.3 6 5 -.2999999 -.08333354 .14305563 6 3.8 6 6 -.2999999 -.1388891 .14305563 6 3.3 6 7 .1000001 -.08333354 .14305563 6 3.3 6 8 -.2999999 -.02777799 .14305563 6 3.8 7 1 -.0999999 .06250016 .14305563 7 3.5 7 2 -.0999999 .3402779 .14305563 7 3.3 7 3 -.0999999 .11805572 .14305563 7 2.8 7 4 -.0999999 -.1041665 .14305563 7 3 7 5 .3000001 -.1041665 .14305563 7 3 7 6 .1000001 -.1041665 .14305563 7 3.2 7 7 .1000001 -.1041665 .14305563 7 3 7 8 -.0999999 -.1041665 .14305563 7 . 8 1 . . .4430556 8 . 8 2 . . .4430556 8 4.5 8 3 .033333253 -.111111 .4430556 8 3.6 8 4 -.3666667 -.111111 .4430556 8 3.4 8 5 .033333253 -.05555545 .4430556 8 3.1 8 6 -.16666675 .22222233 .4430556 8 3.4 8 7 .43333325 .16666678 .4430556 8 3.6 8 8 .033333253 -.111111 .4430556 8 4.1 9 1 .12500003 .3680556 -.05694437 9 4.6 9 2 .12500003 .47916675 -.05694437 9 4.5 9 3 -.27499998 -2.1319444 -.05694437 9 4 9 4 .12500003 .14583342 -.05694437 9 4.5 9 5 .12500003 .14583342 -.05694437 9 4 9 6 -.27499998 .3680556 -.05694437 9 4.5 9 7 .12500003 .25694454 -.05694437 9 4.5 9 8 -.07499997 .3680556 -.05694437 9 3.4 10 1 .2 0 .8430556 10 3.5 10 2 .2 0 .8430556 10 3.7 10 3 .2 0 .8430556 10 3 10 4 0 0 .8430556 10 3.3 10 5 .2 0 .8430556 10 3.6 10 6 -.8 0 .8430556 10 3 10 7 0 0 .8430556 10 3 10 8 0 0 .8430556 10 3.8 11 1 .8 0 .24305563 11 3.8 11 2 -.20000005 0 .24305563 11 3.4 11 3 -.20000005 0 .24305563 11 3.8 11 4 -.20000005 0 .24305563 11 3.6 11 5 -4.768372e-08 0 .24305563 11 4 11 6 -.20000005 0 .24305563 11 3.6 11 7 -4.768372e-08 0 .24305563 11 3.4 11 8 -4.768372e-08 0 .24305563 11 4.1 12 1 -.75 0 -.05694437 12 3.3 12 2 -.55 0 -.05694437 12 3.8 12 3 .25 0 -.05694437 12 3.9 12 4 -.75 0 -.05694437 12 3.9 12 5 -.35 0 -.05694437 12 3.6 12 6 .65 0 -.05694437 12 3 12 7 .25 0 -.05694437 12 4.3 12 8 1.25 0 -.05694437 12 3 13 1 0 -.009259198 -.6569444 13 3 13 2 0 -.009259198 -.6569444 13 3 13 3 0 -.064814754 -.6569444 13 3 13 4 0 -.009259198 -.6569444 13 end xtset id serial set seed 10000 capture program drop MoDMed program define MoDMed, rclass quietly summarize mod1 return list global m=r(mean) global sd=r(sd) mixed med l.iv mod1 cl.iv#c.mod1 || id: l.iv, var cov(exc) iter(50) return scalar al=(_b[med:L.iv]+_b[med:cL.iv#c.mod1]*($m-$sd)) return scalar ah=(_b[med:L.iv]+_b[med:cL.iv#c.mod1]*($m+$sd)) mixed dv med l.iv || id: l.iv med, var cov(exc) iter(50) return scalar b=(_b[dv:med]) end bootstrap indL=(r(al)*r(b)) indH=(r(ah)*r(b)) diff=(r(ah)*r(b)-r(al)*r(b)), reps(50) reject(e(converged)!=1) cluster(id) idcluster(newid) group(serial): MoDMed estat bootstrap program drop MoDMed
Bootstrap replications (50) ----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 .................................................. 50 Bootstrap results Number of obs = 83 Replications = 50 command: MoDMed indL: r(al)*r(b) indH: r(ah)*r(b) diff: r(ah)*r(b)-r(al)*r(b) (Replications based on 13 clusters in id) ------------------------------------------------------------------------------ | Observed Bootstrap Normal-based | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- indL | -.0221546 .1279645 -0.17 0.863 -.2729605 .2286513 indH | -.0682389 .0730797 -0.93 0.350 -.2114725 .0749947 diff | -.0460843 .1628654 -0.28 0.777 -.3652946 .273126 ------------------------------------------------------------------------------ . estat bootstrap Bootstrap results Number of obs = 83 Replications = 50 command: MoDMed indL: r(al)*r(b) indH: r(ah)*r(b) diff: r(ah)*r(b)-r(al)*r(b) (Replications based on 13 clusters in id) ------------------------------------------------------------------------------ | Observed Bootstrap | Coef. Bias Std. Err. [95% Conf. Interval] -------------+---------------------------------------------------------------- indL | -.02215458 -.0104108 .12796453 -.5778048 .1074098 (BC) indH | -.0682389 -.0272837 .07307971 -.2348322 -.0035386 (BC) diff | -.04608432 -.0168729 .16286538 -.1381491 .5742661 (BC) ------------------------------------------------------------------------------ (BC) bias-corrected confidence interval
clear
version 15.1 //sets the version of the interpreter to use for the entirety of the do file
...
xtset id serial
set seed 10000
capture program drop MoDMed
program define MoDMed, rclass
quietly summarize mod1
return list
global m=r(mean)
global sd=r(sd)
mixed med l.iv mod1 cl.iv#c.mod1 || id: l.iv, var cov(exc) iter(50)
return scalar al=(_b[med:L.iv]+_b[med:cL.iv#c.mod1]*($m-$sd))
return scalar ah=(_b[med:L.iv]+_b[med:cL.iv#c.mod1]*($m+$sd))
mixed dv med l.iv || id: l.iv med, var cov(exc) iter(50)
return scalar b=(_b[dv:med])
end
bootstrap indL=(r(al)*r(b)) indH=(r(ah)*r(b)) diff=(r(ah)*r(b)-r(al)*r(b)), reps(50) reject(e(converged)!=1) cluster(id) idcluster(newid) group(serial): MoDMed
estat bootstrap
program drop MoDMed
(running MoDMed on estimation sample) Bootstrap replications (50) ----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 50 insufficient observations to compute bootstrap standard errors no results will be saved r(2000);
capture program drop MoDMed
program define MoDMed, rclass
sort id serial
quietly summarize mod1
return list
global m=r(mean)
global sd=r(sd)
mixed med l.iv mod1 cl.iv#c.mod1 || id: l.iv, var cov(exc) iter(50)
return scalar al=(_b[med:L.iv]+_b[med:cL.iv#c.mod1]*($m-$sd))
return scalar ah=(_b[med:L.iv]+_b[med:cL.iv#c.mod1]*($m+$sd))
mixed dv med l.iv || id: l.iv med, var cov(exc) iter(50)
return scalar b=(_b[dv:med])
end
bootstrap indL=(r(al)*r(b)) indH=(r(ah)*r(b)) diff=(r(ah)*r(b)-r(al)*r(b)), reps(50) reject(e(converged)!=1) cluster(id) idcluster(newid) group(serial): MoDMed
estat bootstrap
program drop MoDMed
scalar m=r(mean) scalar sd=r(sd) return scalar al=(_b[med:L.iv]+_b[med:cL.iv#c.mod1]*(m-sd)) return scalar ah=(_b[med:L.iv]+_b[med:cL.iv#c.mod1]*(m+sd))
help program properties
program define mybs mixed y L.x || id: end
program define mybs, properties(xtbs) mixed y L.x || id: end
help whatsnew15
xtset id t
bootstrap ... , ... cluster(id) idcluster(newid) group(t)
bootstrap ... , ... cluster(id) idcluster(newid) group(id)
* Toy data: set seed 1234 set obs 50 gen id = _n gen u = rnormal(0.0,0.5) gen v = rnormal(0.0,0.2) expand 10 bys id : gen t = _n gen x1 = runiform() gen e = rnormal() tsset id t gen y = 5 + L.x1 + L.x1*v + u + e drop v u e * Program for bootstrapping: program define mybs, rclass preserve bsample, cluster(id) idcluster(newid) sort newid id, stable qui by newid id: gen newgroup = _n == 1 qui replace newgroup = sum(newgroup) drop id rename newgroup id tsset id t mixed y L.x1 || id: L.x1 return scalar b1 = _b[L.x1] return scalar b0 = _b[_cons] return scalar var_i = exp(_b[lns1_1_2:_cons])^2 return scalar var_c = exp(_b[lns1_1_1:_cons])^2 return scalar var_e = exp(_b[lnsig_e:_cons])^2 end * Point estimates: mixed y L.x1 || id: L.x1 * Bootstrapping standard errors: simulate b1=r(b1) b0=r(b0) var_c=r(var_c) var_i=r(var_i) var_e=r(var_e), /// reps(10) seed(123) : mybs * Show results: foreach v of varlist _all { qui sum `v' di as txt "SE `v'" _col(11) "= " as res `r(sd)' }
* Toy data: clear set seed 1234 set obs 50 gen id = _n gen u = rnormal(0.0,0.5) gen v = rnormal(0.0,0.2) expand 10 bys id : gen t = _n gen x1 = runiform() gen e = rnormal() tsset id t gen y = 5 + L.x1 + L.x1*v + u + e drop v u e * Program for -bootstrap-: program define mybs2, rclass mixed y L.x1 || id: L.x1 return scalar b1 = _b[L.x1] return scalar b0 = _b[_cons] return scalar var_i = exp(_b[lns1_1_2:_cons])^2 return scalar var_c = exp(_b[lns1_1_1:_cons])^2 return scalar var_e = exp(_b[lnsig_e:_cons])^2 end * Run -bootstrap-: bootstrap b1=r(b1) b0=r(b0) var_c=r(var_c) var_i=r(var_i) var_e=r(var_e), /// reps(10) seed(123) cluster(id) idcluster(newid) group(id) : mybs2
Comment