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

Question

我有一个看似简单但非常令人沮丧的问题。当你在R中运行一个带有交互项的模型时，R命名生成的参数“var1：var2”等。不幸的是，这个命名约定阻止我计算需要newdata的预测值和CI，因为“：”不是字符可以包含在列标题中，并且原始数据框中的名称必须与新数据中的名称完全匹配。有没有其他人有这个问题？

这里是我的代码的示例：

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)) 

这产生一个表，用于在指定的电平的每个变量的预测值，但不相互作用。

Answer 1

对于您的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. 创建示例数据预测并预测。请注意，我们只创造x1和x2，做不创建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，含量a，b和c），那么模型y ~ x1 * x2将估计的六个系数，描述：

1. 截距（即当x1为零且观察属于参考水平x2时，预测的y）;

2. x1变化时观察属于x2的参考电平（即，斜率，对于x2参考电平）的效果;

3. 属于第二级别的效果（即由于属于第二级别的截距相对于属于参考级别的变化）;

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

5. x1（即斜率的变化）归属于第二级相对于属于参考级的变化;和

6. 由于属于第三级而相对于属于参考级的x1（即斜率变化）的影响的变化。

如果你想以适应和预测与/描述的相互作用预先计算的数据模型，您可以创建一个数据帧，其中包括列：x1; x2b（二进制，表示观察是否属于b）; x2c（二进制，表示观察是否属于c）; x1x2b（x1和x2b的乘积）;和x1x2c（x1和x2c的产品）。

一个快速的方法来做到这一点是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))