*! Package rori v. 0.12 *! Support: Niels Henrik Bruun, niels.henrik.bruun@gmail.com *! Rewritten due to comments. Thanks to Eric Melse * 2022-05-29 Version bug fixed * 2022-05-23 created *capture program drop rori program define rori, rclass version 12 syntax, /* */Targetpopulation(string) /*ABCD[EF...] or AC[E...] or BC[F...] */Responsepopulation(string) /*abcd[ef...] or ac[e...] or bc[e...] */[ /* */REFerence(integer 1) /* */] capture matrix _t = `targetpopulation' if _rc mata _error("Target is not a matrix or a matrix definition") capture matrix _r = `responsepopulation' if _rc mata _error("Response is not a matrix or a matrix definition") mata: rori(st_matrix("_t"), st_matrix("_r"), `reference') return add end mata mata clear mata real matrix ci(real matrix ln_m_se) { real scalar r, R, m, se, z real matrix out z = invnormal(0.5 + c("level") * 0.005) R = rows(ln_m_se) out = J(R,1,(1,.,.)) for(r=1;r<=R;r++){ if ( ln_m_se[r,.] == (.,.) ) out[r,1] = . se = ln_m_se[r,2] if ( se >= . ) continue m = ln_m_se[r,1] out[r,.] = exp((m, m - z * se, m + z * se)) } return(out) } real matrix ln_mean_se_rsp_tgt(real matrix rsp, real matrix tgt) { real scalar r, R, m, se R = rows(rsp) m_se = J(R,2,.) for(r=1; r<=R; r++){ se = sqrt(rsp[r,2]:^2 - tgt[r,2]:^2) if ( missing(se) ) continue m = rsp[r,1] - tgt[r,1] m_se[r,.] = m, se } return(m_se) } real matrix ln_p_se(real colvector n) { real colvector lnp, se_lnp lnp = log(n :/ colsum(n)) se_lnp = sqrt(1 :/ n) return((lnp, se_lnp)) } real matrix ln_or_se(real matrix m, | real scalar ref) { real scalar r, R real matrix tbl, m_se ref = missing(ref) ? 1 : ref R = rows(m) m_se = J(R, 1, (1,.)) for(r=1; r<=R; r++) { if ( r == ref ) continue tbl = m[(ref,r), .] m_se[r,.] = (1, -1) * ln(tbl) * (1, -1)', sqrt(sum(1 :/ tbl)) } return(m_se) } real scalar validate_target_respons(real matrix tgt, real matrix rsp) { if ( cols(tgt) != cols(rsp) | rows(tgt) != rows(rsp) ) _error("Target and response matrix must have same shape.") if ( cols(tgt) == 0 | cols(tgt) > 2 ) _error("Target matrix must have 1 or 2 columns") if ( cols(rsp) == 0 | cols(rsp) > 2 ) _error("respons matrix must have 1 or 2 columns") if ( rows(tgt) < 2 ) _error("Target matrix must have at least 2 rows") if ( rows(rsp) < 2 ) _error("response matrix must have at least 2 rows") if ( any(tgt :< rsp) ) _error("Target cells must be greater than response cells") return(cols(tgt) == 2) } void rori(real matrix tgt, real matrix rsp, | real scalar ref) { real scalar ror_rrf real matrix t_m_se1, r_m_se1, t_m_se2, r_m_se2, rrf, ror string scalar ci_txt st_rclear() ror_rrf = validate_target_respons(tgt, rsp) ci_txt = sprintf("[%2.0f%%", c("level")) if ( ror_rrf ) { r_m_se1 = ln_p_se(rsp[.,1]) t_m_se1 = ln_p_se(tgt[.,1]) r_m_se2 = ln_p_se(rsp[.,2]) t_m_se2 = ln_p_se(tgt[.,2]) rrf = ci(r_m_se1), ci(t_m_se1), ci(ln_mean_se_rsp_tgt(r_m_se1, t_m_se1)), ci(r_m_se2), ci(t_m_se2), ci(ln_mean_se_rsp_tgt(r_m_se2, t_m_se2)) st_matrix("r(rrf)", rrf) st_matrixcolstripe("r(rrf)", (J(1,3,"rsp(col1)"), J(1,3,"tgt(col1)"), J(1,3,"rsp(col1) vs tgt(col1)"), J(1,3,"rsp(col2)"), J(1,3,"tgt(col2)"), J(1,3,"rsp(col2) vs tgt(col2)") \ J(1,2, ("P", ci_txt, "CI]")), "RRF", ci_txt, "CI]", J(1,2, ("P", ci_txt, "CI]")), "RRF", ci_txt, "CI]")') st_matrixrowstripe("r(rrf)", st_matrixrowstripe("_t")) r_m_se1 = ln_or_se(rsp, ref) t_m_se1 = ln_or_se(tgt, ref) ror = ci(r_m_se1), ci(t_m_se1), ci(ln_mean_se_rsp_tgt(r_m_se1, t_m_se1)) st_matrix("r(ror)", ror) st_matrixcolstripe("r(ror)", (J(1,3,"rsp"), J(1,3,"tgt"), J(1,3,"rsp vs tgt") \ J(1,2, ("OR", ci_txt, "CI]")), "ROR", ci_txt, "CI]")') st_matrixrowstripe("r(ror)", st_matrixrowstripe("_t")) } else { rrf = (ci(ln_p_se(rsp)), ci(ln_p_se(tgt)), ci(ln_mean_se_rsp_tgt(ln_p_se(rsp), ln_p_se(tgt)))) st_matrix("r(rrf)", rrf) st_matrixcolstripe("r(rrf)", (J(1,3,"rsp"), J(1,3,"tgt"), J(1,3,"rsp vs tgt") \ J(1,2, ("P", ci_txt, "CI]")), "RRF", ci_txt, "CI]")') st_matrixrowstripe("r(rrf)", st_matrixrowstripe("_t")) } } end exit /* cls rori, t(22193 \ 17277 \ 9065) r(7719 \ 5297 \ 2291) return list matprint r(rrf) rori, t(30, 401 \ 63, 2807) r(25, 347 \ 45, 2313) return list matprint r(rrf) matprint r(ror) */ mata: real matrix rrfbs( real colvector n_target, real colvector n_response, | real scalar reps ) { real scalar median, ll, ul real colvector p_response, p_target, rvst real matrix bsrrf reps = missing(reps) ? 1000 : reps median = round(reps * 0.5) ll = round(reps * (0.5 - c("level") * 0.005)) ul = round(reps * (0.5 + c("level") * 0.005)) p_target = n_target :/ colsum(n_target) lnvar_p_target = 1 :/ n_target p_response = n_response :/ colsum(n_response) lnvar_p_response = 1 :/ n_response bsrrf = J(3,3,.) for(r=1;r<=3;r++) { rsp = rpoisson(reps, 1, n_response[r]) // rsp < tgt rvst = (sort(rsp, 1) :/ colsum(n_response)) :/ (sort(rpoisson(reps, 1, n_target[r]), 1) :/ colsum(n_target)) bsrrf[r,.] = rvst[(median, ll, ul)]' } return(bsrrf) } end