# Forest plot with Mathematica 9

Does anyone have experience with or a package for the creation of forest plots in Mathematica 9?

For instance for subgroup analyses of cox model or meta-analysis (e.g rmeta package in R)?

-
Can you specify your question? Or is it really useful when someone answers "yes, I have expierience!"? What exactly are you trying to achieve? What have you tried? etc.. – halirutan Dec 4 '12 at 8:02
Any reason that it should be specifically for version 9 (out for about a week now)? – Sjoerd C. de Vries Dec 4 '12 at 8:11
I have been using R for survival analysis. This is now better supported in Mathematica version 9. As part of a planned project, in addition to modelling, forest plots of relative risk would be helpful. I accept, if this is not a priority, that I will use R packages (as I have done) while trying to develop my own functions. I merely was asking a community question. – ubpdqn Dec 4 '12 at 8:54
@halirutan I think it's a valid question. It's not obvious for you that he wants to know if such plots exists? And how can he do that? – Murta Dec 4 '12 at 10:27
@Murta Yes, you are right. I was a bit afraid, that there is too less information to give an exact answer but Vitaliy jumped in. – halirutan Dec 4 '12 at 10:31

The plotting is a little more cumbersome than the computation assuming you want relative risk shown in your forest plot. As you point out, one can use CoxModelFit to compute everything that needs computed.

The larynx cancer data features the variables age and four levels of cancer stage.

larynx = ExampleData[{"Statistics", "LarynxCancer"}];
x = larynx[[All, {1, 3}]];
e = EventData[larynx[[All, 2]], larynx[[All, 5]]];
mod = CoxModelFit[{x, e}, {stage, age}, {stage, age}, NominalVariables -> stage];


Creating the plot is a little ugly if you want a lot of control over the appearance. Here is a first crack at it using graphics primitives that allows the diamonds to scale with the plot.

forestPlot[mod_CoxModel] :=
Block[{rr, rrci, vars, p},
rr = mod["RelativeRisk"];
rrci = mod["RelativeRiskConfidenceIntervals"];
vars = mod["BasisFunctions"] /.
HoldPattern[DiscreteIndicator[s_, i_, _]] :> s[i];
p = Length[vars];

Graphics[Join[Join[{{Dashed, Line[{{1, 0}, {1, p + 1}}]}},
MapThread[{Thick, Line[{{#1[[1]], #2}, {#1[[2]], #2}}]} &,
{rrci,Range[p]}]],
Polygon[{{-.2 + #1, #2}, {#1, .1 + #2}, {.2 + #1, #2},
{#1, -.1 + #2}}]} &, {rr, Range[p]}]],
AspectRatio -> 1, Frame -> True,
FrameTicks -> {Automatic, Transpose[{Range[p], vars}], None, None},
GridLines -> Automatic, GridLinesStyle -> White
]
]


With the larnyx cancer data we get...

forestPlot[mod]


-
Thank you for the feedback and this answer. – ubpdqn Dec 5 '12 at 7:40

You can work with RLink to create native R graphics - this and this examples. I guess forestplot is no different if you install related R packages. I may be far off here, but can't you replicate these plots with built-in Mathemtica functionality? - like these examples on "Compare the distribution of salaries for several departments at a university:"

salaries =  ExampleData[{"Statistics", "UniversitySalaries"}, "DataElements"];
depts = {"Mathematics", "History", "English", "Chemistry", "Law", "Physics", "Statistics"};
data = Table[Cases[salaries, {d, _, salary_, "A"} :> salary], {d, depts}];
all = Cases[salaries, {_, _, salary_, "A"} :> salary];

DistributionChart[data,
ChartLabels -> Placed[{depts, Length /@ data}, {Axis, Center}],
ChartStyle -> 54, GridLines -> {{{Median[all], Gray}}, None}, BarOrigin -> Left]


BoxWhiskerChart[data, {"Outliers", {"MedianMarker", 1, Directive[Thick, White]}},
ChartLabels -> Placed[{depts, Length /@ data}, {Axis, Before}],
ChartStyle -> 10, GridLines -> {{{Median[all], Gray}}, None}, BarOrigin -> Left]


-
Thank you very much. I agree the form of the DistributionChart and BoxWhiskerChart match the forest plot. There are two parts. The first calculational (odds ratio, relative risk or hazard ratios and their respective standard errors). This can be done in Mathematica such as the larygeal cancer example in the Mathematica 9 documentation of CoxModelFit. The second part is the plot. The structure for the plot is simpler than for the distributional charts. – ubpdqn Dec 4 '12 at 11:05
I meant to write that the data structure (after calculation) for the plot is simpler and would not work inputs for the chart functions for distributional representations. – ubpdqn Dec 4 '12 at 11:57

I appreciate that the (somewhat critical) commentary regarding the lack of specificity of my question.

I have written my own function for the plotting (once the calculational aspects from the model or meta-analysis have been collected).

forestplot[x_List, tcks_, opts : OptionsPattern[]] :=
Module[{lg, ln, int, mx, lines}, lg = Length[x]; ln = Range[lg];
mx = Max[x[[All, 3]]];
int = {Transpose[{x[[All, 1]], ln}], Transpose[{x[[All, 3]], ln}]};
lines = MapThread[{Black, Line[{#1, #2}]} &, int];
ListPlot[Transpose[{x[[All, 2]], ln}],
PlotMarkers -> Style[\[DiamondSuit], 20],
PlotRange -> {{-2, 1.1 mx}, {0, Length[x] + 1}},
Epilog -> {{Red, Line[{{1, 0}, {1, lg + 1}}], lines},
MapThread[Text[#1, {-1, #2}] &, {tcks, ln}]},
Axes -> {True, False}, Ticks -> {Table[j, {j, 0, mx, 1}], None},
Evaluate[FilterRules[{opts}, Options[ListPlot]]]]]


The form of x is {lower bound, point estimate, upper bound} and tcks is list of labels. There are formatting issues and I am yet to make a summary statistic and formatting for this. Further, options for plot markers could be customized.

The following utilizes the laryngeal cancer dataset and functions for extracting the point estimates, standard errors and then calculation of 95% confidence intervals.

I would suggest that for new and inexperienced (but interested) users like me asking experts regarding functionalities is within the purview and purpose of this exchange. Further, despite the nominal release of Mathematica 9, I am certain that many contributors to this exchange had access to its functionality (as I saw at IMS2012 in May in London).

Finally, comments and commentary that are personally critical or make assumptions regarding the ignorance or other attributes are at the very least unhelpful. (Having surveyed this site there are a spectrum of users with diverse backgrounds and aims. This diversity would appear to something to be encouraged rather than crushed. I have previously fully accepted the NARQ of a previous question I submitted. In my naive opinion, this question is legitimate. ).

-
Maybe a little late, but welcome to the group! It is quite common that questions need to be edited by the author to provide additional detail. – Gustavo Delfino Dec 4 '12 at 16:43
I just asked why you wanted the solution aimed at version 9. The answer you provide here doesn't use any special feature of v9. So the question still is, why do you explicitly mention version 9? What feature of v9 would you want to be used? Is that so strange to ask? – Sjoerd C. de Vries Dec 4 '12 at 21:22
@SjoerdC.deVries I suspect they say M9 because CoxModelFit makes computing this stuff trivial. It was quite difficult to do before M9. – Andy Ross Dec 4 '12 at 23:13
@SjoerdC.deVries I am sorry for my sensitivity. I find this site an excellent resource. As @Andy Ross comments features such as CoxModelFit make this somewhat easier. Further, for hazard ratios, calculation of standard errors from the covariance matrix (as in the examples in documentation) is required. @Gustavo Delfino I appreciate the need for editing and am grateful for the comments. – ubpdqn Dec 5 '12 at 7:44