关于R语言中混合线性模型summary()结果中交互作用beta值的含义

关于R语言中混合线性模型summary()结果中交互作用beta值的含义

本文以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 

验证一下:

> (11.2 - 1.67) == 9.53
[1] TRUE

此时交互作用的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起来方便,我懒…)。
编辑于 2019-08-03

文章被以下专栏收录