Plot three way interaction with 3d graph

Question

I have fitted a binominal logistic glm with a three-way interaction between sex (male & female), tree cover including a quadratic term (1-100%), and the mean tree cover of an area (1-100%). (case is 1 used and 0 for not used)

glm.winter.3 <- glm(case ~ 
                      sex. * mean95 * poly(tree.cover,2),
                    data = rsf.winter.3, family = binomial (link = "logit"))

I found a nice plot in a paper. I would like to do something similar but I can not find a way to approach it :

My data set is large. So it's hard to share it. Maybe somebody has an idea how to approach it anyway? Thanks

Answer 1

Here's an example with some built-in data. Using the mtcars data, I made a binary variable identifying gas guzzlers gg which is 1 for those cars that get less than 19 MPG and 0 for all others. I then modelled it as a multiplicative function of wt (weight), hp (horsepower) and vs whether or not it is a v-shaped engine. Here are the results:

library(lattice)
library(dplyr)
data(mtcars)
mtcars <- mtcars %>% 
  mutate(fast = as.numeric(qsec < mean(qsec)), 
         gg = as.numeric(mpg< 19), 
         carb2 = ifelse(carb <= 2, 0, 1))

mod <- glm(gg ~ hp* wt * vs, data=mtcars, family=binomial)
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(mod)
#> 
#> Call:
#> glm(formula = gg ~ hp * wt * vs, family = binomial, data = mtcars)
#> 
#> Deviance Residuals: 
#>      Min        1Q    Median        3Q       Max  
#> -2.24867  -0.00002   0.00000   0.28763   1.17741  
#> 
#> Coefficients:
#>               Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -3.299e+01  4.336e+01  -0.761    0.447
#> hp           1.346e-01  2.341e-01   0.575    0.565
#> wt           8.744e+00  1.271e+01   0.688    0.491
#> vs          -1.447e+03  5.600e+05  -0.003    0.998
#> hp:wt       -3.231e-02  6.798e-02  -0.475    0.635
#> hp:vs        9.144e+00  6.074e+03   0.002    0.999
#> wt:vs        4.512e+02  1.701e+05   0.003    0.998
#> hp:wt:vs    -2.906e+00  1.815e+03  -0.002    0.999
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 44.236  on 31  degrees of freedom
#> Residual deviance: 11.506  on 24  degrees of freedom
#> AIC: 27.506
#> 
#> Number of Fisher Scoring iterations: 22

The next step is to make sequences of the two continuous variables that go from their minima to maxima across some number of values - usually 25 is enough to give a smooth looking graph.

hp_seq <- seq(min(mtcars$hp, na.rm=TRUE), 
              max(mtcars$hp, na.rm=TRUE), 
              length=25)

wt_seq <- seq(min(mtcars$wt, na.rm=TRUE), 
              max(mtcars$wt, na.rm=TRUE), 
              length=25)

Next, we make all combinations of the two sequences of values and the binary variable (in this case vs).

eg <- expand.grid(hp = hp_seq, 
                  wt = wt_seq, 
                  vs = c(0,1))

We can then make predictions from those data and save them back into the data object. We also make vs a factor - the levels of the factor will show up in the strip above each plot.

eg$fit <- predict(mod, newdata=eg, type="response")
eg$vs <- factor(eg$vs, labels=c("VS = 0", "VS = 1"))

We use wireframe from the lattice package to make the plots

wireframe(fit ~ hp + wt | vs, data=eg, drape=TRUE, 
          default.scales = list(arrows=FALSE))

^{Created on 2022-05-06 by the reprex package (v2.0.1)}

Note, there is another answer about making 3D surface plots here that demonstrates perp from base R and plot_ly from the plotly package and persp3D from the plot3D package.