sub group analysis in clinical trial and forest plot in stata

[Winfred Mwangi]

Kindly help with the command for a forest plot

Am analysis RCT data
A binary outcome with several baseline factors defining subgroups and I am using interaction terms.

[Carlo Lazzaro]

Winfred:
it is impossible to reply positively with such scant details (plese, read and act on the FAQ on how to post more effectively. Thanks).
As far as forest plot is concerned, see -help meta forestplot-.

[Winfred Mwangi]

Thank you.

My data looks like this:
----------------------- copy starting from the next line -----------------------

* Example generated by -dataex-. To install: ssc install dataex
clear
input float(gap6_binary arm) byte Age_protocol float(FacilityxVL timeDxHIV_cat Pregnancyoutcome_eid_5 GapinCare)
1 1 2 4 1 1 1
1 0 2 2 3 1 0
1 1 1 5 2 1 0
0 0 2 2 2 1 1
0 1 2 5 2 0 1
0 0 1 3 2 1 0
0 0 2 3 2 1 0
0 0 2 2 3 1 0
1 1 2 3 3 1 0
1 1 2 2 3 1 0
end
label values gap6_binary gap6_binary_label
label def gap6_binary_label 0 "No", modify
label def gap6_binary_label 1 "Yes", modify
label values arm arm_label
label def arm_label 0 "Control", modify
label def arm_label 1 "Intervention", modify
label values Age_protocol Age_protocol_label
label def Age_protocol_label 1 "18-29yrs", modify
label def Age_protocol_label 2 "30-39yrs", modify
label values FacilityxVL FacilityxVL_label
label def FacilityxVL_label 2 "MTRH", modify
label def FacilityxVL_label 3 "Kitale", modify
label def FacilityxVL_label 4 "UGDH", modify
label def FacilityxVL_label 5 "Huruma", modify
label values timeDxHIV_cat timeDxHIV_cat_label
label def timeDxHIV_cat_label 1 "<1 year", modify
label def timeDxHIV_cat_label 2 "1-5 years", modify
label def timeDxHIV_cat_label 3 ">5 years", modify
label values Pregnancyoutcome_eid_5 Pregnancyoutcome_eid_5_label
label def Pregnancyoutcome_eid_5_label 0 "Pregnancy loss/Unknown", modify
label def Pregnancyoutcome_eid_5_label 1 "Livebirth", modify
label values GapinCare GapinCare_label
label def GapinCare_label 0 "NO", modify
label def GapinCare_label 1 "Yes", modify

------------------ copy up to and including the previous line ------------------

[Winfred Mwangi]

DV=gap6_binary and is binary
Treatment=arm and is binary 0=control 1=intervention
Objective is to determine effect of Treatment on ga6_binary

other IVs= timeDxHIV_cat FacilityxVL education_level Marital_status2 Pregnancyoutcome_eid_5 GapinCare phone_own2 hours_cat etc
Though randomization was effective. We are adjusting for 3 categorical variables: Age_protocol (maternal age -3 levels) ///
FacilityxVL (5 different clinic sites) timeDxHIV_cat (times since diagnosis- 3 levels)
The effect of treatment with or without controlling for above variables is negative
So I would like to do subgroup analysis (was pre-defined)
subgroups: Inaddition to controlling variables am using above I have: education_level, Marital_status2, ///
Pregnancyoutcome_eid_5 hours_cat, -all categorical with 3 levels
plus phone_own2 and GapinCare binary variables with 2 levels

[Winfred Mwangi]

Am using this command:

eststo:poisson gap6_binary ib(first).(arm FacilityxVL Age_protocol timeDxHIV_cat education_level Gest_Age_Enrol_cat Pregnancyoutcome_eid_5 ///
GapinCare Marital_status2 hours_cat years_res3 phone_own2 ) arm##(arm FacilityxVL Age_protocol timeDxHIV_cat Pregnancyoutcome_eid_5 ///
GapinCare Marital_status2 hours_cat education_level years_res3 phone_own2 ) , baselevels vce (robust) irr
esttab using "G:\Manuscript one_delivery VL\Results_14092021\Table 1\VL test2b.rtf",replace wide ci label pr2 b(%9.3f) eform

[Winfred Mwangi]

Thank you for your help. Sorry if these are too long

[Winfred Mwangi]

I would like the output to come out as a forest plot

[Weiwen Ng]

When people say forest plot, I think they're usually referring to meta-analysis. That is, you have the results from several studies. You use meta-analysis to pool the average effects of the studies and you get a forest plot. You get one marker for each study showing the average treatment effect and the standard error around that average. You also get an overall marker for
the pooled average (and its SE).

