机器算法验证 - 如何通过 R 中的连续交互获得多个连续水平的斜率和标准误差？ - 吾爱随笔录

如何通过 R 中的连续交互获得多个连续水平的斜率和标准误差？

机器算法验证 r 回归数据可视化相互作用

2022-04-17 16:17:51

我将几个不同响应变量（DV；代表不同人群）的斜率与一组预测变量（IV）进行比较。对于某些 DV，支持 2 向交互（连续通过连续）。为了便于比较 IV 系数，我想在单个图表上绘制斜率估计值和 95% CI（每个 IV 的单独图表），对于具有交互作用的 DV，我想将斜率绘制在 ~3连续调节变量的值（例如，下图中的“DV 1”）。

在此处输入图像描述

我确信有多种方法可以获取这些值，但我希望有人可以为我指出一些简单的代码或一个可以帮助我自动化这个过程的包。我还应该注意我的模型来自 lme4。

“效果”包可以方便地计算用户指定的调节变量级别的预测值，但据我所知没有提供斜率或 SE（虽然我可以从预测值中算出这些，但我希望有更多的流内衬法）。

这是一些玩具数据，虽然它不会产生我在图中显示的交互；

set.seed(50)
x1 <- rnorm(100,2,10)
x2 <- rnorm(100,2,10)
y1 <- x1+x2+x1*x2+rnorm(100,0,100)

model1<-lm(y1 ~ x1*x2)

这是从“效果”绘制的预测值，但我想要这些线的斜率和 SE...

library(effects)
model1.eff<-effect("x1*x2",model1,xlevels=3)
plot(model1.eff,multiline=T,ci.style="bands")
as.data.frame(model1.eff)

2个回答

为了检查一个连续变量的不同水平的简单斜率，您可以简单地将另一个连续变量居中以关注感兴趣的斜率。在具有连续通过连续交互的模型中，如下所示：

y = β_{0} + β_{1} x_{1} + β_{2} x_{2} + β_{3} x_{1} * x_{2}

$y = \beta_0 + \beta_1x_1 + \beta_2x_2 + \beta_3x_1*x_2$ 两个单预测系数（

β_{1}

$\beta_1$ 和

β_{2}

$\beta_2$ ) 是当另一个预测器（不管它是居中的）等于 0 时预测器的简单斜率。

所以，如果我在上面运行你的练习代码，我会得到以下输出：

Call:
lm(formula = y1 ~ x1 * x2)

Residuals:
     Min       1Q   Median       3Q      Max 
-281.996  -70.148   -3.702   70.190  209.182 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  17.7519    10.8121   1.642    0.104    
x1            1.4175     1.0151   1.397    0.166    
x2            0.8222     1.0614   0.775    0.440    
x1:x2         0.8911     0.1295   6.882 6.04e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 100.6 on 96 degrees of freedom
Multiple R-squared:  0.4283,    Adjusted R-squared:  0.4105 
F-statistic: 23.98 on 3 and 96 DF,  p-value: 1.15e-11

x1 输出为我们提供了 x2 = 0 时 x1 斜率的测试。因此，我们得到了一个斜率、标准误差和（作为奖励）与 0 比较的参数估计值的测试。如果我们想得到简单的斜率当 x2 = 6 时，x1（以及标准误差和 sig.test），我们简单地使用线性变换使 x2 上的值 6 成为 0 点：

x2.6<- x2-6

通过查看汇总统计信息，我们可以看到这与之前完全相同的变量，但它在数轴上向下移动了 6 个单位：

summary(x2)
summary(x2.6)

 > summary(x2)
   Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-31.0400  -5.9520   1.3430   0.8396   8.0090  22.3800 

 > summary(x2.6)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-37.040 -11.950  -4.657  -5.160   2.009  16.380

现在，如果我们重新运行相同的模型，但用 x2 代替我们新居中的变量 x2.6，我们会得到：

model1.6<- lm(y1~x1*x2.6)
summary(model1.6)


Call:
lm(formula = y1 ~ x1 * x2.6)

Residuals:  
     Min       1Q   Median       3Q      Max 
-281.996  -70.148   -3.702   70.190  209.182 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  22.6853    12.6384   1.795   0.0758 .  
x1            6.7639     1.2346   5.479 3.44e-07 ***
x2.6          0.8222     1.0614   0.775   0.4404    
x1:x2.6       0.8911     0.1295   6.882 6.04e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 100.6 on 96 degrees of freedom
Multiple R-squared:  0.4283,    Adjusted R-squared:  0.4105 
F-statistic: 23.98 on 3 and 96 DF,  p-value: 1.15e-11

如果我们将此输出与旧输出进行比较，我们可以看到综合 F 仍然是 23.98，交互作用 t 仍然是 6.882，x2.6 的斜率仍然是 0.822（并且不显着）。然而，我们的 x1 系数现在更大且显着。当 x2 等于 6（或当 x2.6 = 0 时）时，这个斜率现在是 x1 的简单斜率。通过以几个不同的变量为中心，我们可以测试几个不同的简单效果（并获得斜率和标准误差）而无需太多工作。通过使用（在 R 社区中很可怕的）for 循环来遍历列表，我们可以非常有效地测试几种不同的简单效果：

