Custom subgroup means in meta forest

[Jean-Michel Galarneau]

Hello all, I am trying to include custom diamonds for each subgroup in my meta analysis that uses the option subgroup(). I know of the customoverall() option but cannot find one for each group.
Any ideas on how I could achieve this using the meta suite?

Thank you

[David Fisher]

Hi Jean-Michel,

This can certainly be achieved (as can most things! ) with user-written metan. I'm not sure whether you want custom values, or a custom presentation style, so in the example below I use both.

I will use Example 12 of the meta forestplot documentation PDF ("Adding custom columns and overall effect sizes"), in which repeated calls to customoverall() are used to show "custom" diamonds representing predicted (marginal) effects from a meta-regression. The same dataset (a meta-analysis of BCG vaccine efficacy for TB) was used in the 2008 Stata Journal article on metan (although note that metan itself has evolved significantly since then -- I really need to submit an updated article to Stata Journal when I have time!). To run this code you will need the latest version of metan (version 4.08.1 12jul2024) available from the SSC archive.

I hope this is useful!
BW,
David.

use "https://www.stata-press.com/data/r18/bcgset", clear

* Form study subgroups based on latitude
recode latitude (min/23.5 = 1 "Tropical, < 23.5 latitude") (23.5/40 = 2 "23.5-40 latitude") (40/max = 3 "Northern, > 40 latitude"), gen(lat_cat)
assert !missing(latitude)
label var latitude "Latitude"
assert !missing(lat_cat)
label var lat_cat "Latitude region"

* Perform the meta-regression (on latitude as a continuous covariate) and save the desired marginal results in a table (see text of Example 18 for further details)
metareg _meta_es latitude_c, wsse(_meta_se) z
margins, at(latitude_c = (-18.5 -5.5 16.5))
mat marginstable = r(table)

* Now here is the tricky part. We run the meta-analysis, but instead of immediately displaying the forest plot
* we instead replace the data in memory with a "forestplot results set", into which we insert our custom values
preserve
metan _meta_es _meta_se, lcols(studylbl latitude) by(lat_cat) re(reml) rr nograph clear
forvalues i=1/3 {
replace _ES = marginstable["b",`i'] if _USE==3 & _BY==`i'
replace _LCI = marginstable["ll",`i'] if _USE==3 & _BY==`i'
replace _UCI = marginstable["ul",`i'] if _USE==3 & _BY==`i'
}
replace _WT = . if _USE==3
replace _LABELS = "Predicted effect at latitude = 15" if _USE==3 & _BY==1
replace _LABELS = "Predicted effect at latitude = 28" if _USE==3 & _BY==2
replace _LABELS = "Predicted effect at latitude = 50" if _USE==3 & _BY==3
replace _LABELS = "{bf:" + _LABELS + "}" if inlist(_USE, 3, 5)
format %-1s _LABELS

* Finally, we display the forest plot. I have also demonstrated how you can change the look of the diamonds.
forestplot, useopts xlabel(.125 .25 .5 1 2 4) astext(60) plotid(_BY) diam1opts(lcolor(red)) diam2opts(lcolor(green) fcolor(green%50)) diam3opts(lcolor(blue) fcolor(blue))
restore

[Jean-Michel Galarneau]

David, thank you very much for this. I had just about given up.

That is an excellent example.

[Tiago Pereira]

Extremely useful, David! Thank you!

[Jean-Michel Galarneau]

David, is it possible that there is an issue with hlineopts and xtick? Everything else worked beautifully. For the moment I fix it by playing a recording but I was wondering if there was an issue.

In the hlineopts, lwidth throws an error. If i remove that specification at least lpattern works but i cannot edit the line thickness.
xtick just does nothing but does not throw up an error.

My code:

* Finally, we display the forest plot. I have also demonstrated how you can change the look of the diamonds.
forestplot, useopts xlabel(-3(1)3, nogrid) astext(60) plotid(_BY) diam1opts(lcolor(stc1%0) fcolor(stc2)) diam2opts(lcolor(stc2%0) fcolor(stc2)) diam3opts(lcolor(stc2%0) fcolor(stc2)) diam4opts(lcolor(stc2%0) fcolor(stc2)) diam5opts(lcolor(stc3%0) fcolor(stc3)) nobox note("Note: Weights and between-subgroup heterogeneity test are from random-effects model",size(tiny)) effect("Hedge's g") ciopts(lwidth(vthin) lcolor(stc1)) pointopts(mcolor(stc1) lcolor(stc1%0) msymbol(square)) null(0) nlineopts(lpattern(shortdash) lwidth(vthin)) olineopts(lcolor(white%0)) hlineopts(lpattern(solid) lwidth(vthin)) xscale(lwidth(vthin)) xtick(,tlwidth(vthin))

option lwidth() not allowed
r(198);

Stata 18.5
forestplot version
*! version 4.08 David Fisher 17jun2024

[David Fisher]

Hi Jean-Michel,

Thanks for bringing this to my attention. The short answer is: xtick is working as intended, but hlineopts is not -- and I shall fix this as soon as possible.

The longer answers are:

It is perhaps not sufficiently clear in the Stata documentation, but xtick is intended for use with additional, unlabelled ticks. To modify ticks corresponding to labelled values, we use suboptions to xlabel. Hence, instead of:

xlabel(-3(1)3, nogrid) xtick(,tlwidth(vthin))

...you should specify:

xlabel(-3(1)3, nogrid tlwidth(vthin))

The issue with hlineopts is more difficult to explain, but it is a consequence of my use of the yline option to construct the line, which doesn't like "repeated options". That is, I have already "hard-coded" a default value for lwidth, so that when you add your own value, it results in an error. This behaviour does not apply to other lines (e.g. olineopts) because I am instead using twoway function to construct those lines. Hopefully, this will be an easy fix.

Best wishes,

David.

[Jean-Michel Galarneau]

This works, thank you very much.