***Epidemiology 265, Assignment 8
**Stephen Chang

clear

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

capture log close

set more off

log using assign8_stephenchang.smcl, replace

*Assignment8

/*
Modify your data generation code from last week to match the following causal 
DAG:

Where S1 and S2 are indicators for whether an observation is selected into the 
analysis (ie, when S1=0 or S2=0, the observations are excluded).  

Write out data generating rules to match this causal DAG for at least 2 values 
S1 (e.g., S1 = 1 for 90% of observations vs S1=1 for 50% of observations) and 2 
values for S2.  Please be sure that your data generating model entails that U 
influences X, Y, and S2. 
For each, write a simulation model generating data with N=500 individuals 
following these data generating rules.  Report the structural causal model for 
each variable in the DAG, including the distributions of the error terms.  Use
the statistical method of your choice to estimate the effect of X on Y in each
scenario, assuming that U is unmeasured (ie, your estimate will likely be 
biased because you cannot account for U). 

Now assume the null, that X has no effect on Y.  Under each of your four 
scenarios above, what is:
- The average estimate of the effect of X on Y (ie, if you repeat your 
chosen statistical analysis across repeated iterations, what is the average 
value of the parameter estimate for the association between X and Y)? 
- the type 1 error rate using an alpha threshold of 0.05 (ie, the percent of 
iterations under which you would find a statistically significant association 
between X and Y, using a p<.05 threshold)?*/

if 1 {
clear

set seed 98770
set obs 500
gen errorTerm = rnormal(0, 0.5)
gen U = runiform(0, 1.)
gen X = runiform(0, 1) + 0.2 * U
gen Y = 0 + 0.5*X + 1.*X*(U-0.5) + errorTerm

* generate data
foreach p1 in 0.2 1.1 ///
{
  foreach p2 in 1.0 0.1 ///
  {
  
    display "S1 param = `p2', S2 param = `p1'"

    capture drop S1 S2 Xsel Ysel
  
	gen S2 = runiform()>=`p1'*U
    gen S1 = runiform()<=`p2'*(X+Y)

	gen Xsel = X if S1 == 1 & S2 == 1
	gen Ysel = Y if S1 == 1 & S2 == 1
	
    * check & visualize data
    summarize
    * scatter Xsel Ysel
	
	* simple linear regression of Y on X
	regress Ysel Xsel
  }
}  
}

/* The average estimate of the effect of X on Y (ie, if you repeat your chosen 
   statistical analysis across repeated iterations, what is the average value 
   of the parameter estimate for the association between X and Y)? */

if 1 {
clear
set obs 50000

capture drop errorTerm U X Y
gen errorTerm = rnormal(0, 0.5)
gen U = runiform(0, 1.)
gen X = runiform(0, 1) + 0.2 * U
gen Y = 0 + 0.5*X + 1.*X*(U-0.5) + errorTerm   
  
* generate data
foreach p1 in 0.2 1.1 ///
{
  foreach p2 in 1.0 0.1 ///
  {
  
    display "S1 param = `p2', S2 param = `p1'"

    capture drop S1 S2 Xsel Ysel
  
	gen S2 = runiform()>=`p1'*U
    gen S1 = runiform()<=`p2'*(X+Y)

	gen Xsel = X if S1 == 1 & S2 == 1
	gen Ysel = Y if S1 == 1 & S2 == 1
	
	* simple linear regression of Y on X
	bootstrap, reps(100) size(500): regress Ysel Xsel
  }
}
}
  
/* What is the type 1 error rate using an alpha threshold of 0.05 (ie, the percent of 
   iterations under which you would find a statistically significant association 
   between X and Y, using a p<.05 threshold)? */
     
if 1 {	 
clear
set seed 98770
foreach p1 in 0.2 1.1 ///
{
  foreach p2 in 1.0 0.1 ///
  {

    display "S1 param = `p2', S2 param = `p1'"
  
	postfile results p using sim_`p1'_`p2', replace
	forvalues n=1/100 {
	   set obs 500
	   
	   capture drop errorTerm U X Y
       gen errorTerm = rnormal(0, 0.5)
       gen U = runiform(0, 1.)
       gen X = runiform(0, 1) + 0.2 * U
       gen Y = 0 + 0.5*X + 1.*X*(U-0.5) + errorTerm
  
       capture drop S1 S2 Xsel Ysel
	   gen S2 = runiform()>=`p1'*U
       gen S1 = runiform()<=`p2'*(X+Y)

	   gen Xsel = X if S1 == 1 & S2 == 1
	   gen Ysel = Y if S1 == 1 & S2 == 1
	
	   quietly regress Ysel Xsel
	   mat b = r(table)
	   local p = b[4,1]
	   post results (`p')
	}
	postclose results
  }
}

foreach p1 in 0.2 1.1 ///
{
  foreach p2 in 1.0 0.1 ///
  {
    clear
    use sim_`p1'_`p2'
	
	display "S1 param = `p2', S2 param = `p1'"
	
	* show total number of p values
	summarize p
	
	* show number of p values with p < 0.05
	summarize p if p < 0.05
  }
}
}