log using z:documents/teach/courses/b536/hw2Stata.log, replace
* The log file will echo all commands and then store the tabular output
* Graphical output will need to be separately copied to, for instance,
* a Word document

* Note that I can input the file by reading directly from the course web pages.
* I use the "quietly" prefix to suppress the many warnings that arise from
* the use of "NA" as an indicator for missing data. (Stata would not have
* given a warning if I had used the Stata code of "." for missing data.)
quietly: infile id site age male bkrace smoker estrogen prevdis diab2 bmi ///
     systBP aai cholest crp fib ttodth death cvddth using                 ///
	 http://www.emersonstatistics.com/datasets/inflamm.txt

* First line in datafile had variable names
list in 1
drop in 1

* Restrict all analyses to females
drop if male==1

* Remove cases with missing values for estrogen
* (This is for the purposes of the homework. In real life, we would want to
* carefully think through what we would do with cases having missing data.)
tabulate estrogen, missing
drop if estrogen==.

* Create indicator of death attributed to cardiovascular disease within 4 years
* (Again, this is for the purposes of the homework. In real life, we would
* definitely think through what it means to die of other causes prior to
* dying from cardiovascular disease-- a situation termed a "competing risk".)
* This code relies on my prior knowledge that no subjects were missing data
* for either ttodth or cvddth and that no subjects were censored prior to 4 yrs
* If I did not know those facts, I would have to decide how to handle any
* censored data and add code to check for missing data.
g cvddeath4=0
replace cvddeath4=1 if ttodth<= 4*365.25 & cvddth==1

* Formatting variables for descriptive statistics output
format age %8.2f
format estrogen prevdis cvddeath4 %8.3f

* Create a variable modeling interaction between estrogen and prevdis
g estr_prev= estrogen * prevdis

* Descriptive statistics for key variables within groups defined by estrogen
bysort estrogen: tabstat age prevdis cvddeath4, col(stat) ///
      stat(n mean sd min q max) format

* PROBLEM #1: analyses based on RD
*
* OVERVIEW:
* Problem 1a: unadjusted analyses of cvddeath4 - estrogen association
*		Fisher's exact test: tabulate; cs
*		chi-square test: tabulate; cs
*		t-test: ttest
*			- equal variances
*			- unequal variances
*		ordinary linear regression: regress; glm
*			- presuming homoscedasticity
*			- allowing for possible heteroscedasticity
*		weighted linear regression: binreg; glm
*			- estimated binomial weights with standard SE
*			- estimated binomial weights with robust SE
*
* Problem 1b: evidence of effect modification by prevdis
*		descriptive statistics
*		inferential statistics
*
* Problem 1c: evidence of confounding by prevdis
*		Association of estrogen - prevdis in sample
*		Association of prevdis - cvddeath4 after adjustment for estrogen
*
* Problem 1d: analyses of cvddeath4 - estrogen association adjusted for prevdis
*		Adjustment for main effect
*			ordinary linear regression: regress
*				- presuming homoscedasticity
*				- allowing for possible heteroscedasticity
*			weighted linear regression: binreg; glm
*				- estimated binomial weights with standard SE
*				- estimated binomial weights with robust SE
*		Adjustment for main effect and estrogen - prevdis interaction
*			stratified analysis: cs, ir
*			ordinary linear regression: regress
*				- presuming homoscedasticity
*				- allowing for possible heteroscedasticity
*			weighted linear regression: binreg; glm
*				- estimated binomial weights with standard SE
*				- estimated binomial weights with robust SE
*
* Problem 1e: evidence of further confounding by age
*		Association of age - estrogen beyond adjustment for prevdis
*			based on regression
*			graphical assessment
*		Association of age - cvddeath4 after adjustment for estrogen, prevdis
*
* Problem 1f: aassociation of cvddeath4 - estrogen adjusted for prevdis, age
*		Adjustment for main effect (unweighted vs weighted, classical vs robust)
*			Dichotomized
*			Dummy variables
*				Quintiles
*				Scientific intervals (5 year)
*			Grouped linear
*			Continuous linear
*			Quadratic
*			Piecewise linear
*			Splines
*		Adjustment for main effect and age - prevdis interaction
*			(Note could mix and match above models for main effect with
*			 alternative models for interactions)
*

* EXECUTION OF CODE FOR PROBLEM 1
* Problem 1a: unadjusted analyses of cvddeath4 - estrogen association
*		Fisher's exact test: tabulate; cs
tabulate estrogen cvddeath4, row exact
cs cvddeath4 estrogen, exact
*		chi-square test: tabulate; cs
tabulate estrogen cvddeath4, row
cs cvddeath4 estrogen
*		t-test: ttest
*			- equal variances
ttest cvddeath4, by(estrogen)
*			- unequal variances
ttest cvddeath4, by(estrogen) unequal
*		ordinary linear regression: regress; glm
*			- presuming homoscedasticity
regress cvddeath4 estrogen
glm cvddeath4 estrogen, family(gaussian) link(identity)
glm cvddeath4 estrogen
*			- allowing for possible heteroscedasticity
regress cvddeath4 estrogen, robust
glm cvddeath4 estrogen, family(gaussian) link(identity) robust
glm cvddeath4 estrogen, robust
*		weighted linear regression: binreg; glm
*			- estimated binomial weights with standard SE
binreg cvddeath4 estrogen, rd
glm cvddeath4 estrogen, family(binomial) link(identity)
*			- estimated binomial weights with robust SE
binreg cvddeath4 estrogen, rd vce(robust)
glm cvddeath4 estrogen, family(binomial) link(identity) vce(robust)

* Problem 1b: evidence of effect modification by prevdis
*		descriptive statistics
bysort prevdis: tabstat cvddeath4, col(stat) stat(n mean) by(estrogen)
*		inferential statistics
regress cvddeath4 estrogen prevdis estr_prev, robust
glm cvddeath4 estrogen prevdis estr_prev, family(binomial) link(identity)
predict fit0
bysort estrogen prevdis: summ cvddeath4 fit0


* Problem 1c: evidence of confounding by prevdis
*		Association of estrogen - prevdis in sample
tabulate prevdis estrogen, row
*		Association of prevdis - cvddeath4 after adjustment for estrogen
bysort estrogen: regress cvddeath4 prevdis, robust


