How would you plot a regression in 3D that has an interaction and its CIs?

First of all this is my data:

df<-data.frame(fpergandei=c(0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1),
mintemp2=c(1.29569224,1.40509105,0.74869815,1.51448987,2.49907923,-1.00168292,
-0.01709357,-0.34529002,-0.56408765,-0.01709357,0.20170407,0.42050170,-1.00168292,
-1.98627227,0.09230525,-0.23589120,-0.78288528,-0.89228410,-0.01709357,1.95208514,
0.31110289,0.31110289,0.09230525,0.96749579,-0.23589120,-0.01709357,-1.22048055,
-0.23589120,-0.45468883,-0.67348647,-1.65807582,-0.45468883),
moist2=c(0.24947803,-1.17491998,0.05498936,-1.58242006,-1.46943140,2.16472842,1.64979649,
1.65535331,-1.00821540,0.40506897,-0.22840843,-0.28027207,0.92370544,1.04039865,
1.36639871,1.32564870,-0.18210160,-0.25248798,-0.81557900,-0.72481762,-1.75282919,
-0.58219259,-0.27286298,-0.07281749,0.07536436,-0.04688566,-0.41363574,-0.54699940,
-0.19691978,0.89036452,0.84035315,-1.03785177))

df
   fpergandei    mintemp2      moist2
1           0  1.29569224  0.24947803
2           0  1.40509105 -1.17491998
3           0  0.74869815  0.05498936
4           0  1.51448987 -1.58242006
5           0  2.49907923 -1.46943140
6           0 -1.00168292  2.16472842
7           0 -0.01709357  1.64979649
8           0 -0.34529002  1.65535331
9           1 -0.56408765 -1.00821540
10          1 -0.01709357  0.40506897
11          1  0.20170407 -0.22840843
12          1  0.42050170 -0.28027207
13          0 -1.00168292  0.92370544
14          0 -1.98627227  1.04039865
15          0  0.09230525  1.36639871
16          0 -0.23589120  1.32564870
17          1 -0.78288528 -0.18210160
18          1 -0.89228410 -0.25248798
19          1 -0.01709357 -0.81557900
20          1  1.95208514 -0.72481762
21          1  0.31110289 -1.75282919
22          1  0.31110289 -0.58219259
23          1  0.09230525 -0.27286298
24          1  0.96749579 -0.07281749
25          1 -0.23589120  0.07536436
26          1 -0.01709357 -0.04688566
27          1 -1.22048055 -0.41363574
28          1 -0.23589120 -0.54699940
29          1 -0.45468883 -0.19691978
30          1 -0.67348647  0.89036452
31          1 -1.65807582  0.84035315
32          1 -0.45468883 -1.03785177

note i factorized "fpergandei" to make it two levels of 1 and 0

So I recently ran a GLM binomial with two continuous explanatory variables. The results yielded in a significant interaction between the two explanatory variables. I ended up plotting the model using the persp() function

mylogit<- glm(fpergandei~mintemp2*moist2,data=fedelog,family="binomial")
press_grid <- seq(-2.2, 2.2, by = 0.1)
v_grid <- seq(-2.2, 2.2, by = 0.1)
newdat <- data.frame(mintemp2 = rep(press_grid, times = length(v_grid)), 
moist2 = rep(v_grid, each = length(press_grid)))
pred <- predict.glm(mylogit, newdata = newdat, type="response")
z <- matrix(pred, length(press_grid))
res<-persp(press_grid, v_grid, z, xlab = "Min. Temperature", ylab = 
"Moisture", zlab = "Predicted Probability", main = "Plot Name", theta = 60, 
phi = 27, col = mycol)

I've also figured out how to use the persp function to plot confidence intervals.

pred2 <- predict.glm(mylogit, newdata = newdat, type="response", se.fit = 
TRUE)#with confinterval
pred$se.fit #standard errors for all predicted values
CIlow <- exp(pred-1.96*pred2$se.fit)/(1+exp(pred-1.96*pred2$se.fit)) 
#calculating lower confidence interval 
CIup <- exp(pred+1.96*pred2$se.fit)/(1+exp(pred+1.96*pred2$se.fit))  
#calculating upper confidence intervals

