/*==============================================================================
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).

EXPOSURE: EDUCATION (COLLEGE OR MORE)
OUTCOME: FASTING GLUCOSE LEVEL
EFFECT MODIFIER: FOOD INSECURE YES/NO 
==============================================================================*/

clear all
set seed 518653
set obs 1000

gen edu = runiform()<0.6
gen food = runiform()<0.2

gen gluc = 120 - 20*edu + 40*food + 10*edu*food + rnormal() 


/*==============================================================================
DRAW A SIMPLE RANDOM SAMPLE OF 100 INDIVDIUALS FROM THE POPULATION AND ESTIMATE 
THE POPULATION AVERAGE EXPOSURE-OUTCOME ASSOCIATION & THE ASSOCIATON STRATIFIED 
BY YOUR MODIFIER OF INTEREST WITHIN THIS SUBSET. REPEAT THIS 10 TIMES AND WRITE 
THE PARAMETER ESTIMATES AND CI EACH TIME 
==============================================================================*/

*SIMULATION USING PROGRAM AND SIMULATE COMMAND
program define sample100
	preserve
	sample 100, count
	
	regress gluc edu
	lincom edu
	scalar coef1 = r(estimate)
	scalar pval1 = r(p)
	scalar lower1 = r(lb)
	scalar upper1 = r(ub)
	
	regress gluc edu if food==0
	lincom edu
	scalar coef2 = r(estimate)
	scalar pval2 = r(p)
	scalar lower2 = r(lb)
	scalar upper2 = r(ub)
	
	regress gluc edu if food==1
	lincom edu
	scalar coef3 = r(estimate)
	scalar pval3 = r(p)
	scalar lower3 = r(lb)
	scalar upper3 = r(ub)
	
	restore
end

set seed 6541328
simulate effect_all = coef1 pval_all = pval1 lower_all = lower1 upper_all = upper1 ///
		 effect_fd0 = coef2 pval_fd0 = pval2 lower_fd0 = lower2 upper_fd0 = upper2 ///
		 effect_fd1 = coef3 pval_fd1 = pval3 lower_fd1 = lower3 upper_fd1 = upper3, ///
		 reps(10): sample100

list effect_all pval_all lower_all upper_all //POPULATION AVERAGE
list effect_fd0 pval_fd0 lower_fd0 upper_fd0 //FOOD SECURE
list effect_fd1 pval_fd1 lower_fd1 upper_fd1 //FOOD INSECURE


/*==============================================================================
REPEAT THE DATA SET CONSTRUCTION, SETTING THE CAUSAL EFFECT TO THE NULL. AGAIN
REPEAT THIS 10 TIMES AND WRITE THE PARAMETER ESTIMATES 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 THE P-VALUES)
==============================================================================*/

clear all
set seed 6465165
set obs 1000

gen edu = runiform()<0.6
gen food = runiform()<0.2

gen gluc = 120 + 40*food + rnormal() 

program define sample100_2
	preserve
	sample 100, count
	
	regress gluc edu
	lincom edu
	scalar coef1_2 = r(estimate)
	scalar pval1_2 = r(p)
	scalar lower1_2 = r(lb)
	scalar upper1_2 = r(ub)
	
	regress gluc edu if food==0
	lincom edu
	scalar coef2_2 = r(estimate)
	scalar pval2_2 = r(p)
	scalar lower2_2 = r(lb)
	scalar upper2_2 = r(ub)
	
	regress gluc edu if food==1
	lincom edu
	scalar coef3_2 = r(estimate)
	scalar pval3_2 = r(p)
	scalar lower3_2 = r(lb)
	scalar upper3_2 = r(ub)
	
	restore
end

set seed 6546598
simulate effect_all_2 = coef1_2 pval_all_2 = pval1_2 lower_all_2 = lower1_2 upper_all_2 = upper1_2 ///
		 effect_fd0_2 = coef2_2 pval_fd0_2 = pval2_2 lower_fd0_2 = lower2_2 upper_fd0_2 = upper2_2 ///
		 effect_fd1_2 = coef3_2 pval_fd1_2 = pval3_2 lower_fd1_2 = lower3_2 upper_fd1_2 = upper3_2, ///
		 reps(1000): sample100_2

hist effect_all_2, name(hist_all)
hist pval_all_2, name(pval_all)
hist effect_fd0_2, name(hist_fd0)
hist pval_fd0_2, name(pval_fd0)
hist effect_fd1_2, name(hist_fd1)
hist pval_fd1_2, name(pval_fd1)


/*==============================================================================
USE YOUR CODE ABOVE AND ALSO A CONNED SOFTWARE COMMAND TO ESTIMATE STATISTICAL
POWER TO DETECT THE DIFFERENCE IN MEANS UNDER THE SETTINGS BELOW (SEE DATA GEN
RULES). FOR EACH OF THE THREE SETTINGS ABOVE, WHAT IS THE POWER TO DETECT
WHETHER THE RATIO OF THE MEANS = 1?
==============================================================================*/

power twomeans 0.02 0.12, sd(1) n(100) alpha(0.05)
power twomeans 0.02 0.12, sd(2) n(100) alpha(0.05)
power twomeans 0.30 0.30000001, sd(1) n(500) alpha(0.05)