*! opplot version 1.17 - Biostat Global Consulting - 2023-12-22
*******************************************************************************
* Change log
* 				Updated
*				version
* Date 			number 	Name			What Changed
*-------------------------------------------------------------------------------
* 2015-12-21	1.01	MK Trimner		Added Starting comment and VCP globals
*
*										Put error to the vcqi_log as well as screen
*
* 2016-02-24	1.02	Dale Rhoda		Only set VCP and call vcqi_log_comment
*          								if VCQI is running; allows the possibility
*										to call the program from outside VCQI
*
* 2016-03-12	1.03	Dale Rhoda		Added option to save the plotting dataset
*
* 2016-07-04	1.04	Dale Rhoda		I'm seeing ties in coverage across
*										many strata (often at 100%) and so
* 										added a condition that sorts the plots
*										in descending order of coverage, with
*										ties broken by descending order of
*										weight, and ascending order of
*										cluster ID (clustvar)
*
* 2016-09-10	1.05	Dale Rhoda		Improve option to save plotting dataset
*
* 2017-08-26	1.06	Mary Prier		Added version 14.1 line
*
* 2017-12-11 	1.07	Dale Rhoda		Added passthru twoway option
*
* 2017-12-20	1.08	Dale Rhoda		Cleaned up how the pass_thru local works
*
* 2018-06-22	1.09	Dale Rhoda		Added the equalwidth option
*                                       (which will over-ride the weightvar
*                                       option, if necessary, to yield bars
*                                       of equal width)
*
* 2018-08-14	1.10	Mary Prier		Added line of code that strips out
*                      					double quotes of savedata option,
*										if user supplied filename in double quotes
*
* 2019-01-22	1.11	Dale Rhoda		Added options to add a line indicating
*      									the number of respondents per cluster.
*										(PLOTN NLINEColor NLINEWIdth NLINEPattern
*										 YTITLE2 and YROUND2)
*
* 2019-06-05	1.12	Dale Rhoda		When bars are tied on other quantities,
*										sort in descending order of
*										n_respondents
*
* 2019-12-12	1.13	Dale Rhoda		Cleaned up how the name() option is
*										processed
*
*										Allow the user to specify ytitle(none)
* 										to turn off the ytitle altogether
*										(which shifts the plot slightly
*										 lefttward)
*
* 2020-06-01	1.14	Dale Rhoda		Implement `if' and `in'
*
* 2023-07-28	1.15	MK Trimner		Added globals with Output Strings for multi-lingual purposes
* 2023-11-03	1.16	MK Trimner		Updated barheight to be a double variable
*										Added options for bar colors for Opplots
*										Added notes with Covg, CI, DEFF, ICC, StratumID , Nam and CI type information to Opplot dataset
*										Added STRATUMNAME as an option for notes
*										Added HIGH value for the cutoff for the high end of the bars. This defaults to 100
* 2023-12-22	1.17	Dale Rhoda		Suppress horizontal line at y=lcthreshold
*******************************************************************************