My issue is that when I plot surfaces of the upper and lower confidence intervals they overlap over and under each other:

res<-persp(press_grid, v_grid, ci.low, xlab = "Min. Temperature", ylab = 
"Moisture", zlab = "Predicted Probability", main = "Lower Confidence 
Interval", theta = 60, phi = 27, col = "grey") #surface w/ CI low
par(new=TRUE)
res<-persp(press_grid, v_grid, z, xlab = "Min. Temperature", ylab = 
"Moisture", zlab = "Predicted Probability", main = "Logistic Curve)", theta 
= 60, phi = 27, col = mycol) #surface with pred
par(new=TRUE)
res<-persp(press_grid, v_grid, ci.up, xlab = "Min. Temperature", ylab = 
"Moisture", zlab = "Predicted Probability", main = "Upper Confidence 
Interval", theta = 60, phi = 27, col = "grey") #surface w/ CI low

Is there any way I can make it so that the surface plots end up not overlapping each other and just creating a mess?

Solution

Here is an attempt with plotly. Good thing about this approach is that you can rotate the view till you are satisfied with the perspective

library(plotly)

mylogit <- glm(fpergandei ~ mintemp2 * moist2, 
               data = df,
               family = "binomial")

the following two vectors can be called "mintemp2" and "moist2" I kept the following names because of the OP

press_grid <- seq(-2.2, 2.2, by = 0.1)
v_grid <- seq(-2.2, 2.2, by = 0.1)
newdat <- expand.grid(press_grid, v_grid)  #the grid results in the same values as the newdat in the OP
colnames(newdat) <- c("mintemp2", "moist2")

pred <- predict.glm(mylogit, newdata = newdat, type="link", se=TRUE)

ymin <- mylogit$family$linkinv(pred$fit - 1.96 * pred$se.fit)
ymax <- mylogit$family$linkinv(pred$fit + 1.96 * pred$se.fit)
fit <- mylogit$family$linkinv(pred$fit) 
z <- matrix(fit, length(press_grid))
ci.low <- matrix(ymin, length(press_grid))
ci.up <- matrix(ymax, length(press_grid))


plot_ly(x = press_grid, y = v_grid) %>% 
  add_surface(z = z,
              colorscale = list(c(0,1),c("red","blue"))) %>% 
  add_surface(z = ci.low, opacity = 0.5, showscale = FALSE, colorscale = list(c(0,1),c("grey","grey"))) %>% 
  add_surface(z = ci.up, opacity = 0.5, showscale = FALSE, colorscale = list(c(0,1),c("grey","grey")))

Here is an idea how to plot this in 2 dimensions which I think is more effective

library(tidyverse)

data.frame(pred = fit,
           low = ymin,
           high = ymax,
           newdat) %>%
  filter(moist2 %in%  unique(.$moist2)[c(T, rep(F, 3))]) %>%
  mutate(facet = factor(rep(c("-2.2 - -1.4", "-1 - -0.2", "0.2 - 1", "1.4 - 2.2"), each = length(.$moist2)/4),
                        levels = c("-2.2 - -1.4", "-1 - -0.2", "0.2 - 1", "1.4 - 2.2")),
         moist2 = as.factor(moist2)) %>%
  ggplot()+
  geom_line(aes(x = mintemp2, y = pred, color = moist2))+
  geom_ribbon(aes(x = mintemp2, ymin = low, ymax = high, fill = moist2), alpha = 0.1)+
  facet_wrap(~facet, ncol = 2)+
  theme_bw()+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

In general using fewer values of one variable, and spiting the range into facets. At first look one can see the 0.2 - 1 moist2 range changes the direction of the probability function which I think is quite interesting and emphasizes the interaction of these two variables. Something I failed to see with the 3D plot.