persp add factor group in R

Following the margins vignette https://cran.r-project.org/web/packages/margins/vignettes/Introduction.html#Motivation I would like to know how to plot using persp after a logit containing a triple interaction.

Using only persp and effect only part of the interaction is shown (drat and wt)

x1 <- lm(mpg ~ drat * wt * am, data = mtcars)
head(mtcars)
persp(x1, what = "effect")

However I would like to see the same graph above but at am=0 and am=1. I tried:

persp(x1,"drat","wt", at = list(am = 0:1), what = "effect")

But the same graph is produced. How to see two graphs at am=0 and am=1? or at least two curves representing am=0 and am=1 in the same cube.

Thanks

It doesn't look like you can do it with the persp.glm() function in the margins package. You will probably have to do it "by hand".

data(mtcars)
mtcars$hihp <- as.numeric(mtcars$hp > quantile(mtcars$hp,.5))
x1 <- glm(hihp ~ drat * wt * am + disp + qsec, data = mtcars, family=binomial)
#> Warning: glm.fit: algorithm did not converge
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

drat_s <- with(mtcars, seq(min(drat), max(drat),length=25))
wt_s <- with(mtcars, seq(min(wt), max(wt), length=25))

pred_fun <- function(x,y, am=0){
  tmp <- data.frame(drat = x, wt = y, am=am, 
                    disp = mean(mtcars$disp, na.rm=TRUE), 
                    qsec = mean(mtcars$qsec, na.rm=TRUE))
  predict(x1, newdata=tmp, type="response")
}


p0 <- outer(drat_s, wt_s, pred_fun)
p1 <- outer(drat_s, wt_s, pred_fun, am=1)


persp(drat_s, wt_s, p0, zlim=c(0,1), theta=-80, col=rgb(.75,.75, .75, .75), 
      xlab = "Axle Ratio", 
      ylab="Weight", 
      zlab="Predicted Probability")
par(new=TRUE)
persp(drat_s, wt_s, p1, zlim=c(0,1), theta=-80, col=rgb(1,0,0,.75), xlab="", ylab="", zlab="")

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

Edit: what if you add a factor to the model?

If we turn cyl into a factor and add it to the model, we also have to add it to the tmp object in the predfun() function, however it has to have the same properties that it has in the data, i.e., it has to be a factor (that has a single value) that has the same levels and labels as the one in the data. Here's an example:

data(mtcars)
mtcars$hihp <- as.numeric(mtcars$hp > quantile(mtcars$hp,.5))
mtcars$cyl <- factor(mtcars$cyl)
x1 <- glm(hihp ~ drat * wt * am + disp + qsec + cyl, data = mtcars, family=binomial)
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

drat_s <- with(mtcars, seq(min(drat), max(drat),length=25))
wt_s <- with(mtcars, seq(min(wt), max(wt), length=25))

pred_fun <- function(x,y, am=0){
  tmp <- data.frame(drat = x, wt = y, am=am, 
                    disp = mean(mtcars$disp, na.rm=TRUE), 
                    qsec = mean(mtcars$qsec, na.rm=TRUE), 
                    cyl = factor(2, levels=1:3, labels=levels(mtcars$cyl)))
  predict(x1, newdata=tmp, type="response")
}


p0 <- outer(drat_s, wt_s, pred_fun)
p1 <- outer(drat_s, wt_s, pred_fun, am=1)


persp(drat_s, wt_s, p0, zlim=c(0,1), theta=-80, col=rgb(.75,.75, .75, .75), 
      xlab = "Axle Ratio", 
      ylab="Weight", 
      zlab="Predicted Probability")
par(new=TRUE)
persp(drat_s, wt_s, p1, zlim=c(0,1), theta=-80, col=rgb(1,0,0,.75), xlab="", ylab="", zlab="")

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