* Problem 1d: analyses of cvddeath4 - estrogen association adjusted for prevdis
*		Adjustment for main effect
*			ordinary linear regression: regress
*				- allowing for possible heteroscedasticity
regress cvddeath4 estrogen prevdis, robust
*			weighted linear regression: binreg; glm
*				- estimated binomial weights with standard SE
glm cvddeath4 estrogen prevdis, family(binomial) link(identity)
*		Adjustment for main effect and estrogen - prevdis interaction
*			stratified analysis: cs
cs cvddeath4 estrogen, rd by(prevdis) istandard
cs cvddeath4 estrogen, rd by(prevdis) estandard
g std= 544 / 2899
replace std= 1 - std if prevdis==0
cs cvddeath4 estrogen, rd by(prevdis) standard(std)
drop std
g std= 1 / ( (1 / 2045) + (1 / 310) )
replace std = 1 / ( (1 / 514) + (1 / 30) ) if prevdis==1
cs cvddeath4 estrogen, rd by(prevdis) standard(std)
*			ordinary linear regression: regress
*				- allowing for possible heteroscedasticity
regress cvddeath4 estrogen prevdis estr_prev, robust
test estrogen estr_prev
lincom (269.1932 / (269.1932 + 28.34559)) * estrogen + (28.34559 / (269.1932 + 28.34559)) * (estrogen + estr_prev)

* Problem 1e: evidence of further confounding by age
*		Association of age - estrogen beyond adjustment for prevdis
*			based on regression
bysort prevdis: regress estrogen age, robust
*			graphical assessment
egen p4 = mean(cvddeath4), by(estrogen prevdis age)
egen n4 = count(cvddeath4), by(estrogen prevdis age)
scatter n4 age if estrogen==0 & prevdis==0, col(black) t1("N: estrogen 0 prevdis 0") name(n00)
scatter n4 age if estrogen==0 & prevdis==1, col(blue) t1("N: estrogen 0 prevdis 1") name(n01)
scatter n4 age if estrogen==1 & prevdis==1, col(green) t1("N: estrogen 1 prevdis 1") name(n11)
scatter n4 age if estrogen==1 & prevdis==0, col(red) t1("N: estrogen 1 prevdis 0") name(n10)
graph combine n10 n11 n00 n01, t1("Sample Size by Age within Strata")
scatter p4 age if estrogen==0 & prevdis==0, col(black) t1("Mortality: estrogen 0 prevdis 0") name(p00)
scatter p4 age if estrogen==0 & prevdis==1, col(blue) t1("Mortality: estrogen 0 prevdis 1") name(p01)
scatter p4 age if estrogen==1 & prevdis==1, col(green) t1("Mortality: estrogen 1 prevdis 0") name(p10)
scatter p4 age if estrogen==1 & prevdis==0, col(red) t1("Mortality: estrogen 1 prevdis 1") name(p11)
graph combine p10 p11 p00 p01, t1("Estimated Mortality by Age within Strata")
*		Association of age - cvddeath4 after adjustment for estrogen, prevdis
*			based on regression
*			graphical assessment

* Problem 1f: association of cvddeath4 - estrogen adjusted for prevdis, age
* Create variables to facilitate more flexible modeling of age
*       age squared to be used in a quadratic fit
g agesqr= age ^2
*       splitting at median
xtile ageQ2 = age, nq(2)
*       coding within quintiles to be used as dummny variables
xtile ageQ5 = age, nq(5)
*       coding within (approx) 5 year age groups to be used as dummy variables
g age5yr = age
recode age5yr 65/69=1 70/74=2 75/79=3 80/84=4 85/max=5
*		I prefer to code variables by mean value within interval
tabstat age, by(age5yr) col(stat) stat(n mean sd min q max) format
recode age5yr 1=67 2=72 3=77 4=82 5=88
*       creating linear splines within (approx) 5 year age groups
mkspline age65 69.5 age70 74.5 age75 79.5 age80 84.5 age85 = age

*		Adjustment for main effect (unweighted vs weighted, classical vs robust)
*			Dichotomized
regress cvddeath4 estrogen prevdis ageQ2, robust
predict fit1
* (see do file for code used to produce graphs)
graph drop _all
twoway ///
(scatter fit1 age if estrogen==0 & prevdis==0, msymbol(i) lcol(black) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==0, col(black)), /// 
t1("fit1: estrogen 0 prevdis 0 (N= 2045)") name(f00)
twoway ///
(scatter fit1 age if estrogen==0 & prevdis==1, msymbol(i) lcol(blue) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==1, col(blue)), /// 
t1("fit1: estrogen 0 prevdis 1 (N= 514)") name(f01)
twoway ///
(scatter fit1 age if estrogen==1 & prevdis==1, msymbol(i) lcol(green) connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==1, col(green)), ///
 t1("fit1: estrogen 1 prevdis 1 (N= 30)") name(f11)
 twoway ///
(scatter fit1 age if estrogen==1 & prevdis==0, msymbol(i) lcol(red)  connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==0, col(red)), /// 
t1("fit1: estrogen 1 prevdis 0 (N= 310)") name(f10)
graph combine f10 f11 f00 f01, ycommon t1("fit1: Age Dichotomized at Median") name(panel)
twoway ///
(scatter fit1 age if estrogen==0 & prevdis==0, col(black)) ///
(scatter fit1 age if estrogen==0 & prevdis==1, col(blue)) ///
(scatter fit1 age if estrogen==1 & prevdis==0, col(red)) ///
(scatter fit1 age if estrogen==1 & prevdis==1, col(green)), ///
t1("fit1: Age Dichotomized at Median") ///
legend(label(1 "No estrogen, no CVD") label(2 "No estrogen, prior CVD") label(4 "Estrogen, prior CVD") label(3 "Estrogen, no CVD"))

*			Dummy variables : Quintiles
regress cvddeath4 estrogen prevdis i.ageQ5, robust
predict fit2
* (see do file for code used to produce graphs)
graph drop _all
twoway ///
(scatter fit2 age if estrogen==0 & prevdis==0, msymbol(i) lcol(black) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==0, col(black)), /// 
t1("fit2: estrogen 0 prevdis 0 (N= 2045)") name(f00)
twoway ///
(scatter fit2 age if estrogen==0 & prevdis==1, msymbol(i) lcol(blue) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==1, col(blue)), /// 
t1("fit2: estrogen 0 prevdis 1 (N= 514)") name(f01)
twoway ///
(scatter fit2 age if estrogen==1 & prevdis==1, msymbol(i) lcol(green) connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==1, col(green)), ///
 t1("fit2: estrogen 1 prevdis 1 (N= 30)") name(f11)
 twoway ///
(scatter fit2 age if estrogen==1 & prevdis==0, msymbol(i) lcol(red)  connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==0, col(red)), /// 
t1("fit2: estrogen 1 prevdis 0 (N= 310)") name(f10)
graph combine f10 f11 f00 f01, ycommon t1("fit2: Age Categorized at Quintiles") name(panel)
twoway ///
(scatter fit2 age if estrogen==0 & prevdis==0, col(black)) ///
(scatter fit2 age if estrogen==0 & prevdis==1, col(blue)) ///
(scatter fit2 age if estrogen==1 & prevdis==0, col(red)) ///
(scatter fit2 age if estrogen==1 & prevdis==1, col(green)), ///
t1("fit2: Age Categorized at Quintiles") ///
legend(label(1 "No estrogen, no CVD") label(2 "No estrogen, prior CVD") label(4 "Estrogen, prior CVD") label(3 "Estrogen, no CVD"))

