*================================================================================= * Check if the ideas presented in * http://blog.stata.com/2011/08/22/use-poisson-rather-than-regress-tell-a-friend/ * also apply to -propcnsreg- * not working as well as I would have liked/hoped, especially the standard errors * * also requires -simsum- and -simpplot-, both can be installed using -ssc- *================================================================================= clear all program define sim, rclass drop _all set obs 5000 gen z1 = rnormal() gen z2 = rnormal() gen lat = .25*z1 + .5*z2 gen x1 = rnormal() gen x2 = rnormal() gen lam = 1 + .5*x1 gen y = exp(.25*x1 + .25*x2 + lam*lat + .2*rnormal()) propcnsreg y x1 x2 , constr(z1 z2) lambda(x1) lcons poisson robust return scalar w = e(w_p) return scalar bz1 = [constrained]_b[z1] return scalar bz2 = [constrained]_b[z2] return scalar blx1 = [lambda]_b[x1] return scalar bx1 = [unconstrained]_b[x1] return scalar sez1 = [constrained]_se[z1] return scalar sez2 = [constrained]_se[z2] return scalar selx1 = [lambda]_se[x1] return scalar sex1 = [unconstrained]_se[x1] end set seed 123456 simulate w=r(w) /// bz1 = r(bz1) sez1 = r(sez1) /// bz2 = r(bz2) sez2 = r(sez2) /// blx1 = r(blx1) selx1 = r(selx1) /// bx1 = r(bx1) sex1 = r(sex1), /// reps(5000):sim simsum bz1 , se(sez1) true(.25) mcse simsum bz2 , se(sez2) true(.50) mcse simsum blx1, se(selx1) true(.50) mcse simsum bx1 , se(sex1) true(.25) mcse gen pz1 = 2*normal(-abs(( bz1-.25)/sez1 )) gen pz2 = 2*normal(-abs(( bz2-.50)/sez2 )) gen plx1 = 2*normal(-abs((blx1-.50)/selx1)) gen px1 = 2*normal(-abs(( bx1-.25)/sex1 )) simpplot w pz1 pz2 plx1 px1, scheme(s2color) name(pois_coef, replace) exit