Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

how to fill contour colors and write axes names in RSM (R)

Tags:

r

graphics

I have following data

ct<-structure(list(Conc = c(50L, 100L, 150L, 50L, 100L, 150L, 50L, 
100L, 150L, 100L, 100L, 100L), kGy = c(10L, 10L, 10L, 15L, 15L, 
15L, 20L, 20L, 20L, 15L, 15L, 15L), CT.Y. = c(75L, 65L, 51L, 
87L, 93L, 89L, 81L, 86L, 78L, 92L, 93L, 92L)), .Names = c("Conc", 
"kGy", "CT.Y."), class = "data.frame", row.names = c(NA, -12L))

And i am using following R code for response surface

library(rsm)
ct.rsm<-rsm(CT.Y.~SO(Conc, kGy), data=ct)
persp(ct.rsm, Conc ~ kGy, col=rainbow(50), theta=60,
    phi=0, r = 3, d=1, border = NULL, ltheta = -135, lphi = 0
    , shade = 0.75, zlab="CT",ylab="Concentration %", col.axis=37, font.lab=2,col.lab=33,
    contour=("colors"))

One question is that how can i fill the colors in contours? and other question is about the axes labeling. For label of X and Z axes i can label it, but when i want to include the label of Y-axis, i receive following error.

Error in persp.default(dat$x, dat$y, dat$z, xlab = dat$labs[1], ylab = dat$labs[2],  : 
  formal argument "ylab" matched by multiple actual arguments

Hope that some one can help me in this regard. Thanks in advance. Response Surface Plot

like image 970
Iftikhar Avatar asked Oct 28 '11 19:10

Iftikhar


1 Answers

I pieced together a working example of your data (without your rownames). The object returned from the rsm function is classes "rsm" "lm", so it gets handled by persp.lm. That function has a hard-coded ylab specification and no provision for relabeling. It can be fixed (with a puzzling reversal of x and ylabs). I changed the line function in draw.cont.line to polygon and it illustrates the need for further efforts to link up the endpoints as mentioned in my comment below:

    persp.lm <- 
function (x, form, at, bounds, zlim, zlab, xlabs, col = "white", xlab=xlab,
    contours = NULL, hook, atpos = 3, theta = -25, phi = 20, 
    r = 4, border = NULL, box = TRUE, ticktype = "detailed", ylab,
    ... ) 
{
    draw.cont.line = function(line) {
        if (cont.varycol) {
            cont.col = col
            if (length(col) > 1) 
                cont.col = col[cut(c(line$level, dat$zlim), length(col))][1]
        }
        polygon(trans3d(line$x, line$y, cont.z, transf), col = cont.col, 
            lwd = cont.lwd)
    }
    plot.data = contour.lm(x, form, at, bounds, zlim, xlabs, 
        atpos = atpos, plot.it = FALSE)
    transf = list()
    if (missing(zlab)) 
        zlab = ""
    facet.col = col
    cont = !is.null(contours)
    if (mode(contours) == "logical") 
        cont = contours
    cont.first = cont
    cont.z = cz = plot.data[[1]]$zlim[1]
    cont.col = 1
    cont.varycol = FALSE
    cont.lwd = 1
    if (is.character(contours)) {
        idx = charmatch(contours, c("top", "bottom", "colors"), 
            0)
        if (idx == 1) {
            cont.first = FALSE
            cont.z = plot.data[[1]]$zlim[2]
        }
        else if (idx == 2) {
        }
        else if (idx == 3) {
            cont.varycol = TRUE
            if (length(col) < 2) 
                col = rainbow(40)
        }
        else cont.col = contours
    }
    else if (is.list(contours)) {
        if (!is.null(contours$z)) 
            cz = contours$z
        if (is.numeric(cz)) 
            cont.z = cz
        else if (cz == "top") {
            cont.first = FALSE
            cont.z = plot.data[[1]]$zlim[2]
        }
        if (!is.null(contours$col)) 
            cont.col = contours$col
        if (!is.null(contours$lwd)) 
            cont.lwd = contours$lwd
        if (charmatch(cont.col, "colors", 0) == 1) {
            cont.varycol = TRUE
            if (length(col) < 2) 
                col = rainbow(40)
        }
    }
    for (i in 1:length(plot.data)) {
        dat = plot.data[[i]]
        cont.lines = NULL
        if (!missing(hook)) 
            if (!is.null(hook$pre.plot)) 
                hook$pre.plot(dat$labs)
        if (cont) 
            cont.lines = contourLines(dat$x, dat$y, dat$z)
        if (cont && cont.first) {
            transf = persp(dat$x, dat$y, dat$z, zlim = dat$zlim, xlab=ylab,
                theta = theta, phi = phi, r = r, col = NA, border = NA, 
                box = FALSE)
            lapply(cont.lines, draw.cont.line)
            par(new = TRUE)
        }
        if (length(col) > 1) {
            nrz = nrow(dat$z)
            ncz = ncol(dat$z)
            zfacet = dat$z[-1, -1] + dat$z[-1, -ncz] + dat$z[-nrz, 
                -1] + dat$z[-nrz, -ncz]
            zfacet = c(zfacet/4, dat$zlim)
            facet.col = cut(zfacet, length(col))
            facet.col = col[facet.col]
        }
        transf = persp(dat$x, dat$y, dat$z, xlab = xlab, 
             zlab = zlab, zlim = dat$zlim, ylab=ylab,
            col = facet.col, border = border, box = box, theta = theta, 
            phi = phi, r = r, ticktype = ticktype)
        if (atpos == 3) 
            title(sub = dat$labs[5])
        if (cont && !cont.first) 
            lapply(cont.lines, draw.cont.line)
        if (!missing(hook)) 
            if (!is.null(hook$post.plot)) 
                hook$post.plot(dat$labs)
        plot.data[[i]]$transf = transf
    }
    invisible(plot.data)
}

persp(ct.rsm, Conc ~ kGy, col=rainbow(50), theta=60, xlab="Something",
    phi=0, r = 3, d=1, border = NULL, ltheta = -135, lphi = 0
    , shade = 0.75, zlab="CT",ylab="Concentration %", col.axis=37, font.lab=2,col.lab=33,
    contour=("colors"))

enter image description here

like image 105
IRTFM Avatar answered Oct 21 '22 08:10

IRTFM