*			Dummy variables : Scientific intervals (5 year)
regress cvddeath4 estrogen prevdis i.age5yr, robust
predict fit3
* (see do file for code used to produce graphs)
graph drop _all
twoway ///
(scatter fit3 age if estrogen==0 & prevdis==0, msymbol(i) lcol(black) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==0, col(black)), /// 
t1("fit3: estrogen 0 prevdis 0 (N= 2045)") name(f00)
twoway ///
(scatter fit3 age if estrogen==0 & prevdis==1, msymbol(i) lcol(blue) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==1, col(blue)), /// 
t1("fit3: estrogen 0 prevdis 1 (N= 514)") name(f01)
twoway ///
(scatter fit3 age if estrogen==1 & prevdis==1, msymbol(i) lcol(green) connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==1, col(green)), ///
 t1("fit3: estrogen 1 prevdis 1 (N= 30)") name(f11)
 twoway ///
(scatter fit3 age if estrogen==1 & prevdis==0, msymbol(i) lcol(red)  connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==0, col(red)), /// 
t1("fit3: estrogen 1 prevdis 0 (N= 310)") name(f10)
graph combine f10 f11 f00 f01, ycommon t1("fit3: Age Categorized in 5 Year Intervals") name(panel)
twoway ///
(scatter fit3 age if estrogen==0 & prevdis==0, col(black)) ///
(scatter fit3 age if estrogen==0 & prevdis==1, col(blue)) ///
(scatter fit3 age if estrogen==1 & prevdis==0, col(red)) ///
(scatter fit3 age if estrogen==1 & prevdis==1, col(green)), ///
t1("fit3: Age Categorized in 5 Year Intervals") ///
legend(label(1 "No estrogen, no CVD") label(2 "No estrogen, prior CVD") label(4 "Estrogen, prior CVD") label(3 "Estrogen, no CVD"))

*			Grouped linear
regress cvddeath4 estrogen prevdis age5yr, robust
predict fit4
* (see do file for code used to produce graphs)
graph drop _all
twoway ///
(scatter fit4 age if estrogen==0 & prevdis==0, msymbol(i) lcol(black) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==0, col(black)), /// 
t1("fit4: estrogen 0 prevdis 0 (N= 2045)") name(f00)
twoway ///
(scatter fit4 age if estrogen==0 & prevdis==1, msymbol(i) lcol(blue) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==1, col(blue)), /// 
t1("fit4: estrogen 0 prevdis 1 (N= 514)") name(f01)
twoway ///
(scatter fit4 age if estrogen==1 & prevdis==1, msymbol(i) lcol(green) connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==1, col(green)), ///
 t1("fit4: estrogen 1 prevdis 1 (N= 30)") name(f11)
 twoway ///
(scatter fit4 age if estrogen==1 & prevdis==0, msymbol(i) lcol(red)  connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==0, col(red)), /// 
t1("fit4: estrogen 1 prevdis 0 (N= 310)") name(f10)
graph combine f10 f11 f00 f01, ycommon t1("fit4: Age 5 Year Intervals Grouped Linear") name(panel)
twoway ///
(scatter fit4 age if estrogen==0 & prevdis==0, col(black)) ///
(scatter fit4 age if estrogen==0 & prevdis==1, col(blue)) ///
(scatter fit4 age if estrogen==1 & prevdis==0, col(red)) ///
(scatter fit4 age if estrogen==1 & prevdis==1, col(green)), ///
t1("fit4: Age 5 Year Intervals Grouped Linear") ///
legend(label(1 "No estrogen, no CVD") label(2 "No estrogen, prior CVD") label(4 "Estrogen, prior CVD") label(3 "Estrogen, no CVD"))

*			Continuous linear
regress cvddeath4 estrogen prevdis age, robust
predict fit5
* (see do file for code used to produce graphs)
graph drop _all
twoway ///
(scatter fit5 age if estrogen==0 & prevdis==0, msymbol(i) lcol(black) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==0, col(black)), /// 
t1("fit5: estrogen 0 prevdis 0 (N= 2045)") name(f00)
twoway ///
(scatter fit5 age if estrogen==0 & prevdis==1, msymbol(i) lcol(blue) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==1, col(blue)), /// 
t1("fit5: estrogen 0 prevdis 1 (N= 514)") name(f01)
twoway ///
(scatter fit5 age if estrogen==1 & prevdis==1, msymbol(i) lcol(green) connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==1, col(green)), ///
 t1("fit5: estrogen 1 prevdis 1 (N= 30)") name(f11)
 twoway ///
(scatter fit5 age if estrogen==1 & prevdis==0, msymbol(i) lcol(red)  connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==0, col(red)), /// 
t1("fit5: estrogen 1 prevdis 0 (N= 310)") name(f10)
graph combine f10 f11 f00 f01, ycommon t1("fit5: Age Continuous Linear") name(panel)
twoway ///
(scatter fit5 age if estrogen==0 & prevdis==0, col(black)) ///
(scatter fit5 age if estrogen==0 & prevdis==1, col(blue)) ///
(scatter fit5 age if estrogen==1 & prevdis==0, col(red)) ///
(scatter fit5 age if estrogen==1 & prevdis==1, col(green)), ///
t1("fit5: Age Continuous Linear") ///
legend(label(1 "No estrogen, no CVD") label(2 "No estrogen, prior CVD") label(4 "Estrogen, prior CVD") label(3 "Estrogen, no CVD"))

*			Quadratic
regress cvddeath4 estrogen prevdis age agesqr, robust
predict fit6
* (see do file for code used to produce graphs)
graph drop _all
twoway ///
(scatter fit6 age if estrogen==0 & prevdis==0, msymbol(i) lcol(black) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==0, col(black)), /// 
t1("fit6: estrogen 0 prevdis 0 (N= 2045)") name(f00)
twoway ///
(scatter fit6 age if estrogen==0 & prevdis==1, msymbol(i) lcol(blue) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==1, col(blue)), /// 
t1("fit6: estrogen 0 prevdis 1 (N= 514)") name(f01)
twoway ///
(scatter fit6 age if estrogen==1 & prevdis==1, msymbol(i) lcol(green) connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==1, col(green)), ///
 t1("fit6: estrogen 1 prevdis 1 (N= 30)") name(f11)
 twoway ///
(scatter fit6 age if estrogen==1 & prevdis==0, msymbol(i) lcol(red)  connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==0, col(red)), /// 
t1("fit6: estrogen 1 prevdis 0 (N= 310)") name(f10)
graph combine f10 f11 f00 f01, ycommon t1("fit6: Age Quadratic") name(panel)
twoway ///
(scatter fit6 age if estrogen==0 & prevdis==0, col(black)) ///
(scatter fit6 age if estrogen==0 & prevdis==1, col(blue)) ///
(scatter fit6 age if estrogen==1 & prevdis==0, col(red)) ///
(scatter fit6 age if estrogen==1 & prevdis==1, col(green)), ///
t1("fit6: Age Quadratic") ///
legend(label(1 "No estrogen, no CVD") label(2 "No estrogen, prior CVD") label(4 "Estrogen, prior CVD") label(3 "Estrogen, no CVD"))

