线性(混合)模型中如何对因子变量事先生成虚拟变量

线性(混合)模型中如何对因子变量事先生成虚拟变量

线性模型要求预测变量(或称为自变量)为数值型变量。因此当预测变量为因子变量(或分类变量)时,模型会根据其对比方式或对比矩阵将其重新编码为数值变量。在我之前的两篇文章中已经介绍了部分与之有关的内容:

张光耀:线性模型中无序因子变量的对比方式与回归系数的关系:原理,困境及其解决zhuanlan.zhihu.com图标张光耀:线性(混合)模型中如何设计事先对比(priori contrasts)矩阵zhuanlan.zhihu.com图标

这里以treatment codingsum coding为例,简单回顾一下:

对于某个三水平的因素FactorA,包含三个水平A1A2A3(其中A1为基线)。其treatment 编码方式为:

> FactorA = factor(c('A1','A2','A3'),levels = c('A1','A2','A3'))
> contrasts(FactorA) = contr.treatment(3)
> contrasts(FactorA)
#    2 3
# A1 0 0
# A2 1 0
# A3 0 1 

这里的这种对比方式对初学者来说较难理解,我们将这两列的列名改变一下:

> colnames(contrasts(FactorA)) = c('Variable.1','Variable.2')
> contrasts(FactorA)
#      Variable.1 Variable.2
# A1          0          0
# A2          1          0
# A3          0          1

这里,我们用Variable.1Variable.2 两个变量对 FactorA 重新编码:

当FactorA=A2时,Variable.1 = 1;

当FactorA\neq A2时,Variable.1 = 0;

当FactorA=A3时,Variable.2= 1;

当FactorA\neq A3时,Variable.2 = 0.

接下来看一下sum的编码方式:

> contrasts(FactorA) = contr.sum(3)
> colnames(contrasts(FactorA)) = c('Variable.1','Variable.2')
> contrasts(FactorA)
#      Variable.1 Variable.2
# A1          1          0
# A2          0          1
# A3         -1         -1

同样,我们用Variable.1Variable.2 两个变量对 FactorA 重新编码:

当Factor=A1时,Variable.1 = 1,Variable.2 = 0;

当Factor=A2时,Variable.1 = 0,Variable.2 = 1;

当Factor=A3时,Variable.1 = -1,Variable.2 = -1;

可以看出:某个因子的对比方式实质上是用若干新的数值型变量代替该因子;对于其中某个数值型变量(本文中我们称之为“虚拟变量”),它的数值与该因子的各水平有一一对应的关系;对于某个含有 k 个水平的因子,需要生成 k-1 个虚拟变量(和因子的自由度有关);我们常说的线性模型中的某个因子的系数(或者称为斜率),实质上是虚拟变量的系数

既然如此,建模过程中,只要定义好了因子的对比方式,直接把因子作为自变量建立的模型,与先生成与之对应的虚拟变量再把虚拟变量作为自变量建立的模型,是一致的。这里我们用一个简单的例子演示一下:

仍然借助 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  55.979828
#2  A1  2  36.333803
#3  A1  3 168.098854
#4  A1  4  -9.557433
#5  A1  5 116.715957

这里模拟了一个单因素三水平(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

这里我们把A1作为基线,分别考察A2A3A1的差异(即用到treatment的编码方式):

> levels(demodata$B_A) = paste0('A',1:3)
> contrasts(demodata$B_A) = contr.treatment(3)
> Model1 = lm(data = demodata, DV~B_A)
> round(coef(summary(Model1)),digits = 3)
            Estimate Std. Error t value Pr(>|t|)
(Intercept)      100      9.129  10.954    0.000
B_A2              50     12.910   3.873    0.000
B_A3             -30     12.910  -2.324    0.022

可以看到,第一个回归系数的斜率为50,为A2A1的均值差;第二回归系数为-30,为A3A1的均值差。

下面,首先生成对应的虚拟变量,这里借助R中的model.matrix()函数:

> Dummy = model.matrix(~B_A,data = demodata)
> head(Dummy,5)
#  (Intercept) B_A2 B_A3
#1           1    0    0
#2           1    0    0
#3           1    0    0
#4           1    0    0
#5           1    0    0

这里生成了三列变量,第一列为截距,第二列开始依次为根据对比矩阵定义的新变量,定义的规则同上文所述的treatment对比矩阵的构建规则。接下来,把第二列开始的所有列添加到demodata中。

> demodata[c('Variable.1','Variable.2')] = Dummy[,2:3]
> head(demodata,5)
#  B_A id         DV Variable.1 Variable.2
#1  A1  1  55.979828          0          0
#2  A1  2  36.333803          0          0
#3  A1  3 168.098854          0          0
#4  A1  4  -9.557433          0          0
#5  A1  5 116.715957          0          0

接下来检验一下B_A中各水平和新添加的虚拟变量是否一一对应:

> demodata %>% group_by(B_A) %>% summarise(Var1 = paste0(unique(Variable.1)),Var2 = paste0(unique(Variable.2)))
# A tibble: 3 x 3
#  B_A   Var1  Var2 
#  <fct> <chr> <chr>
#1 A1    0     0    
#2 A2    1     0    
#3 A3    0     1 

是一一对应的。接下来用虚拟变量来建模:

> Model2 = lm(data = demodata, DV~Variable.1+Variable.2)
> round(coef(summary(Model2)),digits = 3)
            Estimate Std. Error t value Pr(>|t|)
(Intercept)      100      9.129  10.954    0.000
Variable.1        50     12.910   3.873    0.000
Variable.2       -30     12.910  -2.324    0.022

可以看到结果是一致的。

总结一下:本文主要介绍了如何根据因子的对比方式生成对应的虚拟变量,并根据虚拟变量来建模。这种建模方式的好处之一是,极大地有助于线性混合模型的优化,这部分内容在下篇文章中讲解。

参考文献:

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.

编辑于 01-20

文章被以下专栏收录