3

我想知道是否有办法在 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对数据运行标准以估计 s 之间的变化差异score时。timegroup

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

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

我想知道的是如何在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 交互的准确估计。group*time为此，需要根据上述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)

我将其传递给contrast.lsmeans

contrast(lsMeansLM, contrasts)

输出中的F和p值与线性趋势的自动测试和 SPSS 中混合 ANCOVA 生成的线性趋势的组差异的值相匹配。然而，混合 ANCOVA 不会生成估计值。

使用上述程序估计效果，而不是大约。-1，就像在lmer（并匹配我在数据中构建的差异）是大约。-10，这是非常不准确的。

我认为这与我编码对比度系数的方式有关。我知道我是否groupMain通过将所有系数除以 4 来归一化矩阵的系数，从而准确估计所有时间点上平均的组的主要影响。但我不知道如何准确估计跨组的平均线性趋势 ( linTrend)，或准确估计跨组的线性趋势差异 ( linXGroup)。

我不确定这个问题是否更适合这里或交叉验证。我首先想到这里是因为它似乎与软件相关，但我知道可能涉及更深层次的问题。任何帮助将非常感激。

4

2 回答 2

2

这里的问题是这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 均值是在 = 1.5 时对每个组的预测timeNum——通常称为调整后的均值；因此，差异就是这两个调整后的平均值之间的差异。

关于在获得大约 10.7 的线性对比度时声称的差异：线性对比度系数c(-3,-1,1,3)为您提供了线斜率的倍数。要获得斜率，您需要除以sum(c(-3,-1,1,3)^2)-- 并乘以 2，因为对比度系数增加 2。

于 2017-08-23T03:00:08.430 回答
0

感谢@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)

并将该主列表传递给contrastslsmeans 中的函数

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。如果我们将group因子的编码更改为简单编码

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

summary(modelLMER)再次运行，time估计现在将是 3.51。

最后对于组的主效应，即所有时间点上组间得分的平均差异。我们可以跑

pairs(lsmeans(modelLM,"group"))

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

于 2017-08-24T05:17:48.707 回答