persp 在 R 中添加因子组

Question

按照边距小插图 https://cran.r-project.org/web/packages/margins/vignettes/Introduction.html#Motivation 我想知道如何在包含三重交互的 logit 之后使用 persp 进行绘图。

仅使用 persp 和 effect 仅显示部分交互（drat 和 wt）

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

不过，我希望在 am=0 和 am=1 处看到与上面相同的图表。我试过了：

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

但生成了相同的图形。如何查看am=0 和am=1 的两张图？或者在同一立方体中至少有两条曲线表示 am=0 和 am=1。

谢谢

Answer 1

使用 margins 包中的 persp.glm() 函数似乎无法做到这一点。您可能必须“手动”完成。

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="")

^{由 reprex package (v2.0.1)}

于 2022-05-16 创建

persp 在 R 中添加因子组

persp add factor group in R

3d

r

margins

layoutmargins