*			Piecewise linear
regress cvddeath4 estrogen prevdis i.age5yr##c.age, robust
predict fit7
* (see do file for code used to produce graphs)
graph drop _all
twoway ///
(scatter fit7 age if estrogen==0 & prevdis==0, msymbol(i) lcol(black) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==0, col(black)), /// 
t1("fit7: estrogen 0 prevdis 0 (N= 2045)") name(f00)
twoway ///
(scatter fit7 age if estrogen==0 & prevdis==1, msymbol(i) lcol(blue) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==1, col(blue)), /// 
t1("fit7: estrogen 0 prevdis 1 (N= 514)") name(f01)
twoway ///
(scatter fit7 age if estrogen==1 & prevdis==1, msymbol(i) lcol(green) connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==1, col(green)), ///
 t1("fit7: estrogen 1 prevdis 1 (N= 30)") name(f11)
 twoway ///
(scatter fit7 age if estrogen==1 & prevdis==0, msymbol(i) lcol(red)  connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==0, col(red)), /// 
t1("fit7: estrogen 1 prevdis 0 (N= 310)") name(f10)
graph combine f10 f11 f00 f01, ycommon t1("fit7: Age Piecewise Linear") name(panel)
twoway ///
(scatter fit7 age if estrogen==0 & prevdis==0, col(black)) ///
(scatter fit7 age if estrogen==0 & prevdis==1, col(blue)) ///
(scatter fit7 age if estrogen==1 & prevdis==0, col(red)) ///
(scatter fit7 age if estrogen==1 & prevdis==1, col(green)), ///
t1("fit7: Age Piecewise Linear") ///
legend(label(1 "No estrogen, no CVD") label(2 "No estrogen, prior CVD") label(4 "Estrogen, prior CVD") label(3 "Estrogen, no CVD"))

*			Splines
regress cvddeath4 estrogen prevdis age65 age70 age75 age80 age85, robust
scatter age65 age, msymbol(i) lcol(black) connect(l) sort(age) t1("age65 vs age") name(g1)
scatter age70 age, msymbol(i) lcol(blue) connect(l) sort(age) t1("age70 vs age") name(g2)
scatter age75 age, msymbol(i) lcol(green) connect(l) sort(age) t1("age75 vs age") name(g3)
scatter age80 age, msymbol(i) lcol(red) connect(l) sort(age) t1("age80 vs age") name(g4)
scatter age85 age, msymbol(i) lcol(orange) connect(l) sort(age) t1("age85 vs age") name(g5)
graph combine g1 g2 g3 g4 g5, t1("Coding of Spline Variables vs Age")

predict fit8
* (see do file for code used to produce graphs)
graph drop _all
twoway ///
(scatter fit8 age if estrogen==0 & prevdis==0, msymbol(i) lcol(black) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==0, col(black)), /// 
t1("fit8: estrogen 0 prevdis 0 (N= 2045)") name(f00)
twoway ///
(scatter fit8 age if estrogen==0 & prevdis==1, msymbol(i) lcol(blue) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==1, col(blue)), /// 
t1("fit8: estrogen 0 prevdis 1 (N= 514)") name(f01)
twoway ///
(scatter fit8 age if estrogen==1 & prevdis==1, msymbol(i) lcol(green) connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==1, col(green)), ///
 t1("fit8: estrogen 1 prevdis 1 (N= 30)") name(f11)
 twoway ///
(scatter fit8 age if estrogen==1 & prevdis==0, msymbol(i) lcol(red)  connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==0, col(red)), /// 
t1("fit8: estrogen 1 prevdis 0 (N= 310)") name(f10)
graph combine f10 f11 f00 f01, ycommon t1("fit8: Age Linear Splines") name(panel)
twoway ///
(scatter fit8 age if estrogen==0 & prevdis==0, col(black)) ///
(scatter fit8 age if estrogen==0 & prevdis==1, col(blue)) ///
(scatter fit8 age if estrogen==1 & prevdis==0, col(red)) ///
(scatter fit8 age if estrogen==1 & prevdis==1, col(green)), ///
t1("fit8: Age Linear Splines") ///
legend(label(1 "No estrogen, no CVD") label(2 "No estrogen, prior CVD") label(4 "Estrogen, prior CVD") label(3 "Estrogen, no CVD"))

*		Adjustment for main effect and age - prevdis interaction
*			(Note could mix and match above models for main effect with
*			 alternative models for interactions)
*			Continuous linear and interaction
regress cvddeath4 estrogen i.prevdis##c.age, robust
predict fit9
* (see do file for code used to produce graphs)
graph drop _all
twoway ///
(scatter fit9 age if estrogen==0 & prevdis==0, msymbol(i) lcol(black) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==0, col(black)), /// 
t1("fit9: estrogen 0 prevdis 0 (N= 2045)") name(f00)
twoway ///
(scatter fit9 age if estrogen==0 & prevdis==1, msymbol(i) lcol(blue) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==1, col(blue)), /// 
t1("fit9: estrogen 0 prevdis 1 (N= 514)") name(f01)
twoway ///
(scatter fit9 age if estrogen==1 & prevdis==1, msymbol(i) lcol(green) connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==1, col(green)), ///
 t1("fit9: estrogen 1 prevdis 1 (N= 30)") name(f11)
 twoway ///
(scatter fit9 age if estrogen==1 & prevdis==0, msymbol(i) lcol(red)  connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==0, col(red)), /// 
t1("fit9: estrogen 1 prevdis 0 (N= 310)") name(f10)
graph combine f10 f11 f00 f01, ycommon t1("fit9: Age Continuous Linear Interaction") name(panel)
twoway ///
(scatter fit9 age if estrogen==0 & prevdis==0, col(black)) ///
(scatter fit9 age if estrogen==0 & prevdis==1, col(blue)) ///
(scatter fit9 age if estrogen==1 & prevdis==0, col(red)) ///
(scatter fit9 age if estrogen==1 & prevdis==1, col(green)), ///
t1("fit9: Age Continuous Linear Interaction") ///
legend(label(1 "No estrogen, no CVD") label(2 "No estrogen, prior CVD") label(4 "Estrogen, prior CVD") label(3 "Estrogen, no CVD"))

