Given the following example:
library(mgcv)
set.seed(2)
dat <- gamSim(2, n = 1000, dist = "normal", scale = 1)
m <- gam(y ~ s(x, z, k = 40), data = dat$data, method = "REML")
fit <- predict(m)
vis.gam(m, se=T, type="response")
Is there any way to find peaks and valleys of the smoothed surface and plot it with label(s)? Something similar to the following:
Let's start by modifying the vis.gam so that it also prints the ouput (note that the only thing I changed was "return(list(m1=m1,m2=m2,z=z))" at the end):
my.gam.out = 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(gettextf("view variables must be one of %s",
paste(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(gettextf("View variables must contain more than one value. view = c(%s,%s).",
view[1], view[2]))
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 <- 4
}
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
}
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("zlab" %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 {
max.z <- max(fv$fit + fv$se.fit * se, na.rm = TRUE)
min.z <- min(fv$fit - fv$se.fit * se, na.rm = TRUE)
zlim <- c(min.z, max.z)
}
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))
}
return(list(m1=m1,m2=m2,z=z))
}
Then, you can define a function that will give you the position of the max in Z. I've made it so that you click on the image when prompted: you select the bottom left and then the top right to let R search for the maximum in Z and return you the coordinate in X and in Y.
isolate.surf = function(dat, #list from the my.gam.out() custom function
mod = NULL,
mod.var = c("x1","x2"),
variables = c("m1","m2","z"),
locator = FALSE,
plot=TRUE) {
# Extract the columns from the list
get.df.from.list = data.frame(x = dat[[variables[1]]], y = dat[[variables[2]]])
get.z = dat[[variables[3]]]
# Instead of providing the values, you can locate them by clicking on the graph
if(locator){
cat("Please, provide 2 points starting from the bottom left corner and top right:")
var = 2
par(mfrow = c(1,1))
vis.gam(mod, se=F,plot.type = "contour", type="response")
my.pts = locator(n = var)
xmin = my.pts$x[1]
xmax = my.pts$x[2]
ymin = my.pts$y[1]
ymax = my.pts$y[2]
}
# Isolate the location with the user specified x- and y-axis limits
find.x = which(get.df.from.list$x>xmin & get.df.from.list$x<xmax)
find.y = which(get.df.from.list$y>ymin & get.df.from.list$y<ymax)
# Subset the vector for each
x = get.df.from.list$x[find.x]
y = get.df.from.list$y[find.y]
z = get.z[find.x,find.y]
# Find the maximum value in the subset
max.z = max(z, na.rm = TRUE)
# Find the position of the maximum Z value
loc.z = which(z == max.z, arr.ind = TRUE)
# Find the X and Y coordinate in the subsetted dataset
x.y.z = c(x[loc.z[1]], y[loc.z[2]])
# More robust: find the location of the max z from the ORIGINAL indicies
get.df.from.list$x[find.x[loc.z[1]]]
get.df.from.list$y[find.y[loc.z[2]]]
get.z[find.x[loc.z[1]],find.y[loc.z[2]]]
# Get the coordinate from the original data
find.local.max = data.frame(
x = get.df.from.list$x[find.x[loc.z[1]]],
y = get.df.from.list$y[find.y[loc.z[2]]],
z = get.z[find.x[loc.z[1]],find.y[loc.z[2]]])
if(!is.null(mod)){
tmp.df = data.frame(x = x.y.z[1], y = x.y.z[2])
colnames(tmp.df) <- mod.var
mod.pred.z = predict(mod, newdata = tmp.df, type = "response")
} else {mod.pred.z = NULL}
# Make a visual verification
if(plot){
vis.gam(mod, se=F,plot.type = "contour", type="response")
points(tmp.df[1],tmp.df[2], cex = 1,pch = 19)
}
# Get the data out of the function
return(list(tmp.df = tmp.df,
mod.pred.z = mod.pred.z) # X-Y-Z coordinate in a dataframe
)
}
Then just use it with you data:
################
library(mgcv)
set.seed(2)
dat <- gamSim(2, n = 1000, dist = "normal", scale = 1)
dat$data$y
m <- gam(y ~ s(x, z, k = 40), data = dat$data, method = "REML")
fit <- predict(m)
vis.gam(m, se=F,plot.type = "contour", type="response")
mygame = my.gam.out(m, se=F,plot.type = "contour", type="response")
# you'll need to click on the graph to 'isolate' a peak
isolate.surf(dat = mygame, mod = m,mod.var = c("x","z"),variables = c("m1","m2","z"),locator = TRUE)
The answer:
Unfortunately, it won't do it for multiple peak. But a small change could be made to let the locator function gather more points.
I don't have the 'valleys' in this function, but it might get you started!