persp add factor group in R-CodePudding

Following the margins vignette

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

CodePudding user response：

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)}