| blob | 0b93219f424ec2429b7bd1c2301e539f2f9fd0f7 |
1 ################################################################################
2 ## Copyright (C) 2008 Roger D. Peng <rpeng@jhsph.edu>
3 ##
4 ## This program is free software; you can redistribute it and/or modify
5 ## it under the terms of the GNU General Public License as published by
6 ## the Free Software Foundation; either version 2 of the License, or
7 ## (at your option) any later version.
8 ##
9 ## This program is distributed in the hope that it will be useful,
10 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
11 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 ## GNU General Public License for more details.
13 ##
14 ## You should have received a copy of the GNU General Public License
15 ## along with this program; if not, write to the Free Software
16 ## Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
17 ## 02110-1301, USA
18 ################################################################################
20 library(splines)
23 mvtsplot <- local({
24 drawImage <- function(cx, pal, nlevels, xlim, cn, xtime, group, gcol) {
25 par(las = 1, cex.axis = 0.6)
26 cn <- colnames(cx)
27 nc <- ncol(cx)
28 side2 <- 0.2
30 ## Setup image plot
31 if(!is.null(cn))
32 side2 <- max(side2, max(strwidth(cn, "inches")) + 0.1)
33 else
34 cn <- as.character(1, nc)
35 par(mai = c(0.4, side2, 0.1, 0.1))
36 image(unclass(xtime), seq_len(nc), cx, col = pal(nlevels),
37 xlim = xlim, xaxt = "n", yaxt = "n", ylab = "", xlab = "")
38 axis(2, at = seq_len(nc), cn)
39 Axis(xtime, side = 1)
40 box()
42 if(!is.null(group)) {
43 usrpar <- par("usr")
44 par(usr = c(usrpar[1:2], 0, 1))
45 tg <- table(group)[-nlevels(group)]
46 abline(h = cumsum(tg) / nc, lwd = 2, col = gcol)
47 }
48 }
50 drawImageMargin <- function(cx, pal, nlevels, xlim, cn, xtime, group,
51 gcol, smooth.df, rowm, nr, bottom.ylim, colm, right.xlim,
52 main) {
53 op <- par(no.readonly = TRUE)
54 on.exit(par(op))
55 par(las = 1, cex.axis = 0.6)
57 cn <- colnames(cx)
58 nc <- ncol(cx)
59 side2 <- 0.55
60 utime <- unclass(xtime)
62 layout(rbind(c(1, 1, 1, 1, 1, 1, 3),
63 c(1, 1, 1, 1, 1, 1, 3),
64 c(1, 1, 1, 1, 1, 1, 3),
65 c(2, 2, 2, 2, 2, 2, 4)))
67 ## Setup image plot
68 if(!is.null(cn))
69 side2 <- max(0.55, max(strwidth(cn, "inches")) + 0.1)
70 else
71 cn <- rep("", nc)
72 par(mai = c(0.05, side2, 0.1, 0.05))
73 image(utime, seq_len(nc), cx, col = pal(nlevels),
74 xlim = xlim, xaxt = "n", yaxt = "n", ylab = "", xlab = "")
75 axis(2, at = seq_len(nc), cn)
76 box()
78 if(!is.null(group)) {
79 usrpar <- par("usr")
80 par(usr = c(usrpar[1:2], 0, 1))
81 tg <- table(group)[-nlevels(group)]
82 abline(h = cumsum(tg) / nc, lwd = 2, col = gcol)
83 }
84 ## Plot bottom
85 if(!is.null(smooth.df)) { ## Smooth row stats
86 xx <- seq_along(rowm)
87 tmp.fit <- lm(rowm ~ ns(xx, smooth.df),
88 na.action = na.exclude)
89 rowm <- predict(tmp.fit)
90 }
91 bottom.ylim <- if(is.null(bottom.ylim))
92 range(rowm, na.rm = TRUE)
93 else
94 bottom.ylim
95 par(mai = c(0.4, side2, 0.05, 0.05))
96 plot(utime, rep(0, nr), type = "n",
97 ylim = bottom.ylim, xaxt = "n", xlab = "", ylab = "Level")
98 par(usr = c(xlim, par("usr")[3:4]))
99 nalines(utime, rowm)
100 Axis(xtime, side = 1)
102 ## Plot right side
103 right.xlim <- if(is.null(right.xlim))
104 range(colm, na.rm = TRUE)
105 else
106 right.xlim
107 par(mai = c(0.05, 0.05, 0.1, 0.1))
108 plot(colm[3,], 1:nc, type = "n", ylab = "", yaxt = "n",
109 xlab = "", xlim = right.xlim)
111 usrpar <- par("usr")
112 par(usr = c(usrpar[1:2], 0, 1))
113 ypos <- (1:nc - 1 / 2) / nc
115 segments(colm[1, ], ypos, colm[2, ], ypos, col = gray(0.6))
116 segments(colm[4, ], ypos, colm[5, ], ypos, col = gray(0.6))
117 points(colm[3,], ypos, pch = 19, cex = 0.6)
119 if(!is.null(group))
120 abline(h = cumsum(tg) / nc, lwd = 2, col = gcol)
122 ## Plot lower right
123 blankplot()
124 text(0, 0, main)
126 rval <- list(z = cx, rowm = rowm, colm = colm)
127 invisible(rval)
128 }
130 blankplot <- function() {
131 plot(0, 0, xlab = "", ylab = "", axes = FALSE, type = "n")
132 }
134 ## Discretize columns of a matrix in to categories defined by 'levels'
136 catcols <- function(x, levels, norm) {
137 if(norm == "internal") {
138 ## Each column gets its own set of categories
139 apply(x, 2, function(y) categorize(y, levels))
140 }
141 else {
142 ## Use a 'global' set of categories
143 xv <- as.vector(x)
144 y <- categorize(xv, levels)
145 matrix(y, nrow = nrow(x), ncol = ncol(x))
146 }
147 }
149 ## Discretize a vector 'x' into categories defined by 'levels'.
151 categorize <- function(x, levels, jitter = TRUE) {
152 ## 'x' is a vector; 'levels' is a single integer, or a vector
153 ## of quantiles
154 if(length(levels) == 1)
155 levels <- seq(0, 1, len = levels + 1)
156 qq <- quantile(x, levels, na.rm = TRUE)
157 qqu <- unique(qq)
159 if(length(qqu) != length(qq)) {
160 if(!jitter)
161 return(rep(NA, length(x)))
162 qqu <- try(suppressWarnings({
163 x <- jitter(x)
164 qq <- quantile(x, levels, na.rm = TRUE)
165 unique(qq)
166 }), silent = TRUE)
167 if(inherits(qqu, "try-error") || length(qqu) != length(qq))
168 return(rep(NA, length(x)))
169 }
170 cx <- cut(x, qqu, include.lowest = TRUE)
171 as.numeric(unclass(cx))
172 }
174 smoothX <- function(x, df) {
175 apply(x, 2, function(v) splineFillIn(v, df))
176 }
178 reorderCols <- function(x, group) {
179 if(length(group) != ncol(x))
180 stop("'group' vector should equal 'ncol(x)'")
181 idx.split <- split(seq_len(ncol(x)), group)
182 idx.reordered <- unlist(idx.split)
183 x[, idx.reordered, drop = FALSE]
184 }
186 rangeby <- function(x, f) {
187 use <- complete.cases(x, f)
188 range(x[use])
189 }
191 splineFillIn <- function(x, df) {
192 if(all(is.na(x)))
193 return(x)
194 tt <- seq_along(x)
195 rng <- rangeby(tt, x)
196 fit <- lm(x ~ ns(tt, df), na.action = na.exclude)
197 idx <- seq(rng[1], rng[2])
198 x[idx] <- suppressWarnings({
199 predict(fit, data.frame(tt = idx))
200 })
201 x
202 }
204 sumNA <- function(x) {
205 if(all(is.na(x)))
206 NA
207 else
208 sum(x, na.rm = TRUE)
209 }
211 checkMatrix <- function(x) {
212 if(!is.matrix(x))
213 stop("'x' should be a matrix")
214 if(ncol(x) < 2)
215 stop("'x' should have more than 1 column")
216 if(nrow(x) < 2)
217 stop("'x' should have more than 1 row")
218 TRUE
219 }
221 ## This function is like 'lines' in that it draws lines between
222 ## points. For two points that are consecutive, the line is drawn
223 ## "black". If one or more missing values is between two points, then
224 ## the line is drawn grey (or 'NAcol').
226 nalines <- function(x, y, NAcol = gray(0.6), lattice = FALSE, ...) {
227 use <- complete.cases(x, y)
228 idx <- which(use)
229 n <- length(idx)
231 if(n < 2)
232 return(invisible())
233 for(i in seq_len(n - 1)) {
234 j <- idx[i]
235 k <- idx[i+1]
236 col <- if((k - j) > 1)
237 NAcol
238 else
239 "black"
240 if(lattice)
241 llines(c(x[j], x[k]), c(y[j], y[k]), col = col)
242 else
243 lines(c(x[j], x[k]), c(y[j], y[k]), col = col)
244 }
245 invisible()
246 }
249 function(x, group = NULL, xtime = NULL,
250 norm = c("internal", "global"), levels = 3,
251 smooth.df = NULL, margin = TRUE,
252 sort = NULL,
253 main = "", palette = "PRGn",
254 rowstat = "median",
255 xlim, bottom.ylim = NULL,
256 right.xlim = NULL,
257 gcol = 1) {
258 if(is.data.frame(x))
259 x <- data.matrix(x)
260 checkMatrix(x)
261 norm <- match.arg(norm)
263 if(!is.null(sort))
264 sort <- match.fun(sort)
265 rowstat <- match.fun(rowstat)
267 if(!require(RColorBrewer))
268 stop("'RColorBrewer' package required")
269 if(is.null(xtime)) {
270 xtime <- seq_len(nrow(x))
271 xlim <- c(0, max(xtime))
272 }
273 else
274 xlim <- range(xtime)
275 if(!is.null(group)) {
276 group <- as.factor(group)
277 x <- reorderCols(x, group)
278 }
279 if(!margin && !is.null(sort)) {
280 stat <- apply(x, 2, sort, na.rm = TRUE)
281 x <- x[, order(stat)]
282 }
283 if(margin) {
284 colm <- apply(x, 2, function(x) {
285 grDevices::boxplot.stats(x)$stats
286 })
287 if(!is.null(sort)) {
288 stat <- apply(x, 2, sort, na.rm = TRUE)
289 ord <- order(stat)
290 x <- x[, ord]
291 colm <- colm[, ord]
292 }
293 }
294 if(is.null(smooth.df))
295 cx <- catcols(x, levels, norm)
296 else {
297 x <- smoothX(x, smooth.df)
298 cx <- catcols(x, levels, norm)
299 }
300 if(margin)
301 rowm <- apply(x, 1, rowstat, na.rm = TRUE)
302 colnames(cx) <- colnames(x)
303 empty <- apply(cx, 2, function(x) all(is.na(x)))
305 if(any(empty)) { ## Remove empty columns
306 cx <- cx[, !empty]
308 if(margin)
309 colm <- colm[, !empty]
310 }
311 pal <- colorRampPalette(brewer.pal(4, palette))
313 nlevels <- if(length(levels) == 1)
314 levels
315 else
316 length(levels)
317 if(margin)
318 drawImageMargin(cx, pal, nlevels, xlim, cn, xtime,
319 group, gcol, smooth.df, rowm, nrow(x),
320 bottom.ylim, colm, right.xlim, main)
321 else
322 drawImage(cx, pal, nlevels, xlim, cn, xtime, group, gcol)
323 }
325 })