本文以2*2的实验设计为例，利用lmerTest包在R中进行混合线性模型分析，采用sum的因子编码方式，简单介绍一下在summary的结果中，交互作用的beta值的含义。
数据准备：
library(tidyverse);
library(lmertest)
DF = read_csv('https://raw.githubusercontent.com/usplos/Eye-movement-related/master/DataforShiny.csv')
DF[c('A','B')] = lapply(DF[c('A','B')], factor)
数据中，A和B是自变量（first level），Y是因变量，Sub和Item是群组变量（second level）
采用sum编码方式，有两种操作方式，下面分别介绍：
在option()中直接设置
options(contrasts = c('contr.sum','contr.poly')) # 这里需要同时依次设置无序因子和有序因子的编码方式查看因子A和因子B的编码方式：
> contrasts(DF$A)
[,1]
A1 1
A2 -1
############################
> contrasts(DF$B)
[,1]
B1 1
B2 -1
即A因素中以A2为基线，B因素中以B2为基线。
建立模型：
Model = lmer(data = DF, Y~A*B+(1|Sub)+(1|Item)) # 这里采用最简单的模型为例获取summary中固定效应的结果
> summary(Model) %>% coef(.) %>% round(x = .,digits = 3)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 248.476 10.116 9.671 24.563 0.000
A1 -3.218 2.551 1092.286 -1.262 0.207
B1 -3.906 2.550 1093.941 -1.531 0.126
A1:B1 -2.382 3.028 58.248 -0.787 0.435
即A因素beta值为3.218，B因素为3.906，AB交互作用为2.382。
注意，这里A的beta值不是A1和A2均值差距，而是 (A1-A2)/2
进行A因素的事后检验，验证上面的结论：
> emmeans::emmeans(Model, pairwise~A) %>% .$contr
NOTE: Results may be misleading due to involvement in interactions
contrast estimate SE df t.ratio p.value
A1 - A2 -6.44 5.12 1091 -1.257 0.2091
Results are averaged over the levels of: B可以看到：6.44 = 2 * 3.22.
同理，B的beta值为(B1-B2)/2.（均值差的一半）。
那么交互作用呢？beta值为2.382，代表什么意思？这里进行简单效应分析：
> emmeans::emmeans(Model, pairwise~A|B) %>% .$contr
B = B1:
contrast estimate SE df t.ratio p.value
A1 - A2 -11.20 8.00 146 -1.400 0.1637
B = B2:
contrast estimate SE df t.ratio p.value
A1 - A2 -1.67 7.86 136 -0.213 0.8318对交互作用了解比较深刻的人似乎可以察觉到 2.382 和 11.2、1.67这两个数值之间有一些关系。再细心一点，就可以看到下面的关系：
> (11.2-1.67)/4
[1] 2.3825
诶？11.2–1.67 代表的是条件A在条件B的不同水平上效应的差异，也就是交互作用在数值上的体现。而summary结果中的2.382是这一数值的1/4，为什么是1/4是想必大家也懂了——既然A和B的beta值是其效应量的一半，那么交互作用肯定是其效应量的1/4了。
那么有什么办法可以直接从summary中获取真实的效应量，而不是减半呢？有！
在建模命令中自定义设置为sum编码
重新建立模型：
Model = lmer(data = Df,
Y~A*B+(1|Sub)+(1|Item),
contrasts = list(A = contr.sum(2)/2,
B = contr.sum(2)/2))再来看一下A和B的编码方式：
> contr.sum(2)/2
[,1]
1 0.5
2 -0.5
即从原来的1和-1，变成了0.5和-0.5.
获取summary中固定效应的结果：
> summary(Model) %>% coef(.) %>% round(x = .,digits = 3)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 248.476 10.116 9.671 24.563 0.000
A1 -6.437 5.102 1092.286 -1.262 0.207
B1 -7.811 5.101 1093.941 -1.531 0.126
A1:B1 -9.530 12.110 58.248 -0.787 0.435再来做一下A和B的事后检验：
> emmeans::emmeans(Model, pairwise~A)$contr
contrast estimate SE df t.ratio p.value
A1 - A2 -6.44 5.12 1091 -1.257 0.2091
>
> emmeans::emmeans(Model, pairwise~B)$contr
contrast estimate SE df t.ratio p.value
B1 - B2 -7.81 5.12 1092 -1.526 0.1274
可以看到，此时的beta值等于对应的效应量。
再做一下简单效应：
> emmeans::emmeans(Model, pairwise~A|B)$contr
B = B1:
contrast estimate SE df t.ratio p.value
A1 - A2 -11.20 8.00 146 -1.400 0.1637
B = B2:
contrast estimate SE df t.ratio p.value
A1 - A2 -1.67 7.86 136 -0.213 0.8318验证一下：
此时交互作用的beta值也等于交互作用的效应量。
总结一下：
交互作用其实是看因素A在因素B不同水平上的效应量的差异；
sum编码下，summary中固定效应的效应量不一定是真实的效应，应是编码方式而定；
有两种sum编码的方式：第一种是在环境中声明——option()命令中设置，第二种是建模中自定义设定(通过设置contrasts参数，contr.sum(n)/n)；
第一种设置下，主效应的 beta值 为真实效应量的1/2，交互作用的 beta 值为真实效应量的1/4；
第二种设置下，主效应和交互作用的beta值等于真实效应量；
编码类型的不同影响anova和summary中p值的一致性；
sum编码方式下，矩阵的数值影响summary中固定效应beta值和真实效应量的一致性；
Shiny Performer中默认为第一种设置方式（因为第一种typing codes起来方便，我懒…）。
1671
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
