线性(混合)模型中如何设计事先对比(priori contrasts)矩阵

线性(混合)模型中如何设计事先对比(priori contrasts)矩阵

本文主要参考论文:Schad, D. J., Vasishth, S., Hohenstein, S., & Kliegl, R. (2020). How to capitalize on a priori contrasts in linear (mixed) models: A tutorial.Journal of Memory and Language,110, 104038.

线性模型要求预测变量(或称为自变量)为数值型变量。因此当预测变量为因子变量(或分类变量)时,模型会根据其对比方式或对比矩阵将其重新编码为数值变量。部分相关内容,以及R中常见的两种对比方式(treatment codingsum coding)如何影响模型系数结果的问题,已在以下文章中涉及,这里就不再赘述了。

张光耀:线性模型中无序因子变量的对比方式与回归系数的关系:原理,困境及其解决zhuanlan.zhihu.com图标

这里以表格形式简单介绍一下常见的因子 对比/编码方式:

在实际使用中,用户可调用R中自带的对比方式对因子进行编码,但有时用户并不希望对比矩阵中的对比方式是同类型的。比如,在一个单因素三水平(水平简称为 A,B,C)的数据中,用户希望在一个矩阵中同时检验两个假设:1. 水平A和水平B是否有显著差异,2. 水平A和另外两个水平的均值是否有显著差异。这里的假设1可以借助treatment的对比方式,假设2却要用到helmert的方式。此时,用户不能简单地调用R中自带的对比,而是要事先定义因子的对比方式,即事先对比(priori contrasts)。本文的目的是阐述在线性(混合)模型中如何定义事先对比矩阵。

这里借用Schad等人(2020)编写的函数MixedDesign() (available from OSF)来模拟数据:

> demodata = mixedDesign(B = 3,n = 30,M = matrix(c(100,150,70),ncol = 1),SD = 50,long = T)
> head(demodata,5)
#  B_A id        DV
#1  A1  1 138.54957
#2  A1  2  96.39047
#3  A1  3  99.66439
#4  A1  4  45.67003
#5  A1  5 125.94778

这里模拟了一个单因素三水平(A1, A2, A3)设计的数据,每个水平下有三十个观测值。三个水平的均值分别为100,150,和70,标准差为50:

> demodata %>% group_by(B_A) %>% summarise(Mean = mean(DV), SD = sd(DV))
# A tibble: 3 x 3
#  B_A    Mean    SD
#  <fct> <dbl> <dbl>
#1 A1      100  50. 
#2 A2      150  50. 
#3 A3       70  50.0

要考察两个假设:

  1. A1水平和A2水平是否有显著差异;
  2. A1水平和其余两个水平是否有显著差异;

首先写出零假设:

H_{00} : \frac{y_{A1}+y_{A2}+y_{A3}}{3} = 0 (这里是截距项的假设,定义为检验三水平的均值是否与0有显著差异);

H_{01}: y_{A1} = y_{A2}

H_{02}: y_{A1} = \frac{y_{A2}+y_{A3}}{2}

第二步,将三个式子改写为,等号右边为0,等号左边为某个数 k_{i} 与对应因子相乘的形式:

H_{00}: \frac{1}{3}\times y_{A1}+\frac{1}{3}\times y_{A2}+\frac{1}{3}\times y_{A3} = 0

H_{01}: 1\times y_{A1}+(-1)\times y_{A2}+ 0\times y_{A3} = 0

H_{02}: 1\times y_{A1}+ (-\frac{1}{2})\times y_{A2}+(-\frac{1}{2})\times y_{A3} = 0

第三步,将各个因子前面的系数 k_{i} 提出来,组成假设矩阵:

>  Mat = rbind(c(1/3, 1/3, 1/3), c(1, -1, 0), c(1, -1/2, -1/2)) # 定义假设矩阵
> rownames(Mat) = paste0('H',0:2) # 行名为假设的名称
> colnames(Mat) = paste0('A',1:3) # 列名为水平的名称
> fractions(Mat) # 以分数形式查看假设矩阵
#    A1   A2   A3  
# H0  1/3  1/3  1/3
# H1    1   -1    0
# H2    1 -1/2 -1/2

注意三点:

  1. 假设矩阵的每一行代表一个假设,每一列表示一个水平,矩阵中第 j 行第 k 列的数值表示第 k 列所表示的那个水平在第 j 行所表示的那个假设中的权重;
  2. 假设矩阵的行列可以颠倒,因为其不进入最终的建模;但因对比矩阵是假设矩阵的逆矩阵,且对比矩阵中必须按照“每一列表示一个对比,每一行表示一个水平”,因此在构建假设矩阵时,建议用户按照这里的排版(即行表示假设,列表示水平)方式进行排列。
  3. 假设矩阵中的每一列表示的水平,按照从左到右(这里从左到右依次为A1A2A3),要与因子的水平顺序对应(如不对应,要先调整):
> levels(demodata$B_A)
# [1] "A1" "A2" "A3"

第四步,对假设矩阵求逆矩阵,其逆矩阵即为对比矩阵:

> Con = ginv(Mat) # 求逆
> rownames(Con) = paste0('A',1:3)
> colnames(Con) = paste0('H',0:2)
> fractions(Con)
#      H0   H1   H2  
# A1    1    0  2/3
# A2    1   -1  2/3
# A3    1    1 -4/3

第五步,将Con的第2列到最后一列(第一列为截距项),定义为该因子的对比方式:

> contrasts(demodata$B_A) = Con[,2:3]
> fractions(contrasts(demodata$B_A))
#     H1   H2  
#A1    0  2/3
#A2   -1  2/3
#A3    1 -4/3

建模:

> Model = lm(data = demodata, DV~B_A)
> round(coef(summary(Model)),3) # 提取模型参数部分,并以保留3位小数的格式呈现
#            Estimate Std. Error t value Pr(>|t|)
#(Intercept)  106.667       5.27  20.239    0.000
#B_AH1        -50.000      12.91  -3.873    0.000
#B_AH2        -10.000      11.18  -0.894    0.374

下面简单验证三个参数。

三个水平的均值为: \mu_{A1} = 100;\mu_{A2} = 150;\mu_{A1} = 70

第一个为截距,考察的是三个水平的均值:

> mean(c(100,150,70))
# [1] 106.6667

第二个为A1A2的差异,为-50.

第三个为A1和其余两水平均值的差异:

> 100 - mean(c(150,70))
# [1] -10

以上是定义事先对比方式的一个案例,总结一下步骤:

  1. 写出零假设(第一个是截距项的零假设,后面依次为要考察的问题的零假设;对于含有 n 个水平的因子,可写 n 个零假设);
  2. 将零假设改写为“等号右边为0,等号左边为某个数 k_{i} 乘以对应因子”的格式;
  3. 把改写后的零假设中的 k_{i} 提出来组成假设矩阵,其中每一行表示一个假设,每一列表示一个水平;按照从左到右的顺序,每一列表示的水平要和该因子的水平保持一致;
  4. 对假设矩阵求逆,得到对比矩阵,对比矩阵的行名称与假设矩阵的列名称一致,对比矩阵的列名称与假设矩阵的行名称一致;
  5. 将对比矩阵的第二列到最后一列定义为该因子的对比方式。

多数因子设计的实验往往包含两个或两个以上因子。此时,如果用户想考察的问题不能通过直接调用R中的对比方式实现,建议用户将这些因子合并为一个因子(比如一个 2\ast3 的实验设计,可合并为一个单因素6水平的设计),然后采用事先定义对比矩阵的方法考察。

更多案例和相关内容请见本文开篇索引的论文。

编辑于 01-19

文章被以下专栏收录