3
$\begingroup$

I want to graphically show the change in relative risk along with confidence limits over 17 years.

Can I use a forest plot without meta-analysis because it will have the added advantage to tabulate Z-values and p-values also? Here I want to show the progressive change in the relative risk and NOT meta-analysis.

Or is there a better alternative for the graphic presentation?

$\endgroup$
  • 2
    $\begingroup$ As a reader, I would find that a good way to represent changes in RR. Usually time would be in the horizontal axis, so I would be tempted as a reader to turn the page 90 deg when I'm reading it. $\endgroup$ – Ming-Chih Kao Mar 18 '11 at 15:28
  • $\begingroup$ @Ming-Chin Kao: Thank you for the comment. Is it wrong to use a forest plot? The reason is that journals put a limit of 5 figures/tables. One forest plot will have the matter of 1 table and 1 graph along with confidence limits. In a graph with time on horizontal axis, can you kindly guide me how to draw RR with confidence limits? $\endgroup$ – DrWho Mar 18 '11 at 23:21
  • $\begingroup$ I got an answer to make a bar graph in Excel. ph.ucla.edu/epi/rapidsurveys/RScourse/chartconfinterval.pdfBut Excel will not generate publication quality graphs or eps files. Can any other software generate print quality graph? $\endgroup$ – DrWho Mar 19 '11 at 3:47
  • $\begingroup$ @mbq: Thank you very much for the editing to make the question clearer. $\endgroup$ – DrWho Mar 19 '11 at 3:49
  • $\begingroup$ I got partial answer for Print quality graph. In SPSS, Graph>Legacy dialogs>High-Low will generate this graph $\endgroup$ – DrWho Mar 19 '11 at 5:00
3
$\begingroup$

I've actually want to do the same thing with my data and I've created a forestplot version that allows both for expressions and advanced visualization of multiple methods that I use for displaying my multivariate analysis with 20-30 different variables. This code is in R and based on the package {rmeta}, it is still a work in progress but perhaps it might help:

