线性模型中无序因子变量的对比方式与回归系数的关系:原理,困境及其解决

线性模型中无序因子变量的对比方式与回归系数的关系:原理,困境及其解决

R语言中,进行线性模型分析时,因为模型要求预测变量是数值型,当碰到因子时,它会用一系列与因子水平相对应的数值型的对照变量来代替因子,这些对照变量就是题目中说的对比方式


这里先以二因素的变量为例,以R中的mtcars数据,演示对单因素两水平的变量是如何产生对比方式的。并回答问题:对单个因子变量,对比方式产生的规律是什么

mtcars数据收集了若干种汽车的若干参数,这里想考察汽车的变速器(am, Transmission, 0 = automatic, 1 = manual)对耗油量(mpg, miles per gallon)的影响。

这里把am的变量重新编码一下,方便查看

> mtcarsNew = within(mtcars,{Transmission = NA; Transmission[am == 0] = 'Auto';Transmission[am == 1] = 'Manual'})
> head(mtcarsNew)
                   mpg cyl disp  hp drat    wt  qsec vs am gear carb Transmission
Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4       Manual
Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4       Manual
Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1       Manual
Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1         Auto
Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2         Auto
Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1         Auto
> mtcarsNew$Transmission = factor(mtcarsNew$Transmission)

生成一个新的因子Transmission,标明汽车为手动挡(Manual)还是自动挡(Auto)。

建模:

> Model = lm(data = mtcarsNew, mpg~Transmission)
> summary(Model)

Call:
lm(formula = mpg ~ Transmission, data = mtcarsNew)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.3923 -3.0923 -0.2974  3.2439  9.5077 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)          17.147      1.125  15.247 1.13e-15 ***
TransmissionManual    7.245      1.764   4.106 0.000285 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.902 on 30 degrees of freedom
Multiple R-squared:  0.3598,	Adjusted R-squared:  0.3385 
F-statistic: 16.86 on 1 and 30 DF,  p-value: 0.000285

模型回归系数部分中显示的是当汽车为手动挡(Manual)时的信息,这里和连续变量的结果有区别,是因为模型对自变量生成了一个新变量,设置了其对比方式,可以用contrasts()查看:

> contrasts(mtcarsNew$Transmission)
       Manual
Auto        0
Manual      1

当汽车为自动挡时,该变量为0,当汽车为手动挡时,该变量为1. 之后对两个水平作差,得到两水平间的差异(比如这里的回归系数7.245表示的是手动挡相比于自动挡的耗油量多7.245个单位)。

注意,summary()中的回归系数部分始终是两个水平之间的差值。对任何因子变量,线性模型分析时都会对其产生一种对比方式,对比方式以矩阵的形式存在,比如三水平因子变量的对比方式:

> ThreeLevel = factor(c('A','B','C'))
> contrasts(ThreeLevel)
  B C
A 0 0
B 1 0
C 0 1

四水平的比较方式:

> FourLevel = factor(c('A','B','C','D'))
> contrasts(FourLevel)
  B C D
A 0 0 0
B 1 0 0
C 0 1 0
D 0 0 1

可以看出,因子变量的对比方式是通过如下规律产生的:

  1. 对于某个含有k个水平的因子变量,会产生k-1个对比变量;
  2. 对比变量以一个 k\times(k-1) 的矩阵的形式存在;
  3. 矩阵中每一列表示一个对比变量或对比方式;
  4. 矩阵中的每一行代表某个水平在不同对比方式下的编码;
  5. 默认的对比方式总是将某个水平作为基线(比如上面四水平变量中的),每种对比方式表示用其他的某个水平和基线进行比较。

这一部分要回答的问题是模型是如何利用对比方式来获取回归系数的。以含有四水平的因子变量为例,其对比方式为:

> FourLevel = factor(c('A','B','C','D'))
> contrasts(FourLevel)
  B C D
A 0 0 0
B 1 0 0
C 0 1 0
D 0 0 1

基于这种对比方式,会产生如下模型:

Y = b_{1}X_{B} + b_{2}X_{C}+b_{3}X_{D}+d+\epsilon

最后面两项表示截距和残差。

X_{B} 表示第一列对比方式;X_{C}表示第二列对比方式; X_{D} 表示第三列对比方式。

因此对每个水平模型分别如下:

A水平:Y_{A} = d+\epsilon

B水平:Y_{B} =b_{1}+ d+\epsilon

C水平:Y_{C} =b_{2}+ d+\epsilon

D水平:Y_{A} = b_{3}+d+\epsilon

