*! 1.1 25mar2010 rosa gini // changelog from 1.0 // corrected bug in markupcolor and marklowcolor // moved -twowayopts- at the end of the graph command, so that it might override any other default *! 1.1 31oct2017 minor bug fixes from funnelcompar rosa gini, silvia forni (based on confunnel by tom palmer and on eclplot by roger newson) *! 1.2 fixed [if] [in] & scatteroptsapplied to main scatter even if no markall opt program funnelperform, rclass sortpreserve version 9.1 /*set a maximum conditions for superimposing scatters and specifying colors and legends*/ local maxmarkcond=3 numlist "1(1)`maxmarkcond'" foreach i in "" `r(numlist)' { local markcondi `markcondi' MARKCOND`i'(string asis) local colormarkcondi `colormarkcondi' COLORMARKCOND`i'(string asis) local legendmarkcondi `legendmarkcondi' LEGENDMARKCOND`i'(string asis) local optionsmarkcondi `optionsmarkcondi' OPTIONSMARKCOND`i'(string asis) local labelcondi `labelcondi' LABELCOND`i' } local defaultReal -1.7e+38 /*****************/ /** SYNTAX **/ /*****************/ syntax varlist(min=3 max=4) [if] [in] [iweight/], [BINOMial POISson CONTinuous GAMMA BETA SMR /// ext_stand(real `defaultReal') NOWEIght ext_sd(real `defaultReal') EXACt Contours(numlist) CONSTant(real 1) /// MARKCONtour(real `defaultReal') MARKUP MARKUPCOLor(string) MARKLOW MARKLOWCOLor(string) MARKUnits(string asis) LABELUnits(string asis) MARKAll MARKTEXTOPtions(string) /// MARKSCATTEROPTions(string) MARKCOLor(string) TRUNC0 /// `markcondi' `colormarkcondi' `legendmarkcondi' `optionsmarkcondi' `labelcondi' /// NODRAw VERTical SCATTERCOLor(string) ASPECTratio(string) CONTCOLor(string) LEGENDCONTour UNITLABel(string) EXTRAplot(string) /// FUNCTIONLOWopts(string) FUNCTIONUPPopts(string) LEGENDMORe(string) LEGENDopts(string) ONEsided(string) /// SCATTERopts(string) SHADEDContours SOLIDContours N(integer 100) TWOWAYopts(string) YTITle(passthru) XTITle(passthru) TITle(passthru) LINECOLor(string) /// SAVing(string asis) DISPLAYcommand] marksample touse,strok tokenize `varlist' local value `1' local disp `2' local unit `3' local sdvalue `4' foreach v in ext_stand ext_sd markcontour{ if (``v''==`defaultReal'){ local `v' . } } /* 'value' contains proportions/rates/means -either crude or directly standardised- that are to be plotted; SMRs are considered as rates 'disp' contains the denominators of proportions/rates/means, or the expected number of events 'unit' contains the observation units whose values are to be compared 'sdvalue' contains the standard deviation of the variable 'value'; it is only admitted if value' arises from a continuous distribution and -ext_sd- is not specified */ /*OPTIONS* -BINOMial POISson CONTinuous GAMMA- only one of those options are possible and at least one is mandatory: they specify whether the routine expect to process proportions, rates or means -ext_stand(string) ext_sd(string)- specify external references; if they are not specified references are computed from the data -smr- specifies that rates are to be considered as indirectly standardised -NOWEIght- if references are computed form the data, noweiht tells not to weight the data before computing reference values -CONSTant(string)- contains the multiplicative constant for proportions or rates (eg 100 for percentages) -exact- specifies that the reference curves must be computed with exact methods instead of the defauild normal approximation (only for poisson or binomial) -nodraw- causes the command to generate the data and the command for plotting the funnel plot, but *not* to plot it -VERTical- specifies that the funnel must be vertical (horizontal is the default) -SCATTERCOLor(string)- specifies the color of the scatter points -Contours(string)- specifies significance levels (in percentage) of test 'measure'='ext_stand' at which contours should be plotted; default is .2 and 5 -LEGENDCONTour- puts the significance levels of the contours in the legend -MARKUnits(string asis) MARKTEXTOPtions MARKCOLor(string) MARKSCATTEROPTions()- allows to mark specific points listed in -MARKUnits- with specified color, possible labels and scatter options -MARKAll - marks all points with their labels or values -MARKCOND(string asis) COLORMARKCOND(string asis) LEGENDMARKCOND(string asis)- marks points accoprding to some conditions on data (up to 6 conditions are available of the form markcond(), markcond1(), markcond2()... ) -MARKCONtour(real) MARKUP MARKUPCOLor(string) MARKLOW MARKLOWCOLor(string)- marks point upper/lower the contour at -markcontour- confidence level // deprecated - mirror twoway function n option instead -STEP- set the granularity at which the contour lines are computed; a default is computed from the data -SAVing(string asis)- saves the dataset that generates the curves -DISPLAYcommand- makes the routine display the command that generates the graph -LINECOLor(string)- is the color of the reference line -graph options- SCATTERopts(string) ASPECTratio(string) CONTCOLor(string) UNITLABel(string) EXTRAplot(string) FUNCTIONLOWopts(string) FUNCTIONUPPopts(string) LEGENDMORe(string) LEGENDopts(string) ONEsided(string) SHADEDContours SOLIDContours YTITle(string) XTITle(string) TITle(string) TWOWAYopts(string) the option must be inserted in compound double quotes `"..."' in particular the string in -twowayopts- is placed at the end of the graph command, so that it might override any previous option */ /*****************/ /** ERRORS **/ /*****************/ local namedistr = trim("`poisson' `binomial' `gamma' `continuous' `beta'") if "`namedistr'"==""{ di as error "You must specify -binomial- or -poisson- or -continuous- or -gamma- or -beta-." di as error "If -binomial- is specified, then data are considered to be proportions arising from a binomial distribution" di as error "If -poisson- is specified, then data are considered to be rates arising from a poisson distribution or SMR" di as error "If -continuous- is specified, then data considered to be means arising from a normal distribution" error 198 } if (wordcount("`namedistr'")>1){ di as error "Only one option between -binomial-, -poisson-, -gamma-, -beta- and -continuous- may be specified" error 198 } if "`sdvalue'"!="" &"`continuous'"==""{ di as error "The fourth variable of the command '`sdvalue'' should only be specified when -continuous- option is also specified" di as error "In this case '`sdvalue'' is assumed to contain the standard errors of the means contained in '`value''" di as error "Recall that -continuous- must specified if data contained in '`value'' are means arising from a normal distribution" error 198 } if "`smr'"!="" & !inlist("`namedistr'","poisson","gamma","beta"){ di as error "Option -smr- implies that distribution -poisson-, -gamma- or -beta- must be specified " error 198 } if `ext_stand'==. &"`noweight'"!=""{ di as error "Option -ext_stand- is incompatible with -noweight- option," di as error "since -ext_stand- assumes that the data are to be compared against and external standard " di as error "whereas -noweight- assumes that the data are to be compared against the internal noweighted mean" error 198 } if `ext_sd'!=. & "`continuous'"=="" { di as error "Option -ext_sd- only might be specified if the -continuous- option is specified," di as error "and only in case the data are considered to be means arising from a normal distribution " di as error "and one wants to compare them against a standand mean with external standard deviation -ext_sd-" error 198 } if "`sdvalue'"=="" & `ext_sd'==. & "`continuous'"!=""{ di as error "When the -continuous- option is specified either a fourth variable or the external standard deviation -ext_sd- option should be specified" di as error "in order to let the program decide which reference standard deviation should be assumed" error 198 } if "`sdvalue'"!="" & `ext_sd'!=.{ di as error "You specified both variable '`sdvalue'' and the external standard deviation -ext_sd- option" di as error "You should decide whether the reference standard deviation should be computed from the data contained in '`sdvalue''" di as error "or the program should assume the number `ext_sd' as reference standard deviation" error 198 } if "`smr'"!="" & (`ext_stand'!=. |"`noweight'"!="") { di as error "Option 'smr' only might be specified in case of indirectly standardised rates " di as error "and in this case neither 'ext_stand' nor 'noweight' might be specified" di as error "since the external standard is fixed to the expected rate" error 198 } if inlist("`namedistr'","gamma") & "`smr'"==""{ di as error "Currently the gamma (incomplete gamma (ratio)) option can only be apply to standardised (observed/expected) rates. Please specify SMR" error 198 } if "`gamma'"!="" & !inlist(`ext_stand',.,1){ di as error "The external standard must be 1 if the gamma option is specified" error 198 } if "`beta'"=="" & "`weight'"!="" { di as error "an iweight can only be applied if the beta option is specified" error 198 } if "`beta'"!="" & (("`weight'"=="") != ("`smr'"=="")) { di as error "if using beta option, either SMR & iweight options must both be supplied, or both must be missing" error 198 } if "`weight'"!="" { confirm numeric variable `exp' capture assert `exp' > `disp' if `touse' if _rc != 0 { di as error "the iweight (population at risk) must be greater than the expected number of events" error 198 } } if "`exact'"!="" & !inlist("`namedistr'","poisson","binomial"){ di as error "Option 'exact' can only be applied to the poisson and binomial distributions" error 198 } if "`onesided'" != "" & "`onesided'" != "lower" & "`onesided'" != "upper" { di as error "onesided() must be unspecified, lower or upper" error 198 } if (`"`markunits'"'!=""&"`markall'"!=""){ di as error "only one option between -markunits- and -markall- must be specified" di as error "-markunits- lists the units that must be labelled"" di as error "-markall- specifies that all units must be labelled with their value labels (or with their values if no label is available)" error 198 } if (`"`markunits'"'==""&"`markscatteroptions'"!=""){ di as error "-markscatteroptions- can only be specified if -markunits- is specified" di as error "-markunits- lists the units that must be labelled" di as error "-markscatteroptions- specifies that they should be also marked with some option" error 198 } forv i1=0(1)`maxmarkcond' { local i=cond(`i1'==0,"","`i1'") if `"`markcond`i''"'!=""{ capture count if `markcond`i'' if _rc!=0{ di as error `"Option -markcond`i'- contains the string "`markcond`i''" which is not a valid condition on the dataset"' error 198 } } } if "`contours'" == "" { // default significance contours local contours "95.45 99.73" } if (`"`markup'"'!=""|`"`marklow'"'!="") & `markcontour'!=.{ //checks that `markcontour' is one of the contours if !inlist(`markcontour',`=subinstr("`contours'"," ",",",.)'){ di as error `"Option -markcontour- contains the value "`markcontour'" which is not one of the contours `contours' that you asked to compute"' error 198 } } confirm numeric variable `value' confirm numeric variable `disp' if ("`sdvalue'"!=""){ confirm numeric variable `sdvalue' } /*****************/ /** DEFAULTS **/ /*****************/ if (inlist("`namedistr'","gamma","beta")){ if ("`smr'"!="") { local ext_stand 1 } } if ("`beta'" != ""){ //a bit of a 'hack - uses similar loggic to the exact methods, even though it is not one local exact exact } if "`unitlabel'"==""{ local labelunit:variable label `unit' if "`labelunit'"!=""{ local unitlabel "`labelunit'" } else{ local unitlabel "Units" } } if "`scattercolor'" == "" { // default significance contours local scattercolor "black" } local ncontours = wordcount("`contours'") if (`"`markup'"'!=""|`"`marklow'"'!="") & `markcontour'==.{ //default for contour whose external should be marked local markcontour:word 1 of `contours' } /*a string with a dot cannot be the name of a local macro; hence in the names `c' of the confidence levels dots are substituted by the string -dot-, and the corresponding numbers are stored in local macros number`c'*/ local contours: subinstr local contours "." "_",all local markcontour: subinstr local markcontour "." "_",all foreach c of local contours{ local number`c' : subinstr local c "_" ".",all } if "`vertical'"!=""{ //reverse x and y axis if vertical is not specified (the default is horizontal) local x "x" local y "y" } else{ local x "y" local y "x" } if "`solidcontours'" == "solidcontours" { // use dashed or solid contours forvalues m = 1/`ncontours' { local linepatt `"`linepatt' solid"' } } else if "`solidcontours'" == "" { local linepatt "longdash dash shortdash dot shortdash_dot dash_dot longdash_dot" // line pattern styles for the contours } local j 0 foreach lp in `linepatt' { local lp`++j' `lp' } forvalues j = 1/`ncontours' { // shades for shaded contours if "`contcolor'" == "" { local shadedcontcol black } else { local shadedcontcol `contcolor' } if "`shadedcontours'" == "shadedcontours" { local lc`j' "`shadedcontcol'*`=1 - `j'/(`ncontours'*1.25)'" } else if "`shadedcontours'" == "" & "`contcolor'" == "" { local lc`j' "gs8" } else if "`shadedcontours'" == "" & "`contcolor'" != "" { local lc`j' "`contcolor'" } } if "`aspectratio'" != "" { // default aspectratio local aspectratio "aspectratio(`aspectratio')" } if "`linecolor'" == "" { // default color for reference line local linecolor red } if `"`markup'"'!="" & "`markupcolor'"==""{ //default for marking points upper the contour local markupcolor orange } if `"`marklow'"'!="" & "`marklowcolor'"==""{ //default for marking points lower the contour local marklowcolor blue } if "`markcolor'"==""{ local markcolor black } if ("`colormarkcond'"=="") { local colormarkcond red } /*AUXILIARY LOCAL MACROS*/ local plusub + local pluslb - /*****************/ /** ANALYSIS **/ /*****************/ /*set the value of the standard according to the options*/ local obs=_N /*generate the dataset of the reference curves*/ tempvar _funnel `y'var gen byte `_funnel'=`touse'==1 capture noisily{ qui{ if `ext_stand'==.{ if "`smr'"!=""{ local ext_stand 1 } else{ if "`noweight'"==""{ qui su `value' [fw= `disp'] if `touse',meanonly } else{ qui su `value' if `touse',meanonly } local ext_stand=r(mean) } } if "`binomial'"!=""{ local ext_sd=sqrt(`ext_stand'*(1-`ext_stand')) } if "`poisson'"!=""{ local ext_sd=sqrt(`ext_stand') } if `ext_sd'==. & "`continuous'"!=""{ tempvar variance gen float `variance'=`sdvalue'^2 if "`noweight'"==""{ qui su `variance' [fw= `disp'] if `touse',meanonly } else{ qui su `variance' if `touse',meanonly } local ext_sd=sqrt(r(mean)) } if "`exact'"==""{ foreach c in `contours' { tempvar ub`c' lb`c' if ("`gamma'"!=""){ gen float `ub`c''= invgammap(`disp'+1,(100+`number`c'')/200)/`disp' if `touse' gen float `lb`c''= invgammap(`disp',(100-`number`c'')/200)/`disp' if `touse' } else { local z = invnorm((100+`number`c'')/200) gen float `ub`c''=`ext_stand' + `z'*`ext_sd'/sqrt(`disp') if `touse' gen float `lb`c''=`ext_stand' - `z'*`ext_sd'/sqrt(`disp') if `touse' } } } else{ if ("`beta'"=="" | "`weight'"==""){ su `disp' if `touse', meanonly local ymin `r(min)' local ymax `r(max)' local steps = `ymax'-`ymin' + 1 local step = ceil(`steps'/`n') local n = ceil(`steps'/`step') set obs `=`obs' + `n'' replace `touse' = 1 in `=`obs'+1'/l replace `disp'=(_n-`obs')*`step'+`ymin' in `=`obs'+1'/l replace `_funnel'=2 in `=`obs'+1'/l if ("`beta'" != "") { tempvar estd_rpt gen float `estd_rpt' = `ext_stand' } } else { tempvar proportion gen float `proportion' = . } sort `disp' foreach c in `contours' { tempvar ub`c' lb`c' gen float `lb`c'' = . gen float `ub`c'' = . if "`binomial'"!=""{ mata: getInvbinom2("`disp'","`touse'", "`lb`c''", "`ub`c''", `ext_stand',`number`c'') } else if "`poisson'"!=""{ mata: getInvpoisson2("`disp'","`touse'", "`lb`c''", "`ub`c''", `ext_stand',`number`c'') } else if ("`beta'"!=""){ if ("`weight'"!="") { replace `proportion' = `disp'/`exp' if `touse' mata: getBeta2("`exp'","`proportion'","`touse'", "`lb`c''", "`ub`c''", `number`c'') replace `lb`c'' = `lb`c''/`disp' replace `ub`c'' = `ub`c''/`disp' } else { mata: getBeta2("`disp'","`estd_rpt'","`touse'", "`lb`c''", "`ub`c''", `number`c'') replace `lb`c'' = `lb`c''/`disp' replace `ub`c'' = `ub`c''/`disp' } } foreach lim in ub lb{ replace ``lim'`c''= ``lim'`c''*`constant' label variable ``lim'`c'' `"`lim' for variable `value' at confidence level `lev' assuming `namedistr' distribution"' } } } if `"`saving'"'!=""{ preserve local namefile=cond(strmatch(`"`saving'"',"*.dta")==1,`"`saving'"',`"`saving'.dta"') if ("`binomial'"!=""| "`poisson'"!=""){ local methodcurves=cond("`exact'"=="","normal approximation","exact method") } else{ local methodcurves "" } label data "Dataset for funnel plot (see characheristics with -char list-)" label variable `_funnel' "Instruction for funnel" label define funnel_label 0 "Not to use" 1 "Points for scatter" 2 "Points for curves" label values `_funnel' funnel_label local target_val `ext_stand' local target_sd `ext_sd' foreach char in namedistr target_val target_sd methodcurves{ char define _dta[`char'] `"``char''"' } save `namefile',replace noi di in ye `"File for funnel plot generation was saved in `namefile'"' restore } /**********************/ /** GRAPH GENERATION **/ /**********************/ noi di in ye "Starting plotting the graph" local val_label:variable label `value' if ("`val_label'"=="") { local val_label `value' } if (`constant'!=1){ tempvar scaled gen float `scaled'=`value'*`constant' //what to do with format - could copy format from value, but depending on scale may not be relevant local value `scaled' } //arguments in the correct order according to graph to be plotted vertical or horizontal local scatterargs=cond("`vertical'"!="","`disp' `value'","`value' `disp' " ) //CONTOURS // legend for contours local i 1 foreach c in `contours' { local i = `i' + 1 local h = `i' - 1 if "`onesided'" == "lower" { local lub`c' lc(none) local llb`c' lc(`lc`h'') lp(`lp`h'') lw(thin) `functionlowopts' local contourlabels`c' `"`=2*`h' - 1' "Sign. `number`c''%""' } else if "`onesided'" == "upper" { local lub`c' lc(`lc`h'') lp(`lp`h'') lw(thin) `functionupopts' local llb`c' lc(none) local contourlabels`c' `"`=2*`h' - 1' "Sign. `number`c''%""' } else{ local lub`c' lc(`lc`h'') lp(`lp`h'') lw(thin) `functionupopts' local llb`c' lc(`lc`h'') lp(`lp`h'') lw(thin) `functionlowopts' local contourlabels`c' `"`=2*`h'' "Sign. `number`c''%""' } } //command for contours sum `disp' if `touse', meanonly local range `r(min)' `r(max)' foreach c in `contours'{ local invnorm = invnorm((100+`number`c'')/200) foreach lim in ub lb{ if "`exact'"==""{ local functionArgs range(`range') n(`n') `l`lim'`c'' `contourargs' if ("`gamma'" == ""){ local contourargs=cond("`vertical'"!=""," horizontal","") local singlecontour `constant'*(`ext_stand'`plus`lim''`invnorm'*`ext_sd'/sqrt(x)) if ("`trunc0'"!="" && "`lim'"=="lb"){ local singlecontour max(`singlecontour',0) } local function `function' function `singlecontour', `functionArgs' || } else if ("`lim'"=="ub") { local function `function' function `constant'*invgammap(x+1,(100+`number`c'')/200)/x, `functionArgs' || function `constant'*invgammap(x,(100-`number`c'')/200)/x, `functionArgs' || } } else{ local contourargs=cond("`vertical'"!="","`disp' ``lim'`c''","``lim'`c'' `disp'" ) local function `function' line `contourargs', `l`lim'`c'' || } } local contourlabels `"`contourlabels' `contourlabels`c''"' } //POINTS //marks units specified in -markunits- or -markall- options tempvar marker gen byte `marker' = 1 if (`"`labelunits'"'!="") { local labelwords:word count `labelunits' capture confirm numeric variable `unit' if (_rc==0) { local unitvallabel:value label `unit' tempname newunitlabel if ("`unitvallabel'"==""){ label define `newunitlabel' `labelunits' } else { label copy `unitvallabel' `newunitlabel' label define `newunitlabel' `labelunits', modify } label values `unit' `newunitlabel' forvalues i=1(2)`labelwords'{ local v:word `i' of `labelunits' local markunits `markunits' `v' } } else { tempvar unitclone gen `unitclone' = `unit' forvalues i=1(2)`labelwords'{ local v:word `i' of `labelunits' local r:word `=`i'+1' of `labelunits' replace `unitclone' = `"`r'"' if `unitclone' == "`v'" local markunits "`markunits'" "`r'" } local markunits `""`markunits'""' local unit `unitclone' } } if (`"`markunits'"'!="") { capture confirm numeric variable `unit' //some of what follows would be far better managed with a regex, but as far as I am awatre the regex functions onlty replace the FIRST instance if (_rc==0) { //newer regex replace all command capture local markunits = ustrregexra("`markunits'", " +", ",") if (_rc!=0) { while (regexm("`markunits'", " +")){ local markunits = regexr("`markunits'", " +", ",") } } } else { //if version>=14 ustrregexra(`"`markunits'"', `"" +""', `"",""') capture local markunits = ustrregexra(`"`markunits'"', `"" +""', `"",""') if (_rc!=0) { while (regexm(`"`markunits'"', `"" +""')){ local markunits = regexr(`"`markunits'"', `"" +""', `"",""') } } } local ++maxmarkcond local markcond`maxmarkcond' inlist(`unit',`markunits') local labelcond`maxmarkcond' 1 } if("`markall'"!="") { local mainscatteropts mlabel(`unit') } local mainscatteropts `mainscatteropts' `scatteropts' numlist "1(1)`maxmarkcond'" foreach i in "" `r(numlist)' { if `"`markcond`i''"'!=""{ if "`optionsmarkcond`i''"==""{ local optionsmarkcond`i' `"`scatteropts' `optionsmarkcond'"' } if ("`colormarkcond`i''"==""){ local colormarkcond`i' `colormarkcond' } local markval = `i' +4 replace `marker' = `markval' if `_funnel'==1 & `markcond`i'' local markconditions `markconditions' scatter `scatterargs' if `marker' == `markval' & `touse', mc(`colormarkcond`i'') `optionsmarkcond`i'' if ("`labelcond`i''"!="") { local markconditions `markconditions' mlabel(`unit') mlabcolor(`colormarkcond`i'') } local markconditions `markconditions' || } } //marks units upper and/or lower the contour foreach lev in up low{ if `"`mark`lev''"'!=""{ if ("`lev'"=="up"){ replace `marker' = 2 if `_funnel'==1 & `value' > `ub`markcontour'' * `constant' local ifstate `marker' == 2 } else { replace `marker' = 3 if `_funnel'==1 & `value' < `lb`markcontour'' * `constant' local ifstate `marker' == 3 } if "`markscatter`lev'options'"==""{ local markscatter`lev'options `scatteropts' } local mark`lev'scatter `mark`lev'scatter' scatter `scatterargs' if `ifstate' & `touse', mc(`mark`lev'color') `markscatter`lev'options' mlabel(`unit') mlabcolor(`mark`lev'color') || } } } //end qui // LEGEND if "`legendopts'" != "off" { // position and other options if "`legendopts'" == "" { // default legend options local legendopts "ring(0) pos(2) size(small) symxsize(*.4) cols(1)" } //legend of units local legendtot `"order( 1 `"`unitlabel'"' "' //legend of conditions forv i1=0(1)`=`maxmarkcond'-1' { local i=cond(`i1'==0,"","`i1'") local ip=`i1'+1 if `"`markcond`i''"'!="" &`"`legendmarkcond`i''"'!=""{ local legendtot `"`legendtot' `=2*`ncontours' +1+`ip'' `"`legendmarkcond`i''"' "' } } //legend of contours if "`legendcontour'"!=""{ local legendtot `"`legendtot' `contourlabels' "' } local legendtot `"`legendtot' `legendmore') `legendopts'"' } else{ local legendtot `"off"' } // ACTUAL GENERATION OF GRAPH local graph_command `"twoway scatter `scatterargs' if `marker'==1 & `touse' , mc(`scattercolor') `mainscatteropts' || `function' `markupscatter' `marklowscatter' `markconditions' `extraplot' , `x'line(`=`ext_stand'*`constant'',lcolor(`linecolor') ) ytitle(`"`val_label'"') `aspectratio' `y'scale(`reverse') ylabel(, angle(horizontal)) `xtitle' `ytitle' `title' legend(`legendtot') `marktext' `markuptext' `marklowtext' `twowayopts'"' if "`displaycommand'"!=""{ di in ye `"`graph_command'"' } if "`nodraw'"==""{ `graph_command' } } //end capture //clean up qui drop if `_funnel'==2 if ("`newunitlabel'"!=""){ if ("`unitvallabel'"!=""){ label values `unit' `unitvallabel' } else { label values `unit' . } label drop `newunitlabel' } if _rc!=0 { qui error _rc } // RETURNED RESULTS return local target_val "`ext_stand'" return local target_sd "`ext_sd'" return local graph_command `"`graph_command'"' end /***************************/ /*MATA AUXILIARY FUNCTIONS */ /***************************/ mata: mata set matastrict on class interpolResult { public real matrix counts, weights } class interpolResult stepSearch(real scalar theta, real scalar ci, real colvector pop, pointer(real scalar function) delegate){ real scalar x,n,len,F,found,target,priorF,i; class interpolResult scalar returnVar; len = rows(pop); returnVar.counts=J(len,2,.); returnVar.weights=J(len,2,.); target = (100-ci)/200; x=0; //find lower bounds of first part of graph for(i=1;i<=2;i++) { if(i==2){ x = returnVar.counts[1,1]; target = 1-target; } for(n=1;n<=len;n++){ priorF=.; found=0; do { F = (*delegate)(pop[n], x,theta); if (F>=target){ if (x==0){ returnVar.counts[n,i] = x = 0; returnVar.weights[n,i] = 1; } else { returnVar.counts[n,i]= x; if (priorF==.){ priorF = (*delegate)(pop[n], x-1,theta); } returnVar.weights[n,i] = (F-target)/(F-priorF); } found=-1; } else { x++; priorF=F; } } while (!found) } } return(returnVar); } //for smr 1-ibetatail(x+1,n-x, p) //& ub (x,n-x+1, p) real scalar mybinomial(n,x,theta){ return(binomial(n,x,theta)); } void getInvbinom2(string scalar popVar, string scalar touse, string scalar lb, string scalar ub, real scalar theta, | real scalar ci){ real colvector popView; real matrix transformed; class interpolResult scalar result; if (args() < 6){ ci = 95; } st_view(popView,.,popVar,touse); result = stepSearch(theta, ci, popView, &mybinomial()); transformed = (result.counts:-result.weights):/popView; st_store(.,(lb,ub), touse, transformed); } real scalar mypoisson(n,x,theta){ return(poisson(n*theta,x)); } void getInvpoisson2(string scalar popVar, string scalar touse, string scalar lb, string scalar ub, real scalar theta, | real scalar ci){ real colvector popView; real matrix transformed; class interpolResult scalar result; if (args() < 6){ ci = 95; } st_view(popView,.,popVar,touse); result = stepSearch(theta, ci, popView, &mypoisson()); transformed = (result.counts:-result.weights):/popView; st_store(.,(lb,ub), touse, transformed); } class interpolResult invpoisson2(real scalar theta, real colvector pop){ expected=theta*n; i=0; do { F = poisson(expected,i) if (F >= p) { st_numscalar("r(rp)", i) return(i) } i++ } while (st_isnumfmt("r(rp)")<1) } void getBeta2(string scalar popVar,string scalar thetaVar, string scalar touse, string scalar lb, string scalar ub, | real scalar ci){ real colvector popView, thetaView; real matrix result; if (args() < 6){ ci = 95; } st_view(popView,.,popVar,touse); st_view(thetaView,.,thetaVar,touse); result = GetBetaCis(thetaView, popView, ci); st_store(.,(lb,ub), touse, result); } class searchArgs { public real scalar x,n,p public void toString() } void searchArgs::toString(){ printf("x:%f,n:%f,p:%f\n",this.x,this.n,this.p); } real scalar BinarySearch( real scalar target, real scalar minStart, real scalar maxStart, class searchArgs scalar sa, pointer(real scalar function) delegate, | real scalar approx, real scalar maxIts) { real scalar k, val, argLen, min, max; argLen = args(); if (argLen<6){ approx = 0.00005; } if (argLen < 7){ maxIts = 200; } min = minStart; max = maxStart; for(k=0;k= .) { return(.); } if ((abs(val-target)<=approx)){ return(sa.x) } if (val