/*HW for Uncertainty Class:
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). 

Hypothesis: Each 1 SD increase in depressive symptom scores reduces cognitive score
by 0.3 SD units among individuals with above average education and 0.1 units among individuals with
average education or below. In my causal structure, education has no effect on cognition for people with 
depression=0

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 
xposure-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.*/


clear
set seed 89797
set obs 1000
gen dep=rnormal()
gen educ=runiform()<.5
* Cog_e represents all random variation in cognition and every other source of variation due to unknown factors
* Because I would like cognition to have a variance of 1, I've rescaled cog_e just a little bit so the variance
* of .2*depression (the marginal effect of depression on cognition if the effect for half the sample
* is .1 and for the other half of the sample is .3) plus the variance of cog_e sum to 1.
* This step of multiplying by (1-.2^2)^2 is unnecessary but makes the distributions nicer

gen cog_e=rnormal()*(1-.2^2)^.5
gen cog=-0.1*dep-0.2*dep*educ+cog_e

* If I wanted a binary variable, such as cognitive impairment, two simple options
*(but note they no longer correspond to the effect sizes I specified above):
gen cog_imp=cog<-1.5

gen lnodds_CI=0.1*dep+0.2*dep*educ-3
gen odds_CI=exp(lnodds_CI)
gen prob_CI=odds_CI/(1+odds_CI)
gen CI_e=runiform()  /*If I were really trying to be consistent w/ the definition of the continuous version, I'd use the cog_e above, but it has a weird distribution so I don't want to*/
gen CI=prob_CI>CI_e
sum

* Check the results in the full data
reg cog dep
reg cog_imp dep
reg CI dep
logit cog_imp dep
logit CI dep

* Made the data, now run the samples.  One way:
gen sample1=runiform()<.1
gen sample2=runiform()
gen sample3=runiform()
gen sample4=runiform()
gen sample5=runiform()

reg cog dep if sample1<.1
reg cog dep if sample2<.1
* This works more or less fine but the samples are not precisely n=100, so another option:
sort sample1
reg cog dep if _n<101
sort sample2
reg cog dep if _n<1001
*or do it the way stata would like (probably a good idea)
preserve
sample 100, count
reg cog dep
restore
preserve
sample 100 , count
reg cog dep
restore
preserve
sample 100 , count
reg cog dep
restore

* this also works fine but is boring, so another option is to
* bootstrap, making sure your computer knows where to save the results
cd "C:\Users\mglym\Dropbox\Maria\teaching\Epi265UCSF\epi265_2019\StataExamples"
bootstrap _b[dep] _se[dep] , r(10) si(100) saving(depcoefficients, replace): reg cog dep
clear
use depcoefficients
list

/*  Repeat the data set construction, setting the causal effect to the null. */ 


clear
set seed 8585
set obs 1000
gen dep=rnormal()
gen educ=runiform()<.5
gen cog_e=rnormal()
gen cog=0*dep-0*dep*educ+cog_e
reg cog dep
/*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).*/
bootstrap _b[dep] _se[dep] , r(1000) si(100) saving(depcoefficients, replace): reg cog dep
clear
use depcoefficients
* The bootstrap in stata returns the coefficient and standard error, and the ratio of these 
* follows a t-distribution with 99 degrees of freedom (n-# of coefficients)), the t is 
* pretty similar to the z. I'd like a 2-sided test. 
rename _bs_1 dep_beta
rename _bs_2 dep_se
gen dep_t_stat=dep_beta/dep_se
gen prob_tstat=t(99,dep_t_stat)
gen dep_p_val=2*(1-(t(99,abs(dep_t_stat))))
sum
hist dep_beta
hist dep_p_val

/* 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

clear
set seed 07097
* I could create a huge population and resample from it for n=100, but that is really slow
* in stata's bstrap approach, so I will create a do-file to generate the data and run that 
* many times.  This is what stata's "simulate" command does, though there are other options
* for how to run this.  
simulate beta se , r(5000): do gendatastep1
rename _sim_1 group_beta
rename _sim_2 group_se
* Not essential to this part of the problem, but notice that the SD of the beta estimates across iterations
* of the simulation equals the average of the standard errors across iterations of the simulationsum group_beta group_se
gen group_t_stat=group_beta/group_se
gen prob_tstat=t(99,group_t_stat)
gen group_p_val=2*(1-(t(99,abs(group_t_stat))))
gen stat_sig=group_p_val<.05
* the mean of "stat_sig" is the fraction of times the p-value was <.05, i.e., the power. 
sum stat_sig
power twomeans .02 .12 , n(100) a(.05)

*n=100,μ0=.02, μ1=.12, SD=2, α=.05
clear
set seed 07097
simulate beta se, r(1000): do gendatastep2
rename _sim_1 group_beta
rename _sim_2 group_se
sum group_beta group_se
gen group_t_stat=group_beta/group_se
gen prob_tstat=t(99,group_t_stat)
gen group_p_val=2*(1-(t(99,abs(group_t_stat))))
gen stat_sig=group_p_val<.05
sum stat_sig
power twomeans .02 .12 , n(100) a(.05) sd(2)

*n=500, μ0=.3, μ1=.3, SD=1, α=.05
clear
set seed 07097
simulate beta se, r(5000): do gendatastep3
rename _sim_1 group_beta
rename _sim_2 group_se
gen group_t_stat=group_beta/group_se
gen prob_tstat=t(99,group_t_stat)
gen group_p_val=2*(1-(t(99,abs(group_t_stat))))
gen stat_sig=group_p_val<.05
sum stat_sig
power twomeans .3 .300000001 , n(100) a(.05) sd(1)