各水平与基线对比的结果为:

Y_{B}-Y_{A} = b_{1}

Y_{C}-Y_{A} = b_{2}

Y_{D}-Y_{A} = b_{3}

由此获取了不同水平与基线之间的差异。

这种 [0,1] 的编码方式是R中默认的方式,也称为treatment coding。此外无序因子常用的还有sum coding。下面讲解sum coding中对比方式的设置方式及其对回归系数的影响。

通过contr.sum()函数可以获取sum coding的对比:

> contr.sum(2)
  [,1]
1    1
2   -1
> contr.sum(3)
  [,1] [,2]
1    1    0
2    0    1
3   -1   -1
> contr.sum(4)
  [,1] [,2] [,3]
1    1    0    0
2    0    1    0
3    0    0    1
4   -1   -1   -1

可以看出,该方式采用 [1, -1] 的编码方式。在R中修改模型对比方式的方法有两种,一种是通过options()函数设置:

> options(contrasts = c('contr.sum','contr.poly')) 
# 需要同时依次定义无序因子和有序因子的对比方式,后者为定义有序因子的,不常用,这里不展开说了。

另一种为在模型中设置(以这里的mtcars数据为例):

> ComMatrix = contr.sum(2) # 先定义一个对比矩阵,再传入模型
> ComMatrix
  [,1]
1    1
2   -1
> rownames(ComMatrix) = levels(mtcarsNew$Transmission)
> ComMatrix
       [,1]
Auto      1
Manual   -1
> colnames(ComMatrix) = 'Auto'
> ComMatrix
       Auto
Auto      1
Manual   -1
> Model = lm(data = mtcarsNew, mpg~Transmission, contrasts = list(Transmission = ComMatrix))

查看模型信息:

> summary(Model)

Call:
lm(formula = mpg ~ Transmission, data = mtcarsNew, contrasts = list(Transmission = ComMatrix))

Residuals:
    Min      1Q  Median      3Q     Max 
-9.3923 -3.0923 -0.2974  3.2439  9.5077 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       20.7698     0.8822  23.543  < 2e-16 ***
TransmissionAuto  -3.6225     0.8822  -4.106 0.000285 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.902 on 30 degrees of freedom
Multiple R-squared:  0.3598,	Adjusted R-squared:  0.3385 
F-statistic: 16.86 on 1 and 30 DF,  p-value: 0.000285

注意到什么没有,sum coding下的回归系数变为了treatment coding下的一半。要知道原理,还应从数理统计的角度去看:

Transmission在sum coding下的对比方式为:

> options(contrasts = c('contr.sum','contr.poly'))
> contrasts(mtcarsNew$Transmission)
       [,1]
Auto      1
Manual   -1

同样,会产生如下模型:

Y = b_{1}X_{Transmission}+d+\epsilon

对每个水平模型分别如下:

Manual水平: Y = -b_{1}+d+\epsilon

Auto水平: Y = b_{1}+d+\epsilon

两水平相减:

Y_{Auto}-Y_{Manual} = 2b_{1}

b_{1} = \frac{Y_{Auto}-Y_{Manual}}{2}

即summary()中的回归系数为实际回归系数的一半。

注意:这里只是证明了对于两水平的因子变量,对比方式为sum coding时,回归系数是真实值的一半,如果因子有两个以上水平,用sum coding会出错。

以三水平的因子为例,其对比方式如下(这里以C为基线):

> options(contrasts = c('contr.sum','contr.poly'))
> contrasts(ThreeLevel)
  [,1] [,2]
A    1    0
B    0    1
C   -1   -1

会产生如下模型:

Y = b_{1}X_{A} + b_{2}X_{B}+d+\epsilon

各水平对应的模型为:

A水平: Y = b_{1}+d+\epsilon

B水平: Y = b_{2}+d+\epsilon

C水平: Y = -b_{1}- b_{2}+d+\epsilon

各水平对比如下:

Y_{A}-Y_{C} = 2b_{1}+b_{2}

Y_{B}-Y_{C} = b_{1}+2b_{2}

这里的 b_{1}b_{2} 都不是某两个水平之间的真实差异。因此,当因子水平超过2时,采用sum coding的对比方式是可能出问题的


以上介绍了单因素的情况,下面介绍两因素的情况。这里从简单的2*2的设计考虑。比较treatment coding和sum coding下的回归系数。

一旦涉及到多因素,就要考虑是否有交互作用。这里我们分开来看。

  • 两因素不含交互作用

