在 afex、lsmeans 和 lme4 包中生成类似的交互估计
Generating similar estimates of interactions in afex, lsmeans, and lme4 packages
我想知道是否有办法在 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 的数据寻找 score 相对于 time 之间 group 之间变化差异的估计.
summary(modelLMER <- lmer(score ~ group * timeNum + (timeNum|id), df))
我得到了 -1.07 的 group*time 相互作用的估计值，这意味着在 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)
我在 lsmeans 中传递给 contrast 函数。
contrast(lsMeansLM, contrasts)
输出中的F和p值与线性趋势的自动检验和线性趋势中的组差异相匹配由 SPSS 中的混合 ANCOVA 生成。但是，混合 ANCOVA 不会生成估计值。
使用上述程序估算的效果，而不是近似值。 -1，就像在 lmer 中一样（并匹配我内置到数据中的差异）大约是。 -10，这是非常不准确的。
我认为这与我对对比度系数的编码方式有关。我知道如果我通过将所有系数除以 4 来归一化 groupMain 矩阵的系数，则可以准确估计所有时间点平均组的主要影响。但我不知道如何准确估计跨组平均线性趋势 (linTrend)、 或 准确估计跨组线性趋势差异 (linXGroup).
我不确定这个问题更适合这里还是交叉验证。我首先想到这里是因为它似乎与软件相关，但我知道可能涉及更深层次的问题。任何帮助将不胜感激。
这里的问题是 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.
感谢@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)
'main effect' 趋势（因为缺少更好的描述符），它是四个时间点的分数变化率 两个级别组的平均值 , 为 3.51。如果我们通过
将 group 因子的编码更改为简单编码
contrasts(df$group) <- c(-.5,.5)
和 运行 summary(modelLMER)，time 估计现在是 3.51。
最后是组的主效应，即各时间点平均的组间得分差异。我们可以运行
pairs(lsmeans(modelLM,"group"))
这将是 -1.92。谢谢@rvl。一个很好的答案。使用 afex 和 lsmeans 我们现在强制执行混合 ANCOVA，将重复测量变量视为分类变量，以便我们估计趋势和主要影响的组差异，这些差异与混合效应模型返回的结果相匹配，其中重复测量变量是连续的，and with p- and F-values 匹配 SPSS .
我想知道是否有办法在 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 的数据寻找 score 相对于 time 之间 group 之间变化差异的估计.
summary(modelLMER <- lmer(score ~ group * timeNum + (timeNum|id), df))
我得到了 -1.07 的 group*time 相互作用的估计值，这意味着在 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)
我在 lsmeans 中传递给 contrast 函数。
contrast(lsMeansLM, contrasts)
输出中的F和p值与线性趋势的自动检验和线性趋势中的组差异相匹配由 SPSS 中的混合 ANCOVA 生成。但是，混合 ANCOVA 不会生成估计值。
使用上述程序估算的效果，而不是近似值。 -1，就像在 lmer 中一样（并匹配我内置到数据中的差异）大约是。 -10，这是非常不准确的。
我认为这与我对对比度系数的编码方式有关。我知道如果我通过将所有系数除以 4 来归一化 groupMain 矩阵的系数，则可以准确估计所有时间点平均组的主要影响。但我不知道如何准确估计跨组平均线性趋势 (linTrend)、 或 准确估计跨组线性趋势差异 (linXGroup).
我不确定这个问题更适合这里还是交叉验证。我首先想到这里是因为它似乎与软件相关，但我知道可能涉及更深层次的问题。任何帮助将不胜感激。
这里的问题是 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.
感谢@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)
'main effect' 趋势（因为缺少更好的描述符），它是四个时间点的分数变化率 两个级别组的平均值 , 为 3.51。如果我们通过
将group 因子的编码更改为简单编码
contrasts(df$group) <- c(-.5,.5)
和 运行 summary(modelLMER)，time 估计现在是 3.51。
最后是组的主效应，即各时间点平均的组间得分差异。我们可以运行
pairs(lsmeans(modelLM,"group"))
这将是 -1.92。谢谢@rvl。一个很好的答案。使用 afex 和 lsmeans 我们现在强制执行混合 ANCOVA，将重复测量变量视为分类变量，以便我们估计趋势和主要影响的组差异，这些差异与混合效应模型返回的结果相匹配，其中重复测量变量是连续的，and with p- and F-values 匹配 SPSS .