在afex、lsmeans和lme4软件包中生成类似交互作用的估计

Question

我想知道在afex和lsmeans包中是否可以像在lmer中一样获得交互作用效应的相同估计值。下面是针对两个具有不同截距和斜率的组的玩具数据。

set.seed(1234)
A0 <- rnorm(4,2,1)
B0 <- rnorm(4,2+3,1)
A1 <- rnorm(4,6,1)
B1 <- rnorm(4,6+2,1)
A2 <- rnorm(4,10,1)
B2 <- rnorm(4,10+1,1)
A3 <- rnorm(4,14,1)
B3 <- rnorm(4,14+0,1)
score <- c(A0,B0,A1,B1,A2,B2,A3,B3)
id <- factor(rep(1:8,times = 4, length = 32))
time <- factor(rep(0:3, each = 8, length = 32))
timeNum <- as.numeric(rep(0:3, each = 8, length = 32))
group <- factor(rep(c("A","B"), times =2, each = 4, length = 32))
df <- data.frame(id, group, time, timeNum, score)
df  

这里是情节

(ggplot(df, aes(x = time, y = score, group = group)) + 
    stat_summary(fun.y = "mean", geom = "line", aes(linetype = group)) +
    stat_summary(fun.y = "mean", geom = "point", aes(shape = group), size = 3) +
    coord_cartesian(ylim = c(0,18)))

当我在数据上运行标准的lmer，寻找group之间score变化差异的估计值时。

summary(modelLMER <- lmer(score ~ group * timeNum + (timeNum|id), df))

我得到了一个group*time的相互作用的估计值为-1.07，这意味着时间增加一单位时，在group B 比 group A 得分增加少约1分。这个估计值与我在数据集中设置的差异相匹配。

我想知道如何在afex 和 lsmeans 包中做类似的事情。

library(afex)
library(lsmeans)

首先，我生成了afex模型对象。

modelLM <- aov_ez(id="id", dv="score", data=df, between="group", within="time", 
                 type=3, return="lm")

然后将其传递到lsmeans函数中

lsMeansLM <- lsmeans(modelLM, ~rep.meas:group) 

我的目标是在afex和lsmeans中生成一个准确的组*时间交互估计值。为了实现这一点，需要根据上述lsmeans函数中指定的拆分来指定自定义对比矩阵。

groupMain = list(c(-1,-1,-1,-1,1,1,1,1)) # group main effect
linTrend = list(c(-3,-1,1,3,-3,-1,1,3)) # linear trend
linXGroup = mapply("*", groupMain, linTrend) # group x linear trend interaction

然后我制作了一个总清单。

contrasts <- list(groupMain=groupMain, linTrend=linTrend, linXGroup=linXGroup)

我把它传递到lsmeans中的contrast函数。

contrast(lsMeansLM, contrasts)

输出中的F和p值与在SPSS中生成的混合ANCOVA中的线性趋势自动测试及组别差异的线性趋势相匹配。但是，混合ANCOVA并不生成估计值。

使用上述过程得出的效果估计值，与lmer中的近似-1（与我内置数据的差异匹配）不同，大约为-10，误差极大。

我认为这与我如何编写对比系数有关。我知道如果将groupMain矩阵的系数归一化，即通过将所有系数除以四来获得跨所有时间点平均的组主效应的准确估计值。但我不知道如何获得跨组别平均的线性趋势(linTrend)或跨组别的线性趋势差异(linXGroup)的准确估计值。

我不确定这个问题更适合发布在哪里，是在这里还是Cross Validated。我首先考虑了这里，因为它似乎与软件有关，但我知道可能涉及更深层次的问题。非常感谢任何帮助。

- llewmills

Answer 1

这里的问题在于timeNum是一个数值型预测变量。因此，交互作用是斜率的比较。请注意：

> lstrends(modelLMER, ~group, var = "timeNum")
 group timeNum.trend       SE  df lower.CL upper.CL
 A          4.047168 0.229166 6.2 3.490738 4.603598
 B          2.977761 0.229166 6.2 2.421331 3.534191

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
> pairs(.Last.value)
 contrast estimate        SE  df t.ratio p.value
 A - B    1.069407 0.3240897 6.2     3.3  0.0157

这里是你的1.07 - 符号相反是因为比较方向不同。

