1contourPlot3 = function(x, y, z, response, data = NULL, main, xlab, ylab, zlab, border, form = "linear", col = 1, col.text, cex.axis, axes = TRUE, 2 steps, factors) { 3 DB = FALSE 4 out = list() 5 mdo = data 6 x.c = deparse(substitute(x)) 7 y.c = deparse(substitute(y)) 8 z.c = deparse(substitute(z)) 9 r.c = deparse(substitute(response)) 10 if (missing(col)) 11 col = 1 12 if (missing(col.text)) 13 col.text = 1 14 if (missing(cex.axis)) 15 cex.axis = 1 16 if (missing(main)) 17 main = paste("Response Surface for", r.c) 18 if (missing(ylab)) 19 ylab = y.c 20 if (missing(xlab)) 21 xlab = x.c 22 if (missing(zlab)) 23 zlab = z.c 24 if (missing(border)) 25 border = "white" 26 if (missing(factors)) 27 factors = NULL 28 if (missing(steps)) 29 steps = 100 30 col.axis = par("col.axis") 31 if (!is.function(col)) { 32 if (identical(col, 1)) 33 col = colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000")) 34 if (identical(col, 2)) 35 col = colorRampPalette(c("blue", "white", "red"), space = "Lab") 36 if (identical(col, 3)) 37 col = colorRampPalette(c("blue", "white", "orange")) 38 if (identical(col, 4)) 39 col = colorRampPalette(c("gold", "white", "firebrick")) 40 } 41 nameVec = names(names(mdo)) 42 linStrings = "-1" 43 for (i in seq(along = nameVec)) linStrings = paste(linStrings, "+", nameVec[i]) 44 if (DB) 45 print(linStrings) 46 combList = combn(nameVec, 2, simplify = FALSE) 47 quadStrings = character(length = length(combList)) 48 for (i in seq(along = combList)) if (i == 1) 49 quadStrings[i] = paste(combList[[i]][1], ":", combList[[i]][2]) 50 else quadStrings[i] = paste("+", combList[[i]][1], ":", combList[[i]][2]) 51 quadStrings = paste(quadStrings, collapse = "") 52 if (DB) 53 print(quadStrings) 54 if (identical(form, "linear")) { 55 form = paste(r.c, "~", linStrings) 56 if (DB) 57 print(form) 58 } 59 if (identical(form, "quadratic")) { 60 form = paste(r.c, "~", linStrings, "+", quadStrings) 61 if (DB) 62 print(form) 63 } 64 lm.1 = lm(formula = form, data = mdo) 65 if (DB) 66 print(lm.1) 67 dcList = vector(mode = "list", length = length(names(mdo))) 68 names(dcList) = names(names(mdo)) 69 dcList[1:length(names(mdo))] = 0 70 if (!is.null(factors)) { 71 for (i in names(factors)) dcList[[i]] = factors[[i]][1] 72 } 73 if (DB) 74 print(dcList) 75 help.predict = function(a, b, x.c, y.c, lm.1, ...) { 76 dcList[[x.c]] = 2 * b/sqrt(3) 77 dcList[[y.c]] = 1 - (2 * b/sqrt(3)) - (a - b/sqrt(3)) 78 dcList[[z.c]] = a - b/sqrt(3) 79 temp = do.call(data.frame, dcList) 80 invisible(predict(lm.1, temp)) 81 } 82 a = seq(0, 1, length = steps) 83 b = seq(0, sqrt(3)/2, length = steps) 84 mat = outer(a, b, help.predict, x.c, y.c, lm.1) 85 acc = nrow(mat) 86 v = seq(acc, 1, length = acc) 87 w = c(seq(2, acc, length = acc/2), seq(acc, 2, length = acc/2)) 88 mat[outer(w, v, `+`) <= acc] = NA 89 .mfc(mat, main = main, col = col, axes = FALSE, key.axes = axis(4)) 90 if (axes == TRUE) { 91 segments(0.5, 0, 0.5, 1, col = col.axis) 92 coox1 = rep(0.49, 11) 93 coox2 = rep(0.51, 11) 94 cooy = seq(0, 1, length = 11) 95 for (i in 2:10) { 96 segments(coox1[i], cooy[i], coox2[i], cooy[i], col = col.axis) 97 text(coox2[i] + 0.01, cooy[i], labels = (i - 1)/10, cex = cex.axis, col = col.text) 98 } 99 segments(0, 0, 0.75, 0.5, col = col.axis) 100 coox1 = seq(0.745, -0.005, length = 11) 101 coox2 = seq(0.755, 0.005, length = 11) 102 cooy1 = seq(0.51, 0.01, length = 11) 103 cooy2 = seq(0.49, -0.01, length = 11) 104 for (i in 2:10) { 105 segments(coox1[i], cooy1[i], coox2[i], cooy2[i], col = col.axis) 106 text(coox2[i], cooy1[i] + 0.01, labels = (i - 1)/10, cex = cex.axis, col = col.text) 107 } 108 segments(0.25, 0.5, 1, 0, col = col.axis) 109 coox1 = seq(0.245, 0.995, length = 11) 110 coox2 = seq(0.255, 1.005, length = 11) 111 cooy1 = seq(0.49, -0.01, length = 11) 112 cooy2 = seq(0.51, 0.01, length = 11) 113 for (i in 2:10) { 114 segments(coox1[i], cooy1[i], coox2[i], cooy2[i], col = col.axis) 115 text(coox1[i] + 0.01, cooy2[i] + 0.01, labels = (i - 1)/10, cex = cex.axis, col = col.text) 116 } 117 } 118 segments(-0.005, 0, 0.495, 1, lwd = 5, col = border) 119 segments(0.505, 1, 1.005, 0, lwd = 5, col = border) 120 segments(1, -0.005, 0, -0.005, lwd = 5, col = border) 121 mtext(ylab, 1, at = -0.025, cex = 1.5) 122 mtext(xlab, 3, at = 0.5, cex = 1.5, line = 0.1) 123 mtext(zlab, 1, at = 1.025, cex = 1.5) 124 invisible(mat) 125} 126wirePlot3 = function(x, y, z, response, data = NULL, main, xlab, ylab, zlab, form = "linear", phi, theta, col = 1, steps, factors) { 127 DB = FALSE 128 out = list() 129 mdo = data 130 x.c = deparse(substitute(x)) 131 y.c = deparse(substitute(y)) 132 z.c = deparse(substitute(z)) 133 r.c = deparse(substitute(response)) 134 if (missing(col)) 135 col = 1 136 if (missing(main)) 137 main = paste("Response Surface for", r.c) 138 if (missing(ylab)) 139 ylab = y.c 140 if (missing(xlab)) 141 xlab = x.c 142 if (missing(zlab)) 143 zlab = z.c 144 if (missing(phi)) 145 phi = 30 146 if (missing(theta)) 147 theta = 30 148 if (missing(factors)) 149 factors = NULL 150 if (missing(steps)) 151 steps = 100 152 if (!is.function(col)) { 153 if (identical(col, 1)) 154 col = colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000")) 155 if (identical(col, 2)) 156 col = colorRampPalette(c("blue", "white", "red"), space = "Lab") 157 if (identical(col, 3)) 158 col = colorRampPalette(c("blue", "white", "orange")) 159 if (identical(col, 4)) 160 col = colorRampPalette(c("gold", "white", "firebrick")) 161 } 162 phi = phi%%360 163 .phi = phi 164 theta = theta%%360 165 .theta = theta 166 nameVec = names(names(mdo)) 167 linStrings = "-1" 168 for (i in seq(along = nameVec)) linStrings = paste(linStrings, "+", nameVec[i]) 169 if (DB) 170 print(linStrings) 171 combList = combn(nameVec, 2, simplify = FALSE) 172 quadStrings = character(length = length(combList)) 173 for (i in seq(along = combList)) if (i == 1) 174 quadStrings[i] = paste(combList[[i]][1], ":", combList[[i]][2]) 175 else quadStrings[i] = paste("+", combList[[i]][1], ":", combList[[i]][2]) 176 quadStrings = paste(quadStrings, collapse = "") 177 if (DB) 178 print(quadStrings) 179 if (identical(form, "linear")) { 180 form = paste(r.c, "~", linStrings) 181 if (DB) 182 print(form) 183 } 184 if (identical(form, "quadratic")) { 185 form = paste(r.c, "~", linStrings, "+", quadStrings) 186 if (DB) 187 print(form) 188 } 189 lm.1 = lm(formula = form, data = mdo) 190 if (DB) 191 print(lm.1) 192 dcList = vector(mode = "list", length = length(names(mdo))) 193 names(dcList) = names(names(mdo)) 194 dcList[1:length(names(mdo))] = 0 195 if (!is.null(factors)) { 196 for (i in names(factors)) dcList[[i]] = factors[[i]][1] 197 } 198 if (DB) 199 print(dcList) 200 help.predict = function(a, b, x.c, y.c, lm.1, ...) { 201 dcList[[x.c]] = 2 * b/sqrt(3) 202 dcList[[y.c]] = 1 - (2 * b/sqrt(3)) - (a - b/sqrt(3)) 203 dcList[[z.c]] = a - b/sqrt(3) 204 temp = do.call(data.frame, dcList) 205 invisible(predict(lm.1, temp)) 206 } 207 a = seq(0, 1, length = steps) 208 b = seq(0, sqrt(3)/2, length = steps) 209 mat = outer(a, b, help.predict, x.c, y.c, lm.1) 210 acc = nrow(mat) 211 sca = sin(1/3 * pi) 212 ncmat = ncol(mat) 213 v = seq(acc, 1, length = acc) 214 w = c(seq(2, acc, length = acc/2), seq(acc, 2, length = acc/2)) 215 mat[outer(w, v, `+`) <= acc] = NA 216 if (is.function(col)) { 217 nrMat <- nrow(mat) 218 ncMat <- ncol(mat) 219 nbcol <- 100 220 color <- col(nbcol) 221 matFacet = mat[-1, -1] + mat[-1, -ncmat] + mat[-acc, -1] + mat[-acc, -ncmat] 222 facetcol <- cut(matFacet, nbcol) 223 } 224 else { 225 color = col 226 facetcol = 1 227 } 228 maxim = max(mat, na.rm = TRUE) * acc 229 minim = min(mat, na.rm = TRUE) * acc 230 per = persp(x = seq(0, acc, length = acc), y = seq(0, acc * sca, length = ncmat), mat * acc, phi = .phi, theta = .theta, scale = TRUE, col = "transparent", 231 border = FALSE, box = FALSE, main = main, xlab = xlab, ylab = ylab) 232 lineList = contourLines(x = seq(0, acc, length = acc), y = seq(0, acc * sca, length = ncmat), mat) 233 for (i in seq(along = lineList)) lines(trans3d(lineList[[i]]$x, lineList[[i]]$y, z = minim, pmat = per)) 234 if (.phi < 90) { 235 lines(trans3d(x = seq(0, acc/2, length = 10), y = seq(0, acc * sca, length = 10), z = maxim, pmat = per), lty = 2) 236 lines(trans3d(x = seq(acc, acc/2, length = 10), y = seq(0, acc * sca, length = 10), z = maxim, pmat = per), lty = 2) 237 lines(trans3d(x = 0:acc, y = 0, z = maxim, pmat = per), lty = 2) 238 } 239 if (.theta > 323 || .theta < 37) { 240 lines(trans3d(x = acc/2, y = acc * sca, z = minim:maxim, pmat = per), lty = 2) 241 lines(trans3d(x = 0, y = 0, z = minim:maxim, pmat = per), lty = 2) 242 lines(trans3d(x = acc, y = 0, z = minim:maxim, pmat = per), lty = 2) 243 } 244 if (.theta > 37 && .theta < 156) 245 lines(trans3d(x = 0, y = 0, z = minim:maxim, pmat = per), lty = 2) 246 if (.theta > 156 && .theta < 323) { 247 lines(trans3d(x = acc, y = 0, z = minim:maxim, pmat = per), lty = 2) 248 } 249 lines(trans3d(x = seq(0, acc/2, length = 10), y = seq(0, acc * sca, length = 10), z = minim, pmat = per), lty = 1, lwd = 2) 250 lines(trans3d(x = seq(acc, acc/2, length = 10), y = seq(0, acc * sca, length = 10), z = minim, pmat = per), lty = 1, lwd = 2) 251 lines(trans3d(x = 0:acc, y = 0, z = minim, pmat = per), lty = 1, lwd = 2) 252 text(trans3d(x = acc/2 + acc/50, y = acc * sca + acc * sca/50, z = minim, pmat = per), labels = xlab, lwd = 2) 253 text(trans3d(x = -acc/50, y = -acc * sca/50, z = minim, pmat = per), labels = ylab, lwd = 2) 254 text(trans3d(x = acc + acc/50, 0, z = minim, pmat = per), labels = zlab, cex = 1, lwd = 2) 255 par(new = TRUE) 256 persp(x = seq(0, acc, length = acc), y = seq(0, acc * sca, length = ncmat), mat * acc, phi = .phi, theta = .theta, scale = TRUE, col = color[facetcol], 257 border = FALSE, box = FALSE) 258 if (.phi > 0) { 259 lines(trans3d(x = seq(0, acc/2, length = 10), y = seq(0, acc * sca, length = 10), z = maxim, pmat = per), lty = 2) 260 lines(trans3d(x = seq(acc, acc/2, length = 10), y = seq(0, acc * sca, length = 10), z = maxim, pmat = per), lty = 2) 261 lines(trans3d(x = 0:acc, y = 0, z = maxim, pmat = per), lty = 2) 262 } 263 if (.theta > 37 && .theta < 156) { 264 lines(trans3d(x = acc/2, y = acc * sca, z = minim:maxim, pmat = per), lty = 2) 265 lines(trans3d(x = acc, y = 0, z = minim:maxim, pmat = per), lty = 2) 266 } 267 if (.theta > 156 && .theta < 323) { 268 lines(trans3d(x = acc/2, y = acc * sca, z = minim:maxim, pmat = per), lty = 2) 269 lines(trans3d(x = 0, y = 0, z = minim:maxim, pmat = per), lty = 2) 270 } 271 if (TRUE) { 272 zlim = range(mat, finite = TRUE, na.rm = TRUE) 273 leglevel = pretty(zlim, 6) 274 legcol = col(length(leglevel)) 275 legpretty = as.character(abs(leglevel)) 276 temp = character(length(leglevel)) 277 temp[leglevel > 0] = "+" 278 temp[leglevel < 0] = "-" 279 temp[leglevel == 0] = " " 280 legpretty = paste(temp, legpretty, sep = "") 281 if (.theta <= 180) 282 legend("topright", inset = 0.02, legend = paste(">", legpretty), col = legcol, bg = "white", pt.cex = 1.5, cex = 0.75, pch = 15) 283 if (.theta > 180) 284 legend("topleft", inset = 0.02, legend = paste(">", legpretty), col = legcol, bg = "white", pt.cex = 1.5, cex = 0.75, pch = 15) 285 } 286 invisible(mat) 287} 288