*Set up log file
cap log close
clear
set more off
log using "Reproduce_vis_para_data_figs.log", replace

/*---- STATA code: reproduces figures in the following paper:
	   Presenting parasitological data: the good, 
	   the bad and the error bar ----*/

/*-------------- Small dataset plots --------------*/

*Read in example dataset
insheet using "small_dataset.csv", comma clear

*Create group labels for plotting
gen antigen = "Pf merozoite" if variable=="msp3_od"
replace antigen = "Pf-IE" if variable=="dbl5_std_od"
replace antigen = "Pv merozoite" if variable=="pvama1_std_od"

replace statusf = "Case" if statusf=="case"
replace statusf = "Control" if statusf=="control"

/*Figure 2: Dotplot*/
*Install function stripplot
ssc install stripplot

stripplot value, stack centre vertical width(.01) /// 
 by(antigen, col(3) graphregion(color(white)) note("")) /// 
 over(statusf) ///
 mcolor(black) /// 
 ylabel(,nogrid angle(0)) /// 
 ytitle("Total IgG titre (OD)") /// 
 xtitle(" ") ///
 subtitle(, fcolor(white) lcolor(white)) 

*Save Figure 2
graph export Fig2.tif, width(3900) replace

/*-------------- Large dataset plots --------------*/

*Read in example dataset
insheet using "large_dataset.csv", comma clear

*Title of y-axis
label variable value "Total IgG titre (OD)"
label variable logvalue "Natural-log total IgG titre"

*Create group labels for plotting
gen antigen = "Pf merozoite" if variable=="msp3_od"
replace antigen = "Pf-IE" if variable=="dbl5_std_od"
replace antigen = "Pv merozoite" if variable=="pvama1_std_od"

replace statusf = "Case" if statusf=="case"
replace statusf = "Control" if statusf=="control"

*Create variable required for statistical inference plot
encode statusf, generate(statusf_num)
encode antigen, generate(antigen_num)

gen status_antigen = statusf_num if antigen_num==1
replace status_antigen = statusf_num+3 if antigen_num==2
replace status_antigen = statusf_num+6 if antigen_num==3
table status_antigen

save "large_dataset.dta", replace

/*Figure 3: Histogram*/
histogram value, freq width(0.125) ///
 by(statusf antigen, graphregion(fcolor(white)) note("")) ///
 fcolor(white) bcolor(black) ///
 graphregion(fcolor(white)) ///
 ylabel(, grid glcolor(gs10) glpattern(shortdash) angle(0)) yline(0, lcolor(black)) ///
 subtitle(, fcolor(white) lcolor(white))   

*Save Figure 3
graph export Fig3.tif, width(3900) replace
	
/*Figure 4: Box & whisker plot - untransformed data*/
graph box value, ///
 over(statusf) over(antigen)  ///
 box(1, fcolor(white) bcolor(black)) mark(1, mcolor(black) msymbol(Oh)) ///
 graphregion(color(white)) ///
 ylabel(, grid glcolor(gs10) glpattern(shortdash) angle(0)) 

*Save Figure 4
graph export Fig4.tif, width(3900) replace

/*Figure 5: Plots for discrete variables*/
insheet using "parasitaemia_dataset.csv", comma clear

**Bar chart
*Plot with zero measurements excluded
hist parasit if parasit!=0, discrete freq ///
	graphregion(color(white)) lcolor(black) fcolor(white) lwidth(medthin) ///
	ylabel(,grid glcolor(gs10) glpattern(shortdash)) ///
	ytitle(, height(4)) xtitle("Parasitaemia (number of asexual parasites)", height(4)) 
	
*Save Figure 5
graph export Fig5.tif, width(3900) replace
	
/*Figure 6: Box & whisker plot - log-scale*/
use "large_dataset.dta", clear

graph box logvalue, ///
 over(statusf) over(antigen) ///
 box(1, fcolor(white) bcolor(black)) mark(1, mcolor(black) msymbol(Oh)) ///
 graphregion(color(white)) ///
 ylabel(, grid glcolor(gs10) glpattern(shortdash) angle(0))

*Save Figure 6
graph export Fig6.tif, width(3900) replace

/*Figure 7: Statistical inference plot*/
collapse (mean) meanvalue=logvalue (sd) sdvalue=logvalue (count) n=logvalue, by(status_antigen)

*CI on log-scale
gen ci_ul_ln = meanvalue + 1.96 * (sdvalue / sqrt(n))
gen ci_ll_ln = meanvalue - 1.96 * (sdvalue / sqrt(n))

