#------------------------------------------------------------------------------------- # File name : tas_diamondgraph_figure7.ssc # Date : 12/15/2005 (Version 2) # Authors : Xiuhong Li and Alvaro Muņoz # # Figure 7 in Li X, Buechner J, Tarwater P, Muņoz A. A Diamond-Shaped # Equiponderant Graphical Display of the Effects of Two Categorical Predictors on # Continuous Outcomes. The American Statistician, August 2003, Vol 57, No. 3, 193-199. # # Figure legend: Relative risk of breast cancer by adult weight change and hormone # use among postmenopausal women reported by Huang et al. (1997). # The graph on the left panel depicts the relative risks # (i.e., pi= rri / highest rr) and the one on the right panel depicts # the excess risk (i.e.,pi=(rri - lowest rr)/(highest rr - lowest rr). # # To run code provided, you must have a license to run S-plus. Authors used code to # produce figures in the paper in The American Statistician using S-plus 6.0. # Recipients of code may need to do changes to use in other versions and for # different applications. #------------------------------------------------------------------------------------- #------------------------------------------------------------------------------------- # 1. Obtain the number of categories of the two predictors; # row = number of categories of the predictor on the southwest side of Diamond Graph; # col = number of categories of the predictor on the southeast side of Diamond Graph; #------------------------------------------------------------------------------------- row <- 5; col <- 3; #------------------------------------------------------------------------------------- # 2. Select a square region for plot, and put 2 plots on one page; #------------------------------------------------------------------------------------- par(pty="s",mfrow=c(1,2)); plotfigure7 <- function(panel) # Create a function to plot a panel of figure 7; { #------------------------------------------------------------------------------------- # 3. Assign the coordinate of the lowest point of Diamond Graph to (0,0), and define # the ranges for x and y axes based on the number of categories of two predictors. #------------------------------------------------------------------------------------- plot(0,0, xlim=c(0-row, col), ylim=c(0,row+col),axes=F, type="n", xlab="", ylab=""); #------------------------------------------------------------------------------------- # 4. Obtain the coordinate of the center of each square cell (cx,cy) with diagonal= 2; # i.e., the distance from the center to a vertex of a square cell equals to 1. #------------------------------------------------------------------------------------- for (i in 1:row) { if (i == 1){ cx <- seq(0,(col-1) ) cy <- seq(1,col) } else{ cx <- c(cx, seq( (1-i),(col-i) ) ) cy <- c(cy, seq(i,(col+i-1)) ) } } #------------------------------------------------------------------------------------- # 5. Read in the outcomes, which start from the value in the cell 1 (at bottom) # to cell 2, to cell 3, and finally to cell 4 (at top), as depicted below in # the 2x2 graph. # /4\ <- top # / \/\ # \3/\2/ # \1/ <- bottom # # Empty cells (with no subjects) are assigned -1 and will not be displayed. # Only positive outcomes will be plotted as polygons in square cells. #------------------------------------------------------------------------------------- outcome <- c(2.0, 1.3, 1.7, 1.6, 1.3, 1.7, 1.2, 1.5, 2.1, 1.0, 1.3, 1.8, 1.1, 1.1, 1.9); #------------------------------------------------------------------------------------- # 6.If data are not proportional, convert the outcomes to proportions; # to get relative risks,each value is divided by the highest one within that outcome; # to get excess risks, the difference between each value and the minimum outcome is # divided by the difference between the maximum and the minimum outcome. # If data are proportional, no need to convert. #------------------------------------------------------------------------------------- maxoutcome <- max(outcome) minoutcome <- min(outcome[outcome>=0]) if (panel=="Relative Risk") p<-outcome/maxoutcome #get p for left panel(relative risk) else p<-(outcome-minoutcome)/(maxoutcome-minoutcome)#get p for right panel(excess risk) #------------------------------------------------------------------------------------- # 7. Calculate the 'scale' based on the number of categories of two predictors; # We have found scaling all proportions (p) by 'scale' useful in the # graphical illustration when p is high. # If a page has only one graph, we use scale = 1.005 - 0.005(row+col). # But figure 4 has 2 plots on one page, the actual square cell size is similar to # that of a 5x6 plot. So we use scale = 1.005-0.005*(5+6) here. #------------------------------------------------------------------------------------- scale <- 1.005 - 0.005*(5+6); p[p<1] <- p[p<1] *scale p[p==1] <- p[p==1]*scale+0.005 #------------------------------------------------------------------------------------- # 8. Draw the outline of the square cells. # 1) draw 45 degree lines from southwest to northeast # 2) draw -45 degree lines from southeast to northwest #------------------------------------------------------------------------------------- for (i in 0:row) lines(c(-i,col-i),c(i,col+i),lwd=0.01) for (i in 0:col) lines(c(i, -row+i), c(i, row+i),lwd=0.01) #------------------------------------------------------------------------------------- # 9. Calculate the coordinates of the six vertices of inside polygons; # start clockwise from northwest vertex as (x1, y1). # x1,y1 _____ x2,y2 # / \ # x6,y6 / \ x3,y3 # \ / # \_____/ # x5,y5 x4,y4 #------------------------------------------------------------------------------------- x1<- cx - (1-p)/2 x2<- cx + (1-p)/2 x3<- cx + (1+p)/2 x4<- cx + (1-p)/2 x5<- cx - (1-p)/2 x6<- cx - (1+p)/2 y1<- cy + p y2<- cy + p y3<- cy y4<- cy - p y5<- cy - p y6<- cy #------------------------------------------------------------------------------------- # 10. Plot a polygon inside each square cell, and then print the value of the outcome # at the center in each square cell. # # Here we use cex= 0.7*scale for the labels. Other applications may # need to use a different factor so that labels can be read easily. #------------------------------------------------------------------------------------- n <- row*col for (i in 1:n) { if (p[i]>0) polygon(c(x1[i],x2[i],x3[i],x4[i],x5[i],x6[i]), c(y1[i],y2[i],y3[i],y4[i],y5[i],y6[i]), density=-.1, border=F, col=5) # make the outcome=2.0 shown as 2.0 not as 2; outcome[outcome==2.0] <-'2.0'; if (outcome[i]>=0){ # In Window system, S-plus needs 0.03 to be added to cy so that lable is properly centered text(cx[i], cy[i]+0.03, outcome[i], cex=0.7*scale, lwd=2); } } #------------------------------------------------------------------------------------- # 11. Below are the commands to print labels for two predictors. # It requires trial and error for the proper location; # but users can skip the following steps and add the labels in PowerPoint. #------------------------------------------------------------------------------------- # print labels on the southwest side of Diamond Graph text(c(-4.5,-4.5), c(1,1), adj=0.5, cex=0.85,"Adult") text(c(-4.5,-4.5), c(0.7,0.7), adj=0.5, cex=0.85,"Weight") text(c(-4.5,-4.5), c(0.4,0.4), adj=0.5, cex=0.85,"Change, kg") text(c(-0.7,-0.7), c(.35,.35), adj=1, cex=.7,"Gain >20.0") text(c(-1.7,-1.7), c(1.35,1.35), adj=1, cex=.7,"Gain 10.1-20.0") text(c(-2.7,-2.7), c(2.35,2.35), adj=1, cex=.7,"Gain 2.1-10.0") text(c(-3.7,-3.7), c(3.35,3.35), adj=1, cex=.7,"Loss or Gain 2.0") text(c(-4.7,-4.7), c(4.35,4.35), adj=1, cex=.7,"Loss >2.0") # print labels on the southeast side of Diamond Graph text(c(3, 3), c(0.7,0.7), adj=0.5, cex=0.85,"Hormone") text(c(3, 3), c(0.4,0.4), adj=0.5, cex=0.85,"Use") text(c(.7, .7), c(.35,.35), adj=0, cex=.7,"Never") text(c(1.7, 1.7), c(1.35,1.35), adj=0, cex=.7,"Past") text(c(2.7, 2.7), c(2.35,2.35), adj=0, cex=.7,"Current") title(main=panel) }# end of function plotfigure7; # Call the function to plot relative risk in left panel and excess risk in right panel plotfigure7("Relative Risk") plotfigure7("Excess Risk")