How do you create a forest plot without meta-analysis?

[Matt Savage]

I'm having a problem deciding how best to present my data (correlations), as heterogeneity precludes meta-analysis. One way I'm considering is using forest plots more descriptively. At the minute my code looks like this:

clear
input str12 study str12 time str12 biomech e n p
Azus 6months vGRF 0.285 46 0.08
Azus 6months vGRF 0.051 46 0.76
ErhartHledik 2years vGRF -0.396 28 0.037
ErhartHledik 2years vGRF -0.584 28 0.009
Pietrosimone 6months vGRF 0.42 25 0.001
Pietrosimone 6months vGRF 0.09 25 0.139
Azus 6months KFM -0.237 46 0.15
Azus 6months KFM -0.187 46 0.26
Azus 6months KFM 0.078 46 0.65
Azus 6months KFM 0.137 46 0.42
ErhartHledik 2years KFM -0.57 16 0.026
ErhartHledik 2years KAM -0.74 16 0.002
ErhartHledik 2years IRmoment -0.82 16 0.001
Azus 6months MedialGRF -0.28 46 0.08
Azus 6months MedialGRF -0.286 46 0.08
Azus 6months MedialGRF -0.28 46 0.09
Azus 6months MedialGRF -0.335 46 0.05
Azus 6months MedialGRF -0.35 46 0.03
Azus 6months MedialGRF -0.342 46 0.04
end

* Calculate standard error
gen SE = sqrt((1 - e^2) / (n - 2))

* Display the data including standard error
list


At this point I get that you would normally declare the data as meta and then can use the forestplot command, but is there a way of doing this without running meta-analysis? I have looked online and I'm finding that most suggestions lead me to using the twoway command but I can't seem to get this to work.
Thanks in advance for any help!!

Matt

[David Fisher]

Hi Matt,
It is very possible to make use of forest plots without running a meta-analysis. You just need to suppress (or ignore) all output relating to the statistical analysis.

Here is some example code using metan (user-contributed, from SSC):

* Example generated by -dataex-. For more info, type help dataex
clear
input str12(study time biomech) float(e n p)
"Azus"         "6months" "vGRF"       .285 46  .08
"Azus"         "6months" "vGRF"       .051 46  .76
"ErhartHledik" "2years"  "vGRF"      -.396 28 .037
"ErhartHledik" "2years"  "vGRF"      -.584 28 .009
"Pietrosimone" "6months" "vGRF"        .42 25 .001
"Pietrosimone" "6months" "vGRF"        .09 25 .139
"Azus"         "6months" "KFM"       -.237 46  .15
"Azus"         "6months" "KFM"       -.187 46  .26
"Azus"         "6months" "KFM"        .078 46  .65
"Azus"         "6months" "KFM"        .137 46  .42
"ErhartHledik" "2years"  "KFM"        -.57 16 .026
"ErhartHledik" "2years"  "KAM"        -.74 16 .002
"ErhartHledik" "2years"  "IRmoment"   -.82 16 .001
"Azus"         "6months" "MedialGRF"  -.28 46  .08
"Azus"         "6months" "MedialGRF" -.286 46  .08
"Azus"         "6months" "MedialGRF"  -.28 46  .09
"Azus"         "6months" "MedialGRF" -.335 46  .05
"Azus"         "6months" "MedialGRF"  -.35 46  .03
"Azus"         "6months" "MedialGRF" -.342 46  .04
end

label variable study "Study"
label variable time "Time"
label variable biomech "Biomech"
label variable n "Sample size"
label variable p "p-value"
format p %05.3f
format study time biomech %-12s

gen SE = sqrt((1 - e^2) / (n - 2))
metan e SE, effect(Correlation) lcols(study time biomech) rcols(n p) nooverall nowt astext(60)

...and here is equivalent code using meta set (official Stata, v16+):

meta set e SE, studylabel(study) studysize(n) eslabel(Correlation)
meta forestplot _id time biomech _plot _esci n p, nooverall noohetstats noohomtest noosigtest nonotes

[Matt Savage]

Hi David,

Thanks for your response - this was a great help!

Matt

[Ashani Abayasekara]

Hi David,

I'm trying to generate a forest plot to visualize effects of heterogeneity analysis. In particular, I'm running fixed-effects regression models, regressing job insecurity on unemployment rates. I'm doing some subsample analyses both by mental health status and by demographics such as gender, age, and education, and am looking to plot effects using forest plots. Am I amble to use this command for my case?

My regression code (for mental health and gender) is as follows:

xtreg job_sec_sat unemp_av6 i.sa4 i.year if entry3==0&sex_avg1==1&mh8_qn1==0, fe i(id) cluster(id)
est store m1
xtreg job_sec_sat unemp_av6 i.sa4 i.year if entry3==0&sex_avg1==1&mh8_qn1==1, fe i(id) cluster(id)
est store m2
xtreg job_sec_sat unemp_av6 i.sa4 i.year if entry3==0&sex_avg1==0&mh8_qn1==0, fe i(id) cluster(id)
est store f1
xtreg job_sec_sat unemp_av6 i.sa4 i.year if entry3==0&sex_avg1==0&mh8_qn1==1, fe i(id) cluster(id)
est store f2

I would like to plot the unemployment effects and associated confidence intervals using forest plots.

Below is a sample of my data:

* Example generated by -dataex-. For more info, type help dataex
clear
input long id byte job_sec_sat float unemp_av6 int sa4 float(year sex_avg1 mh8_qn1)
100003 10  8.533333 405 2001 1 0
100003  9  6.816667 405 2002 1 0
100003  2  7.516666 405 2003 1 0
100003  4  5.366667 405 2004 1 0
100003  2       3.2 405 2005 1 0
100003  . 4.6666665 405 2006 1 0
100003  .  5.383333 405 2008 1 0
100003  8       2.5 405 2009 1 0
100003  .  3.383333 405 2011 1 0
100003  .         4 405 2012 1 0
100003  .  6.016667 405 2013 1 0
100003  .  7.266667 403 2014 1 0
100003  .       7.8 403 2015 1 0
100003  .  6.933333 403 2016 1 0
100003  .       6.8 403 2017 1 0
100003  .  6.416667 403 2018 1 0
100010  3       4.5 207 2002 0 0
100010  2  5.116667 207 2003 0 0
100010  1  3.466667 207 2004 0 0
100010  2 4.2166667 207 2005 0 0
100010  .       4.4 207 2006 0 0
100010  .  3.533333 207 2007 0 0
100010  .      4.75 206 2008 0 0
100010  .  5.733334 206 2009 0 0
100010  1 4.7166667 206 2010 0 0
100010  1  3.516667 206 2011 0 0
100010  5 4.4833336 208 2012 0 0
100010  1  4.766667 208 2013 0 0
100010  1      5.05 208 2014 0 0
100010  3  3.583333 208 2015 0 0
100010  .      3.65 208 2016 0 0
100010  2 4.2833333 208 2017 0 0
100010  4       3.2 208 2018 0 0
100010  5 4.0666666 208 2019 0 0
100014  7  9.333334 106 2001 1 1
100014  7  9.666666 106 2002 1 1
100014  9       5.8 106 2003 1 1
100014  7  5.716667 106 2004 1 1
100014  7  5.183333 106 2005 1 1
100014  8  3.866667 106 2006 1 1
100014  6 3.8166666 106 2007 1 1
100014  7 3.9333334 106 2008 1 1
100014  0       6.8 106 2009 1 1
100014  . 4.7333336 106 2010 1 1
100014  .  3.766667 106 2011 1 1
100014  .  2.866667 106 2012 1 1
100014  .  5.316667 106 2013 1 1
100014  .  7.566667 106 2014 1 1
100014  . 11.416667 106 2015 1 1
100014  .  5.066667 106 2016 1 1
100014  .       6.1 106 2018 1 1
100014  .  5.666667 106 2019 1 1
100015  7  9.333334 106 2001 0 0
100015  2  9.666666 106 2002 0 0
100015  8       5.8 106 2003 0 0
100015  9  5.716667 106 2004 0 0
100015 10  5.183333 106 2005 0 0
100015  9  3.866667 106 2006 0 0
100015  9 3.8166666 106 2007 0 0
100015 10       6.8 106 2009 0 0
100015 10 4.7333336 106 2010 0 0
100015 10  3.766667 106 2011 0 0
100015  9  2.866667 106 2012 0 0
100015  8  5.316667 106 2013 0 0
100015  9 11.416667 106 2015 0 0
100015 10  5.066667 106 2016 0 0
100015  8       4.8 106 2017 0 0
100015  8  6.166667 106 2018 0 0
100015  8  5.383333 106 2019 0 0
100016 10      4.85 111 2007 1 0
100016  9 4.7166667 111 2008 1 0
100016  9 4.7166667 111 2009 1 0
100016  9       4.5 111 2010 1 0
100016  9  4.116667 111 2011 1 0
100016  9      5.85 111 2012 1 0
100016  9  6.333333 111 2013 1 0
100016  9       6.7 111 2014 1 0
100016  9  6.766666 111 2015 1 0
100016  9      5.85 111 2016 1 0
100016 10  5.633333 111 2017 1 0
100016 10  6.033333 111 2018 1 0
100016 10  4.366667 111 2019 1 0
100018  .  7.133333 507 2001 0 1
100018  .  5.516666 507 2002 0 1
100018  .  6.766666 507 2003 0 1
100018  .       5.3 507 2004 0 1
100018  .      4.85 507 2005 0 1
100018  2 2.8166666 507 2006 0 1
100018  2  2.966667 507 2007 0 1
100018  3 3.6666665 507 2008 0 1
100018  4      5.45 507 2009 0 1
100018  5  5.433333 507 2010 0 1
100018  6 4.4333334 507 2011 0 1
100018  8  4.133333 507 2012 0 1
100018  7  5.266666 507 2013 0 1
100018  8  5.516666 507 2014 0 1
100018  8       5.7 507 2015 0 1
100018  7  6.316667 507 2016 0 1
100018  7  6.883333 507 2018 0 1
100018  9  7.016666 507 2019 0 1
end
label values job_sec_sat AJBMS
label def AJBMS 0 "[0] Totally dissatisfied", modify
label def AJBMS 10 "[10] Totally satisfied", modify

