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

如何解决绘制将回归系数偏导数与 R、lincom + coefplot 或 plotbeta 中的 CI 相结合的图？

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

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

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

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

系数图（我有的）：

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

我发现 Stata 有一个函数可以完成我正在寻找的功能，称为“plotbeta”。 R 有类似的东西吗？

解决方法

这是一个开始。这定义了一个名为 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,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,"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,xmax=upper)) + 
  geom_vline(xintercept=0,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),position=position_nudge(y=c(-.2,.15))) + 
  scale_shape_manual(values=c(1,16)) + 
  ylab("") + xlab("Estimates Combined") + 
  theme_classic()