交互和聚类的回归系数

Question

我正在运行一个具有交互项和双聚类的线性回归模型，并且很难获得 R 中的组合回归系数和 95% 置信区间。我想知道是否有人知道 R 包或其他解决方法获取此信息。

这是一些虚拟数据和我的基本模型的示例：

dat <- data.frame(a=rnorm(100, 0, 1),
                  b=as.factor(rbinom(100, 1, 0.2)),
                  c=as.factor(rbinom(100, 1, 0.6)),
                  y=rnorm(100, 0, 1),
                  clust1=as.factor(rbinom(100, 10, 0.1)),
                  clust2=as.factor(rbinom(100, 5, 0.5)))

mod <- lm(y ~ a*b*c, data=dat)

在我的数据中，y和a是连续的，而b和c是二进制的。然后我用以下方法解释双聚类：

library(lmtest)
library(sandwich)

coeftest(mod, vcov=vcovCL, cluster= ~clust1 + clust2)

这给了我估计值和标准误差。我需要知道二元预测变量的所有组合的组合回归参数和 95% 置信区间：

当b=0且c=0时；当b=1且c=1时；当b=0且c=1时；当 b=1 且 c=0 时。

我尝试使用biostat3包中的lincom：https://www.rdocumentation.org/packages/biostat3/versions/0.1.9/topics/lincom，它基于汽车的线性假设：https:// www.rdocumentation.org/packages/car/versions/3.1-2/topics/linearHypothesis。此示例适用于 b=0 和 c=1 的情况：

library(biostat3)

lincom(mod, vcov=vcovCL, cluster= ~clust1 + clust2, "a+c+a:c")

这看起来很有希望，但是关于这个命令的文档并不多，所以我不清楚它是否完全符合我的要求。

有谁知道有什么方法可以做到这一点吗？或者，可以确认biostat3包中的lincom是正确的吗？预先感谢您，如果我可以扩展任何内容，请告诉我。

Answer 1

我不太清楚“组合回归参数”是什么意思。在任何情况下，您几乎肯定可以使用

marginaleffects

包计算感兴趣的数量（免责声明：我是维护者）。该网站有超过 30 章的详细教程：https://marginaleffects.com

明显的起点是此处的“入门”页面：https://marginaleffects.com/vignettes/get_started.html

数据与模型

options(width = 10000)
dat <- data.frame(a=rnorm(100, 0, 1),
                  b=as.factor(rbinom(100, 1, 0.2)),
                  c=as.factor(rbinom(100, 1, 0.6)),
                  y=rnorm(100, 0, 1),
                  clust1=as.factor(rbinom(100, 10, 0.1)),
                  clust2=as.factor(rbinom(100, 5, 0.5)))

mod <- lm(y ~ a*b*c, data=dat)

预测

您可以获得模型中每个预测变量组合的拟合值。在

marginaleffects

中，我们称之为“预测”：

https://marginaleffects.com/vignettes/predictions.html

library(marginaleffects)

predictions(
    mod,
    newdata = datagrid(
        a = mean,
        b = unique,
        c = unique)
)
# 
#       a b c Estimate Std. Error       z Pr(>|z|)   S  2.5 % 97.5 %
#  0.0194 0 0 -0.18668      0.157 -1.1863    0.235 2.1 -0.495  0.122
#  0.0194 0 1 -0.00692      0.138 -0.0501    0.960 0.1 -0.278  0.264
#  0.0194 1 0  0.02691      0.415  0.0648    0.948 0.1 -0.787  0.841
#  0.0194 1 1 -0.15562      0.223 -0.6985    0.485 1.0 -0.592  0.281
# 
# Columns: rowid, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, y, a, b, c 
# Type:  response

比较

您可以衡量将

a

改变 1 个单位的“效果”，同时在每一种组合中保持其他预测变量不变。在

marginaleffects

包中，我们称之为“比较”：

https://marginaleffects.com/vignettes/comparisons.html

comparisons(
    mod,
    variables = "a",
    newdata = datagrid(
        b = unique,
        c = unique)
)
# 
#  Term Contrast b c Estimate Std. Error      z Pr(>|z|)   S  2.5 % 97.5 %      a
#     a       +1 0 0  -0.1190      0.153 -0.777    0.437 1.2 -0.419  0.181 0.0194
#     a       +1 0 1  -0.0413      0.114 -0.361    0.718 0.5 -0.265  0.183 0.0194
#     a       +1 1 0   0.1566      0.440  0.356    0.722 0.5 -0.705  1.018 0.0194
#     a       +1 1 1   0.0315      0.219  0.144    0.886 0.2 -0.398  0.461 0.0194
# 
# Columns: rowid, term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, b, c, predicted_lo, predicted_hi, predicted, y, a 
# Type:  response