#------------------------------------------------------------------------------------- # File name : tas_diamondgraph_figure5.ssc # Date : 12/15/2005 (Version 2) # Authors : Xiuhong Li and Alvaro Muņoz # # Figure 5 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:Two-dimensional approach for depicting the age-adjusted rates of end-stage # renal disease due to any cause per 100,000 person-years according # to systolic and diastolic blood pressures. Each shaded area is a # proportion relative to the largest incidence rate with 211.7 being 1 # (full shading). The exact numeric values of the rates are used as the # labels. Cells with no or very few individuals are blank. # # 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 <- 6 col <- 6 #------------------------------------------------------------------------------------- # 2. Select a square region for plot; #------------------------------------------------------------------------------------- par(pty="s") #------------------------------------------------------------------------------------- # 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( 5.3, 6.0, 4.9, 9.8, -1, -1, 6.8, 6.7, 8.4,11.1,40.2,-1, 13.3,12.8,12.4,16.4,25.4,-1, 22.1,27.5,22.6,27.1,35.4,62.3, 90.2,68.4,57.8,62.3,49.7,88.9, 49.6,95.9,211.7,143.5,124.4,205.6) #------------------------------------------------------------------------------------- # 6. If data are not proportional, convert the outcomes to proportions - # each value is divided by the highest one within that outcome. # If data are proportional, no need to convert. #------------------------------------------------------------------------------------- maxoutcome <- max(outcome) p <- outcome/maxoutcome #------------------------------------------------------------------------------------- # 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. #------------------------------------------------------------------------------------- scale <- 1.005 - 0.005*(row+col) 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=6.0 shown as 6.0 not as 6; outcome[outcome==6.0] <-'6.0' if (outcome[i]>=0){ # In Windows 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(-5.5,-5.5), c(1.8,1.8), adj=0.5, cex=1.0,"Systolic (mm Hg)") text(c(-0.7,-0.7), c(.35,.35), adj=1, cex=.8,"<120") text(c(-1.7,-1.7), c(1.35,1.35), adj=1, cex=.8,"120-129") text(c(-2.7,-2.7), c(2.35,2.35), adj=1, cex=.8,"130-139") text(c(-3.7,-3.7), c(3.35,3.35), adj=1, cex=.8,"140-159") text(c(-4.7, -4.7), c(4.35,4.35), adj=1, cex=.8,"160-179") text(c(-5.7, -5.7), c(5.35,5.35), adj=1, cex=.8,">=180") # print labels on the southeast side of Diamond Graph text(c(5.5, 5.5), c(1.8,1.8), adj=0.5, cex=1.0,"Diastolic (mm Hg)") text(c(.7, .7), c(.35,.35), adj=0, cex=.8,"<80") text(c(1.7, 1.7), c(1.35,1.35), adj=0, cex=.8,"80-84") text(c(2.7, 2.7), c(2.35,2.35), adj=0, cex=.8,"85-89") text(c(3.7, 3.7), c(3.35,3.35), adj=0, cex=.8,"90-99") text(c(4.7, 4.7), c(4.35,4.35), adj=0, cex=.8,"100-109") text(c(5.7, 5.7), c(5.35,5.35), adj=0, cex=.8,">=110")