Can you please give me any guidance on whether this is doable and if so how I should proceed?

Thank you very much,
Ashani

[David Fisher]

Hi Ashani,

It looks like you have an "individual-level" dataset (as opposed to summary data). As such, I could recommend two possible approaches: (1) ipdover / forestplot; or (2) coefplot . Both commands are user-contributed, and both are available from the SSC archive.

In keeping with the overall topic of this thread, I'll start with the forestplot -based solution. You will need to install two packages from SSC: one named metan (which contains forestplot), and one named ipdmetan (which contains ipdover).


The ipdover command is designed to produce forest plot -style solutions for subgroup analyses of a single study (as opposed to meta-analysis). Here is some example code, based on your example data above. (Note that I've replaced your xtreg commands with simpler regressions, because there is insufficient data in your example to fit the full models.)

. label define SEX_AVG1 0 "sex_avg1=0" 1 "sex_avg1=1"
. label values sex_avg1 SEX_AVG1
. label define MH8_QN1 0 "mh8_qn1=0" 1 "mh8_qn1=1"
. label values mh8_qn1 MH8_QN1

. ipdover, over(sex_avg1) over(mh8_qn1) : regress job_sec_sat unemp_av6

. egen sex_mh8 = group(sex_avg1 mh8_qn1), label(SEX_MH8)
. ipdover, over(sex_mh8) forestplot(xtitle("Effect of unemp_av6 on job_sec_sat")) : regress job_sec_sat unemp_av6

If you wish to edit the forest plot, you can use the "clear" option to replace the data in memory with the data underlying the plot. If you alter anything in this new dataset, it will be reflected in the plot. For example:

. preserve
. ipdover, over(sex_mh8) clear nograph forestplot(xtitle("Effect of unemp_av6 on job_sec_sat")) : regress job_sec_sat unemp_av6
. replace _LABELS = "foo" in 1
. forestplot, useopts
. restore

coefplot is designed as an all-purpose package for plotting regression results, and is very flexible but has quite complicated syntax. Here is some example code:

. regress job_sec_sat unemp_av6 if sex_avg1==1 & mh8_qn1==0
. est store m1
. regress job_sec_sat unemp_av6 if sex_avg1==1 & mh8_qn1==1
. est store m2
. regress job_sec_sat unemp_av6 if sex_avg1==0 & mh8_qn1==0.
. est store f1
. regress job_sec_sat unemp_av6 if sex_avg1==0 & mh8_qn1==1
. est store f2

. coefplot (m1 m2 f1 f2), drop(_cons) asequation swapnames legend(off) ///
    eqrename(m1="sex_avg1==1 & mh8_qn1==0" m2="sex_avg1==1 & mh8_qn1==1" f1="sex_avg1==0 & mh8_qn1==0" f2="sex_avg1==0 & mh8_qn1==1") ///
    xtitle("Effect of unemp_av6 on job_sec_sat") xline(0)

Documentation for any of these commands is available through Stata once the packages have been installed. I myself am the maintainer of ipdover and forestplot, whilst coefplot is maintained by Ben Jann.
I hope this is helpful.
Best wishes,
David.

[Ashani Abayasekara]

Hi David,

Thank you very much for the detailed explanation and code which is very much appreciated. I was wondering whether there is a way to only display the effect sizes for each group without all the other information in the forest plot? Also, similar to gender in this case, I also have subgroups for age and education; can all of these be combined in one graph? Relatedly, I tried the code you provided for age on its own (similar to gender), but only two sets of effect sizes were produced (i.e. not differentiated by mh8_qn1).

I tried the coefplot command as well which seemed more straightforward to me, but ended up with the attached graph which looks weird.

Thanks again for your help.

Best,
Ashani