library(rmeta)
forestplot2 <- function (labeltext, mean, lower, upper, 
        align = NULL, is.summary = FALSE,
        fontfamily.summary = NULL, fontfamily.labelrow = NULL, 
        clip = c(-Inf, Inf), xlab = "", zero = 0, 
        graphwidth = unit(51, "mm"), col = meta.colors(), 
        xlog = FALSE, xticks = NULL,
        increase.line_height = 1.2,
        lwd.xaxis=NULL,
        lwd.zero=NULL,
        lwd.forestlines=NULL,
        boxsize = NULL, ...) 
{
    require("grid") || stop("`grid' package not found")
    require("rmeta") || stop("`rmeta' package not found")

    # Addon: if two-dimentional mean, lower & upper 
    # the two are printed next to eachother
    if (NCOL(mean) == NCOL(lower) && NCOL(lower) == NCOL(upper) && NCOL(mean) > 0){
        org_mean <- mean
        org_lower <- lower
        org_upper <- upper
        if (NCOL(mean) > 1){
            mean <- as.vector(mean)
            lower <- as.vector(lower)
            upper <- as.vector(upper)
        }   
    }else{
        stop('Mean, lower and upper contain invalid number of columns')
    }
    # A function that is used to draw the different confidence intervals
    drawNormalCI <- function(LL, OR, UL, size, y.offset = 0.5, clr.line, clr.box) {
        size = 0.75 * size

        # Check if line/box outside graph
        clipupper <- convertX(unit(UL, "native"), "npc", valueOnly = TRUE) > 
                1
        cliplower <- convertX(unit(LL, "native"), "npc", valueOnly = TRUE) < 
                0
        box <- convertX(unit(OR, "native"), "npc", valueOnly = TRUE)
        clipbox <- box < 0 || box > 1

        # A version where arrows are added to the part outside 
        # the limits of the graph
        if (clipupper || cliplower) {
            ends <- "both"
            lims <- unit(c(0, 1), c("npc", "npc"))
            if (!clipupper) {
                ends <- "first"
                lims <- unit(c(0, UL), c("npc", "native"))
            }
            if (!cliplower) {
                ends <- "last"
                lims <- unit(c(LL, 1), c("native", "npc"))
            }
            grid.lines(x = lims, y = y.offset, arrow = arrow(ends = ends, 
                            length = unit(0.05, "inches")), gp = gpar(col = clr.line, lwd=lwd.forestlines))
            if (!clipbox) 
                grid.rect(x = unit(OR, "native"), y=y.offset, width = unit(size, 
                                "snpc"), height = unit(size, "snpc"), gp = gpar(fill = clr.box, 
                                col = clr.box))
        }
        else {
            grid.lines(x = unit(c(LL, UL), "native"), y = y.offset, 
                    gp = gpar(col = clr.line, lwd=lwd.forestlines))
            grid.rect(x = unit(OR, "native"), y=y.offset, width = unit(size, 
                            "snpc"), height = unit(size, "snpc"), gp = gpar(fill = clr.box, 
                            col = clr.box))
            if ((convertX(unit(OR, "native") + unit(0.5 * size, 
                                        "lines"), "native", valueOnly = TRUE) > UL) && 
                    (convertX(unit(OR, "native") - unit(0.5 * size, 
                                        "lines"), "native", valueOnly = TRUE) < LL)) 
                grid.lines(x = unit(c(LL, UL), "native"), y = 0.5, 
                        gp = gpar(col = clr.line, lwd=lwd.forestlines))
        }
    }
    drawSummaryCI <- function(LL, OR, UL, size) {
        grid.polygon(x = unit(c(LL, OR, UL, OR), "native"), y = unit(0.5 + 
                                c(0, 0.5 * size, 0, -0.5 * size), "npc"), gp = gpar(fill = col$summary, 
                        col = col$summary))
    }

    # A function used for fetching the text or expression
    # from the supplied labeltext.
    fetchRowLabel <- function(label_type, labeltext, i, j){
        if (label_type=="expression"){
            # Haven't figured out it this is possible with 
            # a multilevel expression
            row_column_text <- labeltext[[i]]
        }
        else if(label_type=="list"){
            # I get annoying warnings with this
            #if (!is.expression(labeltext[[j]][[i]]) && is.na(labeltext[[j]][[i]]))
            #    return(FALSE)
            row_column_text <- labeltext[[j]][[i]]
        }
        else{
            if (is.na(labeltext[i, j])) 
                return(FALSE)
            row_column_text <- labeltext[i, j]
        }
        return(row_column_text)
    }
    validateLabelList <- function(labelList){
        l = length(labelList[[1]])
        if (length(labelList) == 1)
            return(TRUE)

        for(i in 1:length(labelList)){
            # All elements should have the same length
            if (l != length(labelList[[i]]))
                return(FALSE)
        }

        return(TRUE)
    }

    # The previous algorithm failed when I added the expressions
    findWidestGrob <- function (grob.list){
        len <- c()
        for (i in 1:(length(grob.list))){
            if (is.object(grob.list[[i]])){
                # There is a tendency of underestemating grob size
                # when there are expressions
                mm <- convertWidth(grobWidth(grob.list[[i]]), "mm")
                len <- append(len, mm)
            }else{
                len <- append(len, 0)
            }
        }

        return(unit(max(len), "mm"))
    }

    # Workaround for expressions
    if (is.expression(labeltext)){
        widthcolumn <- c(TRUE)
        # Can't figure out multiple levels of expressions
        nc <- 1
        nr <- length(labeltext)
        label_type = "expression"
    }else if(is.list(labeltext)){
        if (validateLabelList(labeltext) == FALSE)
            stop("Invalid labellist, it has to be formed as a matrix m x n elements")

        # Can't figure out multiple levels of expressions
        nc <- length(labeltext)

        # Works but probably not as orig. intended
        widthcolumn <- rep(TRUE, nc)

        nr <- length(labeltext[[1]])
        label_type = "list"
    }else{
        # Original code for matrixes
        widthcolumn <- !apply(is.na(labeltext), 1, any)
        nc <- NCOL(labeltext)
        nr <- NROW(labeltext)
        label_type = "matrix"
    }

    # Prepare the summary and align variables
    if (is.null(align)){
        align <- c("l", rep("r", nc - 1))
    } else {
        align <- rep(align, length = nc)
    } 

    is.summary <- rep(is.summary, length = nr)

    getLabels <- function(){
        labels <- vector("list", nc)

        # Walk through the labeltext
        # Creates a list matrix with 
        # The column part
        for (j in 1:nc) {
            labels[[j]] <- vector("list", nr)

            # The row part
            for (i in 1:nr) {
                txt_out <- fetchRowLabel(label_type, labeltext, i, j)
                if (is.expression(txt_out) || is.character(txt_out) || is.numeric(txt_out)){
                    x <- switch(align[j], l = 0, r = 1, c = 0.5)

                    just <- switch(align[j], 
                            l = "left", 
                            r = "right",
                            c = "center")

                    # Bold the text if this is a summary
                    if (is.summary[i]){
                        if (is.expression(txt_out)){
                            x <- 0.5
                        }else{
                            x <- 0
                        }
                        # Create a textGrob for the summary
                        labels[[j]][[i]] <- textGrob(txt_out, x = x,
                                just = just,
                                hjust = 0,
                                gp = gpar(fontface = "bold",
                                        fontfamily=fontfamily.summary,  
                                        cex = 1.3, 
                                        col = rep(col$text, length = nr)[i]))
                    }else{
                        # Create a textGrob with the current row-cell for the label
                        labels[[j]][[i]] <- textGrob(txt_out, x = x,
                                just = just, 
                                gp = gpar(fontface = "plain",
                                        fontfamily=fontfamily.labelrow, 
                                        col = rep(col$text, length = nr)[i]))
                    }
                }
            }
        }
        return(labels)
    }
    labels <- getLabels()

    # Set the gap between columns
    colgap <- unit(6, "mm")

    # There is always at least one column so grab the widest one
    # and have that as the base for the column widths
    colwidths <- unit.c(findWidestGrob(labels[[1]]), colgap)

    # If multiple row label columns, add the other column widths
    if (nc > 1) {
        for (i in 2:nc){
            colwidths <- unit.c(colwidths, 
                    findWidestGrob(labels[[i]][widthcolumn]), 
                    colgap)  
        } 
    }

    # Add the base grapwh width to the total column width
    # default is 2 inches
    colwidths <- unit.c(colwidths, graphwidth)

    # Create the canvas for the plot
    plot.new()

    # The base viewport, set the increase.line_height paremeter if it seems a little
    # crowded between the lines that might happen when having multiple comparisons
    pushViewport(viewport(layout = grid.layout(nr + 1, length(colwidths), 
                            widths = colwidths,
                            heights = unit(c(rep(1, nr), 0.5)*increase.line_height, "lines"))))

    # Get width of the lines
    cwidth <- (upper - lower)

    # If the borders are smaller than the upper/lower limits 
    # then clip the graph. The line will have arrows indicating
    # that it continues beyond the graph
    # The zero bar has to be on the chart though!
    xrange <- c(min(zero, max(min(lower, na.rm = TRUE), clip[1])), max(min(max(upper, 
                                                                na.rm = TRUE), clip[2]), zero))

    # Create the fourth argument 4 the drawNormalCI() function 

    # If boxsize was provided override the info
    if (!is.null(boxsize)){
        # If matrix is provided this will convert it
        # to a vector but it doesn't matter in this case
        info <- rep(boxsize, length = length(cwidth))
    }else{
        info <- 1/cwidth
        info <- info/max(info[!is.summary], na.rm = TRUE)
        info[is.summary] <- 1
    }

    printLabels <- function(){
        # Output the labels
        # The column
        for (j in 1:nc) {
            # The row
            for (i in 1:nr) {
                if (!is.null(labels[[j]][[i]])) {
                    # The column position is 2 * j - 1 due to the column gap
                    pushViewport(viewport(layout.pos.row = i, 
                                    layout.pos.col = 2 * j - 1))
                    grid.draw(labels[[j]][[i]])
                    popViewport()
                }
            }
        }
    }
    printLabels()

    pushViewport(viewport(layout.pos.col = 2 * nc + 1, xscale = xrange))

    # Print y-axis - the vertical "zero" axis
    grid.lines(x = unit(zero, "native"), y = 0:1, gp = gpar(col = col$zero, lwd=lwd.zero))

    # Print x-axis
    if (xlog) {
        if (is.null(xticks)) {
            ticks <- pretty(exp(xrange))
            ticks <- ticks[ticks > 0]
        }
        else {
            ticks <- xticks
        }
        if (length(ticks)) {
            if (min(lower, na.rm = TRUE) < clip[1]) 
                ticks <- c(exp(clip[1]), ticks)
            if (max(upper, na.rm = TRUE) > clip[2]) 
                ticks <- c(ticks, exp(clip[2]))
            xax <- xaxisGrob(gp = gpar(cex = 0.6, col = col$axes, lwd=lwd.xaxis), 
                    at = log(ticks), name = "xax")
            xax1 <- editGrob(xax, gPath("labels"), label = format(ticks, 
                            digits = 2))
            grid.draw(xax1)
        }
    } else {
        if (is.null(xticks)) {
            grid.xaxis(gp = gpar(cex = 0.6, col = col$axes, lwd=lwd.xaxis))
        }
        else if (length(xticks)) {
            grid.xaxis(at = xticks, gp = gpar(cex = 0.6, col = col$axes, lwd=lwd.xaxis))
        }
    }

    # Write the label for the x-axis
    grid.text(xlab, y = unit(-2, "lines"), gp = gpar(col = col$axes))
    popViewport()

    # Output the different confidence intervals
    for (i in 1:nr) {
        if (is.na(mean[i])) 
            next
        pushViewport(viewport(layout.pos.row = i, layout.pos.col = 2 * 
                                nc + 1, xscale = xrange))

        if (is.matrix(org_mean)){
            low_values <- org_lower[i,] 
            mean_values <- org_mean[i,] 
            up_values <- org_upper[i,]
            info_values <- matrix(info, ncol=length(low_values))[i, ]
        }else{
            low_values <- lower[i] 
            mean_values <- mean[i] 
            up_values <- upper[i]
            info_values <- info[i]
        }

        # The line and box colors may vary
        clr.line <- rep(col$line, length=length(low_values))
        clr.box <- rep(col$box, length=length(low_values))
        if (is.summary[i]) 
            drawSummaryCI(low_values, mean_values, up_values, info_values)
        else{
            if (length(low_values) > 1){
                y.offset_base <- 0.2
                y.offset_increase <- (1-y.offset_base*2)/length(low_values)

                for(j in 1:length(low_values)){
                    drawNormalCI(low_values[j], 
                            mean_values[j], 
                            up_values[j], 
                            info_values[j], 
                            y.offset = y.offset_base + (j-1)*y.offset_increase,
                            clr.line = clr.line[j],
                            clr.box = clr.box[j])
                }
            }else{
                drawNormalCI(low_values, mean_values, up_values, info_values,
                        clr.line = clr.line, clr.box = clr.box)
            }
        }


        popViewport()
    }
    popViewport()
}