*Geometric mean
gen geomean = exp(meanvalue)
gen geomean_ul = exp(ci_ul_ln)
gen geomean_ll = exp(ci_ll_ln)

*Statistical inference plot - log-scale
twoway scatter meanvalue status_antigen if status_antigen==1, /// 
		msymbol(S) msize(medium) mcolor(black) ///
		graphregion(color(white)) || ///
	   scatter meanvalue status_antigen if status_antigen==2, ///
		msymbol(S) msize(medium) mcolor(black) /// 
		graphregion(color(white)) || ///
	   rcap ci_ul_ln ci_ll_ln status_antigen if inlist(status_antigen, 1, 2), ///
		lcolor(black) lwidth(medthin) ///
		ylabel(, grid glcolor(gs10) glpattern(shortdash) angle(0)) ///
		ytitle("Mean natural-log of total IgG titre") ///  
	    xlabel(1 "Case" 2 "Control") ///
		xscale(range(0(1)3))  /// 
		xtitle(" ") ///  
	    legend(off) ///
		name(panela, replace) nodraw

*Statistical inference plot - geometric mean
twoway scatter geomean status_antigen if status_antigen==1, /// 
		msymbol(S) msize(medium) mcolor(black) ///  
		graphregion(color(white)) || ///
	   scatter geomean status_antigen if status_antigen==2, /// 
	    msymbol(S) msize(medium) mcolor(black) /// 
		graphregion(color(white)) || ///
	   rcap geomean_ul geomean_ll status_antigen if inlist(status_antigen, 1, 2), /// 
	    lcolor(black) lwidth(medthin) ///
		ylabel(, grid glcolor(gs10) glpattern(shortdash) angle(0)) ///
		ytitle("Geometric mean total IgG titre (OD)") ///
		xlabel(1 "Case" 2 "Control") ///
	    xscale(range(0(1)3)) ///  
		xtitle(" ")  ///
	    legend(off)  ///
	    name(panelb, replace) nodraw 

*Combine plots onto same graph region
graph combine panela panelb, cols(2) graphregion(color(white))

*Save Figure 7
graph export Fig7.tif, width(3900) replace

/*Figure 8: Detontator plot for descriptive analyses*/ 
use "large_dataset.dta", clear
collapse (mean) meanvalue=value (sd) sdvalue=value (max) maxvalue=value (count) n=value, by(status_antigen)

**Calculate error bars
gen sd_ul = meanvalue + sdvalue
gen sd_ll = meanvalue - sdvalue

*Error bar - sd
twoway bar meanvalue status_antigen if inlist(status_antigen, 2, 5, 8), /// 
		lcolor(black) lwidth(medthin) fcolor(white) /// 
		graphregion(color(white)) || ///
       bar meanvalue status_antigen if inlist(status_antigen, 1, 4, 7), /// 
	    lcolor(black) lwidth(medthin) fcolor(gs10) ///
		graphregion(color(white)) || ///
	   rcap sd_ul meanvalue status_antigen, /// 
	    lcolor(black) lwidth(medthin) ///
	    ylabel(0 0.25 0.5 0.75 1 1.25 1.5, nogrid angle(0)) /// 
		yscale(range(0(0.1)1.5)) ///
		ytitle("Total IgG titre (OD)") ///
		xlabel(1.5 "Pf merozoite" 4.5 "Pf-IE" 7.5 "Pv merozoite") ///
		xtitle(" ") ///
	    legend(row(1) order(2 "Case" 1 "Control" )) ///
	    name(panela, replace) nodraw 

*Error bar - max
twoway bar meanvalue status_antigen if inlist(status_antigen, 2, 5, 8), /// 
		lcolor(black) lwidth(medthin) fcolor(white) /// 
		graphregion(color(white)) || ///
       bar meanvalue status_antigen if inlist(status_antigen, 1, 4, 7), /// 
	    lcolor(black) lwidth(medthin) fcolor(gs10) ///  
		graphregion(color(white)) || ///
	   rcap maxvalue meanvalue status_antigen, /// 
	    lcolor(black) lwidth(medthin) ///
		ylabel(0 0.25 0.5 0.75 1 1.25 1.5, nogrid angle(0)) /// 
		yscale(range(0(0.1)1.5)) ///
		ytitle("Total IgG titre (OD)") ///
		xlabel(1.5 "Pf merozoite" 4.5 "Pf-IE" 7.5 "Pv merozoite") ///
		xtitle(" ") ///
	    legend(row(1) order(2 "Case" 1 "Control" )) ///
	    name(panelb, replace) nodraw 

*Combine plots onto same graph region
graph combine panela panelb, cols(2) graphregion(color(white)) 

