*! version 1.0.3 17may2019 Ben Jann program robbox, eclass version 11 if replay() { if `"`e(cmd)'"'!="robbox" { di as err "last robbox results not found" exit 301 } Replay `0' exit } Estimate `0' ereturn local cmdline `"robbox `0'"' Replay, `table' `graph' `label' `postoutsides' `flag' `generate' `replace' `gropts' end program Replay syntax [, TABle NOGRaph GRaph noLABel POSTOUTsides flag FLAG2(passthru) /// GENerate GENerate2(passthru) replace * ] if "`table'"!="" { Display local nograph nograph } if "`postoutsides'"!="" { Postoutsides local nograph nograph } if `"`flag'`flag2'"'!="" { Flag, vlist(`e(depvar)') `flag2' `replace' local nograph nograph } if `"`generate'`generate2'"'!="" { Generate, `generate2' `replace' `label' local nograph nograph } if "`graph'"!="" | "`nograph'"=="" { Graph, `options' `label' } end program Display // header di as txt _n "`e(title)'" _c di as txt _col(46) "Number of obs = " as res %10.0g e(N) if `"`e(alpha)'"'!="" { di as txt _col(46) "Alpha = " as res %10.0g e(alpha) } if `"`e(bp)'"'!="" { di as txt _col(46) "Breakdown = " as res %10.0g e(bp) } if `"`e(delta)'"'!="" { di as txt _col(46) "Delta = " as res %10.0g e(delta) } // over() legend if `"`e(over)'"'!="" { _svy_summarize_legend if e(N_vars)==1 { local rowti rowtitle(`e(depvar)') } } else di "" // table tempname tmp tmp1 matrix `tmp' = e(b)' mat coln `tmp' = "_:median" mat `tmp1' = e(box)' mat coln `tmp1' = "lower:hinge" "upper:hinge" matrix `tmp' = `tmp', `tmp1' mat `tmp1' = e(whiskers)' mat coln `tmp1' = "lower:whisker" "upper:whisker" matrix `tmp' = `tmp', `tmp1' mat `tmp1' = e(N_out)' mat coleq `tmp1' = "number of outside values" matrix `tmp' = `tmp', `tmp1' mat `tmp1' = e(_N)' mat coln `tmp1' = "mumber:of_obs" matrix `tmp' = `tmp', `tmp1' matlist `tmp', twidth(16) lines(oneline) border(rows) keepcoleq nohalf /// underscore showcoleq(rcombined) noblank `rowti' end program Graph // syntax syntax [, noOUTsides VERTical HORizontal /// sort SORT2(numlist max=1 int >0) DEScending /// BOXWidth(numlist max=1) OVERGap(numlist max=1) /// noLABel LABels(str asis) LLABels(str asis) MEDMarker /// PLOTopts(str) MEDopts(str) BOXopts(str) WHISKopts(str) OUTopts(str) /// PCYCle(numlist max=1 integer >0) addplot(str asis) * ] if "`sort2'"!="" local sort sort if "`boxwidth'"=="" local boxwidth 50 if "`overgap'"=="" local overgap 10 if "`pcycle'"=="" local pcycle 15 local orientation `vertical' `horizontal' if `:list sizeof orientation'>1 { di as err "only one of vertical and horizontal allowed" } if "`orientation'"=="" local orientation vertical // sort results if "`sort'"!="" { tempname ecurrent _estimates hold `ecurrent', copy restore Sort "`descending'" "`sort2'" } // generate boxplot variables local nvars = e(N_vars) local nover = e(N_over) local n = `nvars' * `nover' if `n'>_N { preserve set obs `n' } tempvar gid id med lbox ubox lav uav nobs qui Generate, `label' generate(`gid' `id' `med' `lbox' `ubox' `lav' `uav' `nobs') // prepare plotting positions and collect labels local over `e(over)' if "`over'"=="" | ("`over'"!="" & `nvars'==1) { // no groups local ng `n' local nj 1 local AT `id' local tmp `gid' local gid `id' local id `tmp' local llab off local barwidth = 1/`nj'*(`boxwidth'/100) local range .5 `ng'.5 } else { local gw = 1 - `overgap'/100 local ng = e(k_eq) local nj = `n' / `ng' tempvar AT qui gen `AT' = `gid' - (.5*`gw') + ((`id'-1)*2+1) * ((1*`gw')/(2*`nj')) in 1/`n' local key 2 local llab forv i = 1/`nj' { gettoken lbl llabels : llabels if `"`lbl'"'=="" { local lbl: label `id' `i' } local llab `llab' `key' `"`lbl'"' local key = `key' + 3 } local llab order(`llab') local barwidth = (1*`gw')/`nj'*(`boxwidth'/100) local range `=1-.5*`gw'' `=`ng'+.5*`gw'' } local clab forv i = 1/`ng' { gettoken lbl labels : labels if `"`lbl'"'=="" { local lbl: label `gid' `i' } local clab `clab' `i' `"`lbl'"' } // generate outlier variables if "`outsides'"=="" { local vlist "`e(depvar)'" local flags local flagats forv i=1/`nvars' { tempname flag flagat local flags `flags' `flag' local flagats `flagats' `flagat' } qui Flag, flag(`flags') vlist(`vlist') local i 0 foreach flag of local flags { local ++i local flagat: word `i' of `flagats' if "`over'"=="" { qui gen byte `flagat' = `i' if `flag'==1 } else { qui gen byte `flagat' = . local j 0 foreach o in `e(over_namelist)' { local ++j local at `j' if `nj'>1 { local at = `at' - (.5*`gw') + ((`i'-1)*2+1) * ((1*`gw')/(2*`nj')) } qui replace `flagat' = `at' if `flag'==1 & `over'==`o' } } } } // compile plot command // - collect numbered options forv i = 1/`nj' { Collectnumopts `i', `options' // returns options, plot#opts, med#opts, etc. } // - median, box, and whiskers if "`medmarker'"!="" { if "`orientation'"=="horizontal" local medplot scatter `AT' `med' else local medplot scatter `med' `AT' } else { tempname AT0 AT1 qui gen `AT0' = `AT' - `barwidth'/2 qui gen `AT1' = `AT' + `barwidth'/2 local medplot pcspike `med' `AT0' `med' `AT1' local medorient `orientation' } local barwidth barwidth(`barwidth') local plots forv i = 1/`nj' { local p = mod(`i'-1, `pcycle') + 1 local plots `plots'/* */ (rcap `lav' `uav' `AT' if `id'==`i', pstyle(p`p') /* */ `plotopts' `plot`i'opts' `orientation' `whiskopts' `whisk`i'opts')/* */ (rbar `lbox' `ubox' `AT' if `id'==`i', pstyle(p`p') fintensity(50) /* */ `plotopts' `plot`i'opts' `barwidth' `orientation' `boxopts' `box`i'opts') /* */ (`medplot' if `id'==`i' & `nobs'!=0, pstyle(p`p') /* */ `plotopts' `plot`i'opts' `medorient' `medopts' `med`i'opts') } // - outside values if "`outsides'"=="" { local i 0 local j 1 foreach flag of local flags { local ++i local flagat: word `i' of `flagats' local v: word `i' of `vlist' if "`orientation'"=="horizontal" local xy `flagat' `v' else local xy `v' `flagat' local p = mod(`j'-1, `pcycle') + 1 local plots `plots' (scatter `xy' if `flag'==1, pstyle(p`p') /* */ `plotopts' `plot`j'opts' `outopts' `out`j'opts') if `nj'>1 { local ++j } } } // - addpot if `"`addplot'"'!="" { local plots `plots' || `addplot' || } // - global options local grid forv i = 1/`=`ng'-1' { local grid `grid' `i'.5 } if "`orientation'"=="horizontal" { local axisopts ylabel(`clab', angle(0) nogrid) yti("") /// ytick(`grid', notick) yscale(range(`range') reverse) } else { local axisopts xlabel(`clab', nogrid) xti("") xtick(`grid', notick) /// xscale(range(`range')) } // draw graph twoway `plots', `axisopts' legend(`llab' all) `options' end program Sort, eclass args descending index if "`index'"=="" local index 0 tempname B mata: robbox_sort("`B'", `index', "`descending'"!="") end program Collectnumopts _parse comma i 0 : 0 syntax [, PLOT`i'opts(str) MED`i'opts(str) BOX`i'opts(str) /// WHISK`i'opts(str) OUT`i'opts(str) * ] c_local plot`i'opts `"`plot`i'opts'"' c_local med`i'opts `"`med`i'opts'"' c_local box`i'opts `"`box`i'opts'"' c_local whisk`i'opts `"`whisk`i'opts'"' c_local out`i'opts `"`out`i'opts'"' c_local options `options' end program Generate syntax [, noLABel checkonly GENerate2(str) replace ] local generate `"`generate2'"' // variables names if `"`generate'"'=="" local generate _box_* if `:list sizeof generate'==1 /// & substr(`"`generate'"',-1,1)=="*" /// & strlen(`"`generate'"')>1 { local prefix = substr(`"`generate'"', 1, strlen(`"`generate'"')-1) local generate foreach v in gid id med lbox ubox lav uav N { local generate `generate' `prefix'`v' } } else { if `:list sizeof generate'!=8 { di as err "generate() must either contain 8 names" /// " or a prefix specified as {it:prefix}{bf:*}" exit 198 } local tmp: list uniq generate if `:list sizeof tmp'!=8 { di as err "generate(): variable names must be unique" exit 198 } } confirm names `generate' // check variables if "`replace'"=="" { confirm new variable `generate' } if `"`checkonly'"'!="" exit if "`replace'"!="" { local i 0 foreach v of local generate { local ++i capt confirm new variable `v' if _rc drop `v' if `i'>2 continue capt label drop `v' // gid and id only } } // generate variables gettoken gid tmp : generate gettoken id tmp : tmp gettoken med tmp : tmp gettoken lbox tmp : tmp gettoken ubox tmp : tmp gettoken lav tmp : tmp gettoken uav tmp : tmp gettoken nobs : tmp // - IDs local n = e(N_vars) * e(N_over) local neq = e(k_eq) local nj = `n' / `neq' qui gen `gid'= ceil(_n/`nj') in 1/`n' qui gen `id' = 1 in 1/`n' qui replace `id' = cond(`gid'==`gid'[_n-1], `id'+`id'[_n-1], `id') in 1/`n' // - main data tempname tmp mat `tmp' = e(b)', e(box)', e(whiskers)', e(_N)' mat coln `tmp' = `med' `lbox' `ubox' `lav' `uav' `nobs' svmat `c(type)' `tmp', names(col) // - compress qui compress `generate' // value labels for IDs if (`"`e(over)'"'!="" & e(N_vars)==1) { if "`label'"=="" local lbls `"`e(over_labels)'"' else local lbls `"`e(over_namelist)'"' local i 0 foreach lbl of local lbls { local ++i lab def `id' `i' `"`lbl'"', add } lab def `gid' 1 "_", add } else { local lbls `"`e(depvar)'"' local i 0 foreach lbl of local lbls { local ++i if "`label'"=="" { local vlab: var lab `lbl' if `"`vlab'"'!="" { local lbl `"`vlab'"' } } lab def `id' `i' `"`lbl'"', add } if `"`e(over)'"'!="" { if "`label'"=="" local lbls `"`e(over_labels)'"' else local lbls `"`e(over_namelist)'"' local i 0 foreach lbl of local lbls { local ++i lab def `gid' `i' `"`lbl'"', add } } else { lab def `gid' 1 "_", add } } lab val `gid' `gid', nofix lab val `id' `id', nofix // variable labels and display lab var `gid' "Group (equation) ID" lab var `id' "Within-group (coefficient) ID" lab var `med' "Median" lab var `lbox' "Lower hinge (box)" lab var `ubox' "Upper hinge (box)" lab var `lav' "Lower adjacent value (whisker)" lab var `uav' "Upper adjacent value (whisker)" lab var `nobs' "Number of observations" display _n as txt "Generated boxplot variables:" _c describe `generate' end program Flag syntax [, checkonly vlist(str) FLAG2(str) replace ] local generate `"`flag2'"' local nvars: list sizeof vlist // variables names if `"`generate'"'=="" local generate _flag_* if `:list sizeof generate'==1 /// & substr(`"`generate'"',-1,1)=="*" /// & strlen(`"`generate'"')>1 { local prefix = substr(`"`generate'"', 1, strlen(`"`generate'"')-1) local generate foreach v of local vlist { local generate `generate' `prefix'`v' } } else { if `:list sizeof generate'!=`nvars' { di as err "flag() must either contain `nvars' name(s)" /// " or a prefix specified as {it:prefix}{bf:*}" exit 198 } } confirm names `generate' local tmp: list uniq generate if `:list sizeof tmp'!=`nvars' { di as err "flag(): variable names must be unique" exit 198 } // check variables if "`replace'"=="" { confirm new variable `generate' } if `"`checkonly'"'!="" exit if "`replace'"!="" { foreach v of local generate { capt confirm new variable `v' if _rc drop `v' } } // generate variables local varlist `vlist' foreach vnew of local generate { gettoken v varlist : varlist qui gen byte `vnew' = . lab var `vnew' "`v' outside value" } tempname lav uav local i 0 if `"`e(over)'"'!="" { local over `e(over)' local overvals `e(over_namelist)' foreach o of local overvals { local varlist `vlist' foreach vnew of local generate { gettoken v varlist : varlist local ++i scalar `lav' = el(e(whiskers), 1, `i') scalar `uav' = el(e(whiskers), 2, `i') qui replace `vnew' = (`v'<`lav') | (`v'>`uav') /// if `over'==`o' & `v'<. & e(sample) } } } else { local varlist `vlist' foreach vnew of local generate { gettoken v varlist : varlist local ++i scalar `lav' = el(e(whiskers), 1, `i') scalar `uav' = el(e(whiskers), 2, `i') qui replace `vnew' = (`v'<`lav') | (`v'>`uav') if `v'<. & e(sample) } } display _n as txt "Generated flags:" _c describe `generate' end program Postoutsides, eclass tempname out tmp local vlist "`e(depvar)'" local nvars: list sizeof vlist local flags forv i=1/`nvars' { tempname flag local flags `flags' `flag' } qui Flag, flag(`flags') vlist(`vlist') local over `e(over)' if "`over'"=="" { foreach flag of local flags { gettoken v vlist : vlist capt mkmat `v' if `flag'==1, matrix(`tmp') if _rc==2000 continue else error _rc mat rown `tmp' = "`v'" mat `out' = nullmat(`out'), `tmp'' } } else { local overvals `e(over_namelist)' if `nvars'==1 { foreach o of local overvals { capt mkmat `vlist' if `flag'==1 & `over'==`o', matrix(`tmp') if _rc==2000 continue else error _rc mat rown `tmp' = "`o'" mat `out' = nullmat(`out'), `tmp'' } } else { foreach o of local overvals { local varlist `vlist' foreach flag of local flags { gettoken v varlist : varlist capt mkmat `v' if `flag'==1 & `over'==`o', matrix(`tmp') if _rc==2000 continue else error _rc mat rown `tmp' = "`o':`v'" mat `out' = nullmat(`out'), `tmp'' } } } } capt confirm matrix `out' if _rc { di as txt "no outside values; did not create matrix {bf:e(outsides)}" exit } mat rown `out' = "y1" // same as in e(b) eret matrix outsides = `out' di as txt "outside values added in matrix {bf:e(outsides)}" end program Estimate, eclass // syntax syntax varlist(numeric) [if] [in] [pw iw fw] [, over(varname numeric) cw /// STandard ADJusted general alpha(numlist max=1 >0 <25) /// bp(numlist max=1 >0 <25) delta(numlist max=1 >0) /// noTABle noGRaph noLABel POSTOUTsides /// flag FLAG2(passthru) GENerate GENerate2(passthru) replace * ] // - passthru options if `"`generate'`generate2'"'!="" & "`replace'"=="" { Generate, checkonly `generate2' } if `"`flag'`flag2'"'!="" & "`replace'"=="" { Flag, checkonly vlist(`varlist') `flag2' } if "`table'"=="" c_local table table if "`graph'"=="" c_local graph graph else c_local graph nograph c_local label `label' c_local flag `flag' `flag2' c_local postoutsides `postoutsides' c_local generate `generate' `generate2' c_local replace `replace' c_local gropts `"`options'"' // - other options local method `standard' `adjusted' `general' if `:list sizeof method'>1 { di as err "only one of standard, adjusted, and general allowed" exit 198 } if "`method'"=="" local method "general" if "`method'"=="adjusted" { if "`alpha'"!="" { di as err "{bf:alpha()} not allowed with {bf:`method'}" exit 198 } } else { tempname ALPHA if "`alpha'"!="" { scalar `ALPHA' = `alpha' } else { scalar `ALPHA' = 200 * (1 - normal(invnormal(.75) + /// 1.5 * (invnormal(.75)-invnormal(.25)))) } } if "`method'"!="general" { if "`bp'"!="" { di as err "{bf:bp()} not allowed with {bf:`method'}" exit 198 } if "`delta'"!="" { di as err "{bf:delta()} not allowed with {bf:`method'}" exit 198 } } else { if "`bp'"=="" local bp 10 if "`delta'"=="" local delta 0.1 } // sample if "`cw'"!="" marksample touse else { marksample touse, novarlist tempvar touse1 qui gen byte `touse1' = 0 foreach v of local varlist { qui replace `touse1' = 1 if `touse' & `v'<. } qui replace `touse' = 0 if `touse1'==0 } markout `touse' `over' _nobs `touse' [`weight'`exp'], min(1) local N = r(N) // prepare containers for results local nvars: list sizeof varlist tempname b box whisk whisk0 Nout _N mat `box' = J(2,`nvars', .) mat coln `box' = `varlist' mat `b' = `box'[1,1...] mat `_N' = `b' mat `whisk' = `box' mat `whisk0' = `box' mat `Nout' = `box' \ `b' mat rown `Nout' = "lower" "upper" "total" if "`method'"=="adjusted" { tempname mc mat `mc' = `b' } else if "`method'"=="general" { tempname g h mat `g' = `b' mat `h' = `b' } // handle over() if "`over'"!="" { capt assert ((`over'==floor(`over')) & (`over'>=0)) if `touse' if _rc { di as err "variable in over() must be integer and nonnegative" exit 452 } qui levelsof `over' if `touse', local(overvals) local N_over: list sizeof overvals local over_labels foreach overval of local overvals { local over_labels `"`over_labels' `"`: label (`over') `overval''"'"' } local over_labels: list clean over_labels local rmat b box whisk whisk0 Nout _N if "`method'"=="adjusted" local rmat `rmat' mc if "`method'"=="general" local rmat `rmat' g h tempname tmp foreach mat of local rmat { matrix rename ``mat'' `tmp' foreach overval in `overvals' { mat coleq `tmp' = "`overval'" mat ``mat'' = nullmat(``mat''), `tmp' } mat drop `tmp' } } else { local total total local N_over 1 } // compute results local i = 0 local j = 0 foreach overval in `overvals' `total' { local ++j if `"`overval'"'=="total" local ovr 1 else local ovr `over'==`overval' foreach v of local varlist { local ++i _Estimate `v' if `touse' & `ovr' [`weight'`exp'], /// method(`method') alpha(`ALPHA') bp(`bp') delta(`delta') mat `_N'[1,`i'] = r(N) mat `b'[1,`i'] = r(p50) mat `box'[1,`i'] = r(p25) \ r(p75) mat `whisk'[1,`i'] = r(lav) \ r(uav) mat `whisk0'[1,`i'] = r(lb) \ r(ub) mat `Nout'[1,`i'] = r(lout) \ r(uout) \ r(tout) if "`method'"=="adjusted" { mat `mc'[1,`i'] = r(mc) } else if "`method'"=="general" { mat `g'[1,`i'] = r(g) mat `h'[1,`i'] = r(h) } } } // flip labels if over() and only one variable if "`over'"!="" & `nvars'==1 { local colnm: coleq `b' foreach mat of local rmat { matrix coleq ``mat'' = "" matrix coln ``mat'' = `colnm' } local k_eq = 1 } else { local k_eq = `N_over' } // post results eret post `b' [`weight'`exp'], obs(`N') esample(`touse') eret local cmd "robbox" local title "box plot" if "`method'"=="general" eret local title "Generalized box plot" else if "`method'"=="adjusted" eret local title "Adjusted box plot" else eret local title "Standard box plot" eret local method "`method'" eret local depvar "`varlist'" eret local over_labels `"`over_labels'"' eret local over_namelist `"`overvals'"' eret local over "`over'" eret scalar N_vars = `nvars' eret scalar N_over = `N_over' eret scalar k_eq = `k_eq' eret matrix _N = `_N' if "`method'"=="adjusted" { eret matrix medcouple = `mc' } else if "`method'"=="general" { eret matrix h = `h' eret matrix g = `g' } eret matrix N_out = `Nout' eret matrix whiskers0 = `whisk0' eret matrix whiskers = `whisk' eret matrix box = `box' if "`ALPHA'"!="" eret scalar alpha = `ALPHA' if "`bp'"!="" eret scalar bp = `bp' if "`delta'"!="" eret scalar delta = `delta' Return_clear // clear r() end program Return_clear, rclass local x end program _Estimate, rclass // syntax syntax varname(numeric) if [pw iw fw/] [, method(str) alpha(str) bp(str) delta(str) ] // weights and number of obs if "`weight'"!="" { local wgt [`weight' = `exp'] if "`weight'"=="iweight" local pcwgt [aw = `exp'] else local pcwgt `wgt' } if "`weight'"!="fweight" { qui count `if' & `varlist'<. local N = r(N) } else { su `varlist' `wgt' `if', meanonly local N = r(N) } if `N'==0 { return scalar N = `N' return scalar p50 = 0 // (because e(b) cannot have missings) exit } // compute quartiles if "`method'"=="general" local percentiles 10 25 50 75 90 else local percentiles 25 50 75 _pctile `varlist' `pcwgt' `if', percentiles(`percentiles') local i 0 foreach p of local percentiles { local ++i tempname p`p' scalar `p`p'' = r(r`i') } // compute minimum lower bound and maximum upper bound of whiskers tempname lb ub if "`method'"=="standard" { tempname c scalar `c' = (invnormal(1-`alpha'/200) - invnormal(.75)) /// / (invnormal(.75) - invnormal(.25)) scalar `lb' = `p25' - `c' * (`p75'-`p25') scalar `ub' = `p75' + `c' * (`p75'-`p25') } else if "`method'"=="adjusted" { tempname MC quietly robstat `varlist' `if' `wgt', statistic(MC) nose scalar `MC' = _b[MC] if `MC'>=0 { scalar `lb' = `p25' - 1.5 * exp(-4*`MC') * (`p75'-`p25') scalar `ub' = `p75' + 1.5 * exp(3*`MC') * (`p75'-`p25') } else { scalar `lb' = `p25' - 1.5 * exp(-3*`MC') * (`p75'-`p25') scalar `ub' = `p75' + 1.5 * exp(4*`MC') * (`p75'-`p25') } } else if "`method'"=="general" { tempvar v tempname scale rmin rscale qmed qscale g h lb0 ub0 // step 1: transform the data // - (a) standardize, so that delta in (b) is relative to the scale scalar `scale' = (`p75' - `p25') / (0.6817766 /// + 0.0534282 * abs((`p90' + `p10' - 2*`p50') / (`p90' - `p10')) /// + 0.1794771 * ((`p90' - `p10') / (`p75' - `p25')) /// - 0.0059595 * ((`p90' - `p10') / (`p75' - `p25'))^2) // // note: sk = (p90+p10-2*p50)/(p90-p10), kurt = (p90-p10)/(p75-p25) // why not simply use scale = (p75-p25) as in paper? qui gen `v' = (`varlist' - `p50') / `scale' `if' // - (b) shift to positive domain and fit into (0,1) su `v' `if', mean scalar `rmin' = r(min) scalar `rscale' = 2 * `delta' + r(max) - `rmin' // note: rscale is equal to min+max from shifted data qui replace `v' = (`v' - `rmin' + `delta') / `rscale' `if' // note: resulting data will be in [d, 1-d] with // d = `delta' / (2*`delta' + range/`scale') // were range = max-min of the original data // - (c) transform to quantiles of the standard normal qui replace `v' = invnormal(`v') `if' // - (d) standardize the quantiles _pctile `v' `pcwgt' `if', percentiles(25 50 75) scalar `qmed' = r(r2) scalar `qscale' = ((r(r3) - r(r1)) / (invnormal(0.75)-invnormal(0.25))) qui replace `v' = (`v' - `qmed') / `qscale' // step 2: fit the g and h distribution local plo = `bp' local pup = 100 - `bp' _pctile `v' `pcwgt' `if', percentiles(`plo' `pup') scalar `g' = ln(-r(r2)/r(r1)) / invnormal(`pup'/100) // note: (r(r2)-med)/(med-r(r1)) = -r(r2)/r(r1) because med = 0 scalar `h' = 2 / invnormal(`pup'/100)^2 * ln(-`g' * ((r(r1) * r(r2)) / (r(r1) + r(r2)))) // step 3: get outlier limits (alpha quantiles) of the fitted g-and-h tempname xi_lo xi_up scalar `lb0' = (1/`g') * (exp(`g' * invnormal(`alpha'/200)) - 1) * /// exp(`h' * invnormal(`alpha'/200)^2 / 2) scalar `ub0' = (1/`g') * (exp(`g' * invnormal(1-`alpha'/200)) - 1) * /// exp(`h' * invnormal(1-`alpha'/200)^2 / 2) // step 4: back-transform the limits to the original scale scalar `lb' = `p50' + `scale' * /// ((normal(`lb0'*`qscale' + `qmed') * `rscale') + `rmin' - `delta') scalar `ub' = `p50' + `scale' * /// ((normal(`ub0'*`qscale' + `qmed') * `rscale') + `rmin' - `delta') // step 5: make sure the limits are not inside the box scalar `lb' = min(`p25', `lb') scalar `ub' = max(`p75', `ub') } // compute lower and upper adjacent values (empirical endpoints if whiskers) tempname lav uav summarize `varlist' `if' & `varlist'>=`lb' & `varlist'<=`ub', mean scalar `lav' = min(`p25', r(min)) scalar `uav' = max(`p75', r(max)) // count outside values if "`weight'"!="fweight" { qui count `if' & `varlist'<`lav' & `varlist'<. local lout = r(N) qui count `if' & `varlist'>`uav' & `varlist'<. local uout = r(N) } else { su `varlist' `wgt' `if' & `varlist'<`lav', meanonly local lout = r(N) su `varlist' `wgt' `if' & `varlist'>`uav', meanonly local uout = r(N) } // return results return scalar N = `N' return scalar p25 = `p25' return scalar p50 = `p50' return scalar p75 = `p75' return scalar lb = `lb' return scalar ub = `ub' return scalar lav = `lav' return scalar uav = `uav' return scalar lout = `lout' return scalar uout = `uout' return scalar tout = `lout' + `uout' if "`method'"=="adjusted" { return scalar mc = `MC' } else if "`method'"=="general" { return scalar g = `g' return scalar h = `h' } end version 11 mata: mata set matastrict on // sort results; note that only the results used for graphing will be sorted void robbox_sort(string scalar B, real scalar index, real scalar des) { real scalar over, nvars, nover, i real colvector p, p1, b, b1, N string scalar ename string rowvector tmp string matrix cstripe over = st_global("e(over)")!="" nvars = st_numscalar("e(N_vars)") if (over) { b = st_matrix("e(b)")' if (nvars==1) { p = p1 = order(b, 1) if (des) p = p1 = revorder(p) } else { if (index>nvars) { display("{txt}(specified sort index larger than number of variables;" + " will sort by average)") index = 0 } nover = st_numscalar("e(N_over)") N = st_matrix("e(_N)")' b1 = J(nover, 1, .) for (i=1; i<=nover; i++) { if (index==0) b1[i] = mean(select(b[|(i-1)*nvars+1 \ i*nvars|], N[|(i-1)*nvars+1 \ i*nvars|]:!=0)) else b1[i] = b[(i-1)*nvars + index] } p1 = order(b1, 1) if (des) p1 = revorder(p1) p = J(rows(b), 1, .) for (i=1; i<=nover; i++) { p[|(i-1)*nvars+1 \ i*nvars|] = (p1[i]-1)*nvars+1::p1[i]*nvars } } tmp = tokens(st_global("e(over_namelist)")) st_global("e(over_namelist)", invtokens(tmp[p1])) tmp = ("`"+`"""') :+ tokens(st_global("e(over_labels)")) :+ (`"""'+"'") st_global("e(over_labels)", invtokens(tmp[p1])) } else { b = st_matrix("e(b)")' p = p1 = order(b, 1) if (des) p = p1 = revorder(p) tmp = tokens(st_global("e(depvar)"))[p1] st_global("e(depvar)", invtokens(tmp[p1])) } cstripe = st_matrixcolstripe("e(b)")[p,] st_matrix(B, b[p]') st_matrixcolstripe(B, cstripe) stata("ereturn repost b = " + B + ", rename") tmp = ("box", "whiskers", "_N") for (i=1; i<=length(tmp); i++) { ename = "e(" + tmp[i] + ")" st_replacematrix(ename, st_matrix(ename)[,p]) st_matrixcolstripe(ename, cstripe) } } end exit