*			Quadratic and interaction
regress cvddeath4 estrogen i.prevdis##c.age i.prevdis##c.agesqr, robust
predict fit10
* (see do file for code used to produce graphs)
graph drop _all
twoway ///
(scatter fit10 age if estrogen==0 & prevdis==0, msymbol(i) lcol(black) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==0, col(black)), /// 
t1("fit10: estrogen 0 prevdis 0 (N= 2045)") name(f00)
twoway ///
(scatter fit10 age if estrogen==0 & prevdis==1, msymbol(i) lcol(blue) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==1, col(blue)), /// 
t1("fit10: estrogen 0 prevdis 1 (N= 514)") name(f01)
twoway ///
(scatter fit10 age if estrogen==1 & prevdis==1, msymbol(i) lcol(green) connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==1, col(green)), ///
 t1("fit10: estrogen 1 prevdis 1 (N= 30)") name(f11)
 twoway ///
(scatter fit10 age if estrogen==1 & prevdis==0, msymbol(i) lcol(red)  connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==0, col(red)), /// 
t1("fit10: estrogen 1 prevdis 0 (N= 310)") name(f10)
graph combine f10 f11 f00 f01, ycommon t1("fit10: Age Quadratic Interaction") name(panel)
twoway ///
(scatter fit10 age if estrogen==0 & prevdis==0, col(black)) ///
(scatter fit10 age if estrogen==0 & prevdis==1, col(blue)) ///
(scatter fit10 age if estrogen==1 & prevdis==0, col(red)) ///
(scatter fit10 age if estrogen==1 & prevdis==1, col(green)), ///
t1("fit10: Age Quadratic Interaction") ///
legend(label(1 "No estrogen, no CVD") label(2 "No estrogen, prior CVD") label(4 "Estrogen, prior CVD") label(3 "Estrogen, no CVD"))

*		Adjustment for three-way interaction
regress cvddeath4 i.estrogen##prevdis##c.age i.estrogen##prevdis##c.agesqr, robust
test 1.estrogen 1.estrogen#1.prevdis 1.estrogen#c.age 1.estrogen#c.agesqr 1.estrogen#1.prevdis#c.age 1.estrogen#1.prevdis#c.agesqr
predict fit11
* (see do file for code used to produce graphs)
graph drop _all
twoway ///
(scatter fit11 age if estrogen==0 & prevdis==0, msymbol(i) lcol(black) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==0, col(black)), /// 
t1("fit11: estrogen 0 prevdis 0 (N= 2045)") name(f00)
twoway ///
(scatter fit11 age if estrogen==0 & prevdis==1, msymbol(i) lcol(blue) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==1, col(blue)), /// 
t1("fit11: estrogen 0 prevdis 1 (N= 514)") name(f01)
twoway ///
(scatter fit11 age if estrogen==1 & prevdis==1, msymbol(i) lcol(green) connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==1, col(green)), ///
 t1("fit11: estrogen 1 prevdis 1 (N= 30)") name(f11)
 twoway ///
(scatter fit11 age if estrogen==1 & prevdis==0, msymbol(i) lcol(red)  connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==0, col(red)), /// 
t1("fit11: estrogen 1 prevdis 0 (N= 310)") name(f10)
graph combine f10 f11 f00 f01, ycommon t1("fit11: Three-way Interaction") name(panel)
twoway ///
(scatter fit11 age if estrogen==0 & prevdis==0, col(black)) ///
(scatter fit11 age if estrogen==0 & prevdis==1, col(blue)) ///
(scatter fit11 age if estrogen==1 & prevdis==0, col(red)) ///
(scatter fit11 age if estrogen==1 & prevdis==1, col(green)), ///
t1("fit11: Three-way Interaction") ///
legend(label(1 "No estrogen, no CVD") label(2 "No estrogen, prior CVD") label(4 "Estrogen, prior CVD") label(3 "Estrogen, no CVD"))

drop fit1 fit2 fit3 fit4 fit5 fit6 fit7 fit8 fit9 fit10 fit11

* PROBLEM #2: analyses based on RR
*
* OVERVIEW:
* Problem 2a: unadjusted analyses of cvddeath4 - estrogen association
*		chi-square test: cs
*		generalized linear model: glm
*			- presuming homoscedasticity
*			- allowing for possible heteroscedasticity
*		weighted generalized linear model: glm
*			- estimated binomial weights with standard SE
*			- estimated binomial weights with robust SE
*
* Problem 2b: evidence of effect modification by prevdis
*		descriptive statistics
*		inferential statistics
*
* Problem 2c: evidence of confounding by prevdis
*		Association of estrogen - prevdis in sample
*		Association of prevdis - cvddeath4 after adjustment for estrogen
*
* Problem 2d: analyses of cvddeath4 - estrogen association adjusted for prevdis
*		Adjustment for main effect
*			generalized linear model: glm
*				- presuming homoscedasticity
*				- allowing for possible heteroscedasticity
*			weighted generalized linear model: glm
*				- estimated binomial weights with standard SE
*				- estimated binomial weights with robust SE
*		Adjustment for main effect and estrogen - prevdis interaction
*			stratified analysis: cs
*			generalized linear model: glm
*				- presuming homoscedasticity
*				- allowing for possible heteroscedasticity
*			weighted generalized linear model: glm
*				- estimated binomial weights with standard SE
*				- estimated binomial weights with robust SE
*
* Problem 2e: evidence of further confounding by age
*		Association of age - estrogen beyond adjustment for prevdis
*			based on regression
*			graphical assessment
*		Association of age - cvddeath4 after adjustment for estrogen, prevdis
*
* Problem 2f: association of cvddeath4 - estrogen adjusted for prevdis, age
*		Adjustment for main effect (unweighted vs weighted, classical vs robust)
*			Dichotomized
*			Dummy variables
*				Quintiles
*				Scientific intervals (5 year)
*			Grouped linear
*			Continuous linear
*			Quadratic
*			Piecewise linear
*			Splines
*		Adjustment for main effect and age - prevdis interaction
*			(Note could mix and match above models for main effect with
*			 alternative models for interactions)
*

* EXECUTION OF CODE FOR PROBLEM 2
* Problem 2a: unadjusted analyses of cvddeath4 - estrogen association
*		chi-square test: cs
cs cvddeath4 estrogen
di  (3 / 340)   /    (88 / 2559)
*		generalized linear model: glm
*			- presuming homoscedasticity
glm cvddeath4 estrogen, family(gaussian) link(log)
di   exp(-3.370034), exp( -1.360267)
di   -1.360267 / 1.075538
di   -1.360267  -  1.96 * 1.075538,         -1.360267    +    1.96 * 1.075538
di   exp( -1.360267  -  1.96 * 1.075538),         exp(-1.360267    +    1.96 * 1.075538)
glm cvddeath4 estrogen, family(gaussian) link(log) eform
di .2565924 -  1.96 * .2759748,  .2565924 +  1.96 * .2759748    
glm cvddeath4 estrogen, family(gaussian) link(log) eform irls
*			- allowing for possible heteroscedasticity
glm cvddeath4 estrogen, family(gaussian) link(log) robust
glm cvddeath4 estrogen, family(gaussian) link(log) robust eform
*		weighted generalized linear model: glm
*			- estimated binomial weights with standard SE
binreg cvddeath4 estrogen, rr
glm cvddeath4 estrogen, family(binomial) link(log) eform
*			- estimated binomial weights with robust SE
binreg cvddeath4 estrogen, rr vce(robust)
glm cvddeath4 estrogen, family(binomial) link(log) vce(robust) eform
poisson cvddeath4 estrogen, irr

