***Epidemiology 265, Assignment 7
**Stephen Chang

clear

cd "/Users/stephenchang/Desktop/TICR Program/Spring 2017/EPI 265 Epidemiologic Methods III Research Methods in Chronic Disease Epidemiology/Lecture 7/HW 7"

capture log close

set more off

log using assign7_stephenchang.smcl, replace

/* Assignment Week 7
Specify a hypothesis regarding a particular exposure and outcome and a 
binary effect modifier including specific measures of association 
(specify the magnitudes of that association you anticipate: I suggest making 
everything cross-sectional). 

Using the software of your choice, generate a population with 1000 people under 
a causal structure consistent with this hypothesis. Draw a simple random sample 
100 individuals from this population and estimate the population average 
exposure-outcome association and the association stratified by your modifier of 
interest within this subset.  Repeat this 10 times and write the parameter 
estimates and CI each time.

Repeat the data set construction, setting the causal effect to the null.  Again 
repeat this 10 times and write the parameter estimate and CI each time (if you 
figure out how to automate it, run it 1000 times and post the histogram of the 
parameter estimates and p-values).

Use your code above and also a canned software command to estimate statistical 
power to detect the difference in means under the settings below:
*n=100, μ0=.02, μ1=.12, SD=1, α=.05
*n=100,μ0=.02, μ1=.12, SD=2, α=.05
*n=500, μ0=.3, μ1=.3, SD=1, α=.05
For each of the 3 settings above, what is the power to detect whether the ratio 
of the means=1?
*/

clear

set obs 1000

set seed 98770

/* Model: BMD = beta0 + beta1A*CalciumIntake[Smoking] + beta1B*CalciumIntake[NoSmoking] + beta2*Smoking + ErrorTerm
          BMD and CalcumIntake are continuous,
		  Smoking is a dichotomous effect modifier 
		  
		  if smoking calcium intake less effecive in increasing BMD
		  => beta1A < beta1B
		  
		  pick: 
		  beta0 = 1.0 g/cm^2
		  beta1A = 0.01 g/cm^2 per 1000 mg Calcium / day
		  beta1B = 0.03 g/cm^2 per 1000 mg Calcium / day
		  beta2 = -0.1 g/cm^2  (smokers w/o Calcium intake should have somewhat smaller
		                        BMD values than non-smokers w/o Calcium intake)
		  
		  */

* ---- 1. Data Generation
gen Smoking = runiform()>=0.5  /* effect modifier; total sample has ~half people that are smoking */ 

gen CalciumIntake=rnormal(1000,300) /* Calcium Intake, mean 1000 mg/day, std-dev 300 mg/day */
replace CalciumIntake = 0 if CalciumIntake < 0 /* don't allow negative Intake */

gen CalciumIntakeSmoking = CalciumIntake if Smoking == 1
replace CalciumIntakeSmoking = 0 if Smoking == 0

gen CalciumIntakeNotSmoking = CalciumIntake if Smoking == 0
replace CalciumIntakeNotSmoking = 0 if Smoking == 1

gen ErrorTerm = rnormal(0, 0.1)  /* normally distributed error term */

gen BMD = 1.0 + 0.01 * CalciumIntakeSmoking + 0.03 * CalciumIntakeNotSmoking -0.1 * Smoking + ErrorTerm

* ---- 2. run regression on 10 random samples of 100 patients each

* the population average exposure-outcome association, i.e., don't account for eff. modifier
bootstrap, reps(10) size(100) noisily: regress BMD CalciumIntake

* the association stratified by your modifier of interest
bootstrap, reps(10) size(100) noisily: regress BMD Smoking##c.CalciumIntake

* ---- 3. Setting the causal effect to the null and repeat 1. & 2.
drop BMD
gen BMD = 1.0 + 0.02 * CalciumIntake + ErrorTerm

* the population average exposure-outcome association, i.e., don't account for eff. modifier
bootstrap, reps(10) size(100) noisily: regress BMD CalciumIntake

* run it 1000 times
bootstrap, reps(1000) size(100): regress BMD CalciumIntake

* ---- 4. estimate statistical power to detect the difference in means
power onemean 0.02 0.12, sd(1) n(100)
power onemean 0.02 0.12, sd(2) n(100)
power onemean 0.3 0.3, sd(1) n(500)