bootstrap-table真实交互数据_关于R语言中混合线性模型summary()结果中交互作用beta值的含义...

最新推荐文章于 2022-07-29 22:22:05 发布

weixin_39998906

最新推荐文章于 2022-07-29 22:22:05 发布

阅读量428

点赞数

文章标签： bootstrap-table真实交互数据

本文以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起来方便，我懒…）。