Merge pull request #107 from opensafely/add-visualisations

add plots

.gitignore

@@ -13,3 +13,4 @@ analysis/*
 venv/
 .cache
 *.exe
+!summary_results15Jun2020_ctr.xlsx

analysis/gr_forestplot_asthma.do

@@ -0,0 +1,109 @@
+
+*need to add in analysis column, in between primsec and title
+
+import excel "C:\Users\lsh1300968\Dropbox\OneDrive\LSHTM\Work\COVID\OpenSAFELY\ICS\summary_results15Jun2020_ctr.xlsx", clear firstrow sheet("Asthma")
+
+gen result_label =string(estimate,"%6.2f") + " (" + string(min95,"%6.2f") + "-" + string(max95,"%6.2f") + ")"
+
+egen first_primsec = tag(primsec)
+egen first_analysis = tag(outcome)
+gen order = _n
+
+
+
+
+** Now for the figure
+
+*only keep variable/estimates you need for now
+*keep if parm == "1.t2dm"
+graph set window fontface "Calibri"
+**sort primsec outcome     
+
+cap drop label
+gen label = outcome
+local obs1 = _N + 1
+set obs `obs1'
+replace outcome=0 if outcome == .
+replace order=0 if order == .
+sort  outcome  order 
+
+* set some label values to missing so that we don't replicate information on the graph
+replace primsec="" if first_primsec != 1
+replace analysis="" if first_analysis != 1
+
+
+* create label for outcome column
+cap drop obs
+gen obs = _n
+replace title="Exposure" if obs==1
+replace analysis="Analysis" if obs==1
+replace result_label = "HR (95% CI)" if obs==1
+
+
+*bold face some labels
+gen bftitle = "{bf:" + title + "}" if title == "Exposure"
+replace bftitle = title if bftitle == ""
+gen bfanalysis = "{bf:" + analysis + "}" if analysis == "Analysis"
+replace bfanalysis = analysis if bfanalysis == ""
+gen bf_result = "{bf:" + result_label + "}" if result_label == "HR (95% CI)"
+replace bf_result = result_label if bf_result == ""
+gen bf_primsec = "{bf:" + primsec + "}" 
+
+cap drop x0_*
+gen x0_7 = -6
+gen x0_3=-3.5
+gen x0_1=-1.5
+gen x0_14=7
+
+
+
+cap drop obs
+gen obs = _n 
+cap drop no_obs 
+cap drop obs2
+cap drop obs3
+
+egen no_obs = max(obs) 
+gen obs2 = no_obs - obs + 1
+rename obs obs3 
+rename obs2 obs
+*add space between column titles and data
+replace obs=obs+0.5 if obs==no_obs
+*add space between primary and secondary outcomes
+// replace obs=obs+.5 if outcome <=9
+// replace obs=obs+.5 if outcome <=6
+// replace obs=obs+.5 if outcome <=3
+// *add to no_obs the extra space I added
+// replace no_obs = no_obs+2
+local height=obs
+
+
+graph twoway ///
+|| scatter obs estimate if obs<no_obs, msymbol(circle) msize(vsmall) mcolor(black)		/// data points 
+		xline(1, lp(solid) lw(vthin) lcolor(black))					/// add ref line
+	|| rcap min95 max95 obs if obs<no_obs, horizontal lw(vthin) lcolor(black)					/// add the CIs
+|| scatter obs estimate if obs<no_obs, msymbol(circle) msize(vsmall) mcolor(black)		/// data points 
+		xline(1, lp(solid) lw(vthin) lcolor(black))					/// add ref line
+	|| rcap min95 max95 obs if obs<no_obs, horizontal lw(vthin) lcolor(black)					/// add the CIs
+		/// Primary/Secondary outcomes
+|| scatter obs x0_7 if obs<no_obs, m(i) mlab(bf_primsec) mlabcol(black) mlabsize(vsmall) ///
+|| scatter obs x0_7 if obs>no_obs, m(i) mlab(bf_primsec) mlabcol(black) mlabsize(small) ///		
+		/// Analysis
+|| scatter obs x0_3 if obs<no_obs, m(i) mlab(bfanalysis) mlabcol(black) mlabsize(vsmall) ///
+|| scatter obs x0_3 if obs>no_obs, m(i) mlab(bfanalysis) mlabcol(black) mlabsize(small) ///	
+		/// Exposure
+|| scatter obs x0_1 if obs<no_obs, m(i)  mlab(bftitle) mlabcol(black) mlabsize(vsmall)  ///
+|| scatter obs x0_1 if obs>no_obs, m(i)  mlab(bftitle) mlabcol(black) mlabsize(small)  ///
+		/// add results labels
+|| scatter obs x0_14 if obs<no_obs, m(i)  mlab(bf_result) mlabcol(black) mlabsize(vsmall) mlabposition(0) mlabgap(tiny)   ///
+|| scatter obs x0_14 if obs>no_obs, m(i)  mlab(bf_result) mlabcol(black) mlabsize(small) mlabposition(0) mlabgap(tiny)   ///
+		, legend(off)						/// turn legend off
+		xtitle("Hazard ratio (HR)", size(vsmall) margin(40 0 0 2)) 		/// x-axis title (left right bottom top) - legend off
+		xlab(0(1)6, labsize(vsmall)) /// x-axis tick marks
+		xscale(range(0.01 8))					///	resize x-axis
+		, ylab(none) ytitle("") 	/// y-axis no labels or title
+		yscale(range(1 `height') lcolor(white))					/// resize y-axis
+		graphregion(color(white)) ysize(15) xsize(20) saving("C:\Users\lsh1300968\Dropbox\OneDrive\LSHTM\Work\COVID\OpenSAFELY\ICS\forestplot_asthma", replace)	/// get rid of rubbish grey/blue around graph
+
+graph export "C:\Users\lsh1300968\Dropbox\OneDrive\LSHTM\Work\COVID\OpenSAFELY\ICS\forestplot_asthma.tif"
+** FYI: to rerun the graph, need to rerun the entire do file -- not sure why

analysis/gr_forestplot_copd.do

@@ -0,0 +1,109 @@
+
+*need to add in analysis column, in between primsec and title
+
+import excel "C:\Users\lsh1300968\Dropbox\OneDrive\LSHTM\Work\COVID\OpenSAFELY\ICS\summary_results15Jun2020_ctr.xlsx", clear firstrow sheet("COPD")
+
+gen result_label =string(estimate,"%6.2f") + " (" + string(min95,"%6.2f") + "-" + string(max95,"%6.2f") + ")"
+
+egen first_primsec = tag(primsec)
+egen first_analysis = tag(outcome)
+gen order = _n
+
+
+
+
+** Now for the figure
+
+*only keep variable/estimates you need for now
+*keep if parm == "1.t2dm"
+graph set window fontface "Calibri"
+**sort primsec outcome     
+
+cap drop label
+gen label = outcome
+local obs1 = _N + 1
+set obs `obs1'
+replace outcome=0 if outcome == .
+replace order=0 if order == .
+sort  outcome  order 
+
+* set some label values to missing so that we don't replicate information on the graph
+replace primsec="" if first_primsec != 1
+replace analysis="" if first_analysis != 1
+
+
+* create label for outcome column
+cap drop obs
+gen obs = _n
+replace title="Exposure" if obs==1
+replace analysis="Analysis" if obs==1
+replace result_label = "HR (95% CI)" if obs==1
+
+
+*bold face some labels
+gen bftitle = "{bf:" + title + "}" if title == "Exposure"
+replace bftitle = title if bftitle == ""
+gen bfanalysis = "{bf:" + analysis + "}" if analysis == "Analysis"
+replace bfanalysis = analysis if bfanalysis == ""
+gen bf_result = "{bf:" + result_label + "}" if result_label == "HR (95% CI)"
+replace bf_result = result_label if bf_result == ""
+gen bf_primsec = "{bf:" + primsec + "}" 
+
+cap drop x0_*
+gen x0_7 = -5.5
+gen x0_3=-3.5
+gen x0_1=-1.5
+gen x0_14=4
+
+
+
+cap drop obs
+gen obs = _n 
+cap drop no_obs 
+cap drop obs2
+cap drop obs3
+
+egen no_obs = max(obs) 
+gen obs2 = no_obs - obs + 1
+rename obs obs3 
+rename obs2 obs
+*add space between column titles and data
+replace obs=obs+0.5 if obs==no_obs
+*add space between primary and secondary outcomes
+// replace obs=obs+.5 if outcome <=9
+// replace obs=obs+.5 if outcome <=6
+// replace obs=obs+.5 if outcome <=3
+// *add to no_obs the extra space I added
+// replace no_obs = no_obs+2
+local height=obs
+
+
+graph twoway ///
+|| scatter obs estimate if obs<no_obs, msymbol(circle) msize(vsmall) mcolor(black)		/// data points 
+		xline(1, lp(solid) lw(vthin) lcolor(black))					/// add ref line
+	|| rcap min95 max95 obs if obs<no_obs, horizontal lw(vthin) lcolor(black)					/// add the CIs
+|| scatter obs estimate if obs<no_obs, msymbol(circle) msize(vsmall) mcolor(black)		/// data points 
+		xline(1, lp(solid) lw(vthin) lcolor(black))					/// add ref line
+	|| rcap min95 max95 obs if obs<no_obs, horizontal lw(vthin) lcolor(black)					/// add the CIs
+		/// Primary/Secondary outcomes
+|| scatter obs x0_7 if obs<no_obs, m(i) mlab(bf_primsec) mlabcol(black) mlabsize(vsmall) ///
+|| scatter obs x0_7 if obs>no_obs, m(i) mlab(bf_primsec) mlabcol(black) mlabsize(small) ///		
+		/// Analysis
+|| scatter obs x0_3 if obs<no_obs, m(i) mlab(bfanalysis) mlabcol(black) mlabsize(vsmall) ///
+|| scatter obs x0_3 if obs>no_obs, m(i) mlab(bfanalysis) mlabcol(black) mlabsize(small) ///	
+		/// Exposure
+|| scatter obs x0_1 if obs<no_obs, m(i)  mlab(bftitle) mlabcol(black) mlabsize(vsmall)  ///
+|| scatter obs x0_1 if obs>no_obs, m(i)  mlab(bftitle) mlabcol(black) mlabsize(small)  ///
+		/// add results labels
+|| scatter obs x0_14 if obs<no_obs, m(i)  mlab(bf_result) mlabcol(black) mlabsize(vsmall) mlabposition(0) mlabgap(tiny)   ///
+|| scatter obs x0_14 if obs>no_obs, m(i)  mlab(bf_result) mlabcol(black) mlabsize(small) mlabposition(0) mlabgap(tiny)   ///
+		, legend(off)						/// turn legend off
+		xtitle("Hazard ratio (HR)", size(vsmall) margin(40 0 0 2)) 		/// x-axis title (left right bottom top) - legend off
+		xlab(0(0.5)3, labsize(vsmall)) /// x-axis tick marks
+		xscale(range(0.01 5))					///	resize x-axis
+		, ylab(none) ytitle("") 	/// y-axis no labels or title
+		yscale(range(1 `height') lcolor(white))					/// resize y-axis
+		graphregion(color(white)) ysize(15) xsize(20) saving("C:\Users\lsh1300968\Dropbox\OneDrive\LSHTM\Work\COVID\OpenSAFELY\ICS\forestplot_copd", replace)	/// get rid of rubbish grey/blue around graph
+
+graph export "C:\Users\lsh1300968\Dropbox\OneDrive\LSHTM\Work\COVID\OpenSAFELY\ICS\forestplot_copd.tif"
+** FYI: to rerun the graph, need to rerun the entire do file -- not sure why

analysis/ics-e-value_plot.do

@@ -0,0 +1,44 @@
+**************************
+/* Set output directory */
+**************************
+
+glob OUTPUT_DIR = ""
+
+****************************
+/* Input effect estimates */
+****************************
+
+* COPD
+glob hr_copd = 1.38
+glob lci_copd = 1.08
+glob uci_copd = 1.75
+
+* Asthma low-dose ICS
+glob hr_asthma_ld = 1.10
+glob lci_asthma_ld = 0.82
+glob uci_asthma_ld = 1.49
+
+* Asthma high-dose ICS
+glob hr_asthma_hd = 1.52
+glob lci_asthma_hd = 1.08
+glob uci_asthma_hd = 2.14
+
+**************************
+/* Create e-value plots */
+**************************
+
+* COPD 
+evalue hr $hr_copd, lcl($lci_copd) ucl($uci_copd) true(0.8) figure(scheme("s1color") aspect(1.01) ytitle(,size(small) margin(medium)) xtitle(,size(small) margin(medium)))
+* addplot: 
+graph save "${OUTPUT_DIR}evalue_COPD_protective.gph", replace
+graph export "${OUTPUT_DIR}evalue_COPD_protective.jpg", replace quality(100) width(1000)
+
+* Ashtma low-dose ICS
+evalue hr $hr_asthma_ld, lcl($lci_asthma_ld) ucl($uci_asthma_ld) true(0.8) figure(scheme("s1color") aspect(1.01) ytitle(,size(small) margin(medium)) xtitle(,size(small) margin(medium))) 
+graph save "${OUTPUT_DIR}evalue_Asthma_LD_protective.gph", replace
+graph export "${OUTPUT_DIR}evalue_Asthma_LD_protective.jpg", replace quality(100) width(1000)
+
+* Asthma high-dose ICS
+evalue hr $hr_asthma_hd, lcl($lci_asthma_hd) ucl($uci_asthma_hd) true(0.8) figure(scheme("s1color") aspect(1.01) ytitle(,size(small) margin(medium)) xtitle(,size(small) margin(medium)))
+graph save "${OUTPUT_DIR}evalue_Asthma_HD_protective.gph", replace 
+graph export "${OUTPUT_DIR}evalue_Asthma_HD_protective.jpg", replace quality(100) width(1000)