/* Biostat 200 Lab 8

This lab will implement commands to make informative and nicely formatted and labeled graphics.  A few of these commands were used in previous labs.

First, modify the next line so it gives the path to the directory where you 
have put this do-file and the seizures dataset for this lab.
*/

!gunzip -f *.dta.gz
cd "~/work/mb/atcr/200/labs"

capture log close
log using "lab 8", replace

/* Read in the data for this lab posted on the class web site and briefly familiarize yourself with the variables. */

use seizures_advanced.dta, clear
des
sum

/* First we will look at the distribution of continuous variables including 
total sedation hours (sedhrs), age (age), PRO dose (dose_pro), and length of 
stay (los), using density and QQ plots to examine each distribution. (See 
VGSM page 13-14 for  how to read QQ plots).  Then we will use scatter plots 
to check whether total sedation hours is related to the others. 

You can use the pdf help for the graphics commands to understand the details 
of these commands (and no doubt make them nicer!).  The code to get the labels 
should be familiar from earlier labs. Note that if and else are used to make 
age and PRO dose predictors of sedation hours, then make sedation hours the
predictor of LOS, based on a rough understanding of the likely causal direction 
of any association. */

foreach x in sedhrs age dose_pro los {
	local `x'_label : var label `x'
	twoway (kdensity `x'), ///
		ytitle("{bf:Density}") xtitle("{bf:``x'_label'}") name(kd_`x', replace)
	qnorm `x', name(qq_`x', replace)
	if inlist("`x'","age","dose_pro") {
		twoway ///
			(scatter sedhrs `x', msize(small)) ///
			(lfit sedhrs `x') ///
			(lowess sedhrs `x'), ///
			xtitle("{bf:``x'_label'}") ytitle("{bf:`sedhrs_label'}") ///
			legend(order(2 "linear fit" 3 "lowess fit") ///
				ring(0) pos(11) rows(3) region(lcolor(white))) ///
			graphregion(color(white)) ///
			name(sc_`x'_sedhrs, replace)
		}
	else if "`x'"=="los" {
		twoway ///
			(scatter `x' sedhrs, msize(small)) ///
			(lfit `x' sedhrs) ///
			(lowess `x' sedhrs), ///
			ytitle("{bf:``x'_label'}") xtitle("{bf:`sedhrs_label'}") ///
			legend(order(2 "linear fit" 3 "lowess fit") ///
				ring(0) pos(11) rows(3) region(lcolor(white))) ///
			graphregion(color(white)) ///
			name(sc_`x'_sedhrs, replace)
		}
	}

/* From the regression lines, total sedation hours appears to be related to 
length of stay only, although the discrepancy between the linear and lowess 
regression lines -- the latter does not force the fit to be a straight line 
-- suggests that the relationship may not be linear.  

The density and QQ plots indicate that total sedation hours and length of stay 
are badly right-skewed, but that age and PRO dose look reasonably normal. We 
can use the gladder and qladder command to find optimal transformations for 
the two skewed variables; the transformed distributions are displayed using 
histograms and QQ plots, respectively. */

foreach x in sedhrs los {
	gladder `x', name(gl_`x', replace)
	qladder `x', name(ql_`x', replace)
	}
	
/* My reading of the plots is that the log transformation does the best job 
of normalizing each of these variables. We will make density plots of the 
transformed variables to check. */

foreach x in sedhrs los {
	cap drop log_`x'
	qui gen log_`x' = log(`x')
	local label : var label `x'
	twoway (kdensity log_`x'), ///
		ytitle("{bf:Density}") xtitle("{bf: log `label'}") name(kd_log`x', replace)
	qnorm log_`x', name(qq_log`x', replace)
	}
	
/* Log transformation does well for sedation hours, but does not quite get 
rid of the long right tail for LOS.  That being said, a continuous predictor 
does not need to be normal! Rather, its association with the outcome needs 
to be linear; you will learn more about how to check this in Biostat 208, and 
fix violations, in particular by transforming the continuous predictor.  Now 
we will re-examine the association of sedation hours with LOS, with both 
log-transformed. */

local xlabel : var label sedhrs
local ylabel : var label los
twoway ///
	(scatter log_los log_sedhrs, msize(small)) /// 
	(lfit log_los log_sedhrs) /// 
	(lowess log_los log_sedhrs), ///
	xtitle("{bf:log `xlabel'}") ytitle("{bf:log `ylabel'}") ///
	legend(order(2 "linear fit"  3 "lowess fit") ///
		ring(0) pos(11) rows(3) region(lcolor(white))) ///
	graphregion(color(white)) ///
	name(sc_log`x'_logsedhrs, replace)
	
/* So there is still considerable evidence for departure from linearity here,
as shown by the difference between the linear and lowess fits. 

Now we will examine the distribution of the number of new AEs, a count with 
range 0-7, and assess whether it is associated with either sedation hours or 
LOS. That said, it might be reasonable to use it as continuous. We consider 
a categorization of the count of new AEs, on the grounds that 8 categories 
are too many for a sample of this size, not to mention interpretation. Here 
we hypothesize that sedation hours --> # new AEs --> LOS, determining what 
goes onthe x-axis in the plots. For sedation hours, we'll use the log 
transformation. In Biostat 208 and 209 you will learn about models for 
count and ordinal outcomes, which could be used for new AEs. */

tab n_newae
cap drop n_newaecat
qui recode n_newae (0=0) (1=1) (2/4=2 "2-4") (5/max=3 ">4"), gen(n_newaecat)
local xlabel : var label n_newae
label var n_newaecat "`xlabel'"
* check this
tab n_newae n_newaecat

* QQ plot
qnorm n_newae, name(qq_newae, replace)

global label1 : var label sedhrs
global label2 : var label n_newae
global label3 : var label los

* sedation hours --> new AEs
twoway ///
	(scatter n_newae log_sedhrs, msize(small)) ///
	(lfit  n_newae log_sedhrs) ///
	(lowess  n_newae log_sedhrs), ///
	xtitle("{bf:log $label1}") ytitle("{bf:$label2}") ///
	legend(order(2 "linear fit" 3 "lowess fit") ///
		ring(0) pos(11) rows(3) region(lcolor(white))) ///
	graphregion(color(white)) ///
	name(sc_newae_sedhrs, replace)
	
* new AEs --> LOS
twoway ///
	(scatter los n_newae, msize(small)) ///
	(lfit los n_newae) ///
	(lowess los n_newae), ///
	xtitle("{bf:$label2}") ytitle("{bf:$label3}") ///
	legend(order(2 "linear fit" 3 "lowess fit") ///
		ring(0) pos(11) rows(3) region(lcolor(white))) ///
	graphregion(color(white)) ///
	name(sc_los_newae, replace)
graph box los, over(n_newaecat) title("$label2") name(box_los_newae, replace)

/* There is some evidence that the number of new AEs may increase with sedation 
hours, and the boxplot suggests increased longer LOS amongthose with 4 or more 
new AEs. 

Now we will look at whether discharge to home vs rehab is associated with 
sedation hours or LOS, using boxplots, and test for differences using Wilcoxon 
test (because both these outcomes are right skewed. Although it reverses the 
causal direction, treating discharge location as the predictor makes this 
easier. In Biostat 208, you will learn (more) about estimating the adjusted 
effects of continuous predictors on binary outcomes using logistic models. */

foreach y in sedhrs los {
	local ylabel : var label `y'
	local xlabel : var label home
	graph box `y', over(home) title("`xlabel'") name(box_home_`y', replace)
	tabstat `y', by(home) s(n mean median min max) format(%8.3g)
	kwallis `y', by(home)
	}
	
/* So sedation hours and LOS are both clearly longer among patients discharged 
to rehab.

Finally, we will assess associations of intubation with LOS, as well as 
discharge location and seizure time (which is dichotomized at 48 hours). For 
LOS, we can use the preceding loop.  For the dichotomous outcomes, we will 
use cross tabs with chi-square tests, which we can illustrate using barplots. */

foreach y in los {
	local ylabel : var label `y'
	local xlabel : var label intub
	graph box `y', over(intub) title("`xlabel'") name(box_home_`y', replace)
	tabstat `y', by(intub) s(n mean median min max) format(%8.3g)
	kwallis `y', by(intub)
	}

foreach y in home sz_48 {
	local ylabel : var label `y'
	local xlabel : var label intub
	graph bar `y', asc over(intub) ///
		ytitle("`ylabel'") title("`xlabel'") ///
		name(bar_intub_`y', replace)
	tab intub `y', row nokey chi
	}
	
/* Not much evidence that intubation is associated with any of these outcomes. */
	
log close