*! v 2.0.0 N.Orsini  13sept2013
*! v 1.0.0 N.Orsini  7july2006

capture program drop postrri
program postrri, rclass
version 8.2
syntax , Prior(string) Data(string)  [ Level(int $S_level) LPrior(int $S_level) ///
LData(int $S_level) Format(string) DFormat(string) CHeck GRaph * ]

if `level' <10 | `level'>99 {
di in red "level() invalid"
exit 198
}

if `lprior' <10 | `lprior'>99 {
di in red "lprior() invalid"
exit 198
}

if `ldata' <10 | `ldata'>99 {
di in red "ldata() invalid"
exit 198
}

if "`format'" == "" {
local format = "%3.2f"
}
else {
local format = "`format'"
}

tempname rr_p lb_p ub_p rr_d lb_d ub_d logrr_post rr_post se_post var_post lb_post ub_post levelci levelci_p levelci_d  se_p se_d w_se_p w_se_d

scalar `levelci' = `level' * 0.005 + 0.50
scalar `levelci_p' = `lprior' * 0.005 + 0.50
scalar `levelci_d' = `ldata' * 0.005 + 0.50

if "`prior'" == "" {
	di in red "specify prior LB and UB"
	exit 198
}
else {
	tokenize "`prior'"
	confirm number `1'
	confirm number `2'
	scalar `lb_p' = `1'
	scalar `ub_p' = `2'
	scalar `rr_p' =  exp([log(`lb_p')+log(`ub_p')]/2)
}

if "`data'" == "" {
	di in red "specify data RR, LB and UB"
	exit 198
}
else {
	tokenize "`data'"
	confirm number `1'
	confirm number `2'
	confirm number `3'

	scalar `rr_d' = `1'
	scalar `lb_d' = `2'
	scalar `ub_d' = `3'
}

// check values

if  inrange(`rr_p', 0, .) != 1 {
		di in red "prior RR has to be >= 0"
		exit 198
	}

if  `rr_d' < 0 |  `rr_d' == . {
		di in red "data RR has to be >= 0"
		exit 198
	}

// calculate variances we assume by default 95% CI for both prior and data

scalar `se_p' =  ( log(`ub_p') - log(`lb_p') )/(2*invnorm(`levelci_p'))
scalar `se_d' =  ( log(`ub_d') - log(`lb_d') )/(2*invnorm(`levelci_d'))

// calculate weights

scalar `w_se_p' =  1/`se_p'^2
scalar `w_se_d' =  1/`se_d'^2

// Approximate-bayesian method

// Calculate posterior relative risk and 95% CI

scalar `rr_post' = exp( (log(`rr_p')*`w_se_p' + log(`rr_d')*`w_se_d')/(`w_se_p' +`w_se_d') )
scalar `var_post' = 1/ [1/`se_p'^2 + 1/`se_d'^2]
scalar `logrr_post' = log(`rr_post')
scalar `se_post' = sqrt(`var_post')

scalar `lb_post' = exp( log(`rr_post') - invnorm(`levelci') * `se_post')
scalar `ub_post' = exp( log(`rr_post') + invnorm(`levelci') * `se_post' )

// display results

di _n as text "Posterior median for RR = " in y `format' `rr_post'
di in g "`level'% posterior limits for RR (" `format' in y `lb_post' in g ", " in y `format' `ub_post' in g ")"

// Check compatibility between prior and data

tempname std_diff p_check

scalar `std_diff' = abs(log(`rr_d') - log(`rr_p'))/sqrt(`se_d'^2 + `se_p'^2)
scalar `p_check' = normal(-abs(`std_diff'))*2

if ("`dformat'" == "") {
						local fmt "%4.3fc"
						}
else {
		local fmt "`dformat'"
}

local obs : di "Data Normal(" `fmt'  log(`rr_d') ", " `fmt' `se_d'^2 ")"
local prior : di "Prior Normal(" `fmt' log(`rr_p') ", " `fmt' `se_p'^2 ")"
local posterior : di "Posterior Normal(" `fmt' `logrr_post' ", " `fmt' `se_post'^2 ")"

if "`check'" != "" {
		di "Homogeneity statistic = " as res %5.4fc `std_diff' as txt "    p-value = " as res %5.4fc `p_check'
}

// graph

if "`graph'" != "" {

tempname obs_logrr prior_logrr obs_se obs_se prior_se
scalar `obs_logrr' = log(`rr_d')
scalar `prior_logrr' = log(`rr_p')
scalar `obs_se' = `se_d'
scalar `prior_se' = `se_p'

local dpostrri : display  "Posterior Relative Risk " `format' `rr_post' " (95% Posterior Limits "  `format' `lb_post' ", "  `format' `ub_post' ")"
local dpriorrri : display "Prior        Relative Risk " `format' `rr_p' " (95% Confidence Limits "  `format' `lb_p' ", "  `format' `ub_p' ")"
local ddatarri : display  "Data        Relative Risk " `format' `rr_d' " (95% Confidence Limits "  `format' `lb_d' ", "  `format' `ub_d' ")"

qui twoway  ///
 (function y=normalden(x, `prior_logrr', `prior_se'), range(`=`prior_logrr'-2' `=`prior_logrr'+2') lc(red) lw(thick) lp(l))   ///
   (function y=normalden(x, `obs_logrr', `obs_se'), range(`=`obs_logrr'-2' `=`obs_logrr'+2')  lc(blue) lw(thick) lp(l))   ///
   (function y=normalden(x, `logrr_post', `se_post'), range(`=`logrr_post'-2' `=`logrr_post'+2') lc(green) lw(thick) lp(l))   ///
  ,  ///
  plotregion(style(none)) ///
  yscale(off) ///
  legend( label(1 `"`prior'"') label(2  `"`obs'"')   label(3 `"`posterior'"') ring(1) pos(11) col(1) region(style(none)) ) ///
  scheme(s1mono)  ///
  xlabel(#10) xtitle("log Relative Risk") ///
  note("`dpriorrri'" "`ddatarri'" "`dpostrri'", place(sw)   just(left)  ) ///
  `options'
}


// saved results

return local cmd = "postrri"
return scalar lb_prior = `lb_p'
return scalar ub_prior = `ub_p'
return scalar rr_prior = `rr_p'
return scalar rr_post = `rr_post'
return scalar lb_post = `lb_post'
return scalar ub_post = `ub_post'
return scalar logrr_post = log(`rr_post')
return scalar var_logrr_post =  `se_post'
return scalar var_logrr_prior = `se_p'^2
return scalar var_logrr_data =  `se_d'^2
return scalar h_t = `std_diff'
return scalar p_h = `p_check'
end


exit

postrri , prior(0.25 4) data(3.51 0.8 15.4)  graph