1## Modifications - MF - 1 Dec 2010 2# -- change default colors to more distinguishable values 3# -- allow to work with >3 dimensional arrays 4# -- modified defaults for mfrow/mfcol to give landscape display, nr <= nc, rather than nr >= nc 5 6# Take a 2+D array and return a 3D array, with dimensions 3+ as a single dimension 7# Include as a separate function, since it is useful in other contexts 8array3d <- function(x, sep=':') { 9 if(length(dim(x)) == 2) { 10 x <- if(is.null(dimnames(x))) 11 array(x, c(dim(x), 1)) 12 else 13 array(x, c(dim(x), 1), c(dimnames(x), list(NULL))) 14 return(x) 15 } 16 else if(length(dim(x))==3) return(x) 17 else { 18 x3d <- array(x, c(dim(x)[1:2], prod(dim(x)[-(1:2)]))) 19 if (!is.null(dimnames(x))) { 20 n3d <- paste(names(dimnames(x))[-(1:2)], collapse=sep) 21 d3d <- apply(expand.grid(dimnames(x)[-(1:2)]), 1, paste, collapse=sep) 22 dimnames(x3d) <- c(dimnames(x)[1:2], list(d3d)) 23 names(dimnames(x3d))[3] <- n3d 24 } 25 return(x3d) 26 } 27} 28 29"fourfold" <- 30function(x, 31# color = c("#99CCFF","#6699CC","#FF5050","#6060A0", "#FF0000", "#000080"), 32 color = c("#99CCFF","#6699CC","#FFA0A0","#A0A0FF", "#FF0000", "#000080"), 33 conf_level = 0.95, 34 std = c("margins", "ind.max", "all.max"), margin = c(1, 2), 35 space = 0.2, main = NULL, sub = NULL, mfrow = NULL, mfcol = NULL, extended = TRUE, 36 ticks = 0.15, p_adjust_method = p.adjust.methods, newpage = TRUE, 37 fontsize = 12, default_prefix = c("Row", "Col", "Strata"), sep = ": ", varnames = TRUE, 38 return_grob = FALSE) 39 40{ 41 ## Code for producing fourfold displays. 42 ## Reference: 43 ## Friendly, M. (1994). 44 ## A fourfold display for 2 by 2 by \eqn{k} tables. 45 ## Technical Report 217, York University, Psychology Department. 46 ## http://datavis.ca/papers/4fold/4fold.pdf 47 ## 48 ## Implementation notes: 49 ## 50 ## We need plots with aspect ratio FIXED to 1 and glued together. 51 ## Hence, even if k > 1 we prefer keeping everything in one plot 52 ## region rather than using a multiple figure layout. 53 ## Each 2 by 2 pie is is drawn into a square with x/y coordinates 54 ## between -1 and 1, with row and column labels in [-1-space, -1] 55 ## and [1, 1+space], respectively. If k > 1, strata labels are in 56 ## an area with y coordinates in [1+space, 1+(1+gamma)*space], 57 ## where currently gamma=1.25. The pies are arranged in an nr by 58 ## nc layout, with horizontal and vertical distances between them 59 ## set to space. 60 ## 61 ## The drawing code first computes the complete are of the form 62 ## [0, totalWidth] x [0, totalHeight] 63 ## needed and sets the world coordinates using plot.window(). 64 ## Then, the strata are looped over, and the corresponding pies 65 ## added by filling rows or columns of the layout as specified by 66 ## the mfrow or mfcol arguments. The world coordinates are reset 67 ## in each step by shifting the origin so that we can always plot 68 ## as detailed above. 69 70 if(!is.array(x)) 71 stop("x must be an array") 72 dimx <- dim(x) # save original dimensions for setting default mfrow/mfcol when length(dim(x))>3 73 x <- array3d(x) 74 if(any(dim(x)[1:2] != 2)) 75 stop("table for each stratum must be 2 by 2") 76 dnx <- dimnames(x) 77 if(is.null(dnx)) 78 dnx <- vector("list", 3) 79 for(i in which(sapply(dnx, is.null))) 80 dnx[[i]] <- LETTERS[seq(length = dim(x)[i])] 81 if(is.null(names(dnx))) 82 i <- 1 : 3 83 else 84 i <- which(is.null(names(dnx))) 85 if(any(i > 0)) 86 names(dnx)[i] <- default_prefix[i] 87 dimnames(x) <- dnx 88 k <- dim(x)[3] 89 90 if(!((length(conf_level) == 1) && is.finite(conf_level) && 91 (conf_level >= 0) && (conf_level < 1))) 92 stop("conf_level must be a single number between 0 and 1") 93 if(conf_level == 0) 94 conf_level <- FALSE 95 96 std <- match.arg(std) 97 98 findTableWithOAM <- function(or, tab) { 99 ## Find a 2x2 table with given odds ratio `or' and the margins 100 ## of a given 2x2 table `tab'. 101 m <- rowSums(tab)[1] 102 n <- rowSums(tab)[2] 103 t <- colSums(tab)[1] 104 if(or == 1) 105 x <- t * n / (m + n) 106 else if(or == Inf) 107 x <- max(0, t - m) 108 else { 109 A <- or - 1 110 B <- or * (m - t) + (n + t) 111 C <- - t * n 112 x <- (- B + sqrt(B ^ 2 - 4 * A * C)) / (2 * A) 113 } 114 matrix(c(t - x, x, m - t + x, n - x), nrow = 2) 115 } 116 117 drawPie <- function(r, from, to, n = 500, color = "transparent") { 118 p <- 2 * pi * seq(from, to, length = n) / 360 119 x <- c(cos(p), 0) * r 120 y <- c(sin(p), 0) * r 121 grid.polygon(x, y, gp = gpar(fill = color), default.units = "native") 122 invisible(NULL) 123 } 124 125 stdize <- function(tab, std, x) { 126 ## Standardize the 2 x 2 table `tab'. 127 if(std == "margins") { 128 if(all(sort(margin) == c(1, 2))) { 129 ## standardize to equal row and col margins 130 u <- sqrt(odds(tab)$or) 131 u <- u / (1 + u) 132 y <- matrix(c(u, 1 - u, 1 - u, u), nrow = 2) 133 } 134 else if(margin %in% c(1, 2)) 135 y <- prop.table(tab, margin) 136 else 137 stop("incorrect margin specification") 138 } 139 else if(std == "ind.max") 140 y <- tab / max(tab) 141 else if(std == "all.max") 142 y <- tab / max(x) 143 y 144 } 145 146 odds <- function(x) { 147 ## Given a 2 x 2 or 2 x 2 x k table `x', return a list with 148 ## components `or' and `se' giving the odds ratios and standard 149 ## deviations of the log odds ratios. 150 if(length(dim(x)) == 2) { 151 dim(x) <- c(dim(x), 1) 152 k <- 1 153 } 154 else 155 k <- dim(x)[3] 156 or <- double(k) 157 se <- double(k) 158 for(i in 1 : k) { 159 f <- x[ , , i] 160 if(any(f == 0)) 161 f <- f + 0.5 162 or[i] <- (f[1, 1] * f[2, 2]) / (f[1, 2] * f[2, 1]) 163 se[i] <- sqrt(sum(1 / f)) 164 } 165 list(or = or, se = se) 166 } 167 168 gamma <- 1.25 # Scale factor for strata labels 169 170 angle.f <- c( 90, 180, 0, 270) # `f' for `from' 171 angle.t <- c(180, 270, 90, 360) # `t' for `to' 172 173 byrow <- FALSE 174 if(!is.null(mfrow)) { 175 nr <- mfrow[1] 176 nc <- mfrow[2] 177 } 178 else if(!is.null(mfcol)) { 179 nr <- mfcol[1] 180 nc <- mfcol[2] 181 byrow <- TRUE 182 } 183 else if(length(dimx)>3) { 184 nr <- dimx[3] 185 nc <- prod(dimx[-(1:3)]) 186 } 187 else { 188# nr <- ceiling(sqrt(k)) 189 nr <- round(sqrt(k)) 190 nc <- ceiling(k / nr) 191 } 192 if(nr * nc < k) 193 stop("incorrect geometry specification") 194 if(byrow) 195 indexMatrix <- expand.grid(1 : nc, 1 : nr)[, c(2, 1)] 196 else 197 indexMatrix <- expand.grid(1 : nr, 1 : nc) 198 199 totalWidth <- nc * 2 * (1 + space) + (nc - 1) * space 200 totalHeight <- if(k == 1) 201 2 * (1 + space) 202 else 203 nr * (2 + (2 + gamma) * space) + (nr - 1) * space 204 xlim <- c(0, totalWidth) 205 ylim <- c(0, totalHeight) 206 if (newpage) grid.newpage() 207 if (!is.null(main) || !is.null(sub)) 208 pushViewport(viewport(height = 1 - 0.1 * sum(!is.null(main), !is.null(sub)), 209 width = 0.9, 210 y = 0.5 - 0.05 * sum(!is.null(main), - !is.null(sub)) 211 ) 212 ) 213 214 pushViewport(viewport(xscale = xlim, yscale = ylim, 215 width = unit(min(totalWidth / totalHeight, 1), "snpc"), 216 height = unit(min(totalHeight / totalWidth, 1), "snpc"))) 217 o <- odds(x) 218 219 ## perform logoddsratio-test for each stratum (H0: lor = 0) and adjust p-values 220 if(is.numeric(conf_level) && extended) 221 p.lor.test <- p.adjust(sapply(1 : k, function(i) { 222 u <- abs(log(o$or[i])) / o$se[i] 223 2 * (1 - pnorm(u)) 224 }), 225 method = p_adjust_method 226 ) 227 228 scale <- space / (2 * convertY(unit(1, "strheight", "Ag"), "native", valueOnly = TRUE) ) 229 v <- 0.95 - max(convertX(unit(1, "strwidth", as.character(c(x))), "native", valueOnly = TRUE) ) / 2 230 231 fontsize = fontsize * scale 232 233 for(i in 1 : k) { 234 235 tab <- x[ , , i] 236 237 fit <- stdize(tab, std, x) 238 239 xInd <- indexMatrix[i, 2] 240 xOrig <- 2 * xInd - 1 + (3 * xInd - 2) * space 241 yInd <- indexMatrix[i, 1] 242 yOrig <- if(k == 1) 243 (1 + space) 244 else 245 (totalHeight 246 - (2 * yInd - 1 + ((3 + gamma) * yInd - 2) * space)) 247 pushViewport(viewport(xscale = xlim - xOrig, yscale = ylim - yOrig)) 248 249 ## drawLabels() 250 u <- 1 + space / 2 251 adjCorr <- 0.2 252 grid.text( 253 paste(names(dimnames(x))[1], 254 dimnames(x)[[1]][1], 255 sep = sep), 256 0, u, 257 gp = gpar(fontsize = fontsize), 258 default.units = "native" 259 ) 260 grid.text( 261 paste(names(dimnames(x))[2], 262 dimnames(x)[[2]][1], 263 sep = sep), 264 -u, 0, 265 default.units = "native", 266 gp = gpar(fontsize = fontsize), 267 rot = 90) 268 grid.text( 269 paste(names(dimnames(x))[1], 270 dimnames(x)[[1]][2], 271 sep = sep), 272 0, -u, 273 gp = gpar(fontsize = fontsize), 274 default.units = "native" 275 ) 276 grid.text( 277 paste(names(dimnames(x))[2], 278 dimnames(x)[[2]][2], 279 sep = sep), 280 u, 0, 281 default.units = "native", 282 gp = gpar(fontsize = fontsize), 283 rot = 90) 284 if (k > 1) { 285 grid.text(if (!varnames) 286 dimnames(x)[[3]][i] 287 else 288 paste(names(dimnames(x))[3], 289 dimnames(x)[[3]][i], 290 sep = sep), 291 0, 1 + (1 + gamma / 2) * space, 292 gp = gpar(fontsize = fontsize * gamma), 293 default.units = "native" 294 ) 295 } 296 297 ## drawFrequencies() 298 299 ### in extended plots, emphasize charts with significant logoddsratios 300 emphasize <- if(extended && is.numeric(conf_level)) 301 2 * extended * (1 + (p.lor.test[i] < 1 - conf_level)) 302 else 0 303 304 d <- odds(tab)$or 305 drawPie(sqrt(fit[1,1]), 90, 180, col = color[1 + (d > 1) + emphasize]) 306 drawPie(sqrt(fit[2,1]), 180, 270, col = color[2 - (d > 1) + emphasize]) 307 drawPie(sqrt(fit[1,2]), 0, 90, col = color[2 - (d > 1) + emphasize]) 308 drawPie(sqrt(fit[2,2]), 270, 360, col = color[1 + (d > 1) + emphasize]) 309 310 u <- 1 - space / 2 311 grid.text(as.character(c(tab))[1], 312 -v, u, 313 just = c("left", "top"), 314 gp = gpar(fontsize = fontsize), 315 default.units = "native") 316 grid.text(as.character(c(tab))[2], 317 -v, -u, 318 just = c("left", "bottom"), 319 gp = gpar(fontsize = fontsize), 320 default.units = "native") 321 grid.text(as.character(c(tab))[3], 322 v, u, 323 just = c("right", "top"), 324 gp = gpar(fontsize = fontsize), 325 default.units = "native") 326 grid.text(as.character(c(tab))[4], 327 v, -u, 328 just = c("right", "bottom"), 329 gp = gpar(fontsize = fontsize), 330 default.units = "native") 331 332 ## draw ticks 333 if(extended && ticks) 334 if(d > 1) { 335 grid.lines(c(sqrt(fit[1,1]) * cos(3*pi/4), 336 (sqrt(fit[1,1]) + ticks) * cos(3*pi/4)), 337 c(sqrt(fit[1,1]) * sin(3*pi/4), 338 (sqrt(fit[1,1]) + ticks) * sin(3*pi/4)), 339 gp = gpar(lwd = 1), 340 default.units = "native" 341 ) 342 grid.lines(c(sqrt(fit[2,2]) * cos(-pi/4), 343 (sqrt(fit[2,2]) + ticks) * cos(-pi/4)), 344 c(sqrt(fit[2,2]) * sin(-pi/4), 345 (sqrt(fit[2,2]) + ticks) * sin(-pi/4)), 346 gp = gpar(lwd = 1), 347 default.units = "native" 348 ) 349 } else { 350 grid.lines(c(sqrt(fit[1,2]) * cos(pi/4), 351 (sqrt(fit[1,2]) + ticks) * cos(pi/4)), 352 c(sqrt(fit[1,2]) * sin(pi/4), 353 (sqrt(fit[1,2]) + ticks) * sin(pi/4)), 354 gp = gpar(lwd = 1), 355 default.units = "native" 356 ) 357 grid.lines(c(sqrt(fit[2,1]) * cos(-3*pi/4), 358 (sqrt(fit[2,1]) + ticks) * cos(-3*pi/4)), 359 c(sqrt(fit[2,1]) * sin(-3*pi/4), 360 (sqrt(fit[2,1]) + ticks) * sin(-3*pi/4)), 361 gp = gpar(lwd = 1), 362 default.units = "native" 363 ) 364 } 365 366 ## drawConfBands() 367 if(is.numeric(conf_level)) { 368 or <- o$or[i] 369 se <- o$se[i] 370 ## lower 371 theta <- or * exp(qnorm((1 - conf_level) / 2) * se) 372 tau <- findTableWithOAM(theta, tab) 373 r <- sqrt(c(stdize(tau, std, x))) 374 for(j in 1 : 4) 375 drawPie(r[j], angle.f[j], angle.t[j]) 376 ## upper 377 theta <- or * exp(qnorm((1 + conf_level) / 2) * se) 378 tau <- findTableWithOAM(theta, tab) 379 r <- sqrt(c(stdize(tau, std, x))) 380 for(j in 1 : 4) 381 drawPie(r[j], angle.f[j], angle.t[j]) 382 } 383 384 ## drawBoxes() 385 grid.polygon(c(-1, 1, 1, -1), 386 c(-1, -1, 1, 1), 387 default.units = "native", 388 gp = gpar(fill = "transparent") 389 ) 390 grid.lines(c(-1, 1), c(0, 0), default.units = "native") 391 for(j in seq(from = -0.8, to = 0.8, by = 0.2)) 392 grid.lines(c(j, j), c(-0.02, 0.02), default.units = "native") 393 for(j in seq(from = -0.9, to = 0.9, by = 0.2)) 394 grid.lines(c(j, j), c(-0.01, 0.01), default.units = "native") 395 grid.lines(c(0, 0), c(-1, 1), default.units = "native") 396 for(j in seq(from = -0.8, to = 0.8, by = 0.2)) 397 grid.lines(c(-0.02, 0.02), c(j, j), default.units = "native") 398 for(j in seq(from = -0.9, to = 0.9, by = 0.2)) 399 grid.lines(c(-0.01, 0.01), c(j, j), default.units = "native") 400 401 popViewport(1) 402 403 } 404 405 if(!is.null(main) || !is.null(sub)) { 406 if (!is.null(main)) 407 grid.text(main, 408 y = unit(1, "npc") + unit(1, "lines"), 409 gp = gpar(fontsize = 20, fontface = 2)) 410 if (!is.null(sub)) 411 grid.text(sub, 412 y = unit(0, "npc") - unit(1, "lines"), 413 gp = gpar(fontsize = 20, fontface = 2)) 414 popViewport(1) 415 } 416 417 popViewport(1) 418 419 if (return_grob) 420 return(invisible(grid.grab())) 421 else 422 return(invisible(NULL)) 423} 424