How to combine a 3d persp plot with a contour plot in R

Question

I am analyzing difference scores with polynomial regression in R. Based on [Edwards and Parry's (1993)][1] recommendations I have been trying to combine a persp() plot with a contour() plot. I would also need to plot the first two principal axes on the contour plot. My attempts so far have only provided me with each individual plot, but I don't know how to combine them. An example for the end-result is :

Edwards & Parry (1993) example difference score visualisation

I manage to get the persp() plot just fine. I have also obtained the contour plot. I can't seem to find any way to combine the two. I have managed to make the plot in plotly using the add_surface() option in the pipeline. My problem with the output is that the surface is smooth, and the contourplot lacks the values in the plot. Basically: persp() and contour() are visualised in a way that is extremely similar to the look I'm aiming for, per the example in the source.

My current attempt (in minimalistic code) is as follows:

surface <- function(e, i){
  y <- .2*e + .14*i + .08*e^2 + + .1*e*i + .2*i^2
}

e <- i <- seq(-3, 3, length= 20)
y <- outer(e, i, surface)

persp(e, i, y,
      xlab = 'Explicit',
      ylab = 'Implicit',
      zlab = 'Depression',
      theta = 45)

contour(e,i,y)

So basically my question is: how can I make a plot like Edwards and Parry (1993) make, with a similar visual style, in R. It does not have to be base-R, I'm happy with any method. I've been stuck on this problem for a week now.

My attempt in plotly (to compare it to my desired end-result) is:

if(!"plotly" %in% installed.packages){install.packages('plotly')}
library(plotly)

plot_ly(z = ~y) %>% add_surface(x = ~e, y= ~i, z= ~y,
    contours = list(
      z = list(
        show=TRUE,
        usecolormap=FALSE,
        highlightcolor="#ff0000",
        project=list(z=TRUE)
      )
    )
) %>%
  layout(
    scene=list(
      xaxis = list(title = "Explicit"),
      yaxis = list(title = "Implicit"),
      zaxis = list(title = "Depression")
    )
  )

[1]: Edwards, J. R., & Parry, M. E. (1993). On the use of polynomial regression as an alternative to difference scores. Academy of Management Journal, 36(6), 1577–1613. https://doi.org/10.2307/256822

Answer 1

I have found an answer and I will share it here. It seems it cannot be done in base-R. But the RSM-package allows for the addition of contour lines to the base of the plot.

In this answer I will give a minimal example of:

• the persp() plot
• the contour lines in the base
• addition of x=y and x=-y axis
• calculation and addition of the first and second principal axis

The only thing I could not solve is that the lines now are drawn over the surface. I don't know how to solve it.

library(rsm)

x <- seq(-3,3,by=0.25) 
y <- seq(-3,3,by=0.25) 
d <- expand.grid(x=x,y=y)
z <- c(data=NA,1089)
b0 = .140; b1 = -.441; b2 = -.154; b3 = .161 ; b4 =-.106; b5 = .168

k=1
for (i in 1:25) {
  for (j in 1:25) {
    z[k]=b0+b1*x[i]+b2*y[j]+b3*x[i]*x[i]+b4*x[i]*y[j]+ b5*y[j]*y[j]
    k=k+1
  }  }

data.lm <- lm(z~poly(x,y,degree=2),data=d)

res1 <- persp(data.lm,x~y, 
              zlim=c(-2,max(z)),
              xlabs = c('X','Y'),
              zlab = 'Z',
              contour=list(z="bottom"),
              theta=55,
              phi=25)    

# draw x=y line (lightly dotted)
xy_pos <- matrix(c(-3,-3,3,3),ncol=2,byrow = T)
lines(trans3d(xy_pos[,2], xy_pos[,1], z=-2, pmat = res1$`y ~ x`$transf), 
      lty = 3,
      col = 'darkgrey')

# draw x=-y line (lightly dotted)
xy_neg <- matrix(c(-3,3,3,-3),ncol=2,byrow = T)
lines(trans3d(xy_neg[,2], xy_neg[,1], z=-2, pmat = res1$`y ~ x`$transf), 
      lty = 3,
      col = 'darkgrey')

# Find stationary points:
X0 <- (b2*b4 - 2*b1*b5) / (4*b3*b5 - b4^2)
Y0 <- (b1*b4 - 2*b2*b3) / (4*b3*b5 - b4^2)

# First Principal Axis
p11 = (b5-b3+sqrt((b3-b5)^2+b4^2))/b4
p10 = Y0 - p11*X0
Ypaf1 = p10 + p11*x

# plot first principal axis (full line)
xypaf1 <- matrix(c(Ypaf1[1], -3, Ypaf1[25], 3),ncol=2, byrow=T)
lines(trans3d(xypaf1[,2], xypaf1[,1], z=-2, pmat = res1$`y ~ x`$transf), 
      lty = 1,
      col = 'black')

# Second Principal Axis
p21 = (b5-b3-sqrt((b3-b5)^2+b4^2))/b4
p20 = Y0 - p21*X0
Ypaf2 = p20 + p21*x

# plot second principal axis (dashed line)
xypaf2 <- matrix(c(Ypaf2[1], -3, Ypaf2[25], 3),ncol=2, byrow=T)
lines(trans3d(xypaf2[,2], xypaf2[,1], z=-2, pmat = res1$`y ~ x`$transf), 
      lty = 2,
      col = 'black')