treatment coding对比下的模型为:

Y = b_{1}X_{A}+b_{2}X_{B}+d+\epsilon

X_{A} = 0, A=A_{1};  X_{A} = 1, A=A_{2}

X_{B} = -1, B=B_{1};  X_{B} = 1, B=B_{2}

A的效应为:

Y_{A_{1}} = b_{2}+d+\epsilon

Y_{A_{2}} = b_{1}+b_{2}+d+\epsilon

Y_{A_{2}}-Y_{A_{1}} = b_{1}

同理B的效应为:

Y_{B_{2}}-Y_{B_{1}} = b_{2}

即采用treatment coding编码的对比下,各因子的回归系数等于其真实效应的差异。

sum coding对比下的模型为:

Y = b_{1}X_{A}+b_{2}X_{B}+d+\epsilon

X_{A} = -1, A=A_{1};  X_{A} = 1, A=A_{2}

X_{B} = -1, B=B_{1};  X_{B} = 1, B=B_{2}

A的效应为:

Y_{A_{1}} = -b_{1}+d+\epsilon

Y_{A_{2}} = b_{1}+d+\epsilon

Y_{A_{2}}-Y_{A_{1}} = 2b_{1}b_{1} = \frac{Y_{A_{2}}-Y_{A_{1}} }{2}

同理B的效应为:

Y_{B_{1}} = -b_{2}+d+\epsilon

Y_{B_{2}} = b_{2}+d+\epsilon

Y_{B_{2}}-Y_{B_{1}} = 2b_{2}b_{2} = \frac{Y_{B_{2}}-Y_{B_{1}} }{2}

即采用sum coding编码的对比下,各因子的回归系数等于其真实效应的一半。

  • 两因素含交互作用

treatment coding对比下的模型为:

Y = b_{1}X_{A}+b_{2}X_{B}+b_{3}X_{A}X_{B}+d+\epsilon

X_{A} = 0, A=A_{1};  X_{A} = 1, A=A_{2}

X_{B} = 0, B=B_{1};  X_{B} = 1, B=B_{2}

A的效应为:

Y_{A_{1}} = b_{2}+d+\epsilon

Y_{A_{2}} = b_{1}+b_{2}+b_{3}+d+\epsilon

Y_{A_{2}}-Y_{A_{1}} = b_{1}+b_{3}

同理B的效应为:

Y_{B_{2}}-Y_{B_{1}} = b_{2}+b_{3}

实验中共有四个条件,各个条件下的效应为:

Y_{A_{1}B_{1}} = d+\epsilon

Y_{A_{1}B_{2}} = b_{2}+d+\epsilon

Y_{A_{2}B_{1}} = b_{1}+d+\epsilon

Y_{A_{2}B_{2}} = b_{1}+b_{2}+b_{3}+d+\epsilon

交互作用的效应为:

(Y_{A_{2}B_{2}} -Y_{A_{2}B_{1}} )-(Y_{A_{1}B_{2}} -Y_{A_{1}B_{1}} ) = b_{3}

即treatment coding的对比下,交互作用的回归系数等于真实的效应,主效应的回归系数不是真实的效应

sum coding对比下的模型为:

Y = b_{1}X_{A}+b_{2}X_{B}+b_{3}X_{A}X_{B}+d+\epsilon

X_{A} = -1, A=A_{1};  X_{A} = 1, A=A_{2}

X_{B} = -1, B=B_{1};  X_{B} = 1, B=B_{2}

A的效应为:

Y_{A_{1}} = -b_{1}-b_{2}+d+\epsilon

Y_{A_{2}} = b_{1}-b_{2}+d+\epsilon

Y_{A_{2}}-Y_{A_{1}} = 2b_{1}b_{1} = \frac{Y_{A_{2}}-Y_{A_{1}} }{2}

同理B的效应为:

Y_{B_{2}}-Y_{B_{1}} = 2b_{2}b_{2} = \frac{Y_{B_{2}}-Y_{B_{1}} }{2}

四个条件下的效应为:

Y_{A_{1}B_{1}} = -b_{1}-b_{2}+b_{3}+d+\epsilon

Y_{A_{1}B_{2}} = -b_{1}+b_{2}-b_{3}+d+\epsilon

Y_{A_{2}B_{1}} =b_{1}-b_{2}+b_{3}+d+\epsilon

Y_{A_{2}B_{2}} = b_{1}+b_{2}+b_{3}+d+\epsilon

交互作用的效应为:

(Y_{A_{2}B_{2}} -Y_{A_{2}B_{1}} )-(Y_{A_{1}B_{2}} -Y_{A_{1}B_{1}} ) = 4b_{3}

b_{3}=\frac{(Y_{A_{2}B_{2}} -Y_{A_{2}B_{1}} )-(Y_{A_{1}B_{2}} -Y_{A_{1}B_{1}} )}{4}

sum coding的对比下,交互作用的回归系数等于真实的效应的 \frac{1}{4} ,主效应的回归系数等于真实的效应的 \frac{1}{2}

事实上,将sum coding的对比矩阵减半后(从 [-1,1] 编码改为 [-0.5, 0.5] 编码),交互作用和主效应的回归系数就等于其各自的真实效应了,证明略


以上,我们比较了两种对比方式(treatment coding VS sum coding)在单因素两水平,单因素多水平,两因素无交互,两因素有交互的情况下回归系数与其对应效应的关系,如图 1所示:

图 1. 模型中回归系数与真实效应的关系

可以证明,只要因素包含两个以上的水平,sum coding的对比就可能出现问题。


既然有问题就要解决呀!问题是什么: k 因素 n 水平的情况下( k\geq1 ; n\geq2 ),如果保证回归系数等于真实的效应值呢?其实很简单,采用simple coding (名称取自JAMOVI)的对比就可以了!but……R里面似乎没有这个函数,怎么破?DIY!

source('https://raw.githubusercontent.com/usplos/Eye-movement-related/master/contr.simple.R')

创建后我们执行一下:

> contr.simple(2)
  [,1]
1  0.5
2 -0.5
> contr.simple(3)
           [,1]       [,2]
[1,] -0.3333333 -0.3333333
[2,]  0.6666667 -0.3333333
[3,] -0.3333333  0.6666667
> contr.simple(4)
      [,1]  [,2]  [,3]
[1,] -0.25 -0.25 -0.25
[2,]  0.75 -0.25 -0.25
[3,] -0.25  0.75 -0.25
[4,] -0.25 -0.25  0.75

即每个因素产生 n\times(n-1) 的矩阵,其中每一列代表一个比较,每一行代表一个水平。第一个水平为基线,其编码为 -\frac{1}{n} ,每一列对应基线之外的某个水平 L 与基线的对比, L 所在的行编码为 1-\frac{1}{n} ,其他行编码为 -\frac{1}{n}

以R中的ToothGrowth为例,考察服用剂量(dose,3水平)和营养液类别(supp,2水平)对老鼠牙齿长度(len)的影响:

> Model = lm(data = ToothGrowth, len~supp*dose, contrasts = list(supp=contr.simple(n = 2), dose = contr.simple(3)))
> summary(Model)

Call:
lm(formula = len ~ supp * dose, data = ToothGrowth, contrasts = list(supp = contr.simple(n = 2), 
    dose = contr.simple(3)))

Residuals:
   Min     1Q Median     3Q    Max 
 -8.20  -2.72  -0.27   2.65   8.27 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  18.8133     0.4688  40.130  < 2e-16 ***
supp1         3.7000     0.9376   3.946 0.000231 ***
dose1         9.1300     1.1484   7.951 1.19e-10 ***
dose2        15.4950     1.1484  13.493  < 2e-16 ***
supp1:dose1   0.6800     2.2967   0.296 0.768308    
supp1:dose2  -5.3300     2.2967  -2.321 0.024108 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.631 on 54 degrees of freedom
Multiple R-squared:  0.7937,	Adjusted R-squared:  0.7746 
F-statistic: 41.56 on 5 and 54 DF,  p-value: < 2.2e-16

将其与jamovi中采用simple coding的结果对比,结果是一致的(出现相反数是因为基线不同,并无大碍)

图 2. JAMOVI中模型的回归系数

这里还是不放心,和jamovi中事后检验和简单效应分析的结果比较一下:

图 3. JAMOVI中事后检验部分
图 4. JAMOVI中简单效应分析部分

结果是一致的。OVER!


  • 最后的总结

我们要回答的问题是:

  1. 线性模型中,对因子变量为何要产生对比方式?
  2. 无序因子的对比矩阵内是如何编码的?
  3. 模型是如何利用对比方式来获取回归系数的?
  4. treatment coding和sum coding两种编码方式下对比有何不同?
  5. treatment coding和sum coding两种对比下,模型回归系数和真实效应的关系?
  6. 如何通过设置对比矩阵,保证回归系数等于真实效应?
编辑于 2019-10-27

文章被以下专栏收录