I'm not sure if forest plot is used a bit more widely than that. You can plot the coefficients in a forest plot like presentation, although you don't get an overall ATE because you're just plotting the coefficients from one study. You can use the add-on command coefplot for that (written by Ben Jann, available on SSC). I am not sure if there's a standard term for this. Ben Jann just says "coefficient plot." I have also heard "dot and whisker plot."

You can also use margins and marginsplot. If you do that, the plots are no longer on the same scale as your coefficients (or on the incidence rate ratio scale), they're on the same scale as the raw counts. (By default, at least.)

So, what is it you want to plot? The effects of all the coefficients? Do you want to plot on the IRR scale, or some other scale?

[Winfred Mwangi]

Thank you for your response Weiwen Ng. I was learning this from youtube and examples of forest plots are with meta-analysis.

I guess I meant coefplot. What is recommended or advisable in this case? All I want is a graphical presentation of the interaction terms not in separate models but in a multivariable model.

I would like to show the exponentiated coefficients and the absolute numbers as well if this is possible.

Thank you!

[Andrew Musau]

The following uses coefplot from SSC, displaying the incidence-rate ratios and reporting the coefficients on the right-hand side margin.

capt program drop coefplot_mlbl2
*! version 1.0.0  10jun2021  Ben Jann
program coefplot_mlbl2, sclass
    _parse comma plots 0 : 0
    syntax [, MLabel(passthru) * ]
    if `"`mlabel'"'=="" local mlabel mlabel(string(@b))
    preserve
    qui coefplot `plots', `options' `mlabel' generate replace nodraw
    sreturn clear
    tempvar touse
    qui gen byte `touse' = 0
    local nplots = r(n_plots)
    forv i = 1/`nplots' {
        qui replace `touse' = __plot==`i' & __at<.
        mata: st_global("s(mlbl`i')", ///
            invtokens((strofreal(st_data(.,"__at","`touse'")) :+ " " :+ ///
            "`" :+ `"""' :+ st_sdata(.,"__mlbl","`touse'") :+ `"""' :+ "'")'))
    }
    sreturn local plots `"`plots'"'
    sreturn local options `"`options'"'
end

webuse dollhill3, clear
poisson deaths smokes i.agecat, exposure(pyears)
lab variable smokes "Smoker"
est sto P
set scheme s1color
coefplot_mlbl2 P, drop(_cons) mlabel(string(@b,"%9.3f"))
coefplot `s(plots)', eform `s(options)' ///
    ymlabel(`s(mlbl1)', angle(0) notick labgap(-10) axis(2) add custom labcolor(navy)) ///
    yaxis(1 2) yscale(alt) yscale(axis(2) alt) xsc(r(. 65)) ///
    ylabel(none, axis(2)) yti("", axis(2)) ///
    note(Note: Incidence-rate ratios displayed. Coefficients reported on the right-hand side margin.)

[Winfred Mwangi]

Andrew Musau Thank you. How do I modify above to show comtrol and intervention side by side?

[Weiwen Ng]

Andrew Musau Thank you. How do I modify above to show comtrol and intervention side by side?

poisson gap6_binary arm##(FacilityxVL Age_protocol timeDxHIV_cat Pregnancyoutcome_eid_5 GapinCare Marital_status2 hours_cat education_level years_res3 phone_own2 , baselevels vce (robust) irr

This is a simplification of the code you intended to run, extracted from post #5. You seem to be interacting the treatment group against all the variables in the model. And actually, in post #4, I now see that you initially said you were running subgroup analyses.

In traditional randomized trials, investigators may decide to run a subgroup analysis. Usually this means subsetting the data, so it's more like (simplifying to make this easier to read):

poisson gap6_binary arm if education == 1, vce(robust) irr
estimates store subgroup_education

In the original specification, you don't really have results for the intervention and control arms per se. The main effect for the interaction arm variable is the incidence rate ratio attributable to the intervention after controlling for everything else. That is, it's the ratio of the adjusted incidence rate in the intervention arm divided by the incidence rate in the control arm.

By interacting arm with the explanatory variables, the interactions will tell you if the effect of the explanatory variables is different in the intervention arm. Is that what you want to do?

If this were an RCT and you were running 3 subgroup analyses, you'd just use estimates store to store the estimates from each regression in memory, then you would tell coefplot to give you a plot of all the 4 analyses simultaneously.

Last, in post #1, you say that you are doing an analysis of a binary outcome in an RCT. A couple of points of confusion here. 1) Why are you using a Poisson model? That is for count data, not binary data. 2) If this is an RCT, why are you controlling for so many variables? It's not traditional to do so, and I think it's not required if the arms were balanced on characteristics at the start and you don't have too much differential attrition. That doesn't mean don't do it, but were you aware that you may not have to?