*! v.1.0.0 N.Orsini 9aug2007 capture program drop episens_mcsa_all program episens_mcsa_all, rclass version 9.2 syntax [anything] [ , obs(integer 5000) apprr(string) applb(string) appub(string) studytype(string) seed(string) GRaph DStat /// sseca(string) sspca(string) ssenc(string) sspnc(string) scorrsens(string) scorrspec(string) /// spscex(string) spscun(string) spsnex(string) spsnun(string) ssbfactor(string) /// spexp(string) spunexp(string) srrcd(string) sorce(string) scorrprev(string) /// ] drop _all qui set obs `obs' if "`seed'" != "" set seed `seed' * get the cell counts gettoken a 0 : 0, parse(" ,") gettoken b 0 : 0, parse(" ,") gettoken c 0 : 0, parse(" ,") gettoken d 0 : 0, parse(" ,") * di "`a' `b' `c' `d'" local m1 = `a' + `b' local m0 = `c' + `d' * get the observed RR tempname arrdx albdx aubdx scalar `arrdx' = `apprr' scalar `albdx' = `applb' scalar `aubdx' = `appub' if ("`sseca'"!="") & ("`sspca'"!="") & ("`ssenc'"!="") & ("`sspnc'"!="") { * Generate a pairs of 0,1 random variables with user-specified correlation tempvar v1 v2 v3 v4 v5 v6 se1 se0 sp1 sp0 tempname rsens rspec forv i = 1/6 { gen `v`i'' = logit(uniform()) } if "`scorrsens'" != "" scalar `rsens' = `scorrsens' else scalar `rsens' = 1 if "`scorrspec'" != "" scalar `rspec' = `scorrspec' else scalar `rspec' = 1 gen `se1' = invlogit( sqrt(`rsens')*`v1'+ sqrt(1-`rsens')*`v2') gen `se0' = invlogit( sqrt(`rsens')*`v1'+ sqrt(1-`rsens')*`v3') gen `sp1' = invlogit( sqrt(`rspec')*`v4'+ sqrt(1-`rspec')*`v5') gen `sp0' = invlogit( sqrt(`rspec')*`v4'+ sqrt(1-`rspec')*`v6') if "`sseca'" != "" { tempvar rn_sseca psseca rpsseca local wn : word count `sseca' tokenize "`sseca'" local dist = "`1'" if "`dist'" == "Constant" qui gen double `rn_sseca' = `2' if "`dist'" == "Uniform" { local min = `2' local max = `3' qui gen double `rn_sseca' = `max'-(`max'-`min')*`se1' } if ("`dist'" == "Triangular") | ("`dist'" == "Trapezoidal") { local min = `2' local mod1 = `3' local mod2 = `4' local max = `5' qui gen double `rn_sseca' = (`se1'*(`max'+`mod2'-`min'-`mod1') + (`min' + `mod1') )/2 qui replace `rn_sseca' = `min' + sqrt((`mod1'-`min')*(2*`rn_sseca'-`min'-`mod1')) if `rn_sseca' < `mod1' qui replace `rn_sseca' = `max' - sqrt( 2*(`max'-`mod2')*(`rn_sseca'-`mod2')) if `rn_sseca' > `mod2' } if "`dist'" == "Logit-Logistic" { local m = `2' local s = `3' local bl = `4' local bu = `5' gen double `psseca' = `m' + logit(`se1')*`s' gen double `rn_sseca' = `bl' + (`bu'-`bl')*invlogit(`psseca') } if "`dist'" == "Logit-Normal" { local m = `2' local s = `3' local bl = `4' local bu = `5' gen double `psseca' = `m' + invnorm(`se1')*`s' gen double `rn_sseca' = `bl' + (`bu'-`bl')*invlogit(`psseca') } } if "`sspca'" != "" { tempvar rn_sspca psspca rpsspca local wn : word count `sspca' tokenize "`sspca'" local dist = "`1'" if "`dist'" == "Constant" qui gen double `rn_sspca' = `2' if "`dist'" == "Uniform" { local min = `2' local max = `3' qui gen double `rn_sspca' = `max'-(`max'-`min')*`sp1' } if ("`dist'" == "Triangular") | ("`dist'" == "Trapezoidal") { local min = `2' local mod1 = `3' local mod2 = `4' local max = `5' qui gen double `rn_sspca' = (`sp1'*(`max'+`mod2'-`min'-`mod1') + (`min' + `mod1') )/2 qui replace `rn_sspca' = `min' + sqrt((`mod1'-`min')*(2*`rn_sspca'-`min'-`mod1')) if `rn_sspca' < `mod1' qui replace `rn_sspca' = `max' - sqrt( 2*(`max'-`mod2')*(`rn_sspca'-`mod2')) if `rn_sspca' > `mod2' } if "`dist'" == "Logit-Logistic" { local m = `2' local s = `3' local bl = `4' local bu = `5' gen double `psspca' = `m' + logit(`sp1')*`s' gen double `rn_sspca' = `bl' + (`bu'-`bl')*invlogit(`psspca') } if "`dist'" == "Logit-Normal" { local m = `2' local s = `3' local bl = `4' local bu = `5' gen double `psspca' = `m' + invnorm(`sp1')*`s' gen double `rn_sspca' = `bl' + (`bu'-`bl')*invlogit(`psspca') } } if "`ssenc'" != "" { tempvar rn_ssenc pssenc rpssenc local wn : word count `ssenc' tokenize "`ssenc'" local distssenc = "`1'" if "`distssenc'" == "Constant" qui gen double `rn_ssenc' = `2' if "`distssenc'" == "Uniform" { local min = `2' local max = `3' qui gen double `rn_ssenc' = `max'-(`max'-`min')*`se0' } if ("`distssenc'" == "Triangular") | ("`distssenc'" == "Trapezoidal") { local min = `2' local mod1 = `3' local mod2 = `4' local max = `5' qui gen double `rn_ssenc' = (`se0'*(`max'+`mod2'-`min'-`mod1') + (`min' + `mod1') )/2 qui replace `rn_ssenc' = `min' + sqrt((`mod1'-`min')*(2*`rn_ssenc'-`min'-`mod1')) if `rn_ssenc' < `mod1' qui replace `rn_ssenc' = `max' - sqrt( 2*(`max'-`mod2')*(`rn_ssenc'-`mod2')) if `rn_ssenc' > `mod2' } if "`distssenc'" == "Logit-Logistic" { local m = `2' local s = `3' local bl = `4' local bu = `5' gen double `pssenc' = `m' + logit(`se0')*`s' gen double `rn_ssenc' = `bl' + (`bu'-`bl')*invlogit(`pssenc') } if "`distssenc'" == "Logit-Normal" { local m = `2' local s = `3' local bl = `4' local bu = `5' gen double `pssenc' = `m' + invnorm(`se0')*`s' gen double `rn_ssenc' = `bl' + (`bu'-`bl')*invlogit(`pssenc') } } if "`sspnc'" != "" { tempvar rn_sspnc psspnc rpsspnc local wn : word count `sspnc' tokenize "`sspnc'" local distsspnc = "`1'" if "`distsspnc'" == "Constant" qui gen double `rn_sspnc' = `2' if "`distsspnc'" == "Uniform" { local min = `2' local max = `3' qui gen double `rn_sspnc' = `max'-(`max'-`min')*`sp0' } if ("`distsspnc'" == "Triangular") | ("`distsspnc'" == "Trapezoidal") { local min = `2' local mod1 = `3' local mod2 = `4' local max = `5' qui gen double `rn_sspnc' = (`sp0'*(`max'+`mod2'-`min'-`mod1') + (`min' + `mod1') )/2 qui replace `rn_sspnc' = `min' + sqrt((`mod1'-`min')*(2*`rn_sspnc'-`min'-`mod1')) if `rn_sspnc' < `mod1' qui replace `rn_sspnc' = `max' - sqrt( 2*(`max'-`mod2')*(`rn_sspnc'-`mod2')) if `rn_sspnc' > `mod2' } if "`distsspnc'" == "Logit-Logistic" { local m = `2' local s = `3' local bl = `4' local bu = `5' gen double `psspnc' = `m' + logit(`sp0')*`s' gen double `rn_sspnc' = `bl' + (`bu'-`bl')*invlogit(`psspnc') } if "`distsspnc'" == "Logit-Normal" { local m = `2' local s = `3' local bl = `4' local bu = `5' gen double `psspnc' = `m' + invnorm(`sp0')*`s' gen double `rn_sspnc' = `bl' + (`bu'-`bl')*invlogit(`psspnc') } } tempname rhose rhosp qui corr `rn_sseca' `rn_ssenc' scalar `rhose' = r(rho) qui corr `rn_sspca' `rn_sspnc' scalar `rhosp' = r(rho) if "`dstat'" != "" { qui corr `rn_sseca' `rn_ssenc' di as text "Correlation sensitivities : " %3.2f r(rho) char `rn_sseca'[varname] "Se|Cases" char `rn_ssenc'[varname] "Se|No-Cases" tabstat `rn_sseca' `rn_ssenc' , stat(min median max ) format(%3.2f) } if "`dstat'" != "" { qui corr `rn_sspca' `rn_sspnc' di as text "Correlation specificities : " %3.2f r(rho) char `rn_sspca'[varname] "Sp|Cases" char `rn_sspnc'[varname] "Sp|No-Cases" tabstat `rn_sspca' `rn_sspnc' , stat(min median max ) format(%3.2f) } if "`graph'" != "" qui hist `rn_sseca' if "`graph'" != "" qui hist `rn_sspca' if "`graph'" != "" qui hist `rn_ssenc' if "`graph'" != "" qui hist `rn_sspnc' } // end misclassification parameters * di as err "`spscex' `spscun' `spsnex' `spsnun' `ssbfactor' " if ( ("`spscex'" != "") & ("`spscun'" != "") & ("`spsnex'" != "") & ("`spsnun'" != "") ) /// | ("`ssbfactor'" != "") { if "`spscex'" != "" { tempvar rn_spscex pspscex rpspscex local wn : word count `spscex' tokenize "`spscex'" local dist = "`1'" if "`dist'" == "Constant" qui gen double `rn_spscex' = `2' if "`dist'" == "Uniform" { local min = `2' local max = `3' qui gen double `rn_spscex' = `max'-(`max'-`min')*uniform() } if ("`dist'" == "Triangular") | ("`dist'" == "Trapezoidal") { local min = `2' local mod1 = `3' local mod2 = `4' local max = `5' qui gen double `rn_spscex' = (uniform()*(`max'+`mod2'-`min'-`mod1') + (`min' + `mod1') )/2 qui replace `rn_spscex' = `min' + sqrt((`mod1'-`min')*(2*`rn_spscex'-`min'-`mod1')) if `rn_spscex' < `mod1' qui replace `rn_spscex' = `max' - sqrt( 2*(`max'-`mod2')*(`rn_spscex'-`mod2')) if `rn_spscex' > `mod2' } if "`dist'" == "Logit-Logistic" { local m = `2' local s = `3' local bl = `4' local bu = `5' gen double `pspscex' = `m' + logit(uniform())*`s' gen double `rn_spscex' = `bl' + (`bu'-`bl')*invlogit(`pspscex') } if "`dist'" == "Logit-Normal" { local m = `2' local s = `3' local bl = `4' local bu = `5' gen double `pspscex' = `m' + invnorm(uniform())*`s' gen double `rn_spscex' = `bl' + (`bu'-`bl')*invlogit(`pspscex') } } if "`spscun'" != "" { tempvar rn_spscun pspscun rpspscun local wn : word count `spscun' tokenize "`spscun'" local dist = "`1'" if "`dist'" == "Constant" qui gen double `rn_spscun' = `2' if "`dist'" == "Uniform" { local min = `2' local max = `3' qui gen double `rn_spscun' = `max'-(`max'-`min')*uniform() } if ("`dist'" == "Triangular") | ("`dist'" == "Trapezoidal") { local min = `2' local mod1 = `3' local mod2 = `4' local max = `5' qui gen double `rn_spscun' = (uniform()*(`max'+`mod2'-`min'-`mod1') + (`min' + `mod1') )/2 qui replace `rn_spscun' = `min' + sqrt((`mod1'-`min')*(2*`rn_spscun'-`min'-`mod1')) if `rn_spscun' < `mod1' qui replace `rn_spscun' = `max' - sqrt( 2*(`max'-`mod2')*(`rn_spscun'-`mod2')) if `rn_spscun' > `mod2' } if "`dist'" == "Logit-Logistic" { local m = `2' local s = `3' local bl = `4' local bu = `5' gen double `pspscun' = `m' + logit(uniform())*`s' gen double `rn_spscun' = `bl' + (`bu'-`bl')*invlogit(`pspscun') } if "`dist'" == "Logit-Normal" { local m = `2' local s = `3' local bl = `4' local bu = `5' gen double `pspscun' = `m' + invnorm(uniform())*`s' gen double `rn_spscun' = `bl' + (`bu'-`bl')*invlogit(`pspscun') } } if "`spsnex'" != "" { tempvar rn_spsnex pspsnex rpspsnex local wn : word count `spsnex' tokenize "`spsnex'" local dist = "`1'" if "`dist'" == "Constant" qui gen double `rn_spsnex' = `2' if "`dist'" == "Uniform" { local min = `2' local max = `3' qui gen double `rn_spsnex' = `max'-(`max'-`min')*uniform() } if ("`dist'" == "Triangular") | ("`dist'" == "Trapezoidal") { local min = `2' local mod1 = `3' local mod2 = `4' local max = `5' qui gen double `rn_spsnex' = (uniform()*(`max'+`mod2'-`min'-`mod1') + (`min' + `mod1') )/2 qui replace `rn_spsnex' = `min' + sqrt((`mod1'-`min')*(2*`rn_spsnex'-`min'-`mod1')) if `rn_spsnex' < `mod1' qui replace `rn_spsnex' = `max' - sqrt( 2*(`max'-`mod2')*(`rn_spsnex'-`mod2')) if `rn_spsnex' > `mod2' } if "`dist'" == "Logit-Logistic" { local m = `2' local s = `3' local bl = `4' local bu = `5' gen double `pspsnex' = `m' + logit(uniform())*`s' gen double `rn_spsnex' = `bl' + (`bu'-`bl')*invlogit(`pspsnex') } if "`dist'" == "Logit-Normal" { local m = `2' local s = `3' local bl = `4' local bu = `5' gen double `pspsnex' = `m' + invnorm(uniform())*`s' gen double `rn_spsnex' = `bl' + (`bu'-`bl')*invlogit(`pspsnex') } } if "`spsnun'" != "" { tempvar rn_spsnun pspsnun rpspsnun local wn : word count `spsnun' tokenize "`spsnun'" local dist = "`1'" if "`dist'" == "Constant" qui gen double `rn_spsnun' = `2' if "`dist'" == "Uniform" { local min = `2' local max = `3' qui gen double `rn_spsnun' = `max'-(`max'-`min')*uniform() } if ("`dist'" == "Triangular") | ("`dist'" == "Trapezoidal") { local min = `2' local mod1 = `3' local mod2 = `4' local max = `5' qui gen double `rn_spsnun' = (uniform()*(`max'+`mod2'-`min'-`mod1') + (`min' + `mod1') )/2 qui replace `rn_spsnun' = `min' + sqrt((`mod1'-`min')*(2*`rn_spsnun'-`min'-`mod1')) if `rn_spsnun' < `mod1' qui replace `rn_spsnun' = `max' - sqrt( 2*(`max'-`mod2')*(`rn_spsnun'-`mod2')) if `rn_spsnun' > `mod2' } if "`dist'" == "Logit-Logistic" { local m = `2' local s = `3' local bl = `4' local bu = `5' gen double `pspsnun' = `m' + logit(uniform())*`s' gen double `rn_spsnun' = `bl' + (`bu'-`bl')*invlogit(`pspsnun') } if "`dist'" == "Logit-Normal" { local m = `2' local s = `3' local bl = `4' local bu = `5' gen double `pspsnun' = `m' + invnorm(uniform())*`s' gen double `rn_spsnun' = `bl' + (`bu'-`bl')*invlogit(`pspsnun') } } if "`ssbfactor'" != "" { tempvar rn_ssbfactor local wn : word count `ssbfactor' tokenize "`ssbfactor'" local dist = "`1'" if "`dist'" == "Constant" qui gen double `rn_ssbfactor' = `2' if "`dist'" == "Uniform" { local bl = `2' local bu = `3' local min = `2' local max = `3' qui gen double `rn_ssbfactor' = `max'-(`max'-`min')*uniform() } if ("`dist'" == "Triangular") | ("`distssbfactor'" == "Trapezoidal") { local min = `2' local mod1 = `3' local mod2 = `4' local max = `5' local bl = `2' local bu = `5' qui gen double `rn_ssbfactor' = (uniform()*(`max'+`mod2'-`min'-`mod1') + (`min' + `mod1') )/2 qui replace `rn_ssbfactor' = `min' + sqrt((`mod1'-`min')*(2*`rn_ssbfactor'-`min'-`mod1')) if `rn_ssbfactor' < `mod1' qui replace `rn_ssbfactor' = `max' - sqrt( 2*(`max'-`mod2')*(`rn_ssbfactor'-`mod2')) if `rn_ssbfactor' > `mod2' } if "`dist'" == "Log-Normal" { local m = `2' local s = `3' gen double `rn_ssbfactor' = exp( `m' + invnorm(uniform())*`s') tempvar log_rn_ssbfactor gen `log_rn_ssbfactor' = log(`rn_ssbfactor') } if "`dist'" == "Log-Logistic" { local m = `2' local s = `3' gen double `rn_ssbfactor' = exp( `m' + logit(uniform())*`s') tempvar log_rn_ssbfactor gen `log_rn_ssbfactor' = log(`rn_ssbfactor') } } if "`graph'" != "" & "`spscex'" != "" qui hist `rn_spscex' if "`graph'" != "" & "`spscun'" != "" qui hist `rn_spscun' if "`graph'" != "" & "`spsnex'" != "" qui hist `rn_spsnex' if "`graph'" != "" & "`spsnun'" != "" qui hist `rn_spsnun' if "`graph'" != "" & "`ssbfactor'" != "" qui hist `rn_ssbfactor' } // end selection bias parameters if ("`spunexp'"!="") & ( ("`spexp'"!= "") | ("`sorce'"!="") ) & ("`srrcd'"!= "") { * Generate a pairs of 0,1 random variables with user-specified correlation tempvar h1 h2 h3 p1 p0 tempname rprev forv i = 1/3 { gen `h`i'' = logit(uniform()) } if "`scorrprev'" != "" scalar `rprev' = `scorrprev' else scalar `rprev' = 0 gen `p1' = invlogit( sqrt(`rprev')*`h1'+ sqrt(1-`rprev')*`h2') gen `p0' = invlogit( sqrt(`rprev')*`h1'+ sqrt(1-`rprev')*`h3') if "`spexp'" != "" { tempvar rn_spexp local wn : word count `sseca' tokenize "`spexp'" local dist = "`1'" if "`dist'" == "Constant" qui gen double `rn_spexp' = `2' if "`dist'" == "Uniform" { local min = `2' local max = `3' qui gen double `rn_spexp' = `max'-(`max'-`min')*`p1' } if ("`dist'" == "Triangular") | ("`dist'" == "Trapezoidal") { local min = `2' local mod1 = `3' local mod2 = `4' local max = `5' qui gen double `rn_spexp' = (`p1'*(`max'+`mod2'-`min'-`mod1') + (`min' + `mod1') )/2 qui replace `rn_spexp' = `min' + sqrt((`mod1'-`min')*(2*`rn_spexp'-`min'-`mod1')) if `rn_spexp' < `mod1' qui replace `rn_spexp' = `max' - sqrt( 2*(`max'-`mod2')*(`rn_spexp'-`mod2')) if `rn_spexp' > `mod2' } if "`dist'" == "Logit-Logistic" { local m = `2' local s = `3' local bl = `4' local bu = `5' gen double `rn_spexp' = `bl' + (`bu'-`bl')*invlogit(`m' + logit(`p1')*`s' ) } if "`dist'" == "Logit-Normal" { local m = `2' local s = `3' local bl = `4' local bu = `5' gen double `rn_spexp' = `bl' + (`bu'-`bl')*invlogit(`m' + invnorm(`p1')*`s' ) } } if "`spunexp'" != "" { tempvar rn_spunexp local wn : word count `spunexp' tokenize "`spunexp'" local dist = "`1'" if "`dist'" == "Constant" qui gen double `rn_spunexp' = `2' if "`dist'" == "Uniform" { local min = `2' local max = `3' qui gen double `rn_spunexp' = `max'-(`max'-`min')*`p0' } if ("`dist'" == "Triangular") | ("`dist'" == "Trapezoidal") { local min = `2' local mod1 = `3' local mod2 = `4' local max = `5' qui gen double `rn_spunexp' = (`p0'*(`max'+`mod2'-`min'-`mod1') + (`min' + `mod1') )/2 qui replace `rn_spunexp' = `min' + sqrt((`mod1'-`min')*(2*`rn_spunexp'-`min'-`mod1')) if `rn_spunexp' < `mod1' qui replace `rn_spunexp' = `max' - sqrt( 2*(`max'-`mod2')*(`rn_spunexp'-`mod2')) if `rn_spunexp' > `mod2' } if "`dist'" == "Logit-Logistic" { local m = `2' local s = `3' local bl = `4' local bu = `5' gen double `rn_spunexp' = `bl' + (`bu'-`bl')*invlogit(`m' + invnorm(`p0')*`s' ) } if "`dist'" == "Logit-Normal" { local m = `2' local s = `3' local bl = `4' local bu = `5' gen double `rn_spunexp' = `bl' + (`bu'-`bl')*invlogit(`m' + invnorm(`p0')*`s' ) } } if "`srrcd'" != "" { tempvar rn_srrcd local wn : word count `srrcd' tokenize "`srrcd'" local distsrrcd = "`1'" if "`distsrrcd'" == "Constant" qui gen double `rn_srrcd' = `2' if "`distsrrcd'" == "Uniform" { local blsrrcd = `2' local busrrcd = `3' local min = `2' local max = `3' qui gen double `rn_srrcd' = `max'-(`max'-`min')*uniform() } if ("`distsrrcd'" == "Triangular") | ("`distsrrcd'" == "Trapezoidal") { local min = `2' local mod1 = `3' local mod2 = `4' local max = `5' local blsrrcd = `2' local busrrcd = `5' qui gen double `rn_srrcd' = (uniform()*(`max'+`mod2'-`min'-`mod1') + (`min' + `mod1') )/2 qui replace `rn_srrcd' = `min' + sqrt((`mod1'-`min')*(2*`rn_srrcd'-`min'-`mod1')) if `rn_srrcd' < `mod1' qui replace `rn_srrcd' = `max' - sqrt( 2*(`max'-`mod2')*(`rn_srrcd'-`mod2')) if `rn_srrcd' > `mod2' } if "`distsrrcd'" == "Log-Normal" { local m = `2' local s = `3' gen double `rn_srrcd' = exp( `m' + invnorm(uniform())*`s') tempvar log_rn_srrcd gen `log_rn_srrcd' = log(`rn_srrcd') } if "`distsrrcd'" == "Log-Logistic" { local m = `2' local s = `3' gen double `rn_srrcd' = exp( `m' + logit(uniform())*`s') tempvar log_rn_srrcd gen `log_rn_srrcd' = log(`rn_srrcd') } } if "`sorce'" != "" { tempvar rn_sorce u gen `u' = uniform() local wn : word count `sorce' tokenize "`sorce'" local distsorce = "`1'" if "`distsorce'" == "Constant" qui gen double `rn_sorce' = `2' if "`distsorce'" == "Uniform" { local blsorce = `2' local busorce = `3' local min = `2' local max = `3' qui gen double `rn_sorce' = `max'-(`max'-`min')*uniform() } if ("`distsorce'" == "Triangular") | ("`distsorce'" == "Trapezoidal") { local min = `2' local mod1 = `3' local mod2 = `4' local max = `5' local blsorce = `2' local busorce = `5' qui gen double `rn_sorce' = (uniform()*(`max'+`mod2'-`min'-`mod1') + (`min' + `mod1') )/2 qui replace `rn_sorce' = `min' + sqrt((`mod1'-`min')*(2*`rn_sorce'-`min'-`mod1')) if `rn_sorce' < `mod1' qui replace `rn_sorce' = `max' - sqrt( 2*(`max'-`mod2')*(`rn_sorce'-`mod2')) if `rn_sorce' > `mod2' } if ("`distsorce'" == "Log-Normal") { local m = `2' local s = `3' gen double `rn_sorce' = exp( `m' + invnormal(uniform())*`s' ) tempvar log_rn_sorce gen `log_rn_sorce' = log(`rn_sorce') } if ("`distsorce'" == "Log-Logistic") { local m = `2' local s = `3' gen double `rn_sorce' = exp( `m' + logit(uniform())*`s' ) tempvar log_rn_sorce gen `log_rn_sorce' = log(`rn_sorce') } } tempname rhoprev qui corr `rn_spexp' `rn_spunexp' scalar `rhoprev' = r(rho) if "`dstat'" != "" { qui corr `rn_spexp' `rn_spunexp' di as text "Correlation prevalences : " %3.2f r(rho) tabstat `rn_spexp' `rn_spunexp' , stat(min median max ) format(%3.2f) } if "`graph'" != "" qui hist `rn_spexp' if "`graph'" != "" qui hist `rn_spunexp' if "`graph'" != "" qui hist `log_rn_srrcd' if "`graph'" != "" & "`sorce'" != "" qui hist `log_rn_sorce' } // end unmeasured confounding * keep a random subset of the bias parameters qui sample 1, count tempname rrstep rrdx scalar `rrstep' = `apprr' if ("`sseca'"!="") & ("`sspca'"!="") & ("`ssenc'"!="") & ("`sspnc'"!="") { * From whatever prior distribution the bias parameter are coming from put the values into scalars for the final calculations tempname seec spec seeo speo scalar `seec' = `rn_sseca'[1] scalar `spec' = `rn_sspca'[1] scalar `seeo' = `rn_ssenc'[1] scalar `speo' = `rn_sspnc'[1] tempname fnec fpec fneo fpeo scalar `fnec' = 1-`seec' scalar `fpec' = 1-`spec' scalar `fneo' = 1-`seeo' scalar `fpeo' = 1-`speo' * calculate the bias-adjusted RR tempname b0 b1 a1 a0 percent_bias adj_factor scalar `b1' = (`speo'*`c'-`fpeo'*`d')/(`seeo'*`speo'-`fneo'*`fpeo') scalar `b0' = (`c'+`d')-`b1' scalar `a1' = (`spec'*`a'-`fpec'*`b')/(`seec'*`spec'-`fnec'*`fpec') scalar `a0' = (`a'+`b')-`a1' if "`studytype'" == "cc" scalar `rrdx' = (`a1'*`b0')/(`b1'*`a0') if "`studytype'" == "ir" scalar `rrdx' = (`a1'/`b1')/(`a0'/`b0') if "`studytype'" == "cs" scalar `rrdx' = (`a1'/(`a1'+`b1'))/(`a0'/(`a0'+`b0')) * Treatment of Negative Adjustment * Discard draws that produce negative cell counts * di as y `a1' " " `a0' " " `b1' " " `b0' * Discard draws that produce counts < 1 if (round(`a1')<1) | (round(`a0')<1) | (round(`b1')<1) | (round(`b0')<1) { scalar `rrdx' = . scalar `a1' = . scalar `a0' = . scalar `b1' = . scalar `b0' = . } * Discard draws that produce counts <=0 if (`seeo' < `c'/`m0') | (`speo' < `d'/`m0') | (`seec' < `a'/`m1') | (`spec' < `b'/`m1') { scalar `rrdx' = . scalar `a1' = . scalar `a0' = . scalar `b1' = . scalar `b0' = . } scalar `rrstep' = `rrdx' di "Adj RR step 1 = " %3.2f `rrstep' return scalar pseca = `seec' return scalar psenc = `seeo' return scalar pspca = `spec' return scalar pspnc = `speo' return scalar rhose = `rhose' return scalar rhosp = `rhosp' } // end misclassification if [ ("`spscex'" != "") & ("`spscun'" != "") & ("`spsnex'" != "") & ("`spsnun'" != "") ] /// | ("`ssbfactor'" != "") { * From whatever prior distribution the bias parameter are coming from put the values into scalars for the final calculations tempname pscex pscun psnex psnun sbfact sel_bias_factor percent_bias adj_factor scalar `pscex' = `rn_spscex'[1] scalar `pscun' = `rn_spscun'[1] scalar `psnex' = `rn_spsnex'[1] scalar `psnun' = `rn_spsnun'[1] scalar `sbfact' = `rn_ssbfactor'[1] /* di "pscex = " `pscex' di "pscun = " `pscun' di "psnex =" `psnex' di "psnun =" `psnun' di "sbfact =" `sbfact' */ * calculate the bias-adjusted RR if "`ssbfactor'" == "" scalar `sel_bias_factor' = (`pscex'*`psnun')/(`pscun'*`psnex') else scalar `sel_bias_factor' = `sbfact' scalar `rrdx' = `rrstep'/`sel_bias_factor' scalar `percent_bias' = (`apprr'-`rrdx')/(`rrdx')*100 scalar `rrstep' = `rrdx' di "Adj RR step 2 = " %3.2f `rrstep' if ("`ssbfactor'" == "") { return scalar pscex = `pscex' return scalar pscun = `pscun' return scalar psnex = `psnex' return scalar psnun = `psnun' return scalar sbf = 1 } else { return scalar pscex = 1 return scalar pscun = 1 return scalar psnex = 1 return scalar psnun = 1 return scalar sbf = `sbfact' } } // end selection bias if ("`spunexp'"!="") & ( ("`spexp'"!= "") | ("`sorce'"!="") ) & ("`srrcd'"!= "") { * From whatever prior distribution the bias parameter are coming from put the values into scalars for the final calculations tempname prz1 prz0 rrdz orze rrdx rrxz percent_bias adj_factor scalar `prz1' = `rn_spexp'[1] scalar `prz0' = `rn_spunexp'[1] scalar `rrdz' = `rn_srrcd'[1] scalar `orze' = `rn_sorce'[1] if ("`sorce'" == "") { scalar `rrdx' = `rrstep' / [ (`prz1'*(`rrdz'-1)+1 )/(`prz0'*(`rrdz'-1)+1 ) ] } if ("`sorce'" != "") { scalar `rrdx' = `rrstep'/[ (`rrdz'*`orze'*`prz0'+1-`prz0') / ( (`rrdz'*`prz0'+1-`prz0')*( `orze'*`prz0'+1-`prz0') ) ] } scalar `percent_bias' = (`arrdx'-`rrdx')/(`rrdx')*100 scalar `rrstep' = `rrdx' di "Adj RR step 3 = " %3.2f `rrstep' if ("`sorce'" != "") return scalar orce = `orze' if ("`sorce'" == "") return scalar orce = [(`prz1')*(1-`prz0')]/[(1-`prz1')*(`prz0')] return scalar pc0 = `prz0' return scalar rrcd = `rrdz' return scalar pc1 = `prz1' return scalar rhoprev = `rhoprev' } // end unmeasured confounding // Saved result return scalar adj_rr_all = `rrstep' return scalar perc_bias = (`apprr'-`rrstep')/`rrstep'*100 end /* set trace off set tracedepth 1 episens_mcsa_all 45 94 257 945, obs(5000) /// sseca(Uniform .6 .8) sspca(Triangular .6 .8 .8 .9) /// ssenc(Trapezoidal .5 .6 .8 .9) sspnc(Logit-Normal 0 1 0.7 1) /// spscex(Uniform .8 1) spscun(Uniform .8 1) /// spsnex(Constant .8) spsnun(Uniform .9 1) /// spexp(Uniform .4 .7) spunexp(Uniform .4 .7) /// srrcd(Log-Normal 2.159 .280) /// apprr(1.76) applb(1.17) appub(2.61) studytype(cc) episens_mcsa_all 45 94 257 945, obs(5000) /// sseca(Uniform .6 .8) sspca(Triangular .6 .8 .8 .9) ssenc(Trapezoidal .5 .6 .8 .9) sspnc(Logit-Normal 0 1 0.7 1) /// ssbfactor(Log-Normal 0 .207) /// spexp(Uniform .4 .7) spunexp(Uniform .4 .7) srrcd(Log-Normal 2.159 .280) apprr(1.76) applb(1.17) appub(2.61) studytype(cc) */