* Problem 2b: evidence of effect modification by prevdis
*		descriptive statistics
cs cvddeath4 estrogen, by(prevdis)
*		inferential statistics
glm cvddeath4 estrogen prevdis estr_prev, family(gaussian) link(log) robust
predict fit0
bysort estrogen prevdis: summ fit0 cvddeath4
glm cvddeath4 estrogen prevdis estr_prev, family(gaussian) link(log) eform robust



* Problem 2c: evidence of confounding by prevdis
* See problem 1

* Problem 2d: analyses of cvddeath4 - estrogen association adjusted for prevdis
*		Adjustment for main effect
*			ordinary linear regression: regress
*				- allowing for possible heteroscedasticity
glm cvddeath4 estrogen prevdis, family(gaussian) link(log) robust eform
poisson cvddeath4 estrogen prevdis, irr
*		Adjustment for main effect and estrogen - prevdis interaction
*			stratified analysis: cs
cs cvddeath4 estrogen, by(prevdis) istandard
cs cvddeath4 estrogen, by(prevdis) estandard
g std= 544 / 2899
replace std= 1 - std if prevdis==0
cs cvddeath4 estrogen, by(prevdis) standard(std)
*			ordinary linear regression: regress
*				- allowing for possible heteroscedasticity
glm cvddeath4 estrogen prevdis estr_prev, family(gaussian) link(log) robust eform
test estrogen estr_prev
lincom (1 - 544/2899) * estrogen + (544/2899) * (estrogen + estr_prev), eform
test estrogen estr_prev
lincom (1 - 544/2899) * estrogen + (544/2899) * (estrogen + estr_prev), eform
* Weight for no prior history stratum
di (1 - 544/2899)
* Average risk in no estrogen stratum
di .81234909 * exp(-4.012235) + (1-.81234909) * exp(-4.012235+1.701838)
* Average risk in estrogen exposure stratum
di .81234909 * exp(-4.012235+  -1.03119) + (1-.81234909) * exp(-4.012235+1.701838+  -1.03119 + -.0596098)
* Relative average risk
di .011496/.03331683
* Inference using nlcom
nlcom ( (1 - 544/2899) * exp(_b[_cons] + _b[estrogen]) + ///
	(544/2899) * exp(_b[_cons] + _b[prevdis] + _b[estrogen] + _b[estr_prev]) ) ///
	/ ((1 - 544/2899) * exp(_b[_cons]) + (544/2899) * exp(_b[_cons] + _b[prevdis]))
poisson cvddeath4 estrogen prevdis estr_prev, irr

* Problem 2e: evidence of further confounding by age
*		(see problem 1)

