我想知道是否有办法在 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
时。time
group
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
)。
我不确定这个问题是否更适合这里或交叉验证。我首先想到这里是因为它似乎与软件相关,但我知道可能涉及更深层次的问题。任何帮助将非常感激。