Forest Plots to display interaction

[Leo Martinez]

Hello! I have a couple of technical questions! Thanks in advance for any help anyone can provide! I am using Stata 12.1.

1) How can forest plots be created to assess interaction in common studies? I have seen this done a couple of times. Is this done with the same commands as a regular meta-analyses (metan command or a couple other options)?

2) I want to display relative risks that adjust for household clustering in my forest plot. I am unsure how to do this -- is there a way that the forest plot can plot exactly what I list in my data sheet and not necessarily automatically compute it's own crude relative risk/odds ratio and then automatically plot it? Since I am adjusting for household correlation and the interaction p-values I have computed are based on controlling for clustering these crude values that metan program automatically computes wouldn't make much sense. Hopefully this makes sense!

Thanks again!

Leo

[David Fisher]

Hi Leo,

Sorry for the late response, I must have missed this.

1) Could you give an example, e.g. a screenshot/image or a fuller description? e.g. what data do you have (IPD? Subgroups?), and what how would you like it displayed?

2) If you pass an effect size and standard error (rather than crude counts) to metan, it will honour them in the plot. Similarly, you can pass an effect size and confidence limits, although there are subtleties here: in order to compute the overall pooled estimate, metan needs to convert the confidence limits into a standard error and hence an inverse-variance weight. Therefore the confidence limits must be derived from a normal distribution (i.e. effect size +/- 1.96 * se), and must be symmetric about the effect size (i.e. must be given on the linear scale, not the exponential. You can, of course, re-exponentiate using the eform option).

I wouldn't like to comment about "the interaction p-values I have computed" without further information (see my answer to (1) ). But interactions often require individual participant data (IPD), and if you have access to such data, you may find my SSC package ipdmetan useful.

Best wishes,

David.

[Leo Martinez]

Thanks David and definitely no problem. You've helped me much in the past -- thanks for all your work.
After re-reading my question I realized it may not be completely clear what I am asking. Basically I want to create a table/figure such as the one below:


Basically they are showing a forest plot from a single study to show interaction between their different stratified groups.
You can find the figure in the article: Pilgrim, Thomas, et al. "Ultrathin strut biodegradable polymer sirolimus-eluting stent versus durable polymer everolimus-eluting stent for percutaneous coronary revascularisation (BIOSCIENCE): a randomised, single-blind, non-inferiority trial." The Lancet 384.9960 (2014): 2111-2122.

The subject is different than mine but I am showing it just to show an example of what I'm trying to do. Another example can be seen below:



Same idea but with less detail. If I have a choice I want to make the first figure but both are conveying the idea I want to show. Here is the citation that has this second figure:
McWilliams, J. Michael, et al. "Performance Differences in Year 1 of Pioneer Accountable Care Organizations." New England Journal of Medicine (2015).

Currently, I have almost the exact table/figure (different topic as shown) as the first example except for the forest plot type image showing the different relative risks. I have attached my table below so you can view it (the whole table is not shown because it is large and covers multiple pages -- I couldn't fit it all within the image and is not relevant to my question here):

As you can see my table is similar to the first figure except for the forest plot. A couple of snags that I mentioned in the first post. The stratum-specific RR is a modified poisson regression model with robust standard error variance and that adjusts for household clustering (Zou, Guangyong. "A modified poisson regression approach to prospective studies with binary data." American journal of epidemiology 159.7 (2004): 702-706.). The p-value for interaction is calculated from these stratum-specific relative risks. The RR's and first p-value don't adjust for anything other than household clustering. So just to clarify if I make a forest plot with the specific prevalence numbers in the columns (as is usually done with forest plots from meta-analyses) I would get a different RR and p-value than shown to the right -- therefore the forest plot would likely be inaccurate or convey something different at a minimum. Hopefully that makes sense.

So I have two basic questions based on my specific data:

1) Is it possible to remake this table putting the forest plot within the table (using the shown RR's instead of using prevalence numbers) as seen in the first example above in the Lancet?

2) If I can't put it within the table, is it possible to make a forest plot alone using the relative risks calculated with the modified poisson model that adjusts for household clustering instead of using the prevalence numbers to automatically calculate effect measurements as seen in metan and other meta-analysis type programs? This would be more like the second example (Mcwilliams et al, 2015).

