r - 效果与预测函数

Question

我同时尝试理解 R 的 Predict() 函数和“effects”包effect() 函数。本质上，我正在运行回归来测试 DV 上两个二分 IV 的相互作用，同时控制两个连续协变量。在我的实际数据集中，交互作用很重要，所以现在我想绘制交互作用。因为我的模型中有协变量，所以我应该在控制这些其他变量(即 SPSS 中的估计边际均值)后绘制均值。我之前没有在 R 中这样做过，在搜索时我希望我应该能够使用effect() 或predict() 函数获得绘图所需的值。因此，我尝试在随机生成的数据集上对每个数据集进行操作:

> set.seed(100)
> test <- data.frame(iv1 = factor(round(rnorm(200, mean=.5, sd=.25), 0), levels=c(0,1), labels=c("A","B")), iv2 = factor(round(rnorm(200, mean=.5, sd=.25), 0), levels=c(0,1), labels=c("C","D")), cv1 = rnorm(200, mean=4, sd=1), cv2 = rnorm(200, mean=3, sd=1), dv = rnorm(200, mean=5, sd=1))
> mod <- lm(dv ~ cv1 + cv2 + iv1*iv2, data = test)
> new <- with(test, expand.grid(iv1 = levels(iv1), iv2 = levels(iv2), cv1 = mean(cv1), cv2 = mean(cv2)))
> test$pv <- predict(mod, newdata = new)

> tapply(test$pv, list(test$iv1, test$iv2), mean)
         C        D
A 5.076842 5.086218
B 5.025614 5.065399

> effect("iv1:iv2", mod)

 iv1*iv2 effect
   iv2
iv1        C        D
  A 5.019391 5.167275
  B 5.216955 4.855195

因为我得到了不同的结果，所以我将数据导出到 SPSS 并运行方差分析来做同样的事情，并查看估计的边际均值 (EMMEANS)。这些与 R 中的effect()给出的结果相同。

SPSS 语法:

DATASET ACTIVATE DataSet1.
RECODE iv1 iv2 ('A'=-1) ('B'=1) ('C'=-1) ('D'=1) INTO iv1_recode iv2_recode.
EXECUTE.

UNIANOVA dv BY iv1_recode iv2_recode WITH cv1 cv2
  /METHOD=SSTYPE(3)
  /INTERCEPT=INCLUDE
  /EMMEANS=TABLES(OVERALL) WITH(cv1=MEAN cv2=MEAN) 
  /EMMEANS=TABLES(iv1_recode) WITH(cv1=MEAN cv2=MEAN) 
  /EMMEANS=TABLES(iv2_recode) WITH(cv1=MEAN cv2=MEAN) 
  /EMMEANS=TABLES(iv1_recode*iv2_recode) WITH(cv1=MEAN cv2=MEAN) 
  /PRINT=DESCRIPTIVE
  /CRITERIA=ALPHA(.05)
  /DESIGN=cv1 cv2 iv1_recode iv2_recode iv1_recode*iv2_recode.

作为检查，EMMEANS 的 SPSS 输出显示:“模型中出现的协变量按以下值进行评估:cv1 = 3.996208827095569，cv2 = 3.052881951477868。”这些与我在预测中使用的协变量的值相同:

> new
  iv1 iv2      cv1      cv2
1   A   C 3.996209 3.052882
2   B   C 3.996209 3.052882
3   A   D 3.996209 3.052882
4   B   D 3.996209 3.052882

那么我不明白什么？或者我在这里做了一些愚蠢的事情(一种明显的可能性)？这可能是因为我不明白估计的边际平均值是什么。

非常感谢任何帮助!

Answer 1

因此，对于获取模型本身的结果以及将模型应用于观察到的数据，似乎存在一些困惑。这里出现了一个大问题

test$pv <- predict(mod, newdata = new)

在这里，new有 4 行，所以 predict(mod, newdata = new)有值(value)观。运行就可以得到

predict(mod, newdata = new)
#        1        2        3        4 
# 5.019391 5.216955 5.167275 4.855195 

并注意这些值如何与 effect() 的结果匹配.

当您将它们分配给 test$pv 时，长度为 4 的向量被回收，因此它最终沿着 test 重复 50 次数据框。和test实际上包含您观察到的数据，因此将模型的理论预测和观察到的数据混合起来并不是一个好主意。如果您确实想要每个观测值的“真实”预测值，那么 test$pv<-predict(mod)将是正确的选择。然而，总和超过 test ，这又是观测值`，其中

tapply(test$pv, list(test$iv1, test$iv2), mean)

将使用 cv1 的值和cv2这是实际观察到的，而不仅仅是协变量的总体平均值。

我们已经看到 effect()使用协变量的平均值，但您也可以使用显式设置值

effect("iv1:iv2", mod, given.values=c(cv1=3.996209, cv2=3.052882))

如果你愿意的话。