centeringValues<- c(1,2,3,4,5,6) # Creating a vector of values to center around

for(i in 1:length(centeringValues)){     #Making a for loop that iterates through the list
  x<- x2-i         # Creating a predictor that is the newly centered variable
  print(paste0('x.',centeringValues[i])) # printing x.centering value so you can keep track of output
  print(summary(lm(y1~x1*x))[4]) # printing coefficients for the model with the center variable

}

此代码首先创建一个值向量，该向量希望成为不希望斜率的变量（在本例中为 x2）的 0 点。接下来，创建一个循环遍历此列表中的位置（即，如果列表有 3 个项目，for 循环将遍历值 1 到 3）。接下来，创建一个新变量，它是您不希望斜率居中的变量的居中版本（在这种情况下，我们对 x1 的简单斜率感兴趣，因此我们将 x2 居中）。最后，打印模型中的系数，其中包含新居中的变量代替原始变量。这将产生以下输出：

[1] "x.1"
    $coefficients
              Estimate Std. Error   t value     Pr(>|t|)
(Intercept) 18.5741364 10.8815154 1.7069439 9.106513e-02
x1           2.3085985  1.0143100 2.2760286 2.506664e-02
x            0.8222252  1.0613590 0.7746909 4.404262e-01
x1:x         0.8910530  0.1294695 6.8823366 6.041102e-10

[1] "x.2"
$coefficients
              Estimate Std. Error   t value     Pr(>|t|)
(Intercept) 19.3963616 11.0528627 1.7548722 8.247158e-02
x1           3.1996515  1.0299723 3.1065415 2.489385e-03
x            0.8222252  1.0613590 0.7746909 4.404262e-01
x1:x         0.8910530  0.1294695 6.8823366 6.041102e-10

[1] "x.3"
$coefficients
              Estimate Std. Error   t value     Pr(>|t|)
(Intercept) 20.2185867 11.3215341 1.7858522 7.728065e-02
x1           4.0907045  1.0613132 3.8543802 2.096928e-04
x            0.8222252  1.0613590 0.7746909 4.404262e-01
x1:x         0.8910530  0.1294695 6.8823366 6.041102e-10

[1] "x.4"
$coefficients
              Estimate Std. Error   t value     Pr(>|t|)
(Intercept) 21.0408119 11.6808159 1.8013135 7.479290e-02
x1           4.9817575  1.1070019 4.5002249 1.905339e-05
x            0.8222252  1.0613590 0.7746909 4.404262e-01
x1:x         0.8910530  0.1294695 6.8823366 6.041102e-10

[1] "x.5"
$coefficients
              Estimate Std. Error   t value     Pr(>|t|)
(Intercept) 21.8630371 12.1226545 1.8034859 7.444873e-02
x1           5.8728105  1.1653521 5.0395160 2.193149e-06
x            0.8222252  1.0613590 0.7746909 4.404262e-01
x1:x         0.8910530  0.1294695 6.8823366 6.041102e-10

[1] "x.6"
$coefficients
              Estimate Std. Error   t value     Pr(>|t|)
(Intercept) 22.6852623 12.6383944 1.7949481 7.580894e-02
x1           6.7638636  1.2345698 5.4787212 3.439867e-07
x            0.8222252  1.0613590 0.7746909 4.404262e-01
x1:x         0.8910530  0.1294695 6.8823366 6.041102e-10

在这里，您可以看到输出提供了多个测试的系数，但每次唯一变化的是 x1 的斜率。每个输出中 x1 的斜率表示当 x2 等于我们为该迭代分配的任何中心值时 x1 的斜率。希望这可以帮助！

虽然@wools 的答案似乎绰绰有余，但这是另一种选择，它允许从单个模型输出计算给定 x2 的 x1 的边际效应，而无需将 x 变量居中；

根据http://statistics.ats.ucla.edu/stat/r/faq/concon.htm；模型在哪里

y ~ β0 + β1x1 + β2x2+ β3x1∗x2

那么在给定 x2 值处 x1 的斜率是 β1 + β3 * x2

所以我可以选择 x2 的几个值作为；

at.x2<-c(-6, 1, 6)

slopes <- coef(model1)["x1"] + coef(model1)["x1:x2"] * at.x2

根据如何计算交互作用中边际效应的标准误差（稳健回归）？斜率的标准误差 = sqrt(var(b1) + var(b3) x2^2 + 2 x2 * cov(b1,b3) )

estvar<-vcov(model1); model1.vcov<-as.data.frame(as.matrix(estvar))
var.b1<-model1.vcov["x1","x1"]
var.b3<-model1.vcov["x1:x2","x1:x2"]
cov.b1.b3<-model1.vcov["x1","x1:x2"]

SEs <- rep(NA, length(at.x2))
for (i in 1:length(at.x2)){
  j <- at.x2[i]  
  SEs[i] <- sqrt(var.b1 + var.b3 * j^2 + 2*j* cov.b1.b3)
}

cbind(SEs, slopes)

其它你可能感兴趣的问题

上一篇是否可以在给定 95% 置信区间的情况下从优势比中恢复 2x2 表的原始数量？下一篇说“主题是单词的分布”是什么意思？