* Version 1.0.0 Date: 13th August 2026 program drop _all program define msrrd, rclass version 17.0 * Syntax for the command syntax varname /// [ , subgroup(varname) subtype(string) /// rrdp(integer 4) rddpprop(integer 4) rddppct(integer 2) pdp(integer 4) pctdp(integer 1) vceuncond ] * Local macros local v "`varlist'" local subgroup "`subgroup'" local subtype = lower("`subtype'") local y = e(depvar) local subtype_specified = ("`subtype'" != "") local vceuncond "`vceuncond'" * Default option for subtype is to present both RR and RD if ("`subtype'" == "") local subtype "both" * --- VALIDATE SUBTYPE OPTION --- if !inlist("`subtype'", "rr", "rd", "both") { display as error "ERROR: subtype(`subtype') is invalid." display as error "ERROR: Valid subtype options are: rr, rd, or both." exit 198 } * Default decimal places for percentage in n/N (%) is 1 dp if ("`pctdp'" == "") local pctdp = 1 ********************************************************************************************************************************************************* * MODEL CHECKS - The following will ensure command only works if a logit/logistic/melogit or a glm/meglm with family(binomial) and link(logit) was fitted ********************************************************************************************************************************************************* cap local model = e(cmd) cap local model2 = e(cmd2) cap local family = lower(e(varfunct)) cap local family2 = lower(e(family)) cap local linkt = lower(e(linkt)) cap local link = lower(e(link)) local is_logit_family = ("`model'" == "logit" | "`model'" == "logistic" | "`model2'" == "melogit") local is_glm_binlogit = ("`model'" == "glm" & ("`family'" == "binomial" | "`family'" == "bernoulli") & "`linkt'" == "logit") local is_meglm_binlogit = ("`model2'" == "meglm" & ("`family2'" == "binomial" | "`family2'" == "bernoulli") & "`link'" == "logit") if !`is_logit_family' & !`is_glm_binlogit' & !`is_meglm_binlogit' { display as error "ERROR: msrrd requires a logistic model (logit/logistic/melogit) or a glm/meglm model fitted with family(binomial) and link(logit)." display as error "ERROR: Re-fit the model accordingly and re-run command." exit 198 } ***************************************************************************************************************** * VCE CHECK - This will display a warning prompting used that fitted model was without robust/cluster robust SE's ***************************************************************************************************************** local vce = e(vce) if ("`vce'" != "robust" & "`vce'" != "cluster") { display as error "WARNING: Fitted model did not include vce(robust) or vce(cluster) option." display as error "WARNING: When estimating RR/RD using marginal standardisation, it is recommended that the fitted model to have robust/cluster robust SE's." } ******************************************************************************************************************** * CODING CHECKS - This ensures the command will only run if both outcome and intervention variables are coded as 0/1 ******************************************************************************************************************** quietly levelsof `y', local(outlevels) if "`outlevels'" != "0 1" { display as error "ERROR: Outcome must be coded 0/1." exit 198 } quietly levelsof `v', local(trtlevels) if "`trtlevels'" != "0 1" { display as error "ERROR: For command {cmd:{it:mssrd}} to work, variabel {cmd:{it:`v'}} must be coded 0/1, where 0=Comparator group and 1=Intervention group." exit 198 } ************************************************************************************************************** * This will ensure if subtype option is selected without the subgroup option then an error message will appear ************************************************************************************************************** if ("`subgroup'" == "") & `subtype_specified' { display as error "ERROR: subgroup option not selected." exit 198 } ******************************************************************************************* * STORE ORIGINAL MODEL - This will retain model estimates to be reused later on when needed ******************************************************************************************* estimates store MSRRD_MODEL_ORIG ************************************************************************************************************************************************************* **************************************************************************************************** * OVERALL MARGINAL STANDARDISATION - This section of the code is for overall result without subgroup **************************************************************************************************** if ("`subgroup'" == "") & ("`vceuncond'" != "") { display "" display in yellow "Running marginal standardisation with {cmd:{it:vce(unconditional)}} to estimate SEs allowing for sampling of covariates" } if ("`subgroup'" == "") & ("`vceuncond'" == "") { display "" display in yellow "Running marginal standardisation with {cmd:{it:vce(delta)}} estimate SEs using delta method" } * --- MARGINS VCE CONTROL --- if ("`vceuncond'" != "") { quietly margins `v', post vce(unconditional) } if ("`vceuncond'" == "") { quietly margins `v', post } ************************* * EXTRACT MARGINS RESULTS ************************* matrix M = r(b) matrix VM = r(V) local p0 = M[1,1] local p1 = M[1,2] local var_p0 = VM[1,1] local var_p1 = VM[2,2] local cov_p1p0 = VM[2,1] *********************************************************************** * RR esimation using marginal standardisation and SE using delta method *********************************************************************** local rr = `p1'/`p0' local logrr = ln(`rr') local varlogrr = (`var_p1'/(`p1'^2)) + (`var_p0'/(`p0'^2)) - (2*`cov_p1p0'/(`p1'*`p0')) local rr_se = sqrt(`varlogrr') local log_rr_se = `rr'*`rr_se' local rr_lb = exp(`logrr' - invnorm(0.975)*`rr_se') local rr_ub = exp(`logrr' + invnorm(0.975)*`rr_se') local rr_z = `logrr'/`rr_se' local rr_p = 2*normal(-abs(`rr_z')) *********************************************************************** * RD esimation using marginal standardisation and SE using delta method *********************************************************************** local rd = `p1' - `p0' local vardiff = `var_p1' + `var_p0' - 2*`cov_p1p0' local rd_se = sqrt(`vardiff') local rd_lb = `rd' - invnorm(0.975)*`rd_se' local rd_ub = `rd' + invnorm(0.975)*`rd_se' local rd_z = `rd'/`rd_se' local rd_p = 2*normal(-abs(`rd_z')) ********* * FORMATS ********* local fmt_rr "%0`=2+`rrdp''.`rrdp'f" local fmt_rd_prop "%0`=2+`rddpprop''.`rddpprop'f" local fmt_rd_pct "%0`=2+`rddppct''.`rddppct'f" local fmt_p "%0`=2+`pdp''.`pdp'f" local fmt_pct "%0`=2+`pctdp''.`pctdp'f" ********************************************** * OVERALL EVENT COUNTS (with pctdp formatting) ********************************************** quietly count if `v'==0 & `y'!=. local N0 = r(N) quietly count if `v'==0 & `y'==1 local n0 = r(N) quietly count if `v'==1 & `y'!=. local N1 = r(N) quietly count if `v'==1 & `y'==1 local n1 = r(N) local pct0_raw = 100*`n0'/`N0' local pct1_raw = 100*`n1'/`N1' local pct0 : display `fmt_pct' `pct0_raw' local pct1 : display `fmt_pct' `pct1_raw' ************************ * OVERALL COMBINED TABLE ************************ if ("`subgroup'" == "") { * Extract value labels for variable local lb1 : label (`v') 1 local lb0 : label (`v') 0 * Fallback if no labels exist if ("`lbl'" == "") local lbl "`v'=1" if ("`lb0'" == "") local lbl "`v'=0" local col1 = "`n1'/`N1' (" + "`pct1'" + "%)" local col0 = "`n0'/`N0' (" + "`pct0'" + "%)" display "" display "----------------------------------------------------------------------" display "{bf:FOR OUTCOME `y'}" display "----------------------------------------------------------------------" display %12s " " /// %25s "`v'=`lb1'" /// %25s "`v'=`lb0'" display "----------------------------------------------------------------------" display %12s "n/N (%)" /// %25s "`col1'" /// %25s "`col0'" display "----------------------------------------------------------------------" display "" } ******************************* * ORIGINAL OVERALL RR/RD OUTPUT ******************************* if ("`subgroup'" == "") { display "***************************" display "Direction of comparison is: {cmd:{it:`v'}}=`lb1' versus {cmd:{it:`v'}}=`lb0'" display "***************************" display "" display /// "Risk ratio (95% CI); [SE] = " /// `fmt_rr' `rr' " (" `fmt_rr' `rr_lb' ", " `fmt_rr' `rr_ub' "); [" `fmt_rr' `log_rr_se' "]" _n /// "P-value for RR = " `fmt_p' `rr_p' _n /// _n /// "Risk difference as percentage (95% CI); [SE] = " /// `fmt_rd_pct' (100*`rd') "% (" `fmt_rd_pct' (100*`rd_lb') "%, " `fmt_rd_pct' (100*`rd_ub') "%); [" /// `fmt_rd_pct' (100*`rd_se') "]" _n /// "Risk difference as proportion (95% CI); [SE] = " /// `fmt_rd_prop' `rd' " (" `fmt_rd_prop' `rd_lb' ", " `fmt_rd_prop' `rd_ub' "); [" /// `fmt_rd_prop' `rd_se' "]" _n /// "P-value for RD = " `fmt_p' `rd_p' capture estimates restore MSRRD_MODEL_ORIG } if ("`subgroup'" == "") & `subtype_specified' { display as error "ERROR: subgroup option not selected." capture estimates restore MSRRD_MODEL_ORIG } ************************ * RETURN OVERALL RESULTS ************************ return scalar rr = `rr' return scalar rr_lb = `rr_lb' return scalar rr_ub = `rr_ub' return scalar rr_z = `rr_z' return scalar rr_p = `rr_p' return scalar rr_se = `rr_se' return scalar rd = `rd' return scalar rd_lb = `rd_lb' return scalar rd_ub = `rd_ub' return scalar rd_z = `rd_z' return scalar rd_p = `rd_p' return scalar rd_se = `rd_se' if ("`subgroup'" == "") exit ************************************************************************************************************************************************************* ********************************************************************** * SUBGROUP ANALYSIS - This section of the code is for subgroup results ********************************************************************** quietly levelsof `subgroup', local(levels) if ("`subgroup'" != "") & ("`vceuncond'" != "") { display in yellow "" display "Running marginal standardisation with {cmd:{it:vce(unconditional)}} to estimate SEs allowing for sampling of covariates" } if ("`subgroup'" != "") & ("`vceuncond'" == "") { display in yellow "" display "Running marginal standardisation with {cmd:{it:vce(delta)}} estimate SEs using delta method" } display in red _n "SUBGROUP ANALYSIS BY `subgroup' - {bf:For outcome `y'}" **************************************************** * RESTORE ORIGINAL MODEL BEFORE ANY POST-ESTIMATION **************************************************** capture estimates restore MSRRD_MODEL_ORIG ******************************************************************* * Interaction test - This will display the Interaction test p-value ******************************************************************* capture testparm i.`v'#i.`subgroup' if (_rc == 0) { display in red "Interaction test p-value for {cmd:{it:`v'}} by {cmd:{it:`subgroup':}}" display in red "Chi2(" r(df) ") = " %9.3f r(chi2) ", P = " `fmt_p' r(p) display in yellow "" } ****************************************************************************************************************************** * STRICT INTERACTION CHECK - This check will ensure the command will not work if interaction term not included in fitted model ****************************************************************************************************************************** local has_interaction = 0 foreach g2 of local levels { capture local test = _b[1.`v'#`g2'.`subgroup'] if (!_rc) local has_interaction = 1 } if (`has_interaction' == 0) { display as error "ERROR: The fitted model does NOT contain the interaction i.`v'#i.`subgroup' / i.`v'##i.`subgroup'." exit 198 } ****************** * SUBGROUP MARGINS ****************** capture estimates restore MSRRD_MODEL_ORIG * --- SUBGROUP MARGINS VCE CONTROL --- if ("`vceuncond'" != "") { quietly margins `v'#`subgroup', post vce(unconditional) } if ("`vceuncond'" == "") { quietly margins `v'#`subgroup', post } matrix Mg = r(b) matrix VMg = r(V) *************************** * LOOP OVER SUBGROUP LEVELS *************************** foreach g of local levels { display "______________________________________________________________________" display "" display "{bf:`subgroup' = `g' (`: label (`subgroup') `g'')}" local p0_g = Mg[1,"0.`v'#`g'.`subgroup'"] local p1_g = Mg[1,"1.`v'#`g'.`subgroup'"] local var_p0_g = VMg["0.`v'#`g'.`subgroup'","0.`v'#`g'.`subgroup'"] local var_p1_g = VMg["1.`v'#`g'.`subgroup'","1.`v'#`g'.`subgroup'"] local cov_p1p0_g = VMg["1.`v'#`g'.`subgroup'","0.`v'#`g'.`subgroup'"] *********************************************************************** * RR esimation using marginal standardisation and SE using delta method *********************************************************************** local rr_g = `p1_g'/`p0_g' local logrr_g = ln(`rr_g') local varlogrr_g = (`var_p1_g'/(`p1_g'^2)) + (`var_p0_g'/(`p0_g'^2)) - (2*`cov_p1p0_g'/(`p1_g'*`p0_g')) local rr_se_g = sqrt(`varlogrr_g') local rr_lb_g = exp(`logrr_g' - invnorm(0.975)*`rr_se_g') local rr_ub_g = exp(`logrr_g' + invnorm(0.975)*`rr_se_g') local rr_z_g = `logrr_g'/`rr_se_g' local rr_p_g = 2*normal(-abs(`rr_z_g')) local log_rr_se_g = `rr_g'*`rr_se_g' *********************************************************************** * RD esimation using marginal standardisation and SE using delta method *********************************************************************** local rd_g = `p1_g' - `p0_g' local vardiff_g = `var_p1_g' + `var_p0_g' - 2*`cov_p1p0_g' local rd_se_g = sqrt(`vardiff_g') local rd_lb_g = `rd_g' - invnorm(0.975)*`rd_se_g' local rd_ub_g = `rd_g' + invnorm(0.975)*`rd_se_g' local rd_z_g = `rd_g'/`rd_se_g' local rd_p_g = 2*normal(-abs(`rd_z_g')) *********************************************** * SUBGROUP EVENT COUNTS (with pctdp formatting) *********************************************** quietly count if `subgroup'==`g' & `v'==0 & `y'!=. local N0_g = r(N) quietly count if `subgroup'==`g' & `v'==0 & `y'==1 local n0_g = r(N) quietly count if `subgroup'==`g' & `v'==1 & `y'!=. local N1_g = r(N) quietly count if `subgroup'==`g' & `v'==1 & `y'==1 local n1_g = r(N) local pct0_g_raw = 100*`n0_g'/`N0_g' local pct1_g_raw = 100*`n1_g'/`N1_g' local pct0_g : display `fmt_pct' `pct0_g_raw' local pct1_g : display `fmt_pct' `pct1_g_raw' ************************* * SUBGROUP COMBINED TABLE ************************* * Extract value labels for variable local lb1 : label (`v') 1 local lb0 : label (`v') 0 * Fallback if no labels exist if ("`lbl'" == "") local lbl "`v'=1" if ("`lb0'" == "") local lbl "`v'=0" local col1_g = "`n1_g'/`N1_g' (" + "`pct1_g'" + "%)" local col0_g = "`n0_g'/`N0_g' (" + "`pct0_g'" + "%)" display "----------------------------------------------------------------------" display %12s " " /// %25s "`v'=`lb1'" /// %25s "`v'=`lb0'" display "----------------------------------------------------------------------" display %12s "n/N (%)" /// %25s "`col1_g'" /// %25s "`col0_g'" display "----------------------------------------------------------------------" display "" ******************************** * ORIGINAL SUBGROUP RR/RD OUTPUT ******************************** if ("`subtype'" == "rr" | "`subtype'" == "both") { display "***************************" display "Direction of comparison is: {cmd:{it:`v'}}=`lb1' versus {cmd:{it:`v'}}=`lb0'" display "***************************" display "" display /// "Subgroup RR (95% CI); [SE] = " /// `fmt_rr' `rr_g' " (" `fmt_rr' `rr_lb_g' ", " `fmt_rr' `rr_ub_g' "); [" /// `fmt_rr' `log_rr_se_g' "]" *display "RR P-value for (`: label (`subgroup') `g'') subgroup = " `fmt_p' `rr_p_g' display "" } if ("`subtype'" == "rd" | "`subtype'" == "both") { display /// "Subgroup RD as percentage (95% CI); [SE] = " /// `fmt_rd_pct' (100*`rd_g') "% (" `fmt_rd_pct' (100*`rd_lb_g') "%, " /// `fmt_rd_pct' (100*`rd_ub_g') "%); [" `fmt_rd_pct' (100*`rd_se_g') "]" display /// "Subgroup RD as proportion (95% CI); [SE] = " /// `fmt_rd_prop' `rd_g' " (" `fmt_rd_prop' `rd_lb_g' ", " /// `fmt_rd_prop' `rd_ub_g' "); [" `fmt_rd_prop' `rd_se_g' "]" *display "RD P-value for (`: label (`subgroup') `g'') subgroup = " `fmt_p' `rd_p_g' display "" display "______________________________________________________________________" } } capture estimates restore MSRRD_MODEL_ORIG **************** * END OF PROGRAM **************** end