我同时试图理解 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 并运行 ANOVA 做同样的事情并查看估计的边际均值 (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。”这些与我在 predict 中使用的协变量的值相同:
> 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
那我有什么不明白的？还是我在这里做了一些愚蠢的事情(一种明显的可能性)？这可能是我没有理解估计的边际平均值是什么。
任何帮助是极大的赞赏!
最佳答案
因此，关于获取模型本身的结果以及将模型应用于观察到的数据，这里似乎有些困惑。这里出现了一个大问题
test$pv <- predict(mod, newdata = new)
在这里，
new 有 4 行，所以 predict(mod, newdata = new) 有值。运行只是为了给predict(mod, newdata = new)
# 1 2 3 4
# 5.019391 5.216955 5.167275 4.855195
并注意这些值如何与
effect() 的结果匹配。当您将它们分配给
test$pv 时，长度为 4 的向量会被回收，因此它最终会沿着 test data.frame 重复 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))
如果你喜欢。
关于r - 效果与预测功能，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/24897789/