How to visualize the prediction interval in a forest plot with nooverall ((as a comparison to confidence intervals)

[Manuela Neuenschwander]

Hello everybody,

I'm conducting an umbrella meta-analysis. For each comparison of interest I calculated the summary RR with the corresponding 95%-CI and the prediction interval. Now I'd like to visualize these results in a forest plot. Each row represents a comparison with the summary RR and 95%-CI. So far that is no problem.

But I'd also like to draw a second horizontal line that represents the prediction interval for each comparison . I already calculated the upper and lower limit of the PI. I only want to visualize it as a second line through each square as a comparison to the confidence interval.

I calculated the PIs with the command rfdist. But with that command it only shows the PI as a horizontal line through the diamond. But I used the command nooverall because it doesn't make sense to summarize the data and I also want to have the PI visualized for each comparison indiviually.

Does anybody know how I can do that? Is there a command?

Thank you in advance for your help!

[David Fisher]

Hi Manuela,

I think I may be able to help. I've developed a meta-analysis package called ipdmetan, which (amongst other things) allows you to control exactly what data is plotted in the forest plot, and how it is plotted. Unfortunately I never considered this particular scenario, so this solution will be a bit clunky. I'll give this some thought for the next version of the package.

Consider the following example dataset, where eff, se_eff, lci and uci are the effect sizes, standard errors and confidence limits; and predlci, preduci are the limits of the prediction intervals:

* Example generated by -dataex-. To install: ssc install dataex
clear
input double(eff se_eff lci uci predlci preduci) int npts
  .456 .138  .185  .727  .016   .896  325
  .489 .167  .163  .816 -.318 1.297  331
 -.198 .118  -.428 .033  -.572  .176  487
 -.051 .114  -.274 .173  -.414  .312  499
end

Now: currently, the command forestplot (within the ipdmetan package; type ssc install ipdmetan to install) expects prediction intervals to have been calculated by the meta-analysis routine, rather than by the user. Consequently, it expects the prediction interval data to be presented in a particular way. We will have to emulate this to get it to display the intervals.

We need to (1) make a copy of each observation; (2) in our copies, replace the effect size confidence limits with the prediction interval limits and remove the effect sizes themselves; (2) re-order the observations; (3) generate a variable, "use", that tells forestplot what sort of observations these are (specifically, pooled observations to be presented as diamonds, as opposed to individual studies to be presented as squares). Then, finally, we can run forestplot:

expand 2, gen(extra)
replace eff = . if extra
replace lci = predlci if extra
replace uci = preduci if extra

gen order = _n
bysort extra (order) : gen neworder = _n
sort neworder extra

gen labels = "Meta-analysis " + string(neworder)
replace labels = "with estimated prediction interval" if extra
label variable labels "Meta-analysis name"

gen use = 3
forestplot eff lci uci, use(use) labels(labels) rfdist(predlci preduci) nowt leftjustify

I appreciate that this is not ideal, but hopefully it is of some use. Please let me know if you have any questions or comments about this, or about the ipdmetan package in general.

Best wishes,

David.

[Manuela Neuenschwander]

Dear David,

Thank you very much for your answer. Unfortunately, when I execute the command as described above, I get the error: Error in prediction interval. it works when I use your data, but with mine it doesn't work? I attached my table, do you see the difference / mistake?

And I'd also like to have only one row per comparison and the PI as a second line, rather than diamond (I also attached an example). Is that possible?

Thank you very much!

Best wishes,

Manuela
Data dietary pattern.xlsx Forest plot with PI.pdf

[David Fisher]

Dear Manuela,

Thank you for sending me your data. I think the error is due to Stata reading in additional, empty rows from the spreadsheet. Therefore, try using the following code which makes use of the cellrange() option to import excel (in Stata 12 and above):

import excel using "Data dietary pattern.xlsx", clear firstrow cellrange(A1:O9)

expand 2, gen(extra)
replace eff = . if extra
replace lci = predlci if extra
replace uci = preduci if extra

gen order = _n
bysort extra (order) : gen neworder = _n
sort neworder extra

gen labels = Exposure
replace labels = "with estimated prediction interval" if extra
label variable labels "Exposure"

gen use = 3
forestplot eff lci uci, use(use) labels(labels) rfdist(predlci preduci) nowt leftjustify

Unfortunately, there is currently no way to display the PI as a second line as in your attachment. However, in the next update (coming soon hopefully!) it will be possible to replace the diamond with a box and line representing the standard confidence interval, on which the PI line may be overlaid. Then, you could format the two lines differently (colours, dash/solid, use of "caps") to make each stand out. Attached is an example, using your data. Does this help at all?

Thanks,

David.

[Manuela Neuenschwander]

Hello David, thank you very much. It works in general. But if I want to adjust the textsize (using texts()) and the xlabel (using force xlabel()), which normally work, I also get Errors. But I'm going to have to wait for the update anyway . I do like your attached example. That is also a nice display of the data and a good alternative to the second line... Do you have an idea when it will be available?

Thank you very much for your effort!

Manuela

[David Fisher]

Dear Manuela,

The options texts() and force xlabel() are metan options, not forestplot options. With forestplot, the text size is determined by the shape of the plot, and may be indirectly changed by using options such as astext() and spacing() -- see help forestplot.

The force option, meanwhile, should go inside the xlabel() option; e.g. xlabel(-2 0 2, force).

I'm hoping to get the update out in the next few weeks, but it depends on results of testing, and on my other work commitments. I'll let you know!

Thanks,

David.

[Manuela Neuenschwander]

Dear David,

Thank you very much for your assistance!

I'm looking forward to hearing from you.

Manuela

[David Fisher]

Dear Manuela,

The update is now available!

In addition to the plot I showed you before, below is another presentational approach which is now possible.

Thanks,

David.

import excel using "Data dietary pattern.xlsx", clear firstrow cellrange(A1:O9)
gen byte use = 3
forestplot eff lci uci, use(use) labels(Exposure) rfdist(predlci preduci) nowt leftjustify nodiam ppointopts(msize(small)) rfopts(lcolor(gs12) lwidth(vthick) overlay)

[Manuela Neuenschwander]

That is good to know! Thank you very much!

[Tom Hsiung]

Is there any option of -metan- to show the prediction interval?

Tom

[Nobuatsu Aoki]

Hello, my name is Nobuatsu Aoki. I would like to express prediction interval in the same way in .metan. If possible, I want to represent it like the Forest plot in the attached Figure. I would be grateful if someone could lend us your wisdom.

Warmest regard,
Nobu

[David Fisher]

Hello Nobuatsu Aoki.

This is a lovely example, thankyou for sharing it. Although it is entirely possible to recreate this plot (as we shall see), it has also highlighted a couple of bugs and issues that I will try and resolve for the next version of metan.
(Important note: the current version is v4.06 12oct2022, available from SSC)

First, some basic strategy:
-- although metan can fit multiple models at once (e.g. fixed and random), it currently only saves one set of weights. Therefore we need to run metan twice, and save the weights from the first run;
-- metan allows you to save a dataset in a format which makes it easy for forestplot to create the plot. Furthermore, this dataset may be edited, giving us great control over the look of the plot.

And now, some suggested code (some variations might also work):

Load the data into memory

clear
input str27 Study int(Events Total)
"Afonso, 2013 (RJ)"              8   60
"Augusto, 2014 (RJ)"            15  351
"Autran, 2018 (CE)"             10   20
"Becker, 2001 (RS)"            101  867
"Boon, 2001 (RS)"                1   44
"Brito, 2002 (AM)"               3   49
"Caixeta, 2015 (GO)"            99  251
"Campos, 2012 (GO)"             19  173
"Campos, 2014 (MS)"             22  170
"Carestiato, 2006 (RJ)"       1488 5833
"Ceccato Junior, 2016 (MG)"     33  160
"Coser, 2013 (RS)"              83  337
"Costa-Lira, 2017 (AM)"         29  133
"da Silva, 2012 (PR)"           15  418
"de Abreu, 2016 (PR)"           66  614
"de Brot, 2017 (SP)"             6   48
"de Mattos, 2011 (ES)"          44   77
"dos Santos, 2013 (RO)"         74  334
"Fedrizzi, 2008 (SC)"           15  100
"Fiqueiredo Alves, 2013 (GO)"  107  432
"Fonseca, 2015 (AM)"           173  607
"Franco, 1995 (PB)"             61  525
"Girianelli, 2010 (RJ)"        263 2056
"Golfetto, 2018 (SC)"           17  325
"Krambeck, 2008 (SC)"            1   57
"Lima, 2014 (SP)"               41  100
"Lippman, 2010 (SP)"            78  386
"Lopes, 2001 (SP)"              34  226
"Lorenzato, 2000 (PE)"         202  479
"Lugo, 2018 (MS)"               11   67
"Oliveira, 2010 (RJ)"           54  241
"Paesi, 2015 (RS)"              77  225
"Rama, 2010 (SP)"              133  301
"Rodrigues, 2014"               24   84
"Rodrigues, 2018 (AM)"          32  112
"Roteli-Martins, 2011"         112 1509
"Santos-Filho, 2016 (AL)"       51  515
"Silva, 2009 (RJ)"              46  150
"Simões, 2017 (GO)"            20  100
"Smith, 2002 (SP)"              17  173
"Ströher, 2016 (RS)"            1   51
"Tamegão-Lopes, 2014 (PA)"     19  143
"Teixeira, 2016 (RS)"           47  182
"Vieira, 2015 (PA)"             26  265
end

Run the initial models

* Fixed effects model
metan Events Total, study(Study) proportion transform(logit) denom(100) ci(exact) nogr
rename _WT WeightFixed

* Random effects model
metan Events Total, study(Study) proportion transform(logit) denom(100) ci(exact) nogr ///
  model(re \ iv) /* also print the fixed effect result */ ///
  rfdist              /* prediction interval */ ///
  lcols(Events Total) /* save the Events and Total columns within the new dataset */ ///
  rcols(WeightFixed) /* also save the fixed-effects weights that we generated earlier */ ///
  extraline(yes)   /* minor option which decreases white space on the left-hand side */ ///
  hetinfo(isq tausq pvalue) isqparam  /* save heterogeneity statistics (this doesn't quite work properly, I've noticed */ ///
  clear  /* clear the current dataset and replace with the "forestplot" dataset */

Various bits of editing and tidying -- difficult to explain, best to run the code and see for yourself

gen order = _n
sort _USE order
drop in L
replace _USE = 3 in L
replace _rfLCI = . in L
replace _rfUCI = . in L
replace WeightFixed = 100 in L
replace _USE  = 7 if _USE==4
sort _USE order
replace _USE = 4 if _USE==7
replace _LABELS = subinstr(_LABELS, ")", ", p < 0.001)", 1) if _USE==4

replace _LABELS = "{bf:Fixed effect model}" if _USE==3
replace _LABELS = "{bf:Random effects model}" if _USE==5 & !missing(_ES)
replace _LABELS = "{bf:Prediction interval}" if _USE==5 & missing(_ES)
format %-1s _LABELS
replace _USE = 5 if _USE==3

label variable WeightFixed "Weight (fixed)"
label variable _WT "Weight (random)"
rename (_WT WeightFixed) (WeightRand _WT)

Finally, construct the forest plot
(note: useopts is a useful option that carries forwards plot-related options from the previous call to metan)

forestplot _Prop_ES _Prop_LCI _Prop_UCI, useopts effect(Events per 100 observations) rcols(WeightRand) ///
   rfopts(sepline lcolor(maroon)) /* Change the rendering of the prediction interval line*/ ///
  xlabel(0(20)80) range(-20 80) astext(70) boxsca(50) /* x-axis labelling and general aesthetic options */

I hope this is helpful, both to yourself and to other readers. I am aware you have also contacted me via e-mail, and I will respond separately.

Best wishes,

David.

[David Fisher]

...and here is the resulting plot:

[Nobuatsu Aoki]

Dear David,

Amazing! Thank you for your perfect response!
I was able to create the same figure with the data and commands you provided.

I have one more request: is it possible to change the box display of each ES in the Forest plot to the size of the Random effect model?
I feel uncomfortable that the Forest plot is a Fixed ffect model while the Prediction interval is displayed.
Thank you very much for your time.

Sincerely regards,
Nobu

[Masum Ali]

Hi David,
I am conducting a meta-analysis and attempting to create a forest plot similar to the one you provided for Nobuatsu Aoki. However, I am encountering issues, including insufficient data for the fixed-effect model and the absence of a prediction interval.

Furthermore, when I tried to generate a forest plot using the data and commands you provided for Nobuatsu Aoki, I am faced with insufficient data for the fixed-effect model (similar problem). I am a STATA 18 user. Could you please assist me in resolving these problems?. I am sharing forest plot that i created using Nobuatsu Aoki's data. Thank you in advance.

Best,
Masum