有没有办法改变 R 在模型输出中标记交互参数的方式?

Is there a way to change the way R labels the interaction parameter in model output?

我遇到了一个看似简单但非常令人沮丧的问题。当您 运行 模型在 R 中具有交互项时,R 将生成的参数命名为 "var1:var2" 等。不幸的是,这种命名约定使我无法计算预测值和 CI 的新数据所在的位置必需,因为“:”不是可以包含在列 header 中的字符,并且原始数据框中的名称必须与 newdata 中的名称完全匹配。还有其他人遇到过这个问题吗?

这是我的代码示例:

wemedist2.exp = glm(survive/trials ~ sitedist + type + sitedist*type + roaddist, family =          binomial(logexp(wemedata$expos)), data=wemedata)
summary(wemedist2.exp)
wemepredict3 = with(wemedata, data.frame(sitedist=mean(sitedist),roaddist=mean(roaddist),    type=factor(1:2))) 
wemepredict3 = cbind(wemepredict3, predict(wemedist2.exp, newdata = wemepredict3, type = "link", se = TRUE))

这会生成一个 table,其中包含每个变量在指定水平的预测值,但没有交互作用。

对于您的 newdata 数据框,您不应包含交互列。调用 predict.

时,将为您计算交互变量的乘积(并乘以估计系数)

例如:

  1. 创建一些虚拟数据:

    set.seed(1)
    n <- 10000
    X <- data.frame(x1=runif(n), x2=runif(n))
    X$x1x2 <- X$x1 * X$x2
    
    head(X)
    #          x1         x2        x1x2
    # 1 0.2655087 0.06471249 0.017181728
    # 2 0.3721239 0.67661240 0.251783646
    # 3 0.5728534 0.73537169 0.421260147
    # 4 0.9082078 0.11129967 0.101083225
    # 5 0.2016819 0.04665462 0.009409393
    # 6 0.8983897 0.13091031 0.117608474
    
    b <- runif(4)
    y <- b[1] + c(as.matrix(X) %*% b[-1]) + rnorm(n, sd=0.1)
    
  2. 拟合模型并比较估计系数与真实系数:

    M <- lm(y ~ x1 * x2, X)
    summary(M)
    
    # Call:
    # lm(formula = y ~ x1 * x2, data = X)
    # 
    # Residuals:
    #      Min       1Q   Median       3Q      Max 
    # -0.43208 -0.06743 -0.00170  0.06601  0.37197 
    # 
    # Coefficients:
    #             Estimate Std. Error t value Pr(>|t|)    
    # (Intercept) 0.202040   0.003906   51.72   <2e-16 ***
    # x1          0.128237   0.006809   18.83   <2e-16 ***
    # x2          0.156942   0.006763   23.21   <2e-16 ***
    # x1:x2       0.292582   0.011773   24.85   <2e-16 ***
    # ---
    # Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    # 
    # Residual standard error: 0.09906 on 9996 degrees of freedom
    # Multiple R-squared:  0.5997,  Adjusted R-squared:  0.5996 
    # F-statistic:  4992 on 3 and 9996 DF,  p-value: < 2.2e-16
    
    b
    # [1] 0.2106027 0.1147864 0.1453641 0.3099322
    
  3. 创建要预测的示例数据,并进行预测。请注意,我们只创建 x1x2,并且 而不是 创建 x1:x2:

    X.predict <- data.frame(x1=runif(10), x2=runif(10))
    
    head(X.predict)
    #           x1        x2
    # 1 0.26037592 0.7652155
    # 2 0.73988333 0.3352932
    # 3 0.02650689 0.9788743
    # 4 0.84083874 0.1446228
    # 5 0.85052685 0.7674547
    # 6 0.13568509 0.9612156
    
    predict(M, newdata=X.predict)
    
    #         1         2         3         4         5         6         7 
    # 0.4138194 0.4221251 0.3666572 0.3681432 0.6225354 0.4084543 0.4711018 
    #         8         9        10 
    # 0.7092744 0.3401867 0.2320834 
    

或者...

另一种方法是通过计算交互项的乘积将交互包含在模型拟合数据中,然后将其也包含在新数据中。我们已经完成了上面第 1 点中的第一步,我们在其中创建了一个名为 x1x2 的列。

然后我们将拟合模型:lm(y ~ x1 + x2 + x1x2, X)

并对以下数据进行预测:

X.predict <- data.frame(x1=runif(10), x2=runif(10), x1x2=runif(10)

如果交互中涉及分类变量...

当您有涉及分类变量的交互时,模型会估计系数来描述属于每个级别相对于属于参考级别的影响。因此,例如,如果我们有一个连续预测变量 (x1) 和一个分类预测变量 (x2,水平为 abc),则模型 y ~ x1 * x2 将估计六个系数,描述:

  1. 截距(即x1为零时的预测y且观测属于x2的参考水平);

  2. 当观测值属于x2参考水平时改变x1的效果(即斜率,对于x2参考水平);

  3. 属于二级的影响(即属于二级相对于属于参考级截距的变化);

  4. 属于第三级的影响(即由于属于第三级相对于属于参考级截距的变化);

  5. x1的效果因属于第二层级相对于属于参考层级的变化(即斜率的变化);和

  6. 由于属于第三级相对于属于参考级x1的影响变化(即斜率变化)

如果您想拟合和预测模型 with/to 描述交互的预计算数据,您可以创建一个包含列的数据框:x1x2b(二进制,表示观察是否属于级别b); x2c(二进制,表示观察是否属于级别c); x1x2bx1x2b 的乘积);和 x1x2cx1x2c 的乘积)。

一个快速的方法是 model.matrix:

set.seed(1)
n <- 1000
d <- data.frame(x1=runif(n), x2=sample(letters[1:3], n, replace=TRUE))

head(d)
#          x1 x2
# 1 0.2655087  b
# 2 0.3721239  c
# 3 0.5728534  b
# 4 0.9082078  c
# 5 0.2016819  a
# 6 0.8983897  a

X <- model.matrix(~x1*x2, d)

head(X)
#   (Intercept)        x1 x2b x2c    x1:x2b    x1:x2c
# 1           1 0.2655087   1   0 0.2655087 0.0000000
# 2           1 0.3721239   0   1 0.0000000 0.3721239
# 3           1 0.5728534   1   0 0.5728534 0.0000000
# 4           1 0.9082078   0   1 0.0000000 0.9082078
# 5           1 0.2016819   0   0 0.0000000 0.0000000
# 6           1 0.8983897   0   0 0.0000000 0.0000000

b <- rnorm(6) # coefficients
y <- X  %*% b + rnorm(n, sd=0.1)

您可以将 X 的列重命名为您想要的任何名称,只要您在稍后 predict 将模型转换为新数据时使用一致的命名即可。

现在拟合模型。在这里,我告诉 lm 不要计算截距(使用 -1),因为变量 (Intercept) 已经存在于 X 中,并且将为它计算一个系数。我们也可以通过拟合数据 as.data.frame(X[, -1]):

来做到这一点
(M <- lm(y ~ . - 1, as.data.frame(X)))

# Call:
# lm(formula = y ~ . - 1, data = as.data.frame(X))
# 
# Coefficients:
# `(Intercept)`          x1         x2b         x2c    `x1:x2b`    `x1:x2c`  
#       1.14389     1.09168    -0.88879     0.20405     0.09085    -1.63769  

创建一些新数据进行预测,并执行预测:

d.predict <- expand.grid(x1=seq(0, 1, 0.1), x2=letters[1:3])
X.predict <- model.matrix(~x1*x2, d.predict)
y.predict <- predict(M, as.data.frame(X.predict))