And an example:

row_names <- list(
        list("variable = 0", "variable = 1", expression(variable >= 2)),
        list(expression(bar(x)==1.8), expression(bar(x) == 1.4), "some cell data"))
test_data <- data.frame(coef1=c(1.59, 1.3, 1.24),
        coef2=c(1.7, 1.4, 1.04),
        low1=c(1.3, 1.1, 0.99),
        low2=c(1.6, 1.2, 0.7),
        high1=c(1.94, 1.6, 1.55),
        high2=c(1.8, 1.55, 1.33))
attach(test_data)
coef <- cbind(coef1, coef2)
low <- cbind(low1, low2)
high <- cbind(high1, high2)
forestplot2(row_names, coef, low, high, zero = 1, 
        col=meta.colors(box=c("royalblue", "gold"), ,line=c("darkblue", "orange"), summary="royalblue"))
detach(test_data)

Plots this image: enter image description here

I've also created more advanced plots, that I've also used in one of my questions, here's a plot that I've created to compare a three different methods for survival analysis:

enter image description here

I have some additional functions that help me prepare the regression data for the plot if you find it interesting. I'll try to publish it in the future but I'm still working on the project :)

Update

If you want to display change over time perhaps a simple line plot will do the trick? Two lines with confidence intervals

You can find the full example here