Thank you for directing me to your package ipdmetan. This is actually the first time I have found or heard of this package and am happy to have seen it now. This could help me with another project on creating a meta-analysis forest plot using adjusted effect measures. Great news! Thanks for letting me know.
I will take a closer look at this package now and see if it can help with the current question. I need to play around with it.

Thanks so much -- hopefully I have clarified the question with the above explanation. If not, please let me know.
If the ipdmetan package can accomplish this I will find out and if so, I will write back. Thanks again.

Leo

[David Fisher]

Hi Leo,

This type of plot is precisely what my program "ipdover" (part of the "ipdmetan" package) was designed for!
You can specify your adjusted, clustered regression model to ipdover to obtain RRs within each subgroup, and then use the lcols/rcols options to create the columns of data, e.g. "n/N (%)" and the p-values for interaction.

The "n/N (%)" columns are a bit fiddly to create. The easiest way is to use the collapse-like syntax of lcols on the numerators and denominators; save the dataset using the saving() option, load it up, create an appropriate string from the pieces, and then run forestplot with the lcols option on your newly created variable. Something like:

webuse lbw, clear
recode age (min/24=0 "<25") (25/max=1 "25+"), gen(agegp) label(agegp)
label var agegp "Mother's age (grouped)"
gen low1 = low
gen low0 = 1 - low
gen smoke0 = smoke * low0
gen smoke1 = smoke * low1
ipdover, over(race ht agegp) lcols((sum) smoke0 low0 smoke1 low1) nowt or saving(test.dta) : logistic low smoke

use test, clear
gen prev0 = string(smoke0) + "/" + string(low0) + " (" + string(100*smoke0/low0, "%4.1f") + ")" if !missing(smoke0)
gen prev1 = string(smoke1) + "/" + string(low1) + " (" + string(100*smoke1/low1, "%4.1f") + ")" if !missing(smoke1)
label var prev0 "Non-LBW infants"
label var prev1 "LBW infants"
forestplot, or nowt lcols(prev0 prev1) favours("Odds of LBW decrease" "with maternal smoking" # "Odds of LBW increase" "with maternal smoking")

Similarly with the p-values for interaction. At present, ipdover does not perform interaction testing; it simply fits the desired model within each subgroup and creates the forest plot. So you will need to fit your interaction-testing models separately, and save the p-values in scalars or macros. Then, in the ipdover saved dataset (e.g. "test.dta" in the above example) you can create a new variable containing the p-values (e.g. using a similar approach to the above, using string(n, s) -- see help f_string), and include this new variable in the rcols option to forestplot.

I hope you find this useful, and please let me know if you have any problems or suggestions.

Thanks,

David.

[Leo Martinez]

Hi David,

This was great! Your package is very useful and I haven't seen any other packages/programs that can do something similar.
I do have some questions if you have the time.

Currently this is the code I have:

gen hivind1=hivindx
gen hivind0=1-hivindx
gen ppd10interaction1=hivind1*ppd10new
gen ppd10interaction0=hivind0*ppd10new
ipdover, over(agegroupnew sex smoker relativenew bcg tbpast hiv etoh room1 indagenew indsex afbind indsmok coughnew cxrind cavtind housetyp windcat familysize) forest(xlabel(0.4(0.25)2.0, force) texts(120) nooverall astext(70) boxsca(2) xsize(18) ysize(15)) lcols((sum) ppd10interaction1 hivind1 ppd10interaction0 hivind0) nowt rr saving(test7.dta):  glm ppd10new hivindx if hivind<5 & ppd10new<5 & case<3, fam(poisson) link(log) nolog robust vce(cluster idno) eform
use test7, clear
gen prev0 = string(ppd10interaction1) + "/" + string(hivind1) + " (" + string(100*ppd10interaction1/hivind1, "%4.1f") + ")" if !missing(hivind1)
gen prev1 = string(ppd10interaction0) + "/" + string(hivind0) + " (" + string(100*ppd10interaction0/hivind0, "%4.1f") + ")" if !missing(hivind0)
label variable prev0 "Proportion of infected contacts of coinfected index cases (%)"
label variable prev1 "Proportion of infected contacts of monoinfected index cases (%)"
forestplot, rr nowt lcols(prev0 prev1)  xlabel(0.35 0.5 (0.25) 1.0 1.5 2.0, force) texts(90) nooverall astext(60) boxsca(2) xsize(18) ysize(15)

This created the figure below (not every variable is shown since I had a hard time actually snipping the whole figure since there are so many variables):



So this was almost exactly what I wanted and with your help I was able to put the proportion infection (with percentages) into the figure.

However, I am having some trouble putting the p-values into the figure. I only want to put a p-value for interaction variable as seen in my table below (also in my last post):

Since this value is attached to the row of the variable (so in the row "Age group, years" for example you can see the p-value for interaction) instead of attached to each value of every variable I am unsure how to set up my test.dta dataset in order to place the p-value in this location.
Is there a way to do this within the ipdover or ipdmetan package? I tried to make scalars for each p-value after running the commands. However, when I did this the p-value became attached to every value of the variable. So for "Age group, years" for example I would get the p-value 0.3702 for every value (0-4, 5-14, 15-25, etc) instead of for just the variable. Hopefully that makes sense.

Second question, more just details of the forest plot. Is there a way to make the null line dotted? I know in Stata it is difficult to change many of the graph features. But I was just wondering about this since this line is bolded. If not no worries!

Thanks again for all your help! You have made a very useful package!

Leo

[David Fisher]

Hi Leo,

Glad you're finding it useful! To answer your questions:


1) Placement of p-values: What you need is the _USE variable in test.dta. This variable tells forestplot what sort of data is in each line (observation), but it can also be useful when making manual adjustments. In your case, you're interested in the subgroup titles, for which _USE == 0. You'll also need the variable _OVER, which identifies the subgroups.

Say you have your interaction p-values for 10 subgroups stored in locals `pval1' , ..., `pval10'. You can then code something like:

use test7, clear

gen pval = .
forvalues i=1/10 {
     replace pval = `pval`i'' if _USE==0 & _OVER==`i'
}

format %05.3f pval
// alternatively, could use the a variation on the following:
// gen pvalstr = "p=" + string(pval, "%05.3f")

label var pval "p-value for interaction"
forestplot, rr nowt lcols(prev0 prev1) rcols(pval)

When I tested this just now with the low birthweight example in the ipdover help file, it clipped the end of the word "interaction" -- I'll need to take a look at this when I get a moment. If this happens to you, you can correct it by specifying the number of characters the column should use via the format command, e.g.

gen pvalstr = string(pval, "%05.3f")
format %11s pvalstr     // since there are 11 letters in the word "interaction"
label var pvalstr "p-value for interaction"
forestplot, rr nowt lcols(prev0 prev1) rcols(pvalstr)

2) Adjusting the null line: This is possible in the latest update of the ipdmetan package (dated 29jun2015) via the nlineopts() option to forestplot. I think I forgot to document this in the help file, though, oops!

e.g. to get a dashed null line:

forestplot, nlineopts(lpattern(dash))

Best wishes,

David.

[Michael Lanni]

David,

hello there. I made an account just now looking for some help. You sort of addressed the issue i am having in an earlier post but I am very new to stata (i used R before) and I didn't fully understand your solution. In april of this year you made this post for a problem similar to mine:

#####
http://www.statalist.org/forums/foru...-a-forest-plot
#############

Dear Nadia,

As Nick says, the "meta" command is very old and is not compatible with Stata's modern graphics setup. I would recommend installing the popular command "metan" (type "ssc describe metan"), or as an alternative, my own command "admetan" (within the ipdmetan package; type "ssc describe ipdmetan").

Could I ask what edits you are trying to make to your graph? There may be a way of obtaining the result you desire at the command-line, rather than in the Graph Editor.

############

I am using the metan command for my data and when I finish, the results window has a "more" link at the bottom. I click it and the forrest plot pops up but I cant seem to save it our edit it because it "is not a live graph." Im really not sure how to address this problem let alone what the solution is. I am doing a meta analysis of 4 papers to compare patients who had 1 of 2 surgeries. the outcomes of interest include whether they needed a reoperation or had certain complications. Being VERY new to stata and stats in general, can you help me out. Maybe you could even suggest a good stata for beginners book or youtube series or something? I'd really appreciate it as i just entered a fellowship and plan to use stats and stata all year (and most likely the rest of my career).

[David Fisher]

Hi Michael,

If you type help metan at the command line, towards the end of the documentation there are some clickable examples. What happens when you run one of these, and then try to save or edit? I've just tried it myself with the first example ("Risk difference from raw cell counts, random effects model, "label" specification with counts displayed") and it worked fine.

I'm no expert with the Graph Editor I'm afraid, but a quick Google of your (and Nadia's) error message suggests that the problem arises from the creation of graphs under an older version of Stata (version 7 or less). What happens if you type which metan? It should say "*! 3.04 21Sep2010". If it does not, type adoupdate metan to upgrade to the latest version, and try creating your graph again. Note that metan is a user-written program, and that Stata is fairly liberal with allowing user-written programs to circulate without official review (unless you wish to publish in the Stata Journal, of course).

If you are new to Stata, I have to say that the manuals themselves are a surprisingly good read. Start with [U] User's guide, and then (for your current purposes), [G] Graphics. Be careful to read the introduction sections, as otherwise the layout can seem confusing (or maybe that's just me).

Hope that helps!

[Michael Lanni]

Hey David,

Thank you for your response. It is helpful. I did what you said and attached screen shots. I'm not sure about the clickable examples you mentioned. I didn't have the data set to try and nothing looked clickable.

as for "which metan" I think you are onto something because mine is very old, but the "adoupdate" thing didn't seem to fix it for me.Is there another way to update it? I tried going in the way I originally added them and saw your name pop up a few times!

Thanks

[David Fisher]

Hi Michael,

Yes, you definitely have an old version, and that should explain everything else. The next question is why adoupdate failed. It couldn't find the file "metan.ado", but the results of which metan clearly show that it exists in the usual directory. Could you tell me (or copy-and-paste) the results of typing sysdir and adopath (as two separate commands)?

David.

[Michael Lanni]

That is something I can do, lol. Here are the results. What do you think?

Notes:

. sysdir
STATA: C:\Program Files\Stata13\
BASE: C:\Program Files\Stata13\ado\base\
SITE: C:\Program Files\Stata13\ado\site\
PLUS: c:\ado\plus\
PERSONAL: c:\ado\personal\
OLDPLACE: c:\ado\

. adopath
[1] (BASE) "C:\Program Files\Stata13\ado\base/"
[2] (SITE) "C:\Program Files\Stata13\ado\site/"
[3] "."
[4] (PERSONAL) "c:\ado\personal/"
[5] (PLUS) "c:\ado\plus/"
[6] (OLDPLACE) "c:\ado/"

.

[David Fisher]

Well, I'm stumped. It says exactly what I'd expect it to say.

All I can suggest is that you look in each of the directories labelled "PERSONAL" ("c:\ado\personal") , "PLUS" "c:\ado\plus" and "OLDPLACE" ("c:\ado") and delete any files there to do with the "metan" package. (metan.ado, obviously; but also any of the following if they exist: metan.hlp, metan.sthlp, metan7.*, metannt.*, and metan_examples*.*). Also look in your current directory (that's what [3] "." in the adopath list means).

Once you've made sure that all copies of metan are gone, you can re-install by typing ssc install metan at the command line. If successful, typing which metan should return "*! 3.04 21Sep2010", and your graphs should be editable.

Hope that works for you!

Thanks,

David.

[Michael Lanni]

You sir, are the man. I did what you said but it still wasn't working. I read the error i received and took a wild guess. After deleting labbe I tried to install metan again and it work. Thank you so much. ( see pictures for this story). Final question: I deleted labbe.ado, labbe.hlp, petometan.ado, petometan.scml, metan.hlp, metan.ado, ipdmetan.scml, ipdmetan.ado, admetan.scml, and admetan.ado. Do you think any of these files are useful and should I add back? Maybe I'll just leave them alone for now and add back things as I need them. Thanks again man. I am giving you a big, professional hug through the internet.

[David Fisher]

Brilliant -- sorry, I forgot about labbe as it has a completely different name!
You can leave petometan, ipdmetan and admetan; they're part of a different package entirely.
No worries -- glad you've gotten it working!

[Michael Lanni]

Thanks again mate. I'll add back petometan, ipdmetan and admetan (all by David J. Fisher, which I think is you. I am in the presence of royalty!)