rcolorsgam
Changing the colors in a contour plot from vis.gam in mgcv
I am using vis.gam() to plot the results of a model. The code is the following :
vis.gam(model,view=c("CHLA","SST"),xlim=c(0,15),ylim=c(10,30),xlab="CHL-a (mg.m-3)",ylab="SST (°C)",plot.type="contour",type="response",color="terrain",n.grid=200)
I'd like to set new colors from my own palette as mentioned in the following code :
# colors palette definition
jet.colors <-colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan",
"#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))
vis.gam(model,view=c("CHLA","SST"),xlim=c(0,15),ylim=c(10,30),xlab="CHL-a (mg.m-3)",ylab="SST (°C)",plot.type="contour",type="response",color=jet.colors(10),n.grid=200)
Unfortunately, I have this error :
Error in vis.gam(model, view = c("CHLA", "SST"), xlim = c(0, 15), ylim = c(10, : color scheme not recognised
I don't know how to use the jet.colors palette for the function vis.gam...
Thanks for your help !
Solution
You can customize the vis.gam function to use your jet.colors ramp:
Here's a hacked version (see my inline comments):
jet.colors <-colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan",
"#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))
myvis.gam <- function (x, view = NULL, cond = list(), n.grid = 30, too.far = 0,
col = NA, color = "heat", contour.col = NULL, se = -1, type = "link",
plot.type = "persp", zlim = NULL, nCol = 50, ...)
{
fac.seq <- function(fac, n.grid) {
fn <- length(levels(fac))
gn <- n.grid
if (fn > gn)
mf <- factor(levels(fac))[1:gn]
else {
ln <- floor(gn/fn)
mf <- rep(levels(fac)[fn], gn)
mf[1:(ln * fn)] <- rep(levels(fac), rep(ln, fn))
mf <- factor(mf, levels = levels(fac))
}
mf
}
dnm <- names(list(...))
v.names <- names(x$var.summary)
if (is.null(view)) {
k <- 0
view <- rep("", 2)
for (i in 1:length(v.names)) {
ok <- TRUE
if (is.matrix(x$var.summary[[i]]))
ok <- FALSE
else if (is.factor(x$var.summary[[i]])) {
if (length(levels(x$var.summary[[i]])) <= 1)
ok <- FALSE
}
else {
if (length(unique(x$var.summary[[i]])) == 1)
ok <- FALSE
}
if (ok) {
k <- k + 1
view[k] <- v.names[i]
}
if (k == 2)
break
}
if (k < 2)
stop("Model does not seem to have enough terms to do anything useful")
}
else {
if (sum(view %in% v.names) != 2)
stop(paste(c("view variables must be one of", v.names),
collapse = ", "))
for (i in 1:2) if (!inherits(x$var.summary[[view[i]]],
c("numeric", "factor")))
stop("Don't know what to do with parametric terms that are not simple numeric or factor variables")
}
ok <- TRUE
for (i in 1:2) if (is.factor(x$var.summary[[view[i]]])) {
if (length(levels(x$var.summary[[view[i]]])) <= 1)
ok <- FALSE
}
else {
if (length(unique(x$var.summary[[view[i]]])) <= 1)
ok <- FALSE
}
if (!ok)
stop(paste("View variables must contain more than one value. view = c(",
view[1], ",", view[2], ").", sep = ""))
if (is.factor(x$var.summary[[view[1]]]))
m1 <- fac.seq(x$var.summary[[view[1]]], n.grid)
else {
r1 <- range(x$var.summary[[view[1]]])
m1 <- seq(r1[1], r1[2], length = n.grid)
}
if (is.factor(x$var.summary[[view[2]]]))
m2 <- fac.seq(x$var.summary[[view[2]]], n.grid)
else {
r2 <- range(x$var.summary[[view[2]]])
m2 <- seq(r2[1], r2[2], length = n.grid)
}
v1 <- rep(m1, n.grid)
v2 <- rep(m2, rep(n.grid, n.grid))
newd <- data.frame(matrix(0, n.grid * n.grid, 0))
for (i in 1:length(x$var.summary)) {
ma <- cond[[v.names[i]]]
if (is.null(ma)) {
ma <- x$var.summary[[i]]
if (is.numeric(ma))
ma <- ma[2]
}
if (is.matrix(x$var.summary[[i]]))
newd[[i]] <- matrix(ma, n.grid * n.grid, ncol(x$var.summary[[i]]),
byrow = TRUE)
else newd[[i]] <- rep(ma, n.grid * n.grid)
}
names(newd) <- v.names
newd[[view[1]]] <- v1
newd[[view[2]]] <- v2
if (type == "link")
zlab <- paste("linear predictor")
else if (type == "response")
zlab <- type
else stop("type must be \"link\" or \"response\"")
fv <- predict.gam(x, newdata = newd, se.fit = TRUE, type = type)
z <- fv$fit
if (too.far > 0) {
ex.tf <- exclude.too.far(v1, v2, x$model[, view[1]],
x$model[, view[2]], dist = too.far)
fv$se.fit[ex.tf] <- fv$fit[ex.tf] <- NA
}
if (is.factor(m1)) {
m1 <- as.numeric(m1)
m1 <- seq(min(m1) - 0.5, max(m1) + 0.5, length = n.grid)
}
if (is.factor(m2)) {
m2 <- as.numeric(m2)
m2 <- seq(min(m1) - 0.5, max(m2) + 0.5, length = n.grid)
}
if (se <= 0) {
old.warn <- options(warn = -1)
av <- matrix(c(0.5, 0.5, rep(0, n.grid - 1)), n.grid,
n.grid - 1)
options(old.warn)
max.z <- max(z, na.rm = TRUE)
z[is.na(z)] <- max.z * 10000
z <- matrix(z, n.grid, n.grid)
surf.col <- t(av) %*% z %*% av
surf.col[surf.col > max.z * 2] <- NA
if (!is.null(zlim)) {
if (length(zlim) != 2 || zlim[1] >= zlim[2])
stop("Something wrong with zlim")
min.z <- zlim[1]
max.z <- zlim[2]
}
else {
min.z <- min(fv$fit, na.rm = TRUE)
max.z <- max(fv$fit, na.rm = TRUE)
}
surf.col <- surf.col - min.z
surf.col <- surf.col/(max.z - min.z)
surf.col <- round(surf.col * nCol)
con.col <- 1
if (color == "heat") {
pal <- heat.colors(nCol)
con.col <- 3
}
else if (color == "topo") {
pal <- topo.colors(nCol)
con.col <- 2
}
else if (color == "cm") {
pal <- cm.colors(nCol)
con.col <- 1
}
else if (color == "terrain") {
pal <- terrain.colors(nCol)
con.col <- 2
}
else if (color == "gray" || color == "bw") {
pal <- gray(seq(0.1, 0.9, length = nCol))
con.col <- 1
}
### customized here
else if (color == 'jet') {
pal <- jet.colors(nCol)
con.col = 1
}
####
else stop("color scheme not recognised")
if (is.null(contour.col))
contour.col <- con.col
surf.col[surf.col < 1] <- 1
surf.col[surf.col > nCol] <- nCol
if (is.na(col))
col <- pal[as.array(surf.col)]
z <- matrix(fv$fit, n.grid, n.grid)
if (plot.type == "contour") {
stub <- paste(ifelse("xlab" %in% dnm, "", ",xlab=view[1]"),
ifelse("ylab" %in% dnm, "", ",ylab=view[2]"),
ifelse("main" %in% dnm, "", ",main=zlab"), ",...)",
sep = "")
if (color != "bw") {
txt <- paste("image(m1,m2,z,col=pal,zlim=c(min.z,max.z)",
stub, sep = "")
eval(parse(text = txt))
txt <- paste("contour(m1,m2,z,col=contour.col,zlim=c(min.z,max.z)",
ifelse("add" %in% dnm, "", ",add=TRUE"), ",...)",
sep = "")
eval(parse(text = txt))
}
else {
txt <- paste("contour(m1,m2,z,col=1,zlim=c(min.z,max.z)",
stub, sep = "")
eval(parse(text = txt))
}
}
else {
stub <- paste(ifelse("xlab" %in% dnm, "", ",xlab=view[1]"),
ifelse("ylab" %in% dnm, "", ",ylab=view[2]"),
ifelse("main" %in% dnm, "", ",zlab=zlab"), ",...)",
sep = "")
if (color == "bw") {
op <- par(bg = "white")
txt <- paste("persp(m1,m2,z,col=\"white\",zlim=c(min.z,max.z) ",
stub, sep = "")
eval(parse(text = txt))
par(op)
}
else {
txt <- paste("persp(m1,m2,z,col=col,zlim=c(min.z,max.z)",
stub, sep = "")
eval(parse(text = txt))
}
}
}
else {
if (color == "bw" || color == "gray") {
subs <- paste("grey are +/-", se, "s.e.")
lo.col <- "gray"
hi.col <- "gray"
}
else {
subs <- paste("red/green are +/-", se, "s.e.")
lo.col <- "green"
hi.col <- "red"
}
if (!is.null(zlim)) {
if (length(zlim) != 2 || zlim[1] >= zlim[2])
stop("Something wrong with zlim")
min.z <- zlim[1]
max.z <- zlim[2]
}
else {
z.max <- max(fv$fit + fv$se.fit * se, na.rm = TRUE)
z.min <- min(fv$fit - fv$se.fit * se, na.rm = TRUE)
}
zlim <- c(z.min, z.max)
z <- fv$fit - fv$se.fit * se
z <- matrix(z, n.grid, n.grid)
if (plot.type == "contour")
warning("sorry no option for contouring with errors: try plot.gam")
stub <- paste(ifelse("xlab" %in% dnm, "", ",xlab=view[1]"),
ifelse("ylab" %in% dnm, "", ",ylab=view[2]"), ifelse("zlab" %in%
dnm, "", ",zlab=zlab"), ifelse("sub" %in% dnm,
"", ",sub=subs"), ",...)", sep = "")
txt <- paste("persp(m1,m2,z,col=col,zlim=zlim", ifelse("border" %in%
dnm, "", ",border=lo.col"), stub, sep = "")
eval(parse(text = txt))
par(new = TRUE)
z <- fv$fit
z <- matrix(z, n.grid, n.grid)
txt <- paste("persp(m1,m2,z,col=col,zlim=zlim", ifelse("border" %in%
dnm, "", ",border=\"black\""), stub, sep = "")
eval(parse(text = txt))
par(new = TRUE)
z <- fv$fit + se * fv$se.fit
z <- matrix(z, n.grid, n.grid)
txt <- paste("persp(m1,m2,z,col=col,zlim=zlim", ifelse("border" %in%
dnm, "", ",border=hi.col"), stub, sep = "")
eval(parse(text = txt))
}
}
For usage just use color = 'jet':
library(mgcv)
set.seed(0)
n<-200;sig2<-4
x0 <- runif(n, 0, 1);x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1)
y<-x0^2+x1*x2 +runif(n,-0.3,0.3)
g<-gam(y~s(x0,x1,x2))
myvis.gam(g,ticktype="detailed",color='heat',theta=-35)
myvis.gam(g,ticktype="detailed",color='jet',theta=-35)
- Calculate transition probabilities
- How to subset R dataframe based on specific values in several columns?
- Change colon to period in figure caption in Quarto
- Applying different functions to different groups using 2 data frames
- Disable button for 5 seconds with different label and change after it
- How to group columns with an odd index, based on columns with an even index?
- Different plots of marginal effects with interaction term and offset when using different packages
- Animate a MWE GIF in R
- populating a matrix by rows
- How to relevel in ROC curve in R/plotROC?
- R - glmer different results on different machines (non-deterministic)
- Dynamically operate on column (inside :=)
- Random change in obj_addr() output when including the objects into a list and vectorizing over them
- Add to curly-curly argument
- h2o predict error: Test/Validation dataset has no columns in common with the training set
- Is is possible to convert a dataframe object to a tribble constructor?
- encoding issues in R
- How to Scrape NBA stats page using rvest
- Conversion of Decimal Minutes Seconds into Decimal Degrees when cells within the Columns Contains the Factor 'Missed' in R
- Difficulty using case_when() to add column that, conditionally, pastes value from another column
- Flip the matrix
- Creating test for date columns incorrectly throwing error
- Exists() and getSymbols() can't find symbol listed in yahoo
- spearman's rho plot
- Use OpenMP on M2 Mac with R and data.table
- How to find the connected instances from a minimum spanning trees model in R
- Hide a column in shiny datatable but keep it searchable
- How to extract data from an edited datatable in shiny?
- Render Data Table with double header
- filter dataTables in Shiny Dashboad based on selectInput Values