program define opplot
	version 14.1

	syntax varlist(max=1) [if] [in], CLUSTVAR(varname fv) ///
	[STRATVAR(varname fv) STRATUM(string) WEIGHTvar(varname) ///
	 TITLE(string) SUBtitle(string) FOOTnote(string) NOTE(string) ///
	 BARCOLOR1(string) LINECOLOR1(string) ///
	 BARCOLOR2(string) LINECOLOR2(string) EQUALWIDTH ///
	 XTITLE(string) YTITLE(string) XLABEL(string) YLABEL(string) ///
	 EXPORTSTRAtumname EXPORT(string) EXPORTWidth(integer 2000) ///
	 SAVING(string asis) NAME(string) SAVEDATA(string asis) ///
	 XSIZE(real -9) YSIZE(real -9) TWOWAY(string asis) PLOTN ///
	 NLINEColor(string asis) NLINEWidth(string) NLINEPattern(string) ///
	 YTITLE2(string asis) YROUND2(integer 5) ///
	 LCTHRESHold(integer 999) ///
	 BARCOLORHIGH1(string) BARCOLORMID1(string) BARCOLORLOW1(string) ///
	 BARCOLORHIGH2(string) BARCOLORMID2(string) BARCOLORLOW2(string) ///
	 LINECOLORHIGH1(string) LINECOLORMID1(string) LINECOLORLOW1(string) ///
	 LINECOLORHIGH2(string) LINECOLORMID2(string) LINECOLORLOW2(string)  ///
	 FOOTNOTECOVG FOOTNOTEDEFF FOOTNOTEICC STRATUMNAME(string) HIGH(integer 100)]

	* Note: opplot.ado is distributed as part of the World Health Organization's
	* Vaccination Coverage Quality Indicators (VCQI) suite of Stata programs.
	* When opplot is called from within VCQI, it interacts with a VCQI log.
	* When opplot is not called from VCQI, it ignores the few lines of
	* VCQI-related code.  The program knows if it is begin called from VCQI
	* because of the value of the global macro: VCQI_LOGOPEN.

	if "$VCQI_LOGOPEN" == "1" {
		local oldvcp $VCP
		global VCP opplot
		vcqi_log_comment $VCP 5 Flow "Starting"
	}

	tokenize "`varlist'"
	local yvar `1'

	quietly {

		preserve

		if "`if'" != "" keep `if'
		if "`in'" != "" keep `in'

		tempvar wclust wstrat bartop barwidth yweight wsum1 barheight cumulative_barwidth clustvar_copy n_respondents wtvar2 wtsum

		* Make a copy of clustvar to save at the end of the program
		clonevar `clustvar_copy' = `clustvar'

		* if the user doesn't specify a stratum, generate one for the program to use
		if "`stratvar'" == "" {
			tempvar stratumvariable
			gen `stratumvariable' = 1
			local stratvar `stratumvariable'
			local stratum 1
		}

		* if the user doesn't specify a weight variable, generate one
		if "`weightvar'" == "" {
			tempvar weightvariable
			gen `weightvariable' = 1
			local weightvar `weightvariable'
		}

		* if the user requests bars of equal width then
		* set up a new temporary weight variable (`wtvar2') and
		* scale it so that the sum of weights is equal in every cluster
		if "`equalwidth'" == "equalwidth" {
			clonevar `wtvar2' = `weightvar'
			local weightvar `wtvar2'
			tab `clustvar' if !missing(`yvar') // how many bars in the plot
			local wt_per_cluster = `=_N'/r(r)  // total weight per bar
			bysort `clustvar': egen `wtsum' = total(`weightvar')     if !missing(`yvar') // current sum of weights in each bar
			replace `wtvar2' = `wtvar2' * `wt_per_cluster' / `wtsum' if !missing(`yvar') // rescale sum of weights in each bar
		}

		* if y contains anything but 0 or 1 or . then fail
		capture tostring `stratvar', replace force
		capture assert inlist(`yvar',0,1,.) if `stratvar' == "`stratum'"
		if _rc == 9 {
			display as error "opplot Error: `yvar' should only have values of 0 or 1 or . when `stratvar' == `stratum'."

			if "$VCQI_LOGOPEN" == "1" {
				vcqi_log_comment $VCP 1 Error "opplot Error: `yvar' should only have values of 0 or 1 or . when `stratvar' == `stratum'."
				vcqi_halt_immediately
			}
		}


		* If the lcthreshold was not provided, wipe out the default value
		 if "`lcthreshold'" == "999" local lcthreshold

		 * check to make sure the colors are valid
		 * If the color is not valid, wipe out and the default value will be applied
		 if "`lcthreshold'" != "" {
			 forvalues n = 1/2 {
				 foreach v in bacrcolor`n' linecolor`n' barcolorhigh`n' barcolormid`n' barcolorlow`n' linecolorhigh`n' linecolormid`n' linecolorlow`n'  {
				 * Check for valid colors
				 local valid 0
					if "``v''" != "" {
						local first = substr("``v''",1,1)
						if !inlist("`first'","1","2","3","4","5","6","7","8","9") {
							capture findfile color-``v''.style
							local valid _rc
						}
						else {
							foreach col in ``v'' {
								capture assert `col' >=0 & `col' <= 255
								if _rc != 0 local valid 999
							}
						}

						 if `valid' != 0 {
							di as error "Option `v' currently takes value ``v''; it should be a valid color option"
							di as text "`msg'"
							capture vcqi_log_comment $VCP 1 Error "Option `v' currently takes value ``v''; it should be a valid color option"
							local `v'
						}
					}
				}
			}
		 }

		 * If the VCQI_CI_METHOD is not set, set to WILSON
		 if "$VCQI_CI_METHOD" == "" global VCQI_CI_METHOD WILSON

		 * Lets build the syvset command based on variables provided in opplot command
		 local svy `clustvar', singleunit(scaled)
		 if "`stratvar'" != "" local svy `svy' strata(`stratvar')
		 if "`weightvar'" != "" local svy `svy' weight(`weightvar')

		 * Now set the svyset command so we can calculate the coverage below
		 svyset `svy'

		* establish default values if the user doesn't specify them
		if "`ylabel'" == "" local ylabel 0(50)100, angle(h) nogrid
		if "`ytitle'" == "" local ytitle ${OS_316} //Percent of Cluster
		if "`=lower(trim("`ytitle'"))'" == "none" local ytitle
		if "`barcolor1'"  == "" local barcolor1 pink
		if "`linecolor1'" == "" local linecolor1 `barcolor'*1.5
		if "`barcolor2'"  == "" local barcolor2 white
		if "`linecolor2'" == "" local linecolor2 black*0.5

		if "`barcolorhigh1'" == "" local barcolorhigh1 "103 169 207"  //"103 169 207" // dark blue
		if "`linecolorhigh1'" == "" local linecolorhigh1 gs15
		if "`barcolorhigh2'" == "" local barcolorhigh2 gs15 //white
		if "`linecolorhigh2'" == "" local linecolorhigh2 black*0.5

		if "`barcolormid1'" == "" local barcolormid1 "0 0 128" // light blue
		if "`linecolormid1'" == "" local linecolormid1 gs15
		if "`barcolormid2'" == "" local barcolormid2 gs15 //white
		if "`linecolormid2'" == "" local linecolormid2 black*0.5

		if "`barcolorlow1'" == "" local barcolorlow1 "255 91 0" // orange
		if "`linecolorlow1'" == "" local linecolorlow1 gs15
		if "`barcolorlow2'" == "" local barcolorlow2 gs15 //white
		if "`linecolorlow2'" == "" local linecolorlow2 black*0.5


		if "`exportstratumname'" != "" {
			local export `stratum'.png
		}
		if `"`name'"' != "" {
			local stripname `name'
			local stripname = subinstr("`stripname'", ", replace", "", .)
			local stripname = subinstr("`stripname'", ",replace" , "", .)
			local stripname = substr("`stripname'",1,min(32,length("`stripname'")))
			local namestring name(`stripname' , replace)
		}
		if `"`saving'"' != "" local savingstring saving(`saving')
		if `xsize'  != -9 local xsizestring xsize(`xsize')
		if `ysize'  != -9 local ysizestring ysize(`ysize')

		* keep track of the number of respondents per cluster
		bysort `stratvar' `clustvar': gen `n_respondents' = _N

		* calculate sum of survey weights in each cluster and stratum
		bysort `stratvar' `clustvar': egen `wclust' = total(`weightvar')
		bysort `stratvar'           : egen `wstrat' = total(`weightvar')

		* calculate the proportion of the stratum weight in each cluster
		* (this corresponds to the bar width)
		gen `barwidth' = 100 * `wclust' / `wstrat'

		* calculate sum of survey weights for respondents with outcome = 1
		gen `yweight' = `yvar' * `weightvar'
		bysort `stratvar' `clustvar': egen `wsum1' = total(`yweight')

		* the height of each bar is the weighted proportion of respondents
		* with outcome = 1
		gen double `barheight' = (100*`wsum1'/`wclust')
		format %3.1f `barheight'

		* keep only observations in the stratum of interest
		keep if `stratvar' == "`stratum'"

		* Lets add some notes with the coverage information
		* We will first need to run the svypd command to get the values
		svypd `yvar', adjust method($VCQI_CI_METHOD)

		global covg : di %-3.1f 100*`r(svyp)'
		global lb : di %-2.1f 100*`r(lb_alpha)'
		global ub : di %-3.1f 100*`r(ub_alpha)'
		global deff : di %-5.1f `r(deff)'
		global covg_info ${covg}% (${lb}-${ub}%)

		quietly calcicc `yvar' `clustvar'
		global icc : di %-7.4f `r(anova_icc)'

		note: Estimated Coverage = ${covg_info} // do not inlcude CI in PLOT notes

		* Since we will use the estimated coverage later we will also store this in a char
		*char _dta[coverage_info] "${covg_info}"

		note: DEFF = $deff
		note: ICC = $icc
		note: StratumID = `stratum'
		if "`stratumname'" == "" {
			capture confirm var level3name
			if _rc  == 0 local stratumname =level3name[1]
		}
		if "`stratumname'" != "" note: StratumName = `stratumname'

		note: CI Method = $VCQI_CI_METHOD

		* keep only one observation per cluster
		bysort `stratvar' `clustvar': keep if _n == 1

		* sort the bars, left-to-right, in descending order of height and width
		* (and if there are ties...then sort the ties by descending n_respondents
		*  and ascending clusterID)
		gsort -`barheight' -`barwidth' -`n_respondents' `clustvar'
		*clonevar barheight = `barheight'

		* the background bars always have a height of 100%
		gen `bartop' = 100

		* Stata's facility for barcharts with varying widths requires a
		* variable that codes the cumulative barwidth
		gen `cumulative_barwidth' = sum(`barwidth')

		* add an extra row onto the dataset to make the x values work out correctly
		set obs `=_N+1'
		* shift the width up by one observation to make the x values
		* work out correctly
		forvalues i = `=_N'(-1)2 {
			replace `cumulative_barwidth' = `=`cumulative_barwidth'[`=`i'-1']' in `i'
		}
		replace `cumulative_barwidth' = 0 in 1

		* The option footnote is a synonym for note; if both are provided
		* ignore footnote
		if "`note'" == "" & "`footnote'" != "" local note `footnote'

		 // BUILD THE FOOTNOTES based on user specified input
		 // Estimated Covg = XXX.X%
		 // Design Effect = XX.X;
		 //        Intercluster Correlation Coefficient = X.XXXX
		* Create footnotes for opplot
		if "`footnotecovg'" != "" local opplot_footnote1 Estimated Coverage = $covg%
		if "`footnotedeff'" != "" local opplot_footnote2 Design Effect = $deff
		if "`footnoteicc'" != "" local opplot_footnote3 Intercluster Correlation Coefficient = $icc

		if "`opplot_footnote1'" != "" & "`opplot_footnote2'" != "" & "`opplot_footnote3'" != "" local opplot_footnote "`opplot_footnote1'; `opplot_footnote2';" "`opplot_footnote3'"
		if "`opplot_footnote1'" != "" & "`opplot_footnote2'" != "" & "`opplot_footnote3'" == "" local opplot_footnote "`opplot_footnote1'; `opplot_footnote2'"
		if "`opplot_footnote1'" != "" & "`opplot_footnote2'" == "" & "`opplot_footnote3'" != "" local opplot_footnote "`opplot_footnote1'; `opplot_footnote3'"
		if "`opplot_footnote1'" != "" & "`opplot_footnote2'" == "" & "`opplot_footnote3'" == "" local opplot_footnote "`opplot_footnote1'"

		if "`opplot_footnote1'" == "" & "`opplot_footnote2'" != "" & "`opplot_footnote3'" != "" local opplot_footnote "`opplot_footnote2'; `opplot_footnote3'"
		if "`opplot_footnote1'" == "" & "`opplot_footnote2'" != "" & "`opplot_footnote3'" == "" local opplot_footnote "`opplot_footnote2'"

		if "`opplot_footnote1'" == "" & "`opplot_footnote2'" == "" & "`opplot_footnote3'" != "" local opplot_footnote "`opplot_footnote3'"

		if `"`opplot_footnote'"' != "" local note "`opplot_footnote'" "`note'"
		di `"note"'

		foreach o in ylabel title subtitle { // moved note out to be able to include double quotes
			if "``o''" != "" local pass_thru `pass_thru' `o'(``o'')
		}

		if `"`note'"' != "" local pass_thru `pass_thru' note("`note'")
		if "`ytitle'" != "" local pass_thru `pass_thru' ytitle("`ytitle'")

		* If the user asks to plot the number of respondents (N) then
		* set up the twoway line and twoway scatteri syntax to do it.
		if "`plotn'" != "" {
			if "`nlinecolor'"   == "" local nlinecolor gs10
			if "`nlinewidth'"   == "" local nlinewidth *.5
			if "`nlinepattern'" == "" local nlinepattern dash
			if "`ytitle2'"      == "" local ytitle2 ${OS_317} //Number of Respondents
			sum `n_respondents'
			local y2max = `yround2' * ceil(`=r(max)+1'/`yround2')
			local plotnsyntax (line `n_respondents' `cumulative_barwidth', connect(stairstep) lc(`nlinecolor') lw(`nlinewidth') lp(`nlinepattern') yaxis(2)) ///
			                  (scatteri `=`n_respondents'[`=_N-1']' `=`cumulative_barwidth'[`=_N-1']' `=`n_respondents'[`=_N-1']' `=`cumulative_barwidth'[`=_N']' , ///
							   ms(i) connect(direct) lc(`nlinecolor') lw(`nlinewidth') lp(`nlinepattern') yaxis(2))
			local yaxis2title ytitle("`ytitle2'", axis(2))
			local yaxis2label ylabel(0(`yround2')`y2max', axis(2) angle(0))
		}


		if "`lcthreshold'" == "" {
			graph twoway ///
			(bar `bartop' `cumulative_barwidth', bartype(spanning) fcolor(`barcolor2') ///
				lpattern(solid) lcolor(`linecolor2') lwidth(*.1) ) ///
			(bar `barheight' `cumulative_barwidth', bartype(spanning) fcolor(`barcolor1') ///
				lpattern(solid) lcolor(`linecolor1') lwidth(*.1) ) ///
				`plotnsyntax' ,  ///
			graphregion(fcolor(white)) xtitle("`xtitle'") `pass_thru' `twoway' ///
			legend(off)  `namestring' `savingstring' `xsizestring' `ysizestring' ///
			xlabel(none) `yaxis2label' `yaxis2title' `threshold'
		}
		if "`lcthreshold'" != "" {

			* We will make several bar graphs with the different colors for each threshold
			* and add a yline at the threshold value

			local plotit 	(bar `bartop' `cumulative_barwidth' , bartype(spanning) fcolor(`"`barcolorhigh2'"') ///
				lpattern(solid) lcolor(`"`linecolorhigh2'"') lwidth(*.1) ) ///
			(bar `barheight' `cumulative_barwidth', bartype(spanning) fcolor(`"`barcolorhigh1'"') ///
				lpattern(solid) lcolor(`"`linecolorhigh1'"') lwidth(*.1) ) ///
			///	 // Create bars for mid section
			(bar `bartop' `cumulative_barwidth' if `barheight' < `high' | (_n == _N & `barheight'[_n-1] < `high'), bartype(spanning) fcolor(`"`barcolormid2'"') ///
				lpattern(solid) lcolor(`"`linecolormid2'"') lwidth(*.1) ) ///
			(bar `barheight' `cumulative_barwidth' if  `barheight'  < `high' | (_n == _N & `barheight'[_n-1] < `high'), bartype(spanning) fcolor(`"`barcolormid1'"') ///
				lpattern(solid) lcolor(`"`linecolormid1'"') lwidth(*.1) ) ///
			/// // Create bars for below the threshold
			(bar `bartop' `cumulative_barwidth' if `barheight' <= `lcthreshold' | (_n == _N & `barheight'[_n-1] <=`lcthreshold') , bartype(spanning) fcolor(`"`barcolorlow2'"') ///
				lpattern(solid) lcolor(`"`linecolorlow2'"') lwidth(*.1) ) ///
			(bar `barheight' `cumulative_barwidth' if `barheight'  <= `lcthreshold' | (_n == _N & `barheight'[_n-1] <=`lcthreshold'), bartype(spanning) fcolor(`"`barcolorlow1'"') ///
				lpattern(solid) lcolor(`"`linecolorlow1'"') lwidth(*.1) )

			graph twoway `plotit' `plotnsyntax' ///  (scatteri `lcthreshold' 0 `lcthreshold' 100, lc(gs8) lp(dash) connect(direct) ms(i))  // Removed the dashed line for now because the user doesn't have any control over it
			,  ///
			graphregion(fcolor(white)) xtitle("`xtitle'") `pass_thru' `twoway' ///
			legend(off)  `namestring' `savingstring' `xsizestring' `ysizestring' ///
			xlabel(none) `yaxis2label' `yaxis2title' `threshold'

		}

		if "`export'" != "" {
			graph export "`export'", width(`exportwidth') replace
		}

		* If the user has asked for underlying data to be saved, then
		* trim down to a small dataset that summarizes what is shown in
		* the bars; this is to help users identify the clusterid of a
		* particular bar in the figure; the order in which clusterids
		* appear in the saved dataset is the same order they appear in
		* the plot

		* Strip off double quotes if user supplied filename in double quotes
		local savedata = subinstr(`"`savedata'"', `"""',  "", .)

		* Strip off , replace if user supplied it in the savedata option
		local savedata = subinstr(`"`savedata'"', ", replace",  "", .)
		local savedata = subinstr(`"`savedata'"', ",replace" ,  "", .)

		if "`savedata'" != "" {
			drop in `=_N'
			capture drop yvar
			gen yvar = "`yvar'"
			capture drop stratvar
			gen stratvar = "`stratvar'"
			capture drop stratum
			gen stratum = "`stratum'"
			capture drop cluster
			gen cluster = "`clustvar'"
			capture drop `clustvar'
			rename `clustvar_copy' `clustvar'
			capture drop n_respondents
			rename `n_respondents' n_respondents
			capture drop barorder
			gen barorder = _n
			capture drop barwidth
			rename `barwidth' barwidth
			* replace cumulative barwidth with values that start at the width of bar 1
			capture drop cumulative_barwidth
			gen cumulative_barwidth = sum(barwidth)
			capture drop barheight
			rename `barheight' barheight
			keep  yvar stratvar stratum cluster `clustvar' n_respondents barorder barwidth cumulative_barwidth barheight
			order yvar stratvar stratum cluster `clustvar' n_respondents barorder barwidth cumulative_barwidth barheight
			capture save "`savedata'", replace
		}
	}

	if "$VCQI_LOGOPEN" == "1" {
		vcqi_log_comment $VCP 5 Flow "Exiting"
		global VCP `oldvcp'
	}
end