# Visualizing change in risk ratio along with confidence limits

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?

• 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. – Ming-Chih Kao Mar 18 '11 at 15:28
• @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? – DrWho Mar 18 '11 at 23:21
• 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? – DrWho Mar 19 '11 at 3:47
• @mbq: Thank you very much for the editing to make the question clearer. – DrWho Mar 19 '11 at 3:49
• I got partial answer for Print quality graph. In SPSS, Graph>Legacy dialogs>High-Low will generate this graph – DrWho Mar 19 '11 at 5:00

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")

# 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:

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:

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?

You can find the full example here

• (+1) very informative, thanks for the codez. – suncoolsu Sep 4 '11 at 21:23
• @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. – DrWho Sep 5 '11 at 8:47
• 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 – Max Gordon Sep 5 '11 at 9:33
• (+1) for the code and the graphs together – Andrew Sep 6 '11 at 15:50

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