*Save Figure 8
graph export Fig8.tif, width(3900) replace
	   

/*Figure 9: Detontator plot for statistical inference*/

**Calculate error bars
*Standard error
gen se_ul = meanvalue + (sdvalue / sqrt(n))
gen se_ll = meanvalue - (sdvalue / sqrt(n))

*95% confidence interval
gen ci_ul = meanvalue + 1.96 * (sdvalue / sqrt(n))
gen ci_ll = meanvalue - 1.96 * (sdvalue / sqrt(n))

*Error bar - standard error
twoway bar meanvalue status_antigen if status_antigen==2, /// 
		lcolor(black) fcolor(white) lwidth(medthin) graphregion(color(white)) || ///
       bar meanvalue status_antigen if status_antigen==1, /// 
	    lcolor(black) fcolor(gs10) lwidth(medthin) graphregion(color(white)) || ///
	   rcap se_ul meanvalue status_antigen if inlist(status_antigen, 1, 2), /// 
	    lcolor(black) lwidth(medthin) ///
		ylabel(0(0.1)0.5, nogrid angle(0)) /// 
		yscale(range(0(0.1)0.5)) ///
		ytitle("Total IgG titre (OD)")  ///
		xlabel( 1.5 "Pf merozoite") ///
		xscale(range(0(1)3)) ///
		xtitle(" ") ///
	    legend(row(1) order(2 "Case" 1 "Control" )) ///
	    name(panela, replace) nodraw 

*Error bar - 95% confidence interval	   
twoway bar meanvalue status_antigen if status_antigen==2, /// 
		lcolor(black) lwidth(medthin) fcolor(white) ///  
		graphregion(color(white)) || ///
       bar meanvalue status_antigen if status_antigen==1, /// 
	    lcolor(black) lwidth(medthin) fcolor(gs10) ///  
		graphregion(color(white)) || ///
	   rcap ci_ul meanvalue status_antigen if inlist(status_antigen, 1, 2), /// 
	    lcolor(black) lwidth(medthin) ///
	    ylabel(0(0.1)0.5, nogrid angle(0)) /// 
	    yscale(range(0(0.1)0.5)) ///
	    ytitle("Total IgG titre (OD)")  ///
	    xlabel(1.5 "Pf merozoite") ///
	    xscale(range(0(1)3)) ///  
	    xtitle(" ") ///
	    legend(row(1) order(2 "Case" 1 "Control")) ///
	    name(panelb, replace) nodraw 

*Combine plots onto same graph region
graph combine panela panelb, cols(2) graphregion(color(white)) 

*Save Figure 9
graph export Fig9.tif, width(3900) replace
	   

/*Supplementary figure 1: Statistical inference plot for population proportion*/
clear

*The function mycii.ado save the results from the 'cii' function to a Stata dataset
/*-------------- mycii.ado*/
 program def mycii 
  version 13.1
  *cii `1' `2' `3'
  cii `1' `2' 
  qui replace N = r(N) in l
  qui replace mean = r(mean) in l
  qui replace se = r(se) in l
  qui replace lb = r(lb) in l
  qui replace ub = r(ub) in l
  local n = _N + 1
  qui set obs `n'
 end 
 /*-------------------*/

*Initial a Stata dataset to save the results of the cii function to 
set obs 1
gen N = .
gen mean = .
gen se = .
gen lb = .
gen ub = .

**Run the function 'mycii.do'
*Arg `1'=total number of observations 
*Arg `2'=number of success
*For example the no. preg woment with Pf is 1059 and 266 gave birth to a low-birthweight baby
mycii 1059 266 /*Pf*/
mycii 513 97   /*Pv*/
mycii 5434 712 /*No malaria*/

*Remove last row of dataset that consists of missing values
gen row=_n 
drop if row==4

twoway scatter mean row, ///  
		msymbol(S) msize(medium) mcolor(black) /// 
		graphregion(color(white)) || ///
	   rcap ub lb row, /// 
	    lcolor(black) lwidth(medthin) ///
		ylabel(0(0.1)0.3, grid glcolor(gs10) glpattern(shortdash) angle(0)) ///
		yscale(range(0(0.1)0.3)) ///
		ytitle("Proportion of low-birthweight babies") ///
		xlabel(1 "Pf" 2 "Pv" 3 "No malaria") /// 
	    xscale(range(0(1)4)) /// 
		xtitle("Malaria exposure")  ///
	    legend(off)

*Save Supplementary Figure 1
graph export Supp_Fig1.tif, width(3900) replace

program drop mycii
	   
clear

log close 
exit