* Workshop October 12, 2018 PhD Program in Public Health

* Part 1) Random error and systematic error
 
clear
set seed 1234
quietly set obs 100000

local rho = 0.7
scalar b_eo = 0
scalar b_co = 3

// Generating correlated regressors 
generate e = rnormal()
generate u = `rho'*e + rnormal()

// Generating Model

* Weak ED association vs strong UD and EU associations
* The strong crude association is explained by the confounder. 
scalar b_eo = -0.1
scalar b_co = 3

* Mechanism underlying the disease
quietly generate d = 70 + b_eo*e + b_co*u + rnormal()
reg d e , noheader  

*tw (function normalden(x, _b[e], _se[e]), lw(thick) lc(black) range(`=_b[e]-2.2*_se[e]'  `=_b[e]+2.2*_se[e]')) ///
*, ytitle("Sampling distribution") xtitle("Effect") plotregion(style(none))

reg d e , noheader  
scalar getbc = _b[e]
scalar getsebc = _se[e]
reg d e u , noheader  
scalar getba = _b[e]
scalar getseba = _se[e]
scalar b1c = _b[e]

/*
twoway ///
 (function normalden(x, getbc, getsebc),   lc(black) range(`=getbc-4.2*getsebc' `=getbc+4.2*getsebc') ) ///
 (function normalden(x, getba, getseba),  lc(black) range(`=getba-4.2*getseba' `=getba+4.2*getseba') ) ///
, ytitle("Sampling distribution") xtitle("Effect") plotregion(style(none)) ///
 title("Adjusted vs unadjusted")  legend(off) xlabel(-.1 0 2) xline(0, lp(dot)) 
*/


* What is the association between x and Mean(y) adjusted for z (b1)?

* Mean(y|x,z) = b0 + b1*x + b2*z

* We would like to isolate the effect of x on y from the effect of z on x and z on y. 
* This can be achieved by removing the effect of z on both x and y. 

* Mean(y|z) = a0 + a1*z
* Res(y|z) = y - (a0 + a1*z)

* Mean(x|z) = c0 + c1*z
* Res(x|z) = x - (c0 + c1*z)

* d1 represents the change in mean y for every 1 unit increase in x taking into 
* account the relationship between z and y and between z and x.
/*

Mean(Res(y|z)|Res(x|z)) = d0 + d1*Res(x|z)

Mean(  [y-(a0+a1*z)] | [x-(c0+c1*z)]  )  = d0 + d1*[x-(c0+c1*z)]

Mean(  [y-a1*z] | [x-c1*z]  )  = d0 + d1*[x-c1*z]

d0 = b0
d1 = b1
 
*/
reg d u 
scalar a1 = _b[u]
predict res_d , residual

reg e u
scalar c1 = _b[u]
predict res_e, residual

reg res_d res_e
reg d e u
scalar b1a = _b[e]

* Suppose you know the UD and EU confouding associations
* You need to know also the distribution of the confounder U
capture drop da ea
gen da = d-(a1*u)
gen ea = e-(c1*u)

reg da ea
reg d e u

di b1c
di a1
di c1 
di b1a

* Part 2) Traditional sensitivity analysis 

clear all
scalar a1 = 45
scalar a0 = 94
scalar b1 = 257
scalar b0 = 945

scalar pu1 = 0.4
scalar pu0 = 0.3

scalar or_ud = 5
scalar arr = 1.76
scalar B11 = pu1*b1
scalar B01 = pu0*b0

scalar A11 = or_ud*a1*B11/(or_ud*B11+b1-B11)
scalar A01 = or_ud*a0*B01/(or_ud*B01+b0-B01)

scalar e = a1-A11
scalar f = a0-A01
scalar g = b1-B11 
scalar h = b0-B01

scalar rrcd = 5

scalar bias = [(0.4*(5-1)+1)/(.3*(5-1)+1)]
display 1.76 / bias

				  
cci 45 94 257 945
cci `=round(scalar(A11))' `=round(scalar(A01))' `=round(scalar(B11))' `=round(scalar(B01))'

di (A11/B11)/(A01/B01)

clear 
input d e u count
1 1 1 35
1 0 1 64
0 1 1 103
0 0 1 284
1 1 0 10
1 0 0 30
0 1 0 154 
0 0 0 661
end

logistic d e   [freq=count]

cc d e [freq=count] , woolf

cc d u [freq=count] if e == 1
cc d u [freq=count] if e == 0


cc d e [freq=count], by(u) woolf
cc d u [freq=count], by(e) woolf

logistic d e  [freq=count]
logistic d e u [freq=count], cformat(%9.2f)

table d e u [freq=count]
 
episensi 45 94 257 945, dpexp(c(.4)) dpunexp(c(.3)) drrcd(c(5))
episensrri 1.76 , pexp(.4) punexp(.3) rrcd(5) 

scalar arr = (a1/b1)/(a0/b0)
di arr

scalar rr = arr / [ (pu1*(or_ud-1)+1 )/(pu0*(or_ud-1)+1) ]
di rr
 
exit