有沒有辦法改變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))