*! 1.0.0 MLB 13Okt2010 program define scenreg version 11 if replay() { if (`"`e(cmd)'"' != "sensreg") error 301 Replay `0' } else Estimate `0' end program Estimate, eclass syntax varlist [if] [in] /// [fweight pweight], /// sd(string) /// [ /// link(string) /// family(string) /// undocumented, for now only binomial allowed dist(string) /// rho(string) /// eclear /// noLOg noCONStant /// VCE(passthru) /// Robust /// old options CLuster(passthru) /// old options draws(integer 100) /// start(integer 15) /// Level(integer `c(level)') /// OR HR RR /// * /// ] macro drop S_* gettoken y x : varlist // parse vce local vceopt = `:length local vce' | /// `:length local weight' | /// `:length local cluster' | /// `:length local robust' if `vceopt' { _vce_parse, argopt(CLuster) opt(OIM OPG Robust) old /// : [`weight'`exp'], `vce' `robust' `cluster' local vce if "`r(cluster)'" != "" { local clustvar `r(cluster)' local vce vce(cluster `r(cluster)') } else if "`r(robust)'" != "" { local vce vce(robust) } else if "`r(vce)'" != "" { local vce vce(`r(vce)') } } // mark estimation sample marksample touse if `:length local clustvar' { markout `touse' `clustvar', strok } qui count if `touse' if r(N) == 0 { error 2000 } // weights if "`weight'" != "" local wgt `"[`weight'`exp']"' if "`weight'" == "pweight" local sumwgt `"[aweight`exp']"' // parse sd tempvar sdvar capture gen double `sdvar' = `sd' if _rc { di as err "invalid expression specified in option sd()" exit _rc } capture assert `sdvar' >= 0 & `sdvar' <= 20 if _rc { di as err "sd must result in numbers between 0 and 20" exit 198 } global S_sd "`sdvar'" // parse dist gettoken dist pr : dist local l = max(4, length("`dist'")) local dist = lower("`dist'") // makes sense to allow one to specify the Gaussian distribution // instead of the gaussian distribution if "`dist'" == substr("normal", 1, `l') | /// "`dist'" == substr("gaussian", 1, `l') | /// "`dist'" == "" { global S_dist normal } else if "`dist'" == substr("uniform", 1, `l') { global S_dist uniform } else if "`dist'" == substr("discrete", 1, `l') { global S_dist discrete // parse the probabilities and attach the mass points local ok = 1 if `: word count `pr'' < 2 { local ok = 0 } foreach p of local pr { capture confirm number `p' if _rc { local ok = 0 } if `p' < 0 { local ok = 0 } } tempname mpnts // This function creates equally spaced masspoints such that mean(e) = 0 & sd(e) = 1 // Leaves these masspoints behind in $S_mpnts and matrix `mpnts' // Leaves pr behind in $S_pr // Leaves the sum of pr behind in local sum_pr if `ok' { mata masspoints() } if abs(`sum_pr' - 1) > .01 { local ok = 0 } if !`ok' { di as err "when specifying discrete in dist() it must be followed by two or more positive numbers adding up to 1" exit 198 } // make matrix pretty for later display matrix rownames `mpnts' = "mass_point" "proportion" forvalues i = 1/ `:word count `pr'' { local coln = `"`coln' "point_`i'""' } matrix colnames `mpnts' = `coln' local mpntsopt "mpnts(`mpnts')" } else { di as err "only normal, gaussian, uniform, or discrete are allowed in -dist()-" exit 198 } // parse rho if "`rho'" != "" { if `: word count `rho'' != 2 { di as err "the rho option must contain one variable and one correlation" exit 198 } gettoken var rho : rho if !`: list var in x' { di as err "the first element in rho must be one of the explanatory variables" exit 198 } global S_rhovar "`var'" capture confirm number `rho' if _rc { di as err "the second element in rho must be one number between -.9 and .9" exit 198 } if `rho' < -.9 | `rho' > .9 { di as err "the second element in rho must be one number between -.9 and .9" exit 198 } global S_rho `rho' qui sum $S_rhovar if `touse' `sumwgt' global S_sdx = r(sd) global S_mx = r(mean) } else { global S_rho = . global S_sdx = . global S_mx = . } // check draws and start if `draws' <= 0 | `draws' > 1000 { di as err "draws must be an integer between 0 and 1000" exit 198 } global S_draws = `draws' if `start' < 0 { di as err "start must be a positive integer" exit 198 } global S_start = `start' //parse link local l = max(4,length("`link'")) if "`link'" == "" | "`link'" == substr("logit",1,`l') { global S_link logit } else if "`link'" == substr("identity", 1, `l'){ global S_link identity } else if "`link'" == substr("probit",1,`l') { global S_link probit } else if "`link'" == substr("cloglog", 1, `l') { global S_link cloglog } else if "`link'" == substr("loglog",1,`l') { global S_link loglog } else if "`link'" == "log" { global S_link log } else { di as err "unrecognized link: `link'" exit 198 } // parse eform options opts_exclusive "`hr' `or' `rr'" // make sure the name of the eform option makes sense given the link function if "`hr'" != "" & "$S_link" != "cloglog" { di as err "the hr option is only allowed in combination with the cloglog link function" exit 198 } if "`or'" != "" & "$S_link" != "logit" { di as err "the or option is only allowed in combination with the logit link function" exit 198 } if "`rr'" != "" & "$S_link" != "log" { di as err "the rr option is only allowed in combination with the log link function" exit 198 } if "`rr'" != "" { local `rr' "eform(Risk Ratio)" } /* for now only binomial family allowed reason is that the standard errors in the other families did not perform well in simulations if "`family'" == "" { if "$S_link" == "logit" | "$S_link" == "probit" | "$S_link" == "cloglog" | "$S_link" == "loglog" { global S_family binomial } if "$S_link" == "log" { global S_family poisson } if "$S_link" == "identity" { global S_family gaussian } } local l = max(4, length("`family'")) else if "`family'" == substr("gaussian", 1, `l') | /// "`family'" == substr("normal", 1, `l') { global S_family gaussian local sigma /ln_sigma } else if "`family'" == substr("binomial", 1, `l') { global S_family binomial } else if "`family'" == substr("poisson", 1, `l') { global S_family poisson } else { di as err "unrecognized family: `family'" exit 198 } */ global S_family binomial // check and collect -ml- and display options _get_diopts diopts options, `options' local diopts `diopts' level(`level') `or' `hr' `rr' mlopts mlopts, `options' // display/suppress itteration log local log = cond("`log'" == "", "noisily", "quietly") qui count if `touse' local N = r(N) // initial values qui glm `y' `x' `wgt' if `touse', `constant' link($S_link) family($S_family) iter(20) tempname b0 matrix `b0' = e(b) matrix coleq `b0' = "" // make the samples from the "unobserved variable" mata: mk_e(`N', $S_draws, $S_start, "$S_sd", "`touse'", "$S_dist", /// "$S_pr", "$S_mpnts", $S_rho, "$S_rhovar", $S_sdx, $S_mx, "`eclear'") // full model `log' ml model lf scenreg_lf /// (xb : `y' = `x', `constant') /// `sigma' /// will only occur with the normal family, which is now not allowed `wgt' if `touse', /// `mlopts' `vce' maximize init(`b0') search(off) ereturn local title "scenario model, link($S_link) family($S_family)" ereturn local link "$S_link" ereturn local family "$S_family" ereturn local cmd "scenreg" ereturn local sd "`sd'" if "$S_rhovar" != "" { ereturn scalar rho = $S_rho ereturn local rhovar "$S_rhovar" } ereturn local dist "$S_dist" if "$S_dist" == "discrete" { ereturn matrix mpnts `mpnts' } // clean up mata: rmexternal("S_unobserved_variable") Replay , `diopts' end program Replay syntax [, sd(string) mpnts(name) *] ml display, `options' di _n di as txt "Scenario:" di as txt "Unobserved variable is " as result "${S_dist}ly" as txt " distributed" if "$S_dist" == "discrete" { di as txt "The mass points and probabilities of the standardized unobserved variable are:" matlist e(mpnts), format(%7.3g) di _n } di as txt "The effect of the standardized unobserved variable is " as result "`e(sd)'" if "`e(rhovar)'" != "" { di as txt "The corelation between " as result "`e(rhovar)'" as txt " and the unobserved variable is " as result `e(rho)' } end /*----------- computes mass points when the discrete distribution is specified*/ mata: void masspoints() { real matrix pr, mpnts real scalar Np, i string scalar str_mpnts pr = strtoreal(tokens(st_local("pr"))) // initial guesses for mpnts // centered around 0, each category 1 appart Np = length(pr) mpnts = (1..Np):- (Np+1)/2 // pr*mpnts' is the mean mpnts = mpnts :- (pr*mpnts') // pr*(mpnts :* mpnts)' is the variance mpnts = mpnts :/ sqrt(pr*(mpnts :* mpnts)') // return results str_mpnts = strofreal(mpnts[1]) for(i=2; i<= length(mpnts); i++){ str_mpnts = str_mpnts + " " + strofreal(mpnts[i]) } st_global("S_mpnts", str_mpnts) st_global("S_pr", st_local("pr")) st_matrix(st_local("mpnts"), mpnts \ pr) // return sum of pr to check if adds up to 1 st_local("sum_pr", strofreal(sum(pr))) } // ---------------------------- make unobserved variable void mk_e( real scalar N, real scalar draws, real scalar start, string scalar sdname, string scalar mlsamp, string scalar dist, string scalar prstr, string scalar mpntsstr, real scalar rho, string scalar rhovar, real scalar sdx, real scalar mx, string scalar eclear) { pointer(real matrix) scalar e if ( (e = findexternal("S_unobserved_variable")) == NULL) { e = crexternal("S_unobserved_variable") } else if (eclear == "") { _error("external global matrix S_unobserved_variables already exist, specify the eclear option if you want to replace it") } else { rmexternal("S_unobserved_variable") e = crexternal("S_unobserved_variable") } real matrix ee, sd // create standardized distribution if (dist == "uniform") { ee = ( (halton(draws*N, 1, start) :- .5) :/ sqrt(1/12) ) } else if (dist == "normal") { ee = invnormal(halton(draws*N, 1, start)) } else if (dist == "discrete") { real matrix u, pr, mpnts real scalar cumpr, prev_cumpr, i pr = strtoreal(tokens(prstr)) mpnts = strtoreal(tokens(mpntsstr)) u = halton(draws*N, 1, start) ee = J(draws*N,1,0) cumpr = 0 for (i = 1; i <= length(pr); i++){ prev_cumpr = cumpr cumpr = cumpr + pr[i] ee = ee :+ ( (u:=prev_cumpr) ) :* mpnts[i] } } // add correlation ee = colshape(ee,draws) if (rhovar != "") { real matrix x st_view(x,.,rhovar,mlsamp) ee = (rho:*(x :- mx):/sdx :+ sqrt(1-rho^2):*ee) } // add standard deviation st_view(sd, ., sdname, mlsamp) ee = sd :* ee *e = ee } end