* Problem 2f: association of cvddeath4 - estrogen adjusted for prevdis, age
*		Adjustment for main effect (unweighted vs weighted, classical vs robust
*			Continuous linear (GLM Gaussian)
glm cvddeath4 estrogen prevdis age, robust link(log) eform
predict fit5
* (see do file for code used to produce graphs)
graph drop _all
twoway ///
(scatter fit5 age if estrogen==0 & prevdis==0, msymbol(i) lcol(black) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==0, col(black)), /// 
t1("fit5: estrogen 0 prevdis 0 (N= 2045)") name(f00)
twoway ///
(scatter fit5 age if estrogen==0 & prevdis==1, msymbol(i) lcol(blue) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==1, col(blue)), /// 
t1("fit5: estrogen 0 prevdis 1 (N= 514)") name(f01)
twoway ///
(scatter fit5 age if estrogen==1 & prevdis==1, msymbol(i) lcol(green) connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==1, col(green)), ///
 t1("fit5: estrogen 1 prevdis 1 (N= 30)") name(f11)
 twoway ///
(scatter fit5 age if estrogen==1 & prevdis==0, msymbol(i) lcol(red)  connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==0, col(red)), /// 
t1("fit5: estrogen 1 prevdis 0 (N= 310)") name(f10)
graph combine f10 f11 f00 f01, ycommon t1("fit5: Age Continuous Linear") name(panel)
twoway ///
(scatter fit5 age if estrogen==0 & prevdis==0, col(black)) ///
(scatter fit5 age if estrogen==0 & prevdis==1, col(blue)) ///
(scatter fit5 age if estrogen==1 & prevdis==0, col(red)) ///
(scatter fit5 age if estrogen==1 & prevdis==1, col(green)), ///
t1("fit5: Age Continuous Linear") ///
legend(label(1 "No estrogen, no CVD") label(2 "No estrogen, prior CVD") label(4 "Estrogen, prior CVD") label(3 "Estrogen, no CVD"))

*			Continuous linear (Poisson)
poisson cvddeath4 estrogen prevdis age, irr
drop fit5
predict fit5
* (see do file for code used to produce graphs)
graph drop _all
twoway ///
(scatter fit5 age if estrogen==0 & prevdis==0, msymbol(i) lcol(black) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==0, col(black)), /// 
t1("fit5: estrogen 0 prevdis 0 (N= 2045)") name(f00)
twoway ///
(scatter fit5 age if estrogen==0 & prevdis==1, msymbol(i) lcol(blue) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==1, col(blue)), /// 
t1("fit5: estrogen 0 prevdis 1 (N= 514)") name(f01)
twoway ///
(scatter fit5 age if estrogen==1 & prevdis==1, msymbol(i) lcol(green) connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==1, col(green)), ///
 t1("fit5: estrogen 1 prevdis 1 (N= 30)") name(f11)
 twoway ///
(scatter fit5 age if estrogen==1 & prevdis==0, msymbol(i) lcol(red)  connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==0, col(red)), /// 
t1("fit5: estrogen 1 prevdis 0 (N= 310)") name(f10)
graph combine f10 f11 f00 f01, ycommon t1("fit5: Age Continuous Linear- Poisson") name(panel)
twoway ///
(scatter fit5 age if estrogen==0 & prevdis==0, col(black)) ///
(scatter fit5 age if estrogen==0 & prevdis==1, col(blue)) ///
(scatter fit5 age if estrogen==1 & prevdis==0, col(red)) ///
(scatter fit5 age if estrogen==1 & prevdis==1, col(green)), ///
t1("fit5: Age Continuous Linear - Poisson") ///
legend(label(1 "No estrogen, no CVD") label(2 "No estrogen, prior CVD") label(4 "Estrogen, prior CVD") label(3 "Estrogen, no CVD"))

*			Quadratic
glm cvddeath4 estrogen prevdis age agesqr, robust link(log) eform
predict fit6
* (see do file for code used to produce graphs)
graph drop _all
twoway ///
(scatter fit6 age if estrogen==0 & prevdis==0, msymbol(i) lcol(black) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==0, col(black)), /// 
t1("fit6: estrogen 0 prevdis 0 (N= 2045)") name(f00)
twoway ///
(scatter fit6 age if estrogen==0 & prevdis==1, msymbol(i) lcol(blue) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==1, col(blue)), /// 
t1("fit6: estrogen 0 prevdis 1 (N= 514)") name(f01)
twoway ///
(scatter fit6 age if estrogen==1 & prevdis==1, msymbol(i) lcol(green) connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==1, col(green)), ///
 t1("fit6: estrogen 1 prevdis 1 (N= 30)") name(f11)
 twoway ///
(scatter fit6 age if estrogen==1 & prevdis==0, msymbol(i) lcol(red)  connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==0, col(red)), /// 
t1("fit6: estrogen 1 prevdis 0 (N= 310)") name(f10)
graph combine f10 f11 f00 f01, ycommon t1("fit6: Age Quadratic") name(panel)
twoway ///
(scatter fit6 age if estrogen==0 & prevdis==0, col(black)) ///
(scatter fit6 age if estrogen==0 & prevdis==1, col(blue)) ///
(scatter fit6 age if estrogen==1 & prevdis==0, col(red)) ///
(scatter fit6 age if estrogen==1 & prevdis==1, col(green)), ///
t1("fit6: Age Quadratic") ///
legend(label(1 "No estrogen, no CVD") label(2 "No estrogen, prior CVD") label(4 "Estrogen, prior CVD") label(3 "Estrogen, no CVD"))

*			Splines
glm cvddeath4 estrogen prevdis age65 age70 age75 age80 age85, robust link(log) eform
predict fit8
* (see do file for code used to produce graphs)
graph drop _all
twoway ///
(scatter fit8 age if estrogen==0 & prevdis==0, msymbol(i) lcol(black) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==0, col(black)), /// 
t1("fit8: estrogen 0 prevdis 0 (N= 2045)") name(f00)
twoway ///
(scatter fit8 age if estrogen==0 & prevdis==1, msymbol(i) lcol(blue) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==1, col(blue)), /// 
t1("fit8: estrogen 0 prevdis 1 (N= 514)") name(f01)
twoway ///
(scatter fit8 age if estrogen==1 & prevdis==1, msymbol(i) lcol(green) connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==1, col(green)), ///
 t1("fit8: estrogen 1 prevdis 1 (N= 30)") name(f11)
 twoway ///
(scatter fit8 age if estrogen==1 & prevdis==0, msymbol(i) lcol(red)  connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==0, col(red)), /// 
t1("fit8: estrogen 1 prevdis 0 (N= 310)") name(f10)
graph combine f10 f11 f00 f01, ycommon t1("fit8: Age Linear Splines") name(panel)
twoway ///
(scatter fit8 age if estrogen==0 & prevdis==0, col(black)) ///
(scatter fit8 age if estrogen==0 & prevdis==1, col(blue)) ///
(scatter fit8 age if estrogen==1 & prevdis==0, col(red)) ///
(scatter fit8 age if estrogen==1 & prevdis==1, col(green)), ///
t1("fit8: Age Linear Splines") ///
legend(label(1 "No estrogen, no CVD") label(2 "No estrogen, prior CVD") label(4 "Estrogen, prior CVD") label(3 "Estrogen, no CVD"))

*			Splines: Poisson regression
poisson cvddeath4 estrogen prevdis age65 age70 age75 age80 age85, irr
drop fit8
predict fit8
* (see do file for code used to produce graphs)
graph drop _all
twoway ///
(scatter fit8 age if estrogen==0 & prevdis==0, msymbol(i) lcol(black) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==0, col(black)), /// 
t1("fit8: estrogen 0 prevdis 0 (N= 2045)") name(f00)
twoway ///
(scatter fit8 age if estrogen==0 & prevdis==1, msymbol(i) lcol(blue) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==1, col(blue)), /// 
t1("fit8: estrogen 0 prevdis 1 (N= 514)") name(f01)
twoway ///
(scatter fit8 age if estrogen==1 & prevdis==1, msymbol(i) lcol(green) connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==1, col(green)), ///
 t1("fit8: estrogen 1 prevdis 1 (N= 30)") name(f11)
 twoway ///
(scatter fit8 age if estrogen==1 & prevdis==0, msymbol(i) lcol(red)  connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==0, col(red)), /// 
t1("fit8: estrogen 1 prevdis 0 (N= 310)") name(f10)
graph combine f10 f11 f00 f01, ycommon t1("fit8: Age Linear Splines") name(panel)
twoway ///
(scatter fit8 age if estrogen==0 & prevdis==0, col(black)) ///
(scatter fit8 age if estrogen==0 & prevdis==1, col(blue)) ///
(scatter fit8 age if estrogen==1 & prevdis==0, col(red)) ///
(scatter fit8 age if estrogen==1 & prevdis==1, col(green)), ///
t1("fit8: Age Linear Splines Poisson") ///
legend(label(1 "No estrogen, no CVD") label(2 "No estrogen, prior CVD") label(4 "Estrogen, prior CVD") label(3 "Estrogen, no CVD"))


* PROBLEM #3: analyses based on RR
*
* OVERVIEW:
* Problem 2a: unadjusted analyses of cvddeath4 - estrogen association
*		chi-square test: cc
*		logistic regression
*
* Problem 2b: evidence of effect modification by prevdis
*		descriptive statistics
*		inferential statistics
*
* Problem 2c: evidence of confounding by prevdis
*		Association of estrogen - prevdis in sample
*		Association of prevdis - cvddeath4 after adjustment for estrogen
*
* Problem 2d: analyses of cvddeath4 - estrogen association adjusted for prevdis
*		Adjustment for main effect
*			logistic regression
*		Adjustment for main effect and estrogen - prevdis interaction
*			stratified analysis: cc
*			logistic regression
*
* Problem 2e: evidence of further confounding by age
*		Association of age - estrogen beyond adjustment for prevdis
*			based on regression
*			graphical assessment
*		Association of age - cvddeath4 after adjustment for estrogen, prevdis
*
* Problem 2f: association of cvddeath4 - estrogen adjusted for prevdis, age
*		Adjustment for main effect (unweighted vs weighted, classical vs robust)
*			Continuous linear
*			Quadratic
*			Splines
*

drop fit0 fit5 fit6 fit8 std

* EXECUTION OF CODE FOR PROBLEM 3
* Problem 3a: unadjusted analyses of cvddeath4 - estrogen association
*		chi-square test: cc
cc cvddeath4 estrogen
di  (3 / 337)   /    (88 / 2471)
*		logistic regression
logit cvddeath4 estrogen
logit cvddeath4 estrogen, robust
logistic cvddeath4 estrogen

* Problem 3b: evidence of effect modification by prevdis
*		descriptive statistics
cc cvddeath4 estrogen, by(prevdis)
*		inferential statistics
logistic cvddeath4 estrogen prevdis estr_prev
predict fit0
bysort estrogen prevdis: summ fit0 cvddeath4

* Problem 3c: evidence of confounding by prevdis
* See problem 1

* Problem 3d: analyses of cvddeath4 - estrogen association adjusted for prevdis
*		Adjustment for main effect
*			logistic regression
logistic cvddeath4 estrogen prevdis
*		Adjustment for main effect and estrogen - prevdis interaction
*			stratified analysis: cs
cc cvddeath4 estrogen, by(prevdis) istandard
cc cvddeath4 estrogen, by(prevdis) estandard
g std= 544 / 2899
replace std= 1 - std if prevdis==0
cc cvddeath4 estrogen, by(prevdis) standard(std)
*			logistic regression
logistic cvddeath4 estrogen prevdis estr_prev
test estrogen estr_prev
lincom (1 - 544/2899) * estrogen + (544/2899) * (estrogen + estr_prev), eform

* Problem 3e: evidence of further confounding by age
*		(see problem 1)

* Problem 3f: association of cvddeath4 - estrogen adjusted for prevdis, age
*		Adjustment for main effect (unweighted vs weighted, classical vs robust
*			Continuous linear
logistic cvddeath4 estrogen prevdis age
predict fit5
* (see do file for code used to produce graphs)
graph drop _all
twoway ///
(scatter fit5 age if estrogen==0 & prevdis==0, msymbol(i) lcol(black) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==0, col(black)), /// 
t1("fit5: estrogen 0 prevdis 0 (N= 2045)") name(f00)
twoway ///
(scatter fit5 age if estrogen==0 & prevdis==1, msymbol(i) lcol(blue) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==1, col(blue)), /// 
t1("fit5: estrogen 0 prevdis 1 (N= 514)") name(f01)
twoway ///
(scatter fit5 age if estrogen==1 & prevdis==1, msymbol(i) lcol(green) connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==1, col(green)), ///
 t1("fit5: estrogen 1 prevdis 1 (N= 30)") name(f11)
 twoway ///
(scatter fit5 age if estrogen==1 & prevdis==0, msymbol(i) lcol(red)  connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==0, col(red)), /// 
t1("fit5: estrogen 1 prevdis 0 (N= 310)") name(f10)
graph combine f10 f11 f00 f01, ycommon t1("fit5: Age Continuous Linear") name(panel)
twoway ///
(scatter fit5 age if estrogen==0 & prevdis==0, col(black)) ///
(scatter fit5 age if estrogen==0 & prevdis==1, col(blue)) ///
(scatter fit5 age if estrogen==1 & prevdis==0, col(red)) ///
(scatter fit5 age if estrogen==1 & prevdis==1, col(green)), ///
t1("fit5: Age Continuous Linear") ///
legend(label(1 "No estrogen, no CVD") label(2 "No estrogen, prior CVD") label(4 "Estrogen, prior CVD") label(3 "Estrogen, no CVD"))

*			Quadratic
logistic cvddeath4 estrogen prevdis age agesqr
predict fit6
* (see do file for code used to produce graphs)
graph drop _all
twoway ///
(scatter fit6 age if estrogen==0 & prevdis==0, msymbol(i) lcol(black) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==0, col(black)), /// 
t1("fit6: estrogen 0 prevdis 0 (N= 2045)") name(f00)
twoway ///
(scatter fit6 age if estrogen==0 & prevdis==1, msymbol(i) lcol(blue) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==1, col(blue)), /// 
t1("fit6: estrogen 0 prevdis 1 (N= 514)") name(f01)
twoway ///
(scatter fit6 age if estrogen==1 & prevdis==1, msymbol(i) lcol(green) connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==1, col(green)), ///
 t1("fit6: estrogen 1 prevdis 1 (N= 30)") name(f11)
 twoway ///
(scatter fit6 age if estrogen==1 & prevdis==0, msymbol(i) lcol(red)  connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==0, col(red)), /// 
t1("fit6: estrogen 1 prevdis 0 (N= 310)") name(f10)
graph combine f10 f11 f00 f01, ycommon t1("fit6: Age Quadratic") name(panel)
twoway ///
(scatter fit6 age if estrogen==0 & prevdis==0, col(black)) ///
(scatter fit6 age if estrogen==0 & prevdis==1, col(blue)) ///
(scatter fit6 age if estrogen==1 & prevdis==0, col(red)) ///
(scatter fit6 age if estrogen==1 & prevdis==1, col(green)), ///
t1("fit6: Age Quadratic") ///
legend(label(1 "No estrogen, no CVD") label(2 "No estrogen, prior CVD") label(4 "Estrogen, prior CVD") label(3 "Estrogen, no CVD"))

*			Splines
logistic cvddeath4 estrogen prevdis age65 age70 age75 age80 age85
predict fit8
* (see do file for code used to produce graphs)
graph drop _all
twoway ///
(scatter fit8 age if estrogen==0 & prevdis==0, msymbol(i) lcol(black) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==0, col(black)), /// 
t1("fit8: estrogen 0 prevdis 0 (N= 2045)") name(f00)
twoway ///
(scatter fit8 age if estrogen==0 & prevdis==1, msymbol(i) lcol(blue) connect(l) sort(age)) ///
(scatter p4 age if estrogen==0 & prevdis==1, col(blue)), /// 
t1("fit8: estrogen 0 prevdis 1 (N= 514)") name(f01)
twoway ///
(scatter fit8 age if estrogen==1 & prevdis==1, msymbol(i) lcol(green) connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==1, col(green)), ///
 t1("fit8: estrogen 1 prevdis 1 (N= 30)") name(f11)
 twoway ///
(scatter fit8 age if estrogen==1 & prevdis==0, msymbol(i) lcol(red)  connect(l) sort(age)) ///
(scatter p4 age if estrogen==1 & prevdis==0, col(red)), /// 
t1("fit8: estrogen 1 prevdis 0 (N= 310)") name(f10)
graph combine f10 f11 f00 f01, ycommon t1("fit8: Age Linear Splines") name(panel)
twoway ///
(scatter fit8 age if estrogen==0 & prevdis==0, col(black)) ///
(scatter fit8 age if estrogen==0 & prevdis==1, col(blue)) ///
(scatter fit8 age if estrogen==1 & prevdis==0, col(red)) ///
(scatter fit8 age if estrogen==1 & prevdis==1, col(green)), ///
t1("fit8: Age Linear Splines") ///
legend(label(1 "No estrogen, no CVD") label(2 "No estrogen, prior CVD") label(4 "Estrogen, prior CVD") label(3 "Estrogen, no CVD"))