Change spacing and axis labels on 3D partial plot of gbm model results - visreg package

Question

I am using the visreg package and the visreg2d() function to create a 3D partial dependency plot showing the interaction between two gbm model variables. I need to change the x,y,and z axis labels as well as the spacing to make nice plots for the journal article. I tried changing settings in par(), but that only allows me to change the label size.

I am generating the plot with the following:

visreg2d(final,n.trees=1000,xvar="ProbMn50ppb",yvar="DTW60Jurgens",plot.type="persp",type="conditional",
       theta=140,phi=40)

"final" is my model object.

I can remove the axes by setting axes=FALSE within the visreg2d() command. So I am guessing the solution could be to create custom labels with custom spacing. I can also change the labels with xlab=, ylab= zlab=. However, another issue is I cannot use an expression for the z label, which I need.

Here is what my plot currently looks like with the labels and axes:

Answer 1

visreg2d(plot.type = "persp") uses persp() to make a graph. If I remember correctly, persp() can't use expression() as label and don't have a option to move a label (you can move it roughly using \n(break) such as xlab="\nx-label"). So you need to do it manually using text(). Two value, text angle and center coordinates, are needed to do it. In the past I made a function to calculates these values (I referred to R mailing help), I show it.

# top/bot are coordinates(x,y,z) of the axis; pmat is output of persp()
# pos means the label is c(left(-1) or right(1), down(-1) or up(1)) of the axis

persp_lab_f <- function(top, bot, pmat, space, pos = c(-1, -1))
{
  coords_3d <- list(top = top, bot = bot)
  coords_2d <- lapply(coords_3d, function(v, pmat) unlist(trans3d(v[1], v[2], v[3], pmat)), pmat = pmat)
  coords_2d$mid <- (coords_2d$top + coords_2d$bot)/2  # mid is calculated from 2d-coordinates
  # coords_2d$mid <- unlist(trans3d(((top + bot)/2)[1], ((top + bot)/2)[2], ((top + bot)/2)[3], pmat)) if use mid in 3d-coordinates
  tb_diff <- coords_2d$top - coords_2d$bot
  angle <- 180/pi * atan2(tb_diff[2], tb_diff[1])
  names(angle) <- "angle"
  center <-  coords_2d$mid + sqrt(space^2 / sum(tb_diff^2)) * rev(abs(tb_diff)) * pos
  out <- list(angle = angle, center = as.data.frame(t(center)))
  return(out)
}

You didn't show reproducible example, I made it (plase provide it next time).

fit <- lm(Ozone ~ Solar.R + Wind + Temp + I(Wind^2) + I(Temp^2) +
            I(Wind*Temp)+I(Wind*Temp^2) + I(Temp*Wind^2) + I(Temp^2*Wind^2),
          data=airquality)

vis_d <- visreg2d(fit, xvar = "Wind", yvar = "Temp", plot.type="persp", type="conditional", theta = 140, phi = 40)

x <- vis_d$x
y <- vis_d$y
z <- vis_d$z

x_top <- c(max(x), max(y), min(z))
x_bot <- c(min(x), max(y), min(z))
y_top <- c(max(x), max(y), min(z))
y_bot <- c(max(x), min(y), min(z))
z_top <- c(max(x), min(y), max(z))
z_bot <- c(max(x), min(y), min(z))

pmat <- persp(x, y, z, theta = 140, phi = 40)

xlab_param <- persp_lab_f(x_top, x_bot, pmat, 0.1, pos = c(1, -1))
ylab_param <- persp_lab_f(y_top, y_bot, pmat, 0.1)
zlab_param <- persp_lab_f(z_top, z_bot, pmat, 0.1)

par(mar=c(1,1,1,1))
visreg2d(fit, xvar = "Wind", yvar = "Temp", plot.type="persp", type="conditional", theta = 140, phi = 40, 
         xlab = "", ylab = "", zlab = "")

text(xlab_param$center, srt = xlab_param$angle + 180, "xxxxxxxxxx")
text(ylab_param$center, srt = ylab_param$angle, "yyyyyyyyyy")
text(zlab_param$center, srt = zlab_param$angle + 180, labels = bquote(Sigma ~ .("zzzzzzzzzz")))