r - 用“at”指定时如何绘制预测边距？

Question

我们可以得到线性模型的边际效应，margins::margins()并且可以使用 option 选择感兴趣的变量variables。

fit <- lm(mpg ~ factor(vs) + gear:factor(vs) + qsec, mtcars)

library(margins)
marg1 <- margins(fit, variables="vs")

> summary(marg1)
 factor    AME     SE      z      p   lower   upper
    vs1 4.8023 2.6769 1.7940 0.0728 -0.4443 10.0490

这个包有一个实现的方法plot.margins，所以我们可以绘制边际效应

plot(marg1)

at允许我们指定计算边际效应的值：

marg2 <- margins(fit, variables="vs", at=list(gear=c(3, 4, 5)))

> summary(marg2)
 factor   gear    AME     SE      z      p   lower   upper
    vs1 3.0000 2.8606 3.3642 0.8503 0.3952 -3.7332  9.4544
    vs1 4.0000 5.6849 2.6713 2.1282 0.0333  0.4493 10.9206
    vs1 5.0000 8.5093 3.8523 2.2089 0.0272  0.9588 16.0597

但是，尝试绘制这些指定的边距会产生错误：

plot(marg2)

Error in `[.data.frame`(summ, , names(attributes(x)[["at"]]), drop = FALSE) : 
  undefined columns selected

由于该margins软件包声称是“Stata 'margins' 命令的 R 端口”，我希望有一个类似于 Stata 给出的情节：

那么，我们如何绘制用 指定的预测边距at？

编辑：

请注意，这并不是一个普通的交互图，因为

with(mtcars[mtcars$gear %in% c(3, 4, 5), ], 
     interaction.plot(gear, vs, mpg, pch=rep(1, 2), type="b"))

给出不同的输出：

标签： rplotregressionmarginal-effects

Answer 1

该错误来自似乎是plot类对象方法中的错误"margins"，plot.margins.
这是试图纠正它。更改在函数体中，只需执行此操作或将其保存在文件中"plotmargins.R"，然后source("plotmargins.R").

plot.margins <-
function (x, pos = seq_along(marginal_effects(x, with_at = FALSE)), 
    which = colnames(marginal_effects(x, with_at = FALSE)), labels = gsub("^dydx_", 
        "", which), horizontal = FALSE, xlab = "", ylab = "Average Marginal Effect", 
    level = 0.95, pch = 21, points.col = "black", points.bg = "black", 
    las = 1, cex = 1, lwd = 2, zeroline = TRUE, zero.col = "gray", 
    ...) 
{
    pars <- list(...)
    summ <- summary(x, level = level, by_factor = TRUE)
    MEs <- summ[, "AME", drop = TRUE]
    lower <- summ[, ncol(summ) - 1L]
    upper <- summ[, ncol(summ)]
    r <- max(upper) - min(lower)

    #--- changes start here
    nms <- intersect(names(summ), names(attributes(x)[["at"]]))
    at_levels <- unique(summ[, nms, drop = FALSE])
    #--- changes end here

    n_at_levels <- nrow(at_levels)
    if (n_at_levels > 1) {
        pos2 <- rep(pos, each = n_at_levels)
        pos2 <- pos2 + seq(from = -0.2, to = 0.2, length.out = n_at_levels)
    }
    else {
        pos2 <- pos
    }
    if (isTRUE(horizontal)) {
        xlim <- if ("xlim" %in% names(pars)) 
            xlim
        else c(min(lower) - 0.04 * r, max(upper) + 0.04 * r)
        ylim <- if ("ylim" %in% names(pars)) 
            xlim
        else c(min(pos2) - (0.04 * min(pos2)), max(pos2) + (0.04 * 
            max(pos2)))
    }
    else {
        xlim <- if ("xlim" %in% names(pars)) 
            xlim
        else c(min(pos2) - (0.04 * min(pos2)), max(pos2) + (0.04 * 
            max(pos2)))
        ylim <- if ("ylim" %in% names(pars)) 
            xlim
        else c(min(lower) - 0.04 * r, max(upper) + 0.04 * r)
    }
    if (isTRUE(horizontal)) {
        plot(NA, xlim = xlim, ylim = ylim, yaxt = "n", xlab = ylab, 
            ylab = xlab, las = las, ...)
        if (isTRUE(zeroline)) {
            abline(v = 0, col = zero.col)
        }
        points(MEs, pos2, col = points.col, bg = points.bg, pch = pch)
        axis(2, at = pos, labels = as.character(labels), las = las)
        mapply(function(pos, upper, lower, lwd) {
            segments(upper, pos, lower, pos, col = points.col, 
                lwd = lwd)
        }, pos2, upper, lower, seq(max(lwd), 0.25, length.out = length(MEs)))
    }
    else {
        plot(NA, xlim = xlim, ylim = ylim, xaxt = "n", xlab = xlab, 
            ylab = ylab, las = las, ...)
        if (isTRUE(zeroline)) {
            abline(h = 0, col = zero.col)
        }
        points(pos2, MEs, col = points.col, bg = points.bg, pch = pch)
        axis(1, at = pos, labels = as.character(labels), las = las)
        mapply(function(pos, upper, lower, lwd) {
            segments(pos, upper, pos, lower, col = points.col, 
                lwd = lwd)
        }, pos2, upper, lower, seq(max(lwd), 0.25, length.out = length(MEs)))
    }
    invisible(x)
}

现在你的代码和图表。

source("plotmargins.R")

marg2 <- margins(fit, variables = "vs", 
                 at = list(gear = c(3, 4, 5)))

plot(marg2)