我将进一步解释，你在问题中描述的lsmeans结果是两个组均值的比较，而不是交互作用对比。 lsmeans使用参考网格：

> ref.grid(modelLMER)
'ref.grid' object with variables:
    group = A, B
    timeNum = 1.5

正如您所看到的，timeNum 被固定在其平均值为 1.5。LS 均值是每个组在 timeNum = 1.5 时的预测值，通常称为调整均值；因此差异就是这两个调整均值之间的差异。

关于获得线性对比约为 10.7 的差异：线性对比系数 c(-3,-1,1,3) 给出了线条斜率的 倍数。要获得斜率，需要除以 sum(c(-3,-1,1,3)^2)，并且还要乘以 2，因为对比系数递增 2。

- Russ Lenth

---

感谢您的回复@rvl。很高兴学到两个新函数lstrends()和.Last.value，这是我之前不知道的。在modelLM对象中，我使用了一个因子变量time，而不是timeNum变量。我想知道是否有办法生成使用modelLM对象和我在lsmeans中指定的对比矩阵之间的线性趋势差异的准确估计（即而不是使用modelLMER对象）？ - llewmills

lsmeans(lsmeansLM, interaction = c("poly","pairwise") 是怎么样的呢？ - Russ Lenth

@llewmills 注意，我在我的答案末尾添加了线性对比的讨论。 - Russ Lenth

我猜测可能是缺少了一个右括号。否则的话，它应该可以工作。 - Russ Lenth

糟糕，应该是contrast()而不是lsmeans()。 - Russ Lenth

Answer 2

感谢 @rvl 的宝贵帮助，我成功解决了这个问题。以下是代码。

为了生成正确的对比矩阵，我们首先需要对它们进行归一化处理。

(mainMat <- c(-1,-1,-1,-1,1,1,1,1)) # main effects matrix
(trendMat <- c(-3,-1,1,3,-3,-1,1,3) # linear trend contrast coefficients 
(nTimePoints <- 4) # number of timePoints
(mainNorm <- 1/nTimePoints)
(nGroups <- 2) # number of between-Ss groups
(trendIncrem <- 2) # the incremental increase of each new trend contrast coefficient    
(trendNorm <- trendIncrem/(sum(trendMat^2))) # normalising the trend coefficients

现在我们以列表形式创建几个对比矩阵。我们使用上面创建的对象进行归一化处理。

(groupMain = list(mainMat*mainNorm)) # normalised group main effect
(linTrend = list(trendMat*trendNorm)) # normalised linear trend
(linXGroup = list((mainMat*trendMat)*(nGroups*trendNorm))) # group x linear trend interaction

现在将这些矩阵列表传递到主列表中。

contrasts <- list(groupMain=groupMain, linTrend=linTrend, linXGroup=linXGroup)

然后将该主列表传递到lsmeans中的contrasts函数

contrast(lsMeansLM, contrasts)

这是输出

 contrast                                               estimate        SE df t.ratio p.value
 c(-0.25, -0.25, -0.25, -0.25, 0.25, 0.25, 0.25, 0.25)  1.927788 0.2230903  6   8.641  0.0001
 c(-0.15, -0.05, 0.05, 0.15, -0.15, -0.05, 0.05, 0.15)  3.512465 0.1609290  6  21.826  <.0001
 c(0.3, 0.1, -0.1, -0.3, -0.3, -0.1, 0.1, 0.3)         -1.069407 0.3218581  6  -3.323  0.0160

我们如何检查这些估计值是否准确？

首先注意，现在group*time交互的估计值与返回值大致相同。

summary(modelLMER)

“主效应”趋势（暂无更好的描述），即在四个时间点上得分变化率在两个组别水平上平均为3.51。如果我们通过简单编码改变组别因素的编码方式，则：

contrasts(df$group) <- c(-.5,.5)

然后再次运行summary(modelLMER)，此时time的估计值将为3.51。

最后，对于组的主要效应，即在所有时间点上平均分数之间的差异。我们可以运行

pairs(lsmeans(modelLM,"group"))

这将是-1.92。谢谢@rvl。一个很棒的答案。使用afex和lsmeans，我们现在强制进行了一个将重复测量变量视为分类变量的混合ANCOVA，以给出与将重复测量变量视为连续变量的混合效应模型返回的趋势差异和主效应估计相匹配的组别差异估计，并且具有与SPSS相匹配的p值和F值。

- llewmills