4  空白人参上!(施工中)

在我小学的时候,在同伴中曾经流行过一个游戏。这个游戏中文名有一种跟汉语若即若离的拉扯感,叫《小朋友齐打交2》,英文名是 “Little Fighter 2”。刚刚在网络上搜索了一下,这个游戏的竟然还有留有一个疑似官网(https://www.lf2.net/),非常推荐没玩过的读者去试试(不过记得回来看书)。

提到这个格斗游戏的原因是因为,第一,它真的很好玩;第二,它在里面有很多角色,其中有一个角色,我管他(们)叫空白人。下图是我从游戏资源文件夹里直接找到的素材(如有侵权,请联系我删除)。这个空白人跟其他角色的不同点在于,他(们)除了最基础的拳脚之外,搓不出任何其他技能。在色彩缤纷的角色中,这个线条勾勒的角色也显得单调而格格不入。

不过,虽然这么说有点残酷,但在研究设计中,这样一个空白人却是天生的对照组圣体。想想吧,我们现在有机人和袋人,我们想比较他们的魔法天赋,那么我们就需要这么一组人,他们没有自我性别认同,但他们又是人,有基础的其他一切功能和能力。这样,我们就可以单独去看机人和袋人这两种性别自我认同分别对人的魔法天赋有什么影响了。

好,那我宣布我们要在机人和袋人的基础上增加一组空白人进入数据。请看代码:

## 生成空白人的对应数据
set.seed(234)
df.blank <- data.frame(ID=c(101:150),iden=c(rep("空白人",50)))
df.blank$magic_score <- rnorm(n = 50, mean = 95, sd = 10)

## 合并空白人数据至我们之前的机人+袋人的数据
df2.1 <- df2 <- rbind(df1,df.blank)

## 看看加进去了没有
tail(df2)
     ID   iden magic_score
145 145 空白人    90.64078
146 146 空白人   105.51912
147 147 空白人    84.72741
148 148 空白人    95.39913
149 149 空白人    91.72250
150 150 空白人    99.92976
summary(df2)
       ID                   iden     magic_score    
 Min.   :  1.00   某品牌塑料袋:50   Min.   : 64.64  
 1st Qu.: 38.25   武装直升机  :50   1st Qu.: 80.95  
 Median : 75.50   空白人      :50   Median : 90.14  
 Mean   : 75.50                     Mean   : 90.30  
 3rd Qu.:112.75                     3rd Qu.: 99.52  
 Max.   :150.00                     Max.   :115.22  

我们再来看看描述统计。

某品牌塑料袋武装直升机空白人7090110130
iden某品牌塑料袋武装直升机空白人三种自我认同人士的魔法天赋数据自我认同魔法天赋分数

好的,数据我们生成完了,接下来就该分析了。我们的研究假设是什么呢?回顾我们在前两章做的事情,显然一个假设可以是:机人、袋人和空白人三组人的魔法天赋存在差异。即

μμμ \mu_{机人} \ne \mu_{袋人} \ne \mu_{空白人}

同样的,我们可以找到我们的零假设,机人、袋人和空白人三组人的魔法天赋都一样,把上面的不等号改成等号就行了。那么,用什么统计分析进行验证呢?

4.1 三水平因素的方差分析实战

学过方差分析的读者大概根本不用犹豫,这不就是最典型的单因素方差分析吗?是的,我们在上一章中使用了方差分析检验机人和袋人之间是否存在魔法天赋的差异。但是实际的研究和教学中,其实大多数使用的都是三水平以上的自变量作为单因素方差分析的例子,这也是它某种程度上优于t检验的地方——它可以检验多于两组以上的组间差异是否显著。

所以我们现在先做一下方差分析。

## 不要在意我的模型和变量命名有什么规律,命名这种事情看心情的
my.aov2.1 <- aov(magic_score~iden, df2.1)
## 我们之后的方差分析会尽量使用car包里的Anova命令报告结果,原因之后会介绍
car::Anova(my.aov2.1)
Anova Table (Type II tests)

Response: magic_score
           Sum Sq  Df F value    Pr(>F)    
iden       6180.1   2  35.415 2.794e-13 ***
Residuals 12826.3 147                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

结果很是喜人,方差分析发现了一个显著的差异,即机人、袋人和空白人在魔法天赋上存在显著差异。

不过,显然本书的主要目的不是带读者们重走方差分析路的。我们在上一章说到,本质上方差分析和回归分析都从属于一般线性模型。那么我们能不能用线性回归模型的办法检验我们的假设呢?

4.2 编码:把类别变量变为连续变量的魔法

使用线性回归分析三组人魔法分数的差异,显然预测变量是一个由这三组人组成的分类变量。而进行分析最大的障碍显然是,如何把这样一个分类变量作为预测变量放进回归模型。其实这个问题我们在上一章的时候就有所涉及,还专门在 与大家一起思考和讨论过。只不过,在上一章中,我们的自我性别认同还只是一个两水平的变量。但现在,这个分类变量因为空白人的插足,变成了三水平的分类变量。我们该怎么办呢?

我们在上一章中的做法,是给机人和袋人分别赋予两个数值,机人=0,袋人=1。最直观的思路,就是顺着这个逻辑再加一个数,比如,空白人=2。这样,行不行呢?我们试试:

df2.1$iden_n <- c(rep(0,50), rep(1,50), rep(2,50))
my.lm2.1 <- lm(magic_score~iden_n, df2.1)
summary(my.lm2.1)

Call:
lm(formula = magic_score ~ iden_n, data = df2.1)

Residuals:
    Min      1Q  Median      3Q     Max 
-32.228  -7.273  -0.312   6.211  24.918 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  83.7360     1.2865  65.091  < 2e-16 ***
iden_n        6.5653     0.9965   6.589  7.3e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.965 on 148 degrees of freedom
Multiple R-squared:  0.2268,    Adjusted R-squared:  0.2216 
F-statistic: 43.41 on 1 and 148 DF,  p-value: 7.296e-10

显然,不行。上一章中,我们比较过回归模型的模型拟合和方差分析的结果应该一致。但显然上面的回归模型的拟合结果的显著性检验 F-statistic: 43.41 on 1 and 148 DF, p-value: 7.296e-10,与我们在 中得到的结果不一致。那么,问题出在哪儿呢?

在上一章中,当我们在使用这种编码方式的时候,其实我们强调过,两水平变量可以被视为一种特殊的“等距“变量处理,所以它能被赋任意两个不相等的值放入模型。那么,三水平变量改变了什么呢?回到我们的例子,当我们只有机人和袋人两组人时,不管我们怎么赋值,二者都只有一个“距离”,但是当空白人加入这个家庭之后,三者就出现了两两之间的3个“距离”。这三个距离相等吗?显然不相等,或者说,我们根本不知道他们的“距离”是什么,他们在这个时候根本无法被看作是等距变量!所以,机人=0,袋人=1,空白人=2,或者任意什么其他同时给他们仨赋值的其他方式,似乎都在假设他们之间的距离情况,而这种假设是不成立的。

我们在 中分析过,不管如何编码两水平分类变量,最终回归模型预测的结果变量都会是机人和袋人的平均魔法天赋分数。那么上面这个模型呢?

table(predict(my.lm2.1))

83.7360375414107 90.3013813036805 96.8667250659503 
              50               50               50 

跟我们实际数据的平均数比较一下。

group_means <- aggregate(magic_score ~ iden, data = df2.1, FUN = mean)
group_means
          iden magic_score
1 某品牌塑料袋    95.29439
2   武装直升机    81.23953
3       空白人    94.37022

不能说一模一样吧,只能说是毫无关系了。把三组人作为等距变量放进模型,看来导致了模型没能

显然,我们需要一种新的编码方式。这种编码方式,需要给2水平以上的分类变量施加魔法,让它们能像两水平分类变量一样,作为数值型的变量进入回归模型。

4.2.1 哑变量编码 (dummy coding)

从上面的分析我们可以知道,要让回归模型正确地运行,我们首先需要让它的预测值变成每组的。我们不妨再回到两水平变量的情况,重新审视一下它预测的过程。

summary(lm1)

Call:
lm(formula = magic_score ~ iden_n, data = df1.1)

Residuals:
    Min      1Q  Median      3Q     Max 
-18.414  -6.369  -1.198   6.157  25.225 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   81.240      1.301  62.426  < 2e-16 ***
iden_n        14.055      1.840   7.637  1.5e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.202 on 98 degrees of freedom
Multiple R-squared:  0.3731,    Adjusted R-squared:  0.3667 
F-statistic: 58.32 on 1 and 98 DF,  p-value: 1.498e-11

我们的回归模型是 ŷ=81.240+14.055x\hat y=81.240+14.055x 。当x=0时,我们的预测值等于b0就是机人的魔法天赋分数均值。而当x=1时, 我们的预测值等于 b0+b1 (81.240+14.055=95.295)就是袋人的魔法天赋平均分了。

在这个模型的基础上我们缺什么呢?不就是空白人的平均魔法分数吗,我们想想办法怎么能补上呢?

显然,上面用的分别赋值0,1,2的办法不行,因为变量并不等距。分类变量只有在两个水平的时候可以视作等距。那简单啊,我们找到一个变量,让它用空白人跟袋人之间的差值作为它对应的回归系数(b2),也就是 b0+b2 =89.23953,问题不就解决了吗?

怎么做到这一点呢?照猫画虎,看看机人和袋人是怎么搞的?我们之前设置了一个变量,叫它D1 好了,机人=0,所以变成了截距。袋人=1,就是b1了。那简单,就让另一个变量,叫它D2,赋值机人=0,空白人等于1,不就解决问题了吗?

就像下面这张表所列的赋值关系。如果我们分别把D1 和 D2 放入两个回归方程,那其实就等价于我们做了两个独立样本t检验。可问题是我们要检验的不是这两个两组间的差异,我们想检验的是三组间的魔法天赋有没有差异呀。

自我性别认同 D1(袋人-机人) D2(空白人-机人)
机人 0 0
袋人 1 -
空白人 - 1

就像上面分析的,我们需要同时在一个方程中的结果得到,b0+b1 和b0+b2 分别对应袋人和机人的平均魔法天赋数值。所以关键问题是,把D1和D2同时放进方程,D1中的空白人和D2中的袋人应该怎么赋值呢?

我们上面说,我们把两水平的分类变量视作等距变量。现在,我们有了D1 和 D2, 我们如果把它们视作两个两水平的变量,分别代表两组人,以此来达到“等距”的效果,显然这时候的赋值就只能在0和1中选择了。那么,给D1中的空白人和D2中的袋人是选择1还是0呢?

我们不妨回到方程:

ŷ=b0+b1D1+b2D2 \hat y= b_0+b_1D_1+b_2D_2

那么,对于机人的预测值,也就是机人的平均数,根据上面的表,D1=D2=0,代入方程中,显然是:

ŷ=b0(4.1) \hat y_{机人} = b_0 \qquad(4.1)

对于袋人呢,我们不知道D2中如何给它赋值,就先用x表示,那么:

ŷ=b0+b1+b2x \hat y_{袋人} = b_0+b_1+b_2x

想重复出之前两水平变量的效果,显然让上式中的 x=0 就可以了,即:

ŷ=b0+b1(4.2) \hat y_{袋人} = b_0+b_1 \qquad(4.2)

同理,对于空白人来说,在D1中对它赋值为0,就可以得到:

ŷ=b0+b2(4.3) \hat y_{空白人} = b_0+b_2 \qquad(4.3)

看上去我们的目标达成了!新的编码表如下。

自我性别认同 D1(袋人-机人) D2(空白人-机人)
机人 0 0
袋人 1 0
空白人 0 1

赶紧来试试能不能成功。

## 按照编码表生成D1和D2两个变量

df2.1$D1 <- ifelse(df2.1$iden == "某品牌塑料袋", 1, 0)
df2.1$D2 <- ifelse(df2.1$iden == "空白人", 1, 0)

## 加入回归模型分析

my.lm2.2 <- lm(magic_score~D1+D2,df2.1)
summary(my.lm2.2)

Call:
lm(formula = magic_score ~ D1 + D2, data = df2.1)

Residuals:
     Min       1Q   Median       3Q      Max 
-29.7311  -6.3261   0.3617   5.8677  25.2248 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   81.240      1.321  61.498  < 2e-16 ***
D1            14.055      1.868   7.523 4.91e-12 ***
D2            13.131      1.868   7.029 7.26e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.341 on 147 degrees of freedom
Multiple R-squared:  0.3252,    Adjusted R-squared:  0.316 
F-statistic: 35.41 on 2 and 147 DF,  p-value: 2.794e-13

对比一下我们之前方差分析的结果。

car::Anova(my.aov2.1)
Anova Table (Type II tests)

Response: magic_score
           Sum Sq  Df F value    Pr(>F)    
iden       6180.1   2  35.415 2.794e-13 ***
Residuals 12826.3 147                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

巧了不是,这不是对上了?这回两个F检验一模一样了。

我们还可以注意一下模型中的回归系数。对比一下我们机人袋人和空白人的魔法天赋分数。

group_means
          iden magic_score
1 某品牌塑料袋    95.29439
2   武装直升机    81.23953
3       空白人    94.37022

截距的估计值显然是机人的魔法天赋分数。b1和b2,也就是D1和D2对应的回归系数呢?其实整理, 。不难推导出:

$$ b_0=\hat y_{机人}\\ b_1=\hat y_{袋人}-\hat y_{机人}\\ b_2=\hat y_{空白人}-\hat y_{机人} $$

看看我们的回归方程结果输出,数值上没错。那么,它们后面对应的,对于系数显著性的检验,检验的是什么呢?谜底藏在谜面上,就是对回归系数数值与0的比较。那么,如果b1不等于0,意味着什么呢?袋人的魔法天赋平均数不等于机人的。同理,b2的显著性检验,检验的就是空白人和机人之间的魔法天赋存不存在显著差异。

注意这个结果跟我们在第二章中, 的结果,尽管数值看上去很接近,但并不一样。

t.test(magic_score~iden, df1.1, var.equal = TRUE)

    Two Sample t-test

data:  magic_score by iden
t = 7.6368, df = 98, p-value = 1.498e-11
alternative hypothesis: true difference in means between group 某品牌塑料袋 and group 武装直升机 is not equal to 0
95 percent confidence interval:
 10.40264 17.70709
sample estimates:
mean in group 某品牌塑料袋   mean in group 武装直升机 
                  95.29439                   81.23953 

我们上面提到的使用b1的显著性对机人和袋人的比较,是考虑了空白人的数据的前提下,比较机人和袋人的魔法天赋的差异。如果有同学在方差分析中学过对应的知识,就会知道,这实际上就是一种事前比较(planned contrast)检验。

引申:事前比较是对单因素方差分析的结果的进一步检验吗?

【待施工】

那么,事后检验是对单因素方差分析结果的进一步检验吗?

LSD SE=MS(1ni+1nj)SE = \sqrt{MS_{组内}(\frac{1}{n_i} + \frac{1}{n_j})}

DescTools::PostHocTest(my.aov2.1,method="lsd")

  Posthoc multiple comparisons of means : Fisher LSD 
    95% family-wise confidence level

$iden
                               diff     lwr.ci     upr.ci    pval    
武装直升机-某品牌塑料袋 -14.0548617 -17.746848 -10.362875 4.9e-12 ***
空白人-某品牌塑料袋      -0.9241741  -4.616161   2.767812  0.6216    
空白人-武装直升机        13.1306875   9.438701  16.822674 7.3e-11 ***

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Bonferoni

DescTools::PostHocTest(my.aov2.1,method="bonferroni")

  Posthoc multiple comparisons of means : Bonferroni 
    95% family-wise confidence level

$iden
                               diff     lwr.ci    upr.ci    pval    
武装直升机-某品牌塑料袋 -14.0548617 -18.579040 -9.530683 1.5e-11 ***
空白人-某品牌塑料袋      -0.9241741  -5.448353  3.600005  1.0000    
空白人-武装直升机        13.1306875   8.606509 17.654866 2.2e-10 ***

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DescTools::PostHocTest(my.aov2.1,method="lsd")$iden[,"pval"]*3
武装直升机-某品牌塑料袋     空白人-某品牌塑料袋       空白人-武装直升机 
           1.472920e-11            1.864675e+00            2.177378e-10 
DescTools::PostHocTest(my.aov2.1,method="bonferroni")$iden[,"pval"]
武装直升机-某品牌塑料袋     空白人-某品牌塑料袋       空白人-武装直升机 
           1.472920e-11            1.000000e+00            2.177378e-10 

是的,我们通过上面的分析,其实一起推导出了如何对多水平分类变量进行哑变量编码(dummy coding),以及如何解读每个虚拟变量的回归系数显著性检验。总结来说,对于随便一个多水平分类变量,假设它有k个水平,那么需要k-1个变量对它进行编码。而哑变量编码的方法,就是在k个水平中选择一个水平,所有虚拟变量对这个水平的赋值均为0,剩下的其他水平,每个虚拟变量对其中的一个水平赋值为1且不重复。用文字表达还挺难理解的,我们可以看表:

分类变量 D1 D2 Dn
A 0 0 0
B 1 0 0
C 0 1 0
N 0 0 1

其实不难发现,我们在进行哑变量编码时,那个所有虚拟变量都赋值为0的水平,就是那个被我们安排到截距的倒霉蛋。而剩下的每个虚拟变量,都是在检验自己赋值为1的那个水平,跟截距对应的那个水平的平均数差异是否显著。换句话说,截距对应的这个变量,通常来说我们会放的就是变量中适合作为对照组的水平,我们也把这个水平叫做基线水平(baseline level)。

不过视线回到这章的一开头,我们加入了空白人这个先天对照组圣体,显然不是为了让它占据D2。所以接下来,我们把它放到截距的位置发光发热,看看结果如何。

首先,我们再重新画一个编码表。

自我性别认同 D1(袋人-空白人) D2(机人-空白人)
机人 0 1
袋人 1 0
空白人 0 0

然后生成新的虚拟变量,并进行回归分析。

## 按照编码表生成D1n和D2n两个变量

df2.1$D1n <- ifelse(df2.1$iden == "某品牌塑料袋", 1, 0)
df2.1$D2n <- ifelse(df2.1$iden == "武装直升机", 1, 0)

## 加入回归模型分析

my.lm2.3 <- lm(magic_score~D1n+D2n,df2.1)
summary(my.lm2.3)

Call:
lm(formula = magic_score ~ D1n + D2n, data = df2.1)

Residuals:
     Min       1Q   Median       3Q      Max 
-29.7311  -6.3261   0.3617   5.8677  25.2248 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  94.3702     1.3210  71.438  < 2e-16 ***
D1n           0.9242     1.8682   0.495    0.622    
D2n         -13.1307     1.8682  -7.029 7.26e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.341 on 147 degrees of freedom
Multiple R-squared:  0.3252,    Adjusted R-squared:  0.316 
F-statistic: 35.41 on 2 and 147 DF,  p-value: 2.794e-13

结果不出意料,模型的拟合情况一致。D2和D2n两个虚拟变量的结果,除了系数本身的符号因为对照组不同而相反之外,绝对值和显著性检验的结果也一致。D1n的估计值也与袋人和空白人的魔法天赋平均数之差一致。而根据这一结果,考虑到空白人的种族特性,说明自我认同为机人会降低人的魔法天赋,而自我认同为袋人似乎并不会相比基线水平增加人的魔法天赋。虽然是我自己生成的数据,但这是一个很有意思的结果,不是吗?

简而言之,使用哑变量编码进行回归分析,可以在检验单因素方差分析检验的结果的同时,使用回归系数进行变量间水平的两两事前比较。

思考:如果不把D1中的空白人或D2中的袋人赋值为0会发生什么?

【待施工】

自我性别认同 D1 D2
机人 0 0
袋人 1 0
空白人 1 1
df2.1$R1 <- c(rep(0,50),rep(1,50),rep(1,50))
df2.1$R2 <- c(rep(0,50),rep(0,50),rep(1,50))
my.lm2.4 <- lm(magic_score~D1+D2,df2.1)
summary(my.lm2.4)

Call:
lm(formula = magic_score ~ D1 + D2, data = df2.1)

Residuals:
     Min       1Q   Median       3Q      Max 
-29.7311  -6.3261   0.3617   5.8677  25.2248 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   81.240      1.321  61.498  < 2e-16 ***
D1            14.055      1.868   7.523 4.91e-12 ***
D2            13.131      1.868   7.029 7.26e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.341 on 147 degrees of freedom
Multiple R-squared:  0.3252,    Adjusted R-squared:  0.316 
F-statistic: 35.41 on 2 and 147 DF,  p-value: 2.794e-13

4.2.2 效应变量编码 (effect coding)

虽然我们在本章的一开头在我们的魔法世界宇宙中设定了空白人是一种先天对照组圣体。然而,世事不会尽如人意。正如我们会讨论机器人或人工智能在觉醒了意识后,我们应该如何处理他们的个体权利问题一样;魔法世界中的空白人有一天突然觉醒了,他们向世界宣告,他们不认为自己的空白自我性别认同是”性空“。相反,他们认为他们的空白自我性别认同是实在的,具体的,是一种“有”。所谓“菩提本无树,明镜亦非台”,“无”本身就是一种“有”,而不是“空”。

好吧,用大白话来说就是,他们觉得空白不是没有自我性别认同,而是自我性别认同的一种。这些卑鄙的空白人,他们破坏了我们研究的对照组!

虽然以上是开玩笑,但这个例子也想揭示,很多时候,特别是当我们使用的是现实生活中存在却极难修改的变量时,很难找到一个合适的且有足够数量的对照组。那么有没有什么办法,能让我们检验空白人或其他两组人相比于这三种自我性别认同人一般魔法天赋的情况呢?

相信有读者已经想到了,所谓的一般情况,不就是平均情况吗?那比较每组人的平均魔法天赋和三组人的总平均魔法天赋之间的差异不就行了?其实认真想想,在方差分析中,使用的组间变异(各组平均数减总平均数的平方和),不就是这个思路吗?

我们在哑变量编码中,实际上是把自我性别认同对魔法天赋的效应拆成了这样的形式:

b0:空白人平均魔法天赋

b1:空白人平均魔法天赋+(袋人-空白人的魔法天赋)

b2:空白人平均魔法天赋+(机人-空白人的魔法天赋)

那么我们能不能转变思路,把基线水平换成三组人的总平均数,即:

b0:总平均魔法天赋

b1:总平均魔法天赋+(袋人-总平均魔法天赋)

b2:总平均魔法天赋+(机人-总平均魔法天赋)

这样,我们不就可以通过检验回归系数来看袋人和机人相对于三族人的平均水平而言,他们二者的魔法天赋情况了吗? 回顾我们的哑变量编码,是这样的:

$$ b_0=\hat y_{机人}\\ b_1=\hat y_{袋人}-\hat y_{机人}\\ b_2=\hat y_{空白人}-\hat y_{机人} $$

那么我们的新编码情况,就是:

$$ b_0=(\hat y_{机人}+\hat y_{袋人}+\hat y_{空白人})/3\\ b_1=\hat y_{袋人}-(\hat y_{机人}+\hat y_{袋人}+\hat y_{空白人})/3\\ b_2=\hat y_{机人}-(\hat y_{机人}+\hat y_{袋人}+\hat y_{空白人})/3 $$

现在的问题是,我们怎么找到对应的编码表,生成虚拟变量呢?很简单,在哑变量编码中,我们是已知编码,来看我们的斜率的实际意义。把组平均数看成了已知常数,求解三个回归系数的解。现在我们反过来,把三个回归系数当作常数,求解三个组平均数的解即可。由上面的式子,不难得到

$$ \hat y_{袋人}=b_0+b_1\\ \hat y_{机人}=b_0+b_2\\ \hat y_{空白人}=b_0-b_1-b_2 $$

那么,编码表应该长什么样呢?根据上面的式子,假设我们生成的两个虚拟变量叫E1和E2,那么E1对应b1的赋值,E2对应b2的赋值。这样的话,我们就能用这两个虚拟编码在回归模型中预测出对应的平均数。即,如下的编码表。大家也可以用下面的编码表代入回去重新检验一遍是不是能推导出上面的回归系数的式子。

自我性别认同 E1 E2
机人 0 1
袋人 1 0
空白人 -1 -1

上面的编码方式就是我们的效应变量编码(effect coding)。简单总结,这种编码方式同样需要找到一个基线水平,在哑变量编码中,我们把这个水平在所有虚拟变量上都赋值为0,在效应变量编码中,我们把这个赋值改为-1就可以了。

我们再实战检验一下。

## 按照编码表生成E1和E2两个变量

df2.1$E1 <- ifelse(df2.1$iden == "某品牌塑料袋", 1, 
                   ifelse(df2.1$iden == "空白人",-1,0))
df2.1$E2 <- ifelse(df2.1$iden == "武装直升机", 1, 
                   ifelse(df2.1$iden == "空白人",-1,0))

## 加入回归模型分析

my.lm2.3 <- lm(magic_score~E1+E2,df2.1)
summary(my.lm2.3)

Call:
lm(formula = magic_score ~ E1 + E2, data = df2.1)

Residuals:
     Min       1Q   Median       3Q      Max 
-29.7311  -6.3261   0.3617   5.8677  25.2248 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  90.3014     0.7627 118.399  < 2e-16 ***
E1            4.9930     1.0786   4.629 8.00e-06 ***
E2           -9.0618     1.0786  -8.401 3.45e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.341 on 147 degrees of freedom
Multiple R-squared:  0.3252,    Adjusted R-squared:  0.316 
F-statistic: 35.41 on 2 and 147 DF,  p-value: 2.794e-13

4.2.3 对比变量编码 (contrast coding)

orthogonal: 交互作用

my.lm2.3 <- lm(magic_score~D1n*D2n,df2.1)
summary(my.lm2.3)

Call:
lm(formula = magic_score ~ D1n * D2n, data = df2.1)

Residuals:
     Min       1Q   Median       3Q      Max 
-29.7311  -6.3261   0.3617   5.8677  25.2248 

Coefficients: (1 not defined because of singularities)
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  94.3702     1.3210  71.438  < 2e-16 ***
D1n           0.9242     1.8682   0.495    0.622    
D2n         -13.1307     1.8682  -7.029 7.26e-11 ***
D1n:D2n           NA         NA      NA       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.341 on 147 degrees of freedom
Multiple R-squared:  0.3252,    Adjusted R-squared:  0.316 
F-statistic: 35.41 on 2 and 147 DF,  p-value: 2.794e-13
df2.1$R1 <- c(rep(1,50),rep(20,50),rep(300,50))
df2.1$R2 <- c(rep(233,50),rep(666,50),rep(4396,50))
my.lm2.4 <- lm(magic_score~D1+D2,df2.1)
summary(my.lm2.4)

Call:
lm(formula = magic_score ~ D1 + D2, data = df2.1)

Residuals:
     Min       1Q   Median       3Q      Max 
-29.7311  -6.3261   0.3617   5.8677  25.2248 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   81.240      1.321  61.498  < 2e-16 ***
D1            14.055      1.868   7.523 4.91e-12 ***
D2            13.131      1.868   7.029 7.26e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.341 on 147 degrees of freedom
Multiple R-squared:  0.3252,    Adjusted R-squared:  0.316 
F-statistic: 35.41 on 2 and 147 DF,  p-value: 2.794e-13
自我性别认同 R1 R2
机人 1 233
袋人 20 666
空白人 300 4396

  1. 也许有读者会问,我们上面推导的是用两个虚拟变量编码一个3水平的分类变量,也没说不能用三个虚拟变量编码呀?其实我们不妨这么思考,三个水平对应三个组平均数,要分别用方程中的所有回归系数(包括截距)指代。可以看作,所有组平均数已知,而回归系数未知的方程组求解,也就是我们上面列出的:ŷ=b0;ŷ=b0+b1;ŷ=b0+b1\hat y_{机人} = b_0; \hat y_{袋人} = b_0+b_1; \hat y_{空白人} = b_0+b_1 。如果你再加入一个变量,就会新增一个变量的回归系数,变成三个式子求解四个未知数了。换句话说,因为截距已经使用了一个信息,所以只剩下两个信息可以给虚拟变量编码使用了,是不是跟自由度的概念差不多?↩︎

  2. 好吧,我们假定我们把空白人设定为一种特殊的种族,而他们的自我性别认同是空白。↩︎

  3. 最典型的例子就是生理性别、族裔、国籍等等人口统计学变量了。↩︎

  4. 是的,因为空白人搞事情坏了我们的大事,所以我们不带他们玩了,去总平均里充当分母去吧。↩︎

  5. 虽然我调侃过这事。不过我强烈建议读者们在这一章,拿出纸笔自己认真推导一遍我们使用的所有编码,理解平均数和回归系数之间如何在线性模型中对应,并为我们所用检验假设的。甚至,我建议大家如果将来使用这种方法检验假设时,用这个办法推导一遍,保证回归系数检验的确实是你要的假设。 比如我们这次的推导,把第一个式子代入后两个式子,得到:

    $$ b_1=\hat y_{袋人}-b_0\\ b_2=\hat y_{机人}-b_0 $$

    再将机人和袋人平均数代回第一个式子得到,

    b0=(b0+b2+b0+b1+ŷ)/3 b_0=(b_0+b_2+b_0+b_1+\hat y_{空白人})/3

    三个平均数对应回归系数的关系就求解出来了。是不是不难得到呢?↩︎