$\endgroup$
  • $\begingroup$ (+1) very informative, thanks for the codez. $\endgroup$ – suncoolsu Sep 4 '11 at 21:23
  • $\begingroup$ @MaxG: Thank you very much for the detailed explanation. I require sometime to understand your answer. Kindly tell me which software you have used for simple line plot. $\endgroup$ – DrWho Sep 5 '11 at 8:47
  • $\begingroup$ I used a free statistical software R. I've used SPSS previously but since I do a lot of programming I like tweaking all the settings and R allows you to do amazing stuff. It is good to start off easy though and you should probably look into some sites like statmethods.net for som hints. Best thing to do is to look at the code, copy it and start changing step by step until you have what you want. The link has the full code. You can find R at cran.r-project.org and I strongly recommend combining it with RStudio rstudio.org $\endgroup$ – Max Gordon Sep 5 '11 at 9:33
  • $\begingroup$ (+1) for the code and the graphs together $\endgroup$ – Andrew Sep 6 '11 at 15:50
1
$\begingroup$

There is no reason why a forest plot is the exclusive provenance of a meta-analysis. It is admittedly an extremely useful plot for meta-analysis, but the forest plot is merely a very, very good way to show a point estimate and confidence interval over a continuous variable.

I would consider two modifications to the "typical" forest plot. If you're trying to just show the change over time, I would dispense with any graphing of a summary estimate (the usual diamond), and I would consider tilting the graph on its side. Outside meta-analysis, time seems most often shown on the x-axis.

Any basic meta-analysis package should be able to do this. There is new free software to do this called 'Forest Plot Viewer' freely available at http://ntp.niehs.nih.gov/go/tools_forestplotviewer

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.