********************************************************************************
* Data generation process
********************************************************************************

clear all
drop _all	
set more off
local seed = 1987
set seed `seed'

local total_runs = 1000
local total_obs = 4000

tempname dich_conf1
postfile `dich_conf1' double (n_x run exposure_n OLR_crude_OR OLR_adj_OR OLR_adj_lower95 OLR_adj_upper95 OLR_adj_coverage MWS_crude MWS_adj MWS_adj_lower95 MWS_adj_upper95 MWS_adj_coverage) using dich_conf1, replace
	
forvalues run = 1/`total_runs'{
	
	clear 
	set obs `total_obs'

		local n_x = 1
		
		gen x1 = .
		gen x2 = rbinomial(1, 0.5)
	/*	gen x3 = rbinomial(1, 0.5)
		gen x4 = rbinomial(1, 0.5)
		gen x5 = rbinomial(1, 0.5)
		gen x6 = rbinomial(1, 0.5)	*/
		
		local beta2 = log(1.5)
	/*	local beta3 = log(1.5)
		local beta4 = log(1.5)
		local beta5 = log(1.5)
		local beta6 = log(1.5)       */
	
		gen u = runiform()			/* which we use to introduce some error */
				
* replace the exposure based on the presence of caovariates
		replace x1 = uniform() < exp(log(0.5) + (`beta2' * x2) - (`beta2' * 0.5))
		
* every cut-off point is represented by its own latent variable (y*) (effect sizes met "local"_) Note that x1 is not repsonsible for the value of the latents what so ever!
		gen ystar1 = ln(u/(1-u)) + (`beta2' * x2) - (`beta2' * 0.5)   /* all betas are equal for all latent variables (proportionality)*/

		egen ystar1z = std(ystar1)				  /* here we standardize the whole latent vars. all have the same distribution */
		
		gen ystar1z_s = ystar1z + invnormal(1/7)  /* here we add the shift in the distributions of the latents for the different categories*/
		gen ystar2z_s = ystar1z + invnormal(2/7)  /* the shift is represented by (_s)*/
		gen ystar3z_s = ystar1z + invnormal(3/7)  /* we chose equal distributions here*/	
		gen ystar4z_s = ystar1z + invnormal(4/7)
		gen ystar5z_s = ystar1z + invnormal(5/7)
		gen ystar6z_s = ystar1z + invnormal(6/7)
				
		/*kdensity ystar1z_s
		kdensity ystar1z_s, addplot	(kdensity ystar2z_s ||	///
									 kdensity ystar3z_s || 	/// 
									 kdensity ystar4z_s || 	/// 
									 kdensity ystar5z_s || 	/// 
									 kdensity ystar6z_s || 	)*/
		
		gen y = cond(ystar1z_s > 0, 6, ///	
				cond(ystar2z_s > 0, 5, ///
				cond(ystar3z_s > 0, 4, ///
				cond(ystar4z_s > 0, 3, ///
				cond(ystar5z_s > 0, 2, ///
				cond(ystar6z_s > 0, 1, 0))))))
				
				tab y

* OLR
	* Crude
	ologit y x1, or
	local OLR_crude_OR = exp(_b[x1])
	
	* Adjusted
	ologit y x1 x2, or
	local OLR_adj_OR = exp(_b[x1])
	local OLR_adj_lower95 = exp(_b[x1] - (1.959964*_se[x1]))
	local OLR_adj_upper95 = exp(_b[x1] + (1.959964*_se[x1]))
	local OLR_adj_coverage = inrange(1,`OLR_adj_lower95',`OLR_adj_upper95')
	
* MWS
	* Crude
	somersd x1 y, tr(c)
	matrix MWM = e(b)
    local MWS_crude = MWM[1,1]
	
	* Stratified
	somersd x1 y, tr(c) wstrata(x2) cimatrix(test)
	matrix MWM = e(b)
    local MWS_adj = MWM[1,1]
	matrix CI = test
	local MWS_adj_lower95 = CI[1,2]
    local MWS_adj_upper95 = CI[1,3]
	local MWS_adj_coverage = inrange(0.5,`MWS_adj_lower95',`MWS_adj_upper95')
			
	* Propensity score	
	* Not applicable in this scenario
	
	tab x1 if x1 == 1
	matrix exposure = r(N)
	local exposure_n = exposure[1,1]
	
	post `dich_conf1' (`n_x') (`run') (`exposure_n') (`OLR_crude_OR') (`OLR_adj_OR') (`OLR_adj_lower95') (`OLR_adj_upper95') (`OLR_adj_coverage') (`MWS_crude') (`MWS_adj') (`MWS_adj_lower95') (`MWS_adj_upper95') (`MWS_adj_coverage')
      }
//

postclose `dich_conf1'
use dich_conf1, clear
gen MWS_OLR_crude_OR = (OLR_crude_OR/(OLR_crude_OR-1)^2)*((OLR_crude_OR-1)-ln(OLR_crude_OR))
gen MWS_OLR_adj_OR = (OLR_adj_OR/(OLR_adj_OR-1)^2)*((OLR_adj_OR-1)-ln(OLR_adj_OR))
save dich_conf1, replace

tempname dich_conf2
postfile `dich_conf2' double (n_x run exposure_n OLR_crude_OR OLR_adj_OR OLR_adj_lower95 OLR_adj_upper95 OLR_adj_coverage MWS_crude MWS_adj MWS_adj_lower95 MWS_adj_upper95 MWS_adj_coverage MWS_adj_pscore MWS_adj_pscore_lower95 MWS_adj_pscore_upper95 MWS_adj_pscore_coverage) using dich_conf2, replace
	
forvalues run = 1/`total_runs'{
	
	clear
	set obs `total_obs'
		
		local n_x = 2
		
		gen x1 = .
		gen x2 = rbinomial(1, 0.5)
		gen x3 = rbinomial(1, 0.5)
	/*	gen x4 = rbinomial(1, 0.5)
		gen x5 = rbinomial(1, 0.5)
		gen x6 = rbinomial(1, 0.5)	*/
		
		local beta2 = log(1.5)
		local beta3 = log(1.5)
	/*	local beta4 = log(1.5)
		local beta5 = log(1.5)
		local beta6 = log(1.5)   */
	
		gen u = runiform()			/* which we use to introduce some error */
				
* replace the exposure based on the presence of caovariates
		replace x1 = uniform() < exp(log(0.5) + (`beta2' * x2) - (`beta2' * 0.5) + (`beta3' * x3) - (`beta3' * 0.5))
		
* every cut-off point is represented by its own latent variable (y*) (effect sizes met "local"_) Note that x1 is not repsonsible for the value of the latents what so ever!
		gen ystar1 = ln(u/(1-u)) + (`beta2' * x2) - (`beta2' * 0.5) + (`beta3' * x3) - (`beta3' * 0.5)   /* all betas are equal for all latent variables (proportionality)*/

		egen ystar1z = std(ystar1)				  /* here we standardize the whole latent vars. all have the same distribution */
		
		gen ystar1z_s = ystar1z + invnormal(1/7)  /* here we add the shift in the distributions of the latents for the different categories*/
		gen ystar2z_s = ystar1z + invnormal(2/7)  /* the shift is represented by (_s)*/
		gen ystar3z_s = ystar1z + invnormal(3/7)  /* we chose equal distributions here*/	
		gen ystar4z_s = ystar1z + invnormal(4/7)
		gen ystar5z_s = ystar1z + invnormal(5/7)
		gen ystar6z_s = ystar1z + invnormal(6/7)
				
		/*kdensity ystar1z_s
		kdensity ystar1z_s, addplot	(kdensity ystar2z_s ||	///
									 kdensity ystar3z_s || 	/// 
									 kdensity ystar4z_s || 	/// 
									 kdensity ystar5z_s || 	/// 
									 kdensity ystar6z_s || 	)*/
		
		gen y = cond(ystar1z_s > 0, 6, ///	
				cond(ystar2z_s > 0, 5, ///
				cond(ystar3z_s > 0, 4, ///
				cond(ystar4z_s > 0, 3, ///
				cond(ystar5z_s > 0, 2, ///
				cond(ystar6z_s > 0, 1, 0))))))
				
				tab y
	
* OLR
	* Crude
	ologit y x1, or
	local OLR_crude_OR = exp(_b[x1])
	
	* Adjusted
	ologit y x1 x2 x3, or
	local OLR_adj_OR = exp(_b[x1])
	local OLR_adj_lower95 = exp(_b[x1] - (1.959964*_se[x1]))
	local OLR_adj_upper95 = exp(_b[x1] + (1.959964*_se[x1]))
	local OLR_adj_coverage = inrange(1,`OLR_adj_lower95',`OLR_adj_upper95')
	
* MWS
	* Crude
	somersd x1 y, tr(c)
	matrix MWM = e(b)
    local MWS_crude = MWM[1,1]
	
	* Stratified
	somersd x1 y, tr(c) wstrata(x2 x3) cimatrix(test)
	matrix MWM = e(b)
    local MWS_adj = MWM[1,1]
	matrix CI = test
	local MWS_adj_lower95 = CI[1,2]
    local MWS_adj_upper95 = CI[1,3]
	local MWS_adj_coverage = inrange(0.5,`MWS_adj_lower95',`MWS_adj_upper95')
		
	* Propensity score
	logit x1 x2 x3
	predict pscore 	
	sum pscore	
	
	xtile pscore_quartile = pscore, nq(4)
			
	somersd x1 y, tr(c) wstrata(pscore_quartile) cimatrix(test)
	matrix MWM = e(b)
    local MWS_adj_pscore = MWM[1,1]
	matrix CI = test
	local MWS_adj_pscore_lower95 = CI[1,2]
	local MWS_adj_pscore_upper95 = CI[1,3]
	local MWS_adj_pscore_coverage = inrange(0.5,`MWS_adj_pscore_lower95',`MWS_adj_pscore_upper95')
		
	tab x1 if x1 == 1
	matrix exposure = r(N)
	local exposure_n = exposure[1,1]
	
	post `dich_conf2' (`n_x') (`run') (`exposure_n') (`OLR_crude_OR') (`OLR_adj_OR') (`OLR_adj_lower95') (`OLR_adj_upper95') (`OLR_adj_coverage') (`MWS_crude') (`MWS_adj') (`MWS_adj_lower95') (`MWS_adj_upper95') (`MWS_adj_coverage') (`MWS_adj_pscore') (`MWS_adj_pscore_lower95') (`MWS_adj_pscore_upper95') (`MWS_adj_pscore_coverage')
      }
//

postclose `dich_conf2'
use dich_conf2, clear
gen MWS_OLR_crude_OR = (OLR_crude_OR/(OLR_crude_OR-1)^2)*((OLR_crude_OR-1)-ln(OLR_crude_OR))
gen MWS_OLR_adj_OR = (OLR_adj_OR/(OLR_adj_OR-1)^2)*((OLR_adj_OR-1)-ln(OLR_adj_OR))
save dich_conf2, replace

tempname dich_conf3
postfile `dich_conf3' double (n_x run exposure_n OLR_crude_OR OLR_adj_OR OLR_adj_lower95 OLR_adj_upper95 OLR_adj_coverage MWS_crude MWS_adj MWS_adj_lower95 MWS_adj_upper95 MWS_adj_coverage MWS_adj_pscore MWS_adj_pscore_lower95 MWS_adj_pscore_upper95 MWS_adj_pscore_coverage) using dich_conf3, replace
	
forvalues run = 1/`total_runs'{
	
	clear 
	set obs `total_obs'
		
		local n_x = 3
		
		gen x1 = .
		gen x2 = rbinomial(1, 0.5)
		gen x3 = rbinomial(1, 0.5)
		gen x4 = rbinomial(1, 0.5)
	/*	gen x5 = rbinomial(1, 0.5)
		gen x6 = rbinomial(1, 0.5)	*/
		
		local beta2 = log(1.5)
		local beta3 = log(1.5)
		local beta4 = log(1.5)
	/*	local beta5 = log(1.5)
		local beta6 = log(1.5)      */
	
		gen u = runiform()			/* which we use to introduce some error */
				
* replace the exposure based on the presence of caovariates
		replace x1 = uniform() < exp(log(0.5) + (`beta2' * x2) - (`beta2' * 0.5) + (`beta3' * x3) - (`beta3' * 0.5) + (`beta4' * x4) - (`beta4' * 0.5))
		
* every cut-off point is represented by its own latent variable (y*) (effect sizes met "local"_) Note that x1 is not repsonsible for the value of the latents what so ever!
		gen ystar1 = ln(u/(1-u)) + (`beta2' * x2) - (`beta2' * 0.5) + (`beta3' * x3) - (`beta3' * 0.5) + (`beta4' * x4) - (`beta4' * 0.5)  /* all betas are equal for all latent variables (proportionality)*/

		egen ystar1z = std(ystar1)				  /* here we standardize the whole latent vars. all have the same distribution */
		
		gen ystar1z_s = ystar1z + invnormal(1/7)  /* here we add the shift in the distributions of the latents for the different categories*/
		gen ystar2z_s = ystar1z + invnormal(2/7)  /* the shift is represented by (_s)*/
		gen ystar3z_s = ystar1z + invnormal(3/7)  /* we chose equal distributions here*/	
		gen ystar4z_s = ystar1z + invnormal(4/7)
		gen ystar5z_s = ystar1z + invnormal(5/7)
		gen ystar6z_s = ystar1z + invnormal(6/7)
				
		/*kdensity ystar1z_s
		kdensity ystar1z_s, addplot	(kdensity ystar2z_s ||	///
									 kdensity ystar3z_s || 	/// 
									 kdensity ystar4z_s || 	/// 
									 kdensity ystar5z_s || 	/// 
									 kdensity ystar6z_s || 	)*/
		
		gen y = cond(ystar1z_s > 0, 6, ///	
				cond(ystar2z_s > 0, 5, ///
				cond(ystar3z_s > 0, 4, ///
				cond(ystar4z_s > 0, 3, ///
				cond(ystar5z_s > 0, 2, ///
				cond(ystar6z_s > 0, 1, 0))))))
				
				tab y
	
* OLR
	* Crude
	ologit y x1, or	
	matrix coefficient = e(b)
	local OLR_crude_OR = exp(coefficient[1,1])
	
	* Adjusted
	ologit y x1 x2 x3 x4, or
	local OLR_adj_OR = exp(_b[x1])	
	local OLR_adj_lower95 = exp(_b[x1] - (1.959964*_se[x1]))
	local OLR_adj_upper95 = exp(_b[x1] + (1.959964*_se[x1]))
	local OLR_adj_coverage = inrange(1,`OLR_adj_lower95',`OLR_adj_upper95')
	
* MWS
	* Crude
	somersd x1 y, tr(c)
	matrix MWM = e(b)
    local MWS_crude = MWM[1,1]
	
	* Stratified
	somersd x1 y, tr(c) wstrata(x2 x3 x4) cimatrix(test)
	matrix MWM = e(b)
    local MWS_adj = MWM[1,1]
	matrix CI = test
	local MWS_adj_lower95 = CI[1,2]
    local MWS_adj_upper95 = CI[1,3]
	local MWS_adj_coverage = inrange(0.5,`MWS_adj_lower95',`MWS_adj_upper95')
		
	* Propensity score
	logit x1 x2 x3 x4
	predict pscore 	
	sum pscore	
	
	xtile pscore_quartile = pscore, nq(4)
			
	somersd x1 y, tr(c) wstrata(pscore_quartile) cimatrix(test)
	matrix MWM = e(b)
    local MWS_adj_pscore = MWM[1,1]
	matrix CI = test
	local MWS_adj_pscore_lower95 = CI[1,2]
	local MWS_adj_pscore_upper95 = CI[1,3]
	local MWS_adj_pscore_coverage = inrange(0.5,`MWS_adj_pscore_lower95',`MWS_adj_pscore_upper95')
	
	tab x1 if x1 == 1
	matrix exposure = r(N)
	local exposure_n = exposure[1,1]
	
	post `dich_conf3' (`n_x') (`run') (`exposure_n') (`OLR_crude_OR') (`OLR_adj_OR') (`OLR_adj_lower95') (`OLR_adj_upper95') (`OLR_adj_coverage') (`MWS_crude') (`MWS_adj') (`MWS_adj_lower95') (`MWS_adj_upper95') (`MWS_adj_coverage') (`MWS_adj_pscore') (`MWS_adj_pscore_lower95') (`MWS_adj_pscore_upper95') (`MWS_adj_pscore_coverage')
      }
//

postclose `dich_conf3'
use dich_conf3, clear
gen MWS_OLR_crude_OR = (OLR_crude_OR/(OLR_crude_OR-1)^2)*((OLR_crude_OR-1)-ln(OLR_crude_OR))
gen MWS_OLR_adj_OR = (OLR_adj_OR/(OLR_adj_OR-1)^2)*((OLR_adj_OR-1)-ln(OLR_adj_OR))
save dich_conf3, replace

tempname dich_conf4
postfile `dich_conf4' double (n_x run exposure_n OLR_crude_OR OLR_adj_OR OLR_adj_lower95 OLR_adj_upper95 OLR_adj_coverage MWS_crude MWS_adj MWS_adj_lower95 MWS_adj_upper95 MWS_adj_coverage MWS_adj_pscore MWS_adj_pscore_lower95 MWS_adj_pscore_upper95 MWS_adj_pscore_coverage) using dich_conf4, replace
	
forvalues run = 1/`total_runs'{
	
	clear 
	set obs `total_obs'
		
		local n_x = 4
		
		gen x1 = .
		gen x2 = rbinomial(1, 0.5)
		gen x3 = rbinomial(1, 0.5)
		gen x4 = rbinomial(1, 0.5)
		gen x5 = rbinomial(1, 0.5)
	/*	gen x6 = rbinomial(1, 0.5)	*/
		
		local beta2 = log(1.5)
		local beta3 = log(1.5)
		local beta4 = log(1.5)
		local beta5 = log(1.5)
	/*	local beta6 = log(1.5)      */
	
		gen u = runiform()			/* which we use to introduce some error */
				
* replace the exposure based on the presence of caovariates
		replace x1 = uniform() < exp(log(0.5) + (`beta2' * x2) - (`beta2' * 0.5) + (`beta3' * x3) - (`beta3' * 0.5) + (`beta4' * x4) - (`beta4' * 0.5) + (`beta5' * x5) - (`beta5' * 0.5))
		
* every cut-off point is represented by its own latent variable (y*) (effect sizes met "local"_) Note that x1 is not repsonsible for the value of the latents what so ever!
		gen ystar1 = ln(u/(1-u)) + (`beta2' * x2) - (`beta2' * 0.5) + (`beta3' * x3) - (`beta3' * 0.5) + (`beta4' * x4) - (`beta4' * 0.5) + (`beta5' * x5) - (`beta5' * 0.5)  /* all betas are equal for all latent variables (proportionality)*/

		egen ystar1z = std(ystar1)				  /* here we standardize the whole latent vars. all have the same distribution */
		
		gen ystar1z_s = ystar1z + invnormal(1/7)  /* here we add the shift in the distributions of the latents for the different categories*/
		gen ystar2z_s = ystar1z + invnormal(2/7)  /* the shift is represented by (_s)*/
		gen ystar3z_s = ystar1z + invnormal(3/7)  /* we chose equal distributions here*/	
		gen ystar4z_s = ystar1z + invnormal(4/7)
		gen ystar5z_s = ystar1z + invnormal(5/7)
		gen ystar6z_s = ystar1z + invnormal(6/7)
				
		/*kdensity ystar1z_s
		kdensity ystar1z_s, addplot	(kdensity ystar2z_s ||	///
									 kdensity ystar3z_s || 	/// 
									 kdensity ystar4z_s || 	/// 
									 kdensity ystar5z_s || 	/// 
									 kdensity ystar6z_s || 	)*/
		
		gen y = cond(ystar1z_s > 0, 6, ///	
				cond(ystar2z_s > 0, 5, ///
				cond(ystar3z_s > 0, 4, ///
				cond(ystar4z_s > 0, 3, ///
				cond(ystar5z_s > 0, 2, ///
				cond(ystar6z_s > 0, 1, 0))))))
				
				tab y
	
* OLR
	* Crude
	ologit y x1, or	
	local OLR_crude_OR = exp(_b[x1])
	
	* Adjusted
	ologit y x1 x2 x3 x4 x5, or
	local OLR_adj_OR = exp(_b[x1])	
	local OLR_adj_lower95 = exp(_b[x1] - (1.959964*_se[x1]))
	local OLR_adj_upper95 = exp(_b[x1] + (1.959964*_se[x1]))
	local OLR_adj_coverage = inrange(1,`OLR_adj_lower95',`OLR_adj_upper95')
	
* MWS
	* Crude
	somersd x1 y, tr(c)
	matrix MWM = e(b)
    local MWS_crude = MWM[1,1]
	
	* Stratified
	somersd x1 y, tr(c) wstrata(x2 x3 x4 x5) cimatrix(test)
	matrix MWM = e(b)
    local MWS_adj = MWM[1,1]
	matrix CI = test
	local MWS_adj_lower95 = CI[1,2]
    local MWS_adj_upper95 = CI[1,3]
	local MWS_adj_coverage = inrange(0.5,`MWS_adj_lower95',`MWS_adj_upper95')
		
	* Propensity score
	logit x1 x2 x3 x4 x5
	predict pscore 	
	sum pscore	
	
	xtile pscore_quartile = pscore, nq(4)
			
	somersd x1 y, tr(c) wstrata(pscore_quartile) cimatrix(test)
	matrix MWM = e(b)
    local MWS_adj_pscore = MWM[1,1]
	matrix CI = test
	local MWS_adj_pscore_lower95 = CI[1,2]
	local MWS_adj_pscore_upper95 = CI[1,3]
	local MWS_adj_pscore_coverage = inrange(0.5,`MWS_adj_pscore_lower95',`MWS_adj_pscore_upper95')
	
	tab x1 if x1 == 1
	matrix exposure = r(N)
	local exposure_n = exposure[1,1]
	
	post `dich_conf4' (`n_x') (`run') (`exposure_n') (`OLR_crude_OR') (`OLR_adj_OR') (`OLR_adj_lower95') (`OLR_adj_upper95') (`OLR_adj_coverage') (`MWS_crude') (`MWS_adj') (`MWS_adj_lower95') (`MWS_adj_upper95') (`MWS_adj_coverage') (`MWS_adj_pscore') (`MWS_adj_pscore_lower95') (`MWS_adj_pscore_upper95') (`MWS_adj_pscore_coverage')
      }
//

postclose `dich_conf4'
use dich_conf4, clear
gen MWS_OLR_crude_OR = (OLR_crude_OR/(OLR_crude_OR-1)^2)*((OLR_crude_OR-1)-ln(OLR_crude_OR))
gen MWS_OLR_adj_OR = (OLR_adj_OR/(OLR_adj_OR-1)^2)*((OLR_adj_OR-1)-ln(OLR_adj_OR))
save dich_conf4, replace
	
tempname dich_conf5
postfile `dich_conf5' double (n_x run exposure_n OLR_crude_OR OLR_adj_OR OLR_adj_lower95 OLR_adj_upper95 OLR_adj_coverage MWS_crude MWS_adj MWS_adj_lower95 MWS_adj_upper95 MWS_adj_coverage MWS_adj_pscore MWS_adj_pscore_lower95 MWS_adj_pscore_upper95 MWS_adj_pscore_coverage) using dich_conf5, replace
	
forvalues run = 1/`total_runs'{
	
	clear 
	set obs `total_obs'
		
		local n_x = 5
		
		gen x1  = .
		gen x2  = rbinomial(1, 0.5)
		gen x3  = rbinomial(1, 0.5)
		gen x4  = rbinomial(1, 0.5)
		gen x5  = rbinomial(1, 0.5)
		gen x6  = rbinomial(1, 0.5)
		
		local beta2  = log(1.5)
		local beta3  = log(1.5)
		local beta4  = log(1.5)
		local beta5  = log(1.5)
		local beta6  = log(1.5)
		
		gen u = runiform()			/* which we use to introduce some error */
				
* replace the exposure based on the presence of caovariates
		replace x1 = uniform() < exp(log(0.5) + (`beta2' * x2) - (`beta2' * 0.5) + (`beta3' * x3) - (`beta3' * 0.5) + (`beta4' * x4) - (`beta4' * 0.5) + (`beta5' * x5) - (`beta5' * 0.5) + (`beta6' * x6) - (`beta6' * 0.5))
		
* every cut-off point is represented by its own latent variable (y*) (effect sizes met "local"_) Note that x1 is not repsonsible for the value of the latents what so ever!
		gen ystar1 = ln(u/(1-u)) + (`beta2' * x2) - (`beta2' * 0.5) + (`beta3' * x3) - (`beta3' * 0.5) + (`beta4' * x4) - (`beta4' * 0.5) + (`beta5' * x5) - (`beta5' * 0.5) + (`beta6' * x6) - (`beta6' * 0.5)   /* all betas are equal for all latent variables (proportionality)*/

		egen ystar1z = std(ystar1)				  /* here we standardize the whole latent vars. all have the same distribution */
		
		gen ystar1z_s = ystar1z + invnormal(1/7)  /* here we add the shift in the distributions of the latents for the different categories*/
		gen ystar2z_s = ystar1z + invnormal(2/7)  /* the shift is represented by (_s)*/
		gen ystar3z_s = ystar1z + invnormal(3/7)  /* we chose equal distributions here*/	
		gen ystar4z_s = ystar1z + invnormal(4/7)
		gen ystar5z_s = ystar1z + invnormal(5/7)
		gen ystar6z_s = ystar1z + invnormal(6/7)
				
		/*kdensity ystar1z_s
		kdensity ystar1z_s, addplot	(kdensity ystar2z_s ||	///
									 kdensity ystar3z_s || 	/// 
									 kdensity ystar4z_s || 	/// 
									 kdensity ystar5z_s || 	/// 
									 kdensity ystar6z_s || 	)*/
		
		gen y = cond(ystar1z_s > 0, 6, ///	
				cond(ystar2z_s > 0, 5, ///
				cond(ystar3z_s > 0, 4, ///
				cond(ystar4z_s > 0, 3, ///
				cond(ystar5z_s > 0, 2, ///
				cond(ystar6z_s > 0, 1, 0))))))
				
				tab y
	
* OLR
	* Crude
	ologit y x1, or	
	local OLR_crude_OR = exp(_b[x1])
	
	* Adjusted
	ologit y x1 x2 x3 x4 x5 x6, or
	local OLR_adj_OR = exp(_b[x1])	
	local OLR_adj_lower95 = exp(_b[x1] - (1.959964*_se[x1]))
	local OLR_adj_upper95 = exp(_b[x1] + (1.959964*_se[x1]))
	local OLR_adj_coverage = inrange(1,`OLR_adj_lower95',`OLR_adj_upper95')
	
* MWS
	* Crude
	somersd x1 y, tr(c)
	matrix MWM = e(b)
    local MWS_crude = MWM[1,1]
	
	* Stratified
	somersd x1 y, tr(c) wstrata(x2 x3 x4 x5 x6) cimatrix(test)
	matrix MWM = e(b)
    local MWS_adj = MWM[1,1]
	matrix CI = test
	local MWS_adj_lower95 = CI[1,2]
    local MWS_adj_upper95 = CI[1,3]
	local MWS_adj_coverage = inrange(0.5,`MWS_adj_lower95',`MWS_adj_upper95')
			
	* Propensity score
	logit x1 x2 x3 x4 x5 x6
	predict pscore 	
	sum pscore	
	
	xtile pscore_quartile = pscore, nq(4)
		
	somersd x1 y, tr(c) wstrata(pscore_quartile) cimatrix(test)
	matrix MWM = e(b)
    local MWS_adj_pscore = MWM[1,1]
	matrix CI = test
	local MWS_adj_pscore_lower95 = CI[1,2]
	local MWS_adj_pscore_upper95 = CI[1,3]
	local MWS_adj_pscore_coverage = inrange(0.5,`MWS_adj_pscore_lower95',`MWS_adj_pscore_upper95')
	
	tab x1 if x1 == 1
	matrix exposure = r(N)
	local exposure_n = exposure[1,1]
	
	post `dich_conf5' (`n_x') (`run') (`exposure_n') (`OLR_crude_OR') (`OLR_adj_OR') (`OLR_adj_lower95') (`OLR_adj_upper95') (`OLR_adj_coverage') (`MWS_crude') (`MWS_adj') (`MWS_adj_lower95') (`MWS_adj_upper95') (`MWS_adj_coverage') (`MWS_adj_pscore') (`MWS_adj_pscore_lower95') (`MWS_adj_pscore_upper95') (`MWS_adj_pscore_coverage')
      }
//

postclose `dich_conf5'
use dich_conf5, clear
gen MWS_OLR_crude_OR = (OLR_crude_OR/(OLR_crude_OR-1)^2)*((OLR_crude_OR-1)-ln(OLR_crude_OR))
gen MWS_OLR_adj_OR = (OLR_adj_OR/(OLR_adj_OR-1)^2)*((OLR_adj_OR-1)-ln(OLR_adj_OR))
save dich_conf5, replace	

* Combine all datasets
use dich_conf1, clear
append using dich_conf2
append using dich_conf3
append using dich_conf4
append using dich_conf5

* Check number of exposed subjects
bysort n_x: sum exposure_n

* Replace invalid MWS results
replace MWS_adj_lower95 = .a if MWS_adj_lower95 == MWS_adj_upper95
replace MWS_adj_upper95 = .a if MWS_adj_lower95 == .a

replace MWS_adj_upper95 = .b if MWS_adj_lower95 <= 0 & !missing(MWS_adj_lower95)
replace MWS_adj_lower95 = .b if MWS_adj_lower95 <= 0 & !missing(MWS_adj_lower95)

replace MWS_adj_lower95 = .c if MWS_adj_upper95 >= 1 & !missing(MWS_adj_upper95) 
replace MWS_adj_upper95 = .c if MWS_adj_upper95 >= 1 & !missing(MWS_adj_upper95) 

replace MWS_adj_coverage = .a if MWS_adj_upper95 == .a
replace MWS_adj_coverage = .b if MWS_adj_upper95 == .b
replace MWS_adj_coverage = .c if MWS_adj_upper95 == .c

replace MWS_adj = .c if MWS_adj == 0
replace MWS_adj = .c if MWS_adj == 1
replace MWS_adj = .c if MWS_adj == 1/2
replace MWS_adj = .c if MWS_adj == 1/3
replace MWS_adj = .c if MWS_adj == 2/3
replace MWS_adj = .c if MWS_adj > 0.33333333 & MWS_adj < 0.33333334
replace MWS_adj = .c if MWS_adj > 0.66666666 & MWS_adj < 0.66666667
replace MWS_adj = .c if MWS_adj == 1/4
replace MWS_adj = .c if MWS_adj == 3/4
replace MWS_adj = .c if MWS_adj == 1/5
replace MWS_adj = .c if MWS_adj == 2/5
replace MWS_adj = .c if MWS_adj == 3/5
replace MWS_adj = .c if MWS_adj == 4/5
replace MWS_adj = .c if MWS_adj == 1/6
replace MWS_adj = .c if MWS_adj == 5/6
replace MWS_adj = .c if MWS_adj == 1/7
replace MWS_adj = .c if MWS_adj == 2/7
replace MWS_adj = .c if MWS_adj == 3/7
replace MWS_adj = .c if MWS_adj == 4/7
replace MWS_adj = .c if MWS_adj == 5/7
replace MWS_adj = .c if MWS_adj == 6/7
replace MWS_adj = .c if MWS_adj == 1/8
replace MWS_adj = .c if MWS_adj == 3/8
replace MWS_adj = .c if MWS_adj == 5/8
replace MWS_adj = .c if MWS_adj == 7/8
replace MWS_adj = .c if MWS_adj == 1/9
replace MWS_adj = .c if MWS_adj == 2/9
replace MWS_adj = .c if MWS_adj == 4/9
replace MWS_adj = .c if MWS_adj == 5/9
replace MWS_adj = .c if MWS_adj == 7/9
replace MWS_adj = .c if MWS_adj == 8/9
replace MWS_adj = .c if MWS_adj == 1/10
replace MWS_adj = .c if MWS_adj == 3/10
replace MWS_adj = .c if MWS_adj == 7/10
replace MWS_adj = .c if MWS_adj == 9/10

* Calculate our definition of bias out of MWS results
gen bias_OLR   = MWS_OLR_crude_OR - 0.5
gen bias_aOLR  = MWS_OLR_adj_OR - 0.5
gen bias_MWS   = MWS_crude - 0.5
gen bias_sMWS  = MWS_adj - 0.5
gen bias_psMWS = MWS_adj_pscore - 0.5

* Calculate median bias	
by n_x: egen p50_aOLR  = pctile(bias_aOLR), p(50)
by n_x: egen p50_sMWS  = pctile(bias_sMWS), p(50)
by n_x: egen p50_psMWS = pctile(bias_psMWS), p(50)

********************************************************************************
* Other outcome measures			
********************************************************************************
* How many sMWS are validly performed (= non-missing)?
bysort n_x: egen sMWS_valid=count(MWS_adj)

* How many sMWS CI's are generated to calculate the coverage probability?
bysort n_x: egen sMWS_CI_ll_valid=count(MWS_adj_lower95)
bysort n_x: egen sMWS_CI_ul_valid=count(MWS_adj_upper95)

* How many times bias is more than 0,05 (MWS) off?
gen aOLR_5 =  cond(MWS_OLR_adj_OR>0.55, 1, ///	
		      cond(MWS_OLR_adj_OR<0.45, 1, 0))

gen sMWS_5 =  cond(MWS_adj>0.55, 1, ///	
			  cond(MWS_adj<0.45, 1, 0))

replace sMWS_5 = .c if MWS_adj == .c				
				
gen psMWS_5 = cond(MWS_adj_pscore>0.55, 1, ///	
			  cond(MWS_adj_pscore<0.45, 1, 0))

by n_x: egen  aOLR_5_prob = mean(aOLR_5)
by n_x: egen  sMWS_5_prob = mean(sMWS_5)
by n_x: egen psMWS_5_prob = mean(psMWS_5)

* Generate the precision
gen  aOLR_precision = (1 - aOLR_5_prob)
gen  sMWS_precision = (1 - sMWS_5_prob) * (sMWS_valid / 1000)
gen psMWS_precision = (1 - psMWS_5_prob)

* Coverage probability
by n_x: egen  aOLR_coverage_prob = mean(OLR_adj_coverage)
by n_x: egen  sMWS_coverage_prob = mean(MWS_adj_coverage)
by n_x: egen psMWS_coverage_prob = mean(MWS_adj_pscore_coverage)

********************************************************************************
* Graph		
********************************************************************************
label variable bias_OLR "OLR"
label variable bias_aOLR "aOLR"
label variable bias_MWS "MWS"
label variable bias_sMWS "sMWS"
label variable bias_psMWS "psMWS"

set scheme s1color

graph box 	bias_MWS bias_OLR bias_sMWS bias_aOLR bias_psMWS ///
			, over(n_x) /// 
			yscale(range(-.10 .2)) ///
			ylabel(-.10(.05).2) ///
			nooutsides ///
			yline(0, lwidth(thin) lpattern(solid) lcolor(black)) ///
			box(1, color(dkgreen%85*.3)) ///
			box(2, color(dknavy%85*.3)) ///
			box(3, color(dkgreen%85*.7)) ///
			box(4, color(dknavy%85*.7)) ///
			box(5, color(dkgreen%85*1)) ///
			ytitle(Bias) ///
			title(Dichotomous confounders) ///	
			legend(style(column) size(vsmall) symy(*0.6) symx(*0.4) pos(3))