在 R、lincom + coefplot 或 plotbeta 中将回归系数（偏导数）与 CI 相结合的绘图？

Question

大多数时候我们运行一个带有交互项的回归，我们对偏导数感兴趣。例如，考虑下面的模型，

如果我想知道X1对P(Y)的影响，或者X1对P(Y)的偏导数，我需要下面的系数组合：

我可以使用 R 中的 lincom 函数来计算回归参数的线性组合，而不是手动计算。但我不仅想知道这样的计算得出的数字；我想绘制它们。问题是，如果我使用 R 包来绘制系数（例如，coefplot），它会绘制我模型中的系数，但没有系数线性组合的选项。有什么方法可以将 lincom 函数（或其他计算参数组合的函数）与 coefplot（或其他具有此选项的系数绘图包）结合起来？

当然，在上面的例子中，我只考虑了 X1 的导数，如果我绘制它，我将得到一个只有一个点及其置信区间的图，但我想在图中显示系数对于 X1、X2 和 Z 的偏导数，如下例所示。

系数图（我有的那个）：

参数组合或偏导数图（我想得到的那个）：

我发现 Stata 有一个函数可以满足我的需求，称为“plotbeta”。 R有类似的东西吗？

Answer 1

这是一个开始。这定义了一个名为 plotBeta() 的函数，... 是传递给估计文本的 geom_text() 的参数。

plotBeta <- function(mod, confidence_level = .95, include_est=TRUE, which.terms=NULL, plot=TRUE, ...){
  require(glue)
  require(ggplot2)
  b <- coef(mod)
  mains <- grep("^[^:]*$", names(b), value=TRUE)
  mains.ind <- grep("^[^:]*$", names(b))
  if(!is.null(which.terms)){
    if(!(all(which.terms %in% mains)))stop("Not all terms in which.terms are in the model\n")
    ins <- match(which.terms, mains)
    mains <- mains[ins]
    mains.ind <- mains.ind[ins]
  }
  icept <- grep("Intercept", mains)
  if(length(icept) > 0){
    mains <- mains[-icept]
    mains.ind <- mains.ind[-icept]
  }
  if(inherits(mod, "lm") & !inherits(mod, "glm")){
    crit <- qt(1-(1-confidence_level)/2, mod$df.residual)
  }else{
    crit <- qnorm(1-(1-confidence_level)/2)
  }
  out.df <- NULL
  for(i in 1:length(mains)){
    others <- grep(glue("^{mains[i]}:"), names(b))
    others <- c(others, grep(glue(":{mains[i]}:"), names(b)))
    others <- c(others, grep(glue(":{mains[i]}$"), names(b)))
    all.inds <- c(mains.ind[i], others)
    ones <- rep(1, length(all.inds))
    est <- c(b[all.inds] %*% ones)
    se.est <- sqrt(c(ones %*% vcov(mod)[all.inds, all.inds] %*% ones))
    lower <- est - crit*se.est
    upper <- est + crit*se.est
    tmp <- data.frame(var = mains[i],  
                      lab = glue("dy/d{mains[i]} = {paste('B', all.inds, sep='', collapse=' + ')}"), 
                      labfac = i, 
                      est = est, 
                      se.est = se.est, 
                      lower = lower, 
                      upper=upper)
    tmp$est_text <- sprintf("%.2f (%.2f, %.2f)", tmp$est, tmp$lower, tmp$upper)
    out.df <- rbind(out.df, tmp)
  }
  out.df$labfac <- factor(out.df$labfac, labels=out.df$lab)
  if(!plot){
    return(out.df)
  }else{
    g <- ggplot(out.df, aes(x=est, y=labfac, xmin=lower, xmax=upper)) + 
      geom_vline(xintercept=0, lty=2, size=.25, col="gray50") + 
      geom_errorbarh(height=0) + 
      geom_point() + 
      ylab("") + xlab("Estimates Combined") + 
      theme_classic() 
    if(include_est){
      g <- g + geom_text(aes(label=est_text), vjust=0, ...)
    }
    g
  }
}

这是一个包含一些虚构数据的示例：

set.seed(2101)
dat <- data.frame(
  X1 = rnorm(500), 
  X2 = rnorm(500), 
  Z = rnorm(500), 
  W = rnorm(500)
)
dat <- dat %>% 
  mutate(yhat = X1 - X2 + X1*X2 - X1*Z + .5*X2*Z - .75*X1*X2*Z + W, 
         y = yhat + rnorm(500, 0, 1.5))

mod <- lm(y ~ X1*X2*Z + W, data=dat)
plotBeta(mod, position=position_nudge(y=.1), size=3) + xlim(-2.5,2)

编辑：比较两个模型

利用新增的plot=FALSE，我们可以生成数据，然后组合绘图。

mod <- lm(y ~ X1*X2*Z + W, data=dat)
p1 <- plotBeta(mod, plot=FALSE)
mod2 <- lm(y ~ X1*X2 + Z + W, data=dat)
p2 <- plotBeta(mod2, plot=FALSE)
p1 <- p1 %>% mutate(model = factor(1, levels=1:2, 
                                   labels=c("Model 1", "Model 2")))
p2 <- p2 %>% mutate(model = factor(2, levels=1:2, 
                                   labels=c("Model 1", "Model 2")))

p_both <- bind_rows(p1, p2)
p_both <- p_both %>% 
  arrange(var, model) %>% 
  mutate(labfac = factor(1:n(), labels=paste("dy/d", var, sep="")))

ggplot(p_both, aes(x=est, y=labfac, xmin=lower, xmax=upper)) + 
  geom_vline(xintercept=0, lty=2, size=.25, col="gray50") + 
  geom_linerange(position=position_nudge(y=c(-.1, .1))) + 
  geom_point(aes(shape=model), 
             position=position_nudge(y=c(-.1, .1))) + 
  geom_text(aes(label=est_text), vjust=0,
            position=position_nudge(y=c(-.2, .15))) + 
  scale_shape_manual(values=c(1,16)) + 
  ylab("") + xlab("Estimates Combined") + 
  theme_classic()

在 R、lincom + coefplot 或 plotbeta 中将回归系数（偏导数）与 CI 相结合的绘图？

Plot combining regression coefficients (partial derivatives) with CIs in R, lincom + coefplot or plotbeta?

plot

r

derivative

coefficients

coefplot

编辑：比较两个模型