哑变量的基本介绍及R语言设置

1. 哑变量的基本介绍【摘自医咖会】

在构建回归模型时，如果自变量X为连续性变量，回归系数β可以解释为：在其他自变量不变的条件下，X每改变一个单位，所引起的因变量Y的平均变化量；如果自变量X为二分类变量，例如是否饮酒（1=是，0=否），则回归系数β可以解释为：其他自变量不变的条件下，X=1（饮酒者）与X=0（不饮酒者）相比，所引起的因变量Y的平均变化量。

但是，当自变量X为多分类变量时，例如职业、学历、血型、疾病严重程度等等，此时仅用一个回归系数来解释多分类变量之间的变化关系，及其对因变量的影响，就显得太不理想。

此时，我们通常会将原始的多分类变量转化为哑变量，每个哑变量只代表某两个级别或若干个级别间的差异，通过构建回归模型，每一个哑变量都能得出一个估计的回归系数，从而使得回归的结果更易于解释，更具有实际意义。

1.1 什么是哑变量？

哑变量（Dummy Variable），又称为虚拟变量、虚设变量或名义变量，从名称上看就知道，它是人为虚设的变量，通常取值为0或1，来反映某个变量的不同属性。对于有n个分类属性的自变量，通常需要选取1个分类作为参照，因此可以产生n-1个哑变量。

将哑变量引入回归模型，虽然使模型变得较为复杂，但可以更直观地反映出该自变量的不同属性对于因变量的影响，提高了模型的精度和准确度。

举一个例子，如职业因素，假设分为学生、农民、工人、公务员、其他共5个分类，其中以“其他职业”作为参照，此时需要设定4个哑变量X1-X4，如下所示：

X1=1，学生；X1=0，非学生；
X2=1，农民；X2=0，非农民；
X3=1，工人；X3=0，非工人；
X4=1，公务员；X4=0，非公务员；
那么对于每一种职业分类，其赋值就可以转化为以下形式：

1.2 什么情况下需要设置哑变量?

(1) 对于无序多分类变量，引入模型时需要转化为哑变量

举一个例子，如血型，一般分为A、B、O、AB四个类型，为无序多分类变量，通常情况下在录入数据的时候，为了使数据量化，我们常会将其赋值为1、2、3、4。

从数字的角度来看，赋值为1、2、3、4后，它们是具有从小到大一定的顺序关系的，而实际上，四种血型之间并没有这种大小关系存在，它们之间应该是相互平等独立的关系。如果按照1、2、3、4赋值并带入到回归模型中是不合理的，此时我们就需要将其转化为哑变量。

(2) 对于有序多分类变量，引入模型时需要酌情考虑

例如疾病的严重程度，一般分为轻、中、重度，可认为是有序多分类变量，通常情况下我们也常会将其赋值为1、2、3（等距）或1、2、4（等比）等形式，通过由小到大的数字关系，来体现疾病严重程度之间一定的等级关系。

但需要注意的是，一旦赋值为上述等距或等比的数值形式，这在某种程度上是认为疾病的严重程度也呈现类似的等距或等比的关系。而事实上由于疾病在临床上的复杂性，不同的严重程度之间并非是严格的等距或等比关系，因此再赋值为上述形式就显得不太合理，此时可以将其转化为哑变量进行量化。

(3) 对于连续性变量，进行变量转化时可以考虑设定为哑变量

对于连续性变量，很多人认为可以直接将其带入到回归模型中即可，但有时我们还需要结合实际的临床意义，对连续性变量作适当的转换。例如年龄，以连续性变量带入模型时，其解释为年龄每增加一岁时对于因变量的影响。但往往年龄增加一岁，其效应是很微弱的，并没有太大的实际意义。

此时，我们可以将年龄这个连续性变量进行离散化，按照10岁一个年龄段进行划分，如0-10、11-20、21-30、31-40等等，将每一组赋值为1、2、3、4，此时构建模型的回归系数就可以解释为年龄每增加10岁时对因变量的影响。

以上赋值方式是基于一个前提，即年龄与因变量之间存在着一定的线性关系。但有时候可能会出现以下情况，例如在年龄段较低和较高的人群中，某种疾病的死亡率较高，而在中青年人群中，死亡率却相对较低，年龄和死亡结局之间呈现一个U字型的关系，此时再将年龄段赋值为1、2、3、4就显得不太合理了。

因此，当我们无法确定自变量和因变量之间的变化关系，将连续性自变量离散化时，可以考虑进行哑变量转换。

还有一种情况，例如将BMI按照临床诊断标准分为体重过低、正常体重、超重、肥胖等几种分类时，由于不同分类之间划分的切点是不等距的，此时赋值为1、2、3就不太符合实际情况，也可以考虑将其转化为哑变量。

1.3 如何设置哑变量的参照组？

在上面的内容中我们提到，对于有n个分类的自变量，需要产生n-1个哑变量，当所有n-1个哑变量取值都为0的时候，这就是该变量的第n类属性，即我们将这类属性作为参照。

例如上面提到的以职业因素为例，共分为学生、农民、工人、公务员、其他共5个分类，设定了4个哑变量，其中职业因素中“其它”这个属性，每个哑变量的赋值均为0，此时我们就将“其它”这个属性作为参照，在最后进行模型解释时，所有类别哑变量的回归系数，均表示该哑变量与参照相比之后对因变量的影响。

在设定哑变量时，应该选择哪一类作为参照呢？

(1) 一般情况下，可以选择有特定意义的，或者有一定顺序水平的类别作为参照
例如，婚姻状态分为未婚、已婚、离异、丧偶等情况，可以将“未婚”作为参照；或者如学历，分为小学、中学、大学、研究生等类别，存在着一定的顺序，可以将“小学”作为参照，以便于回归系数更容易解释。

(2) 可以选择临床正常水平作为参照
例如，BMI按照临床诊断标准分为体重过低、正常体重、超重、肥胖等类别，此时可以选择“正常体重”作为参照，其他分类都与正常体重进行比较，更具有临床实际意义。

(3) 还可以将研究者所关注的重点类别作为参照
例如血型，分为A、B、O、AB四个类型，研究者更关注O型血的人，因此可以将O型作为参照，来分析其他血型与O型相比后对于结局产生影响的差异。

1.4 设置哑变量时的注意事项

(1) 原则上哑变量在模型中应同进同出，也就是说在一个模型中，如果同一个分类变量的不同哑变量，出现了有些哑变量有统计学显著性，有些无统计学显著性的情况下，为了保证所有哑变量代表含义的正确性，应当在模型中纳入所有的哑变量。

(2) 在如何选择哑变量的参照组时需要注意的是，被选为参照的那一类分组，应该保证有一定的样本量。如果参照组样本量太少，则将会导致其他分类与参照相比时，参数估计的标准误较大，可信区间较大，精度降低，会出现估计参数极大或极小的现象。

2. R语言中哑变量的设置

在R语言中对包括分类变量(factor)的数据建模时，一般会将其自动处理为虚拟变量或哑变量(dummy variable)。但有一些特殊的函数，如neuralnet包中的neuralnet函数就不会预处理。如果直接将原始数据扔进去，会出现”requires numeric/complex matrix/vector arguments”需要数值/复数矩阵/矢量参数错误。

这个时候，除了将这些变量删除，我们只能手动将factor variable转换为取值(0,1)的虚拟变量。所用的函数一般有model.matrix(),nnet package中的class.ind()。

2.1 示例数据

下面以UCI的german credit data为例说明。

# 从UCI网站上下载到german.data数据集
data <- read.table("http://archive.ics.uci.edu/ml/machine-learning-databases/statlog/german/german.data")
head(data)
##  V1 V2  V3  V4   V5  V6  V7 V8  V9  V10 V11  V12 V13  V14  V15 V16  V17 V18  V19  V20 V21
## A11  6 A34 A43 1169 A65 A75  4 A93 A101   4 A121  67 A143 A152   2 A173   1 A192 A201   1
## A12 48 A32 A43 5951 A61 A73  2 A92 A101   2 A121  22 A143 A152   1 A173   1 A191 A201   2
## A14 12 A34 A46 2096 A61 A74  2 A93 A101   3 A121  49 A143 A152   1 A172   2 A191 A201   1
## A11 42 A32 A42 7882 A61 A74  2 A93 A103   4 A122  45 A143 A153   1 A173   2 A191 A201   1
## A11 24 A33 A40 4870 A61 A73  3 A93 A101   4 A124  53 A143 A153   2 A173   2 A191 A201   2
## A14 36 A32 A46 9055 A65 A73  2 A93 A101   4 A124  35 A143 A153   1 A172   2 A192 A201   1
str(data)
## 'data.frame':    1000 obs. of  21 variables:
##  $ V1 : Factor w/ 4 levels "A11","A12","A13",..: 1 2 4 1 1 4 4 2 4 2 ...
##  $ V2 : int  6 48 12 42 24 36 24 36 12 30 ...
##  $ V3 : Factor w/ 5 levels "A30","A31","A32",..: 5 3 5 3 4 3 3 3 3 5 ...
##  $ V4 : Factor w/ 10 levels "A40","A41","A410",..: 5 5 8 4 1 8 4 2 5 1 ...
##  $ V5 : int  1169 5951 2096 7882 4870 9055 2835 6948 3059 5234 ...
##  $ V6 : Factor w/ 5 levels "A61","A62","A63",..: 5 1 1 1 1 5 3 1 4 1 ...
##  $ V7 : Factor w/ 5 levels "A71","A72","A73",..: 5 3 4 4 3 3 5 3 4 1 ...
##  $ V8 : int  4 2 2 2 3 2 3 2 2 4 ...
##  $ V9 : Factor w/ 4 levels "A91","A92","A93",..: 3 2 3 3 3 3 3 3 1 4 ...
##  $ V10: Factor w/ 3 levels "A101","A102",..: 1 1 1 3 1 1 1 1 1 1 ...
##  $ V11: int  4 2 3 4 4 4 4 2 4 2 ...
##  $ V12: Factor w/ 4 levels "A121","A122",..: 1 1 1 2 4 4 2 3 1 3 ...
##  $ V13: int  67 22 49 45 53 35 53 35 61 28 ...
##  $ V14: Factor w/ 3 levels "A141","A142",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ V15: Factor w/ 3 levels "A151","A152",..: 2 2 2 3 3 3 2 1 2 2 ...
##  $ V16: int  2 1 1 1 2 1 1 1 1 2 ...
##  $ V17: Factor w/ 4 levels "A171","A172",..: 3 3 2 3 3 2 3 4 2 4 ...
##  $ V18: int  1 1 2 2 2 2 1 1 1 1 ...
##  $ V19: Factor w/ 2 levels "A191","A192": 2 1 1 1 1 2 1 2 1 1 ...
##  $ V20: Factor w/ 2 levels "A201","A202": 1 1 1 1 1 1 1 1 1 1 ...
##  $ V21: int  1 2 1 1 2 1 1 1 1 2 ...

该数据有21个变量，其中V21为目标变量，V1-V20中包括integer和factor两种类型。下面将用V1分类变量（包含4个level）和V2,V5,V8三个数值型变量作为解释变量建模。
首先加载neuralnet包尝试一下，只用数值型变量建模，没有报错。

library("neuralnet")
NNModelAllNum <- neuralnet(V21 ~ V2 + V5 + V8, data)
NNModelAllNum
## Call: neuralnet(formula = V21 ~ V2 + V5 + V8, data = data)
## 
## 1 repetition was calculated.
## 
##         Error Reached Threshold Steps
## 1 104.9993578    0.005128715177    55

当我们把V1放入解释变量中出现了如下错误：

NNModel <- neuralnet(V21 ~ V1 + V2 + V5 + V8, data)
## Error in neurons[[i]] %*% weights[[i]] : 需要数值/复数矩阵/矢量参数

2.2 哑变量设置的4种方式

# (1) stats::model.matrix()  #Construct Design Matrices
# 将V1转化为三个虚拟变量，V1A12,V1A13,V1A14
dummyV1 <- model.matrix(~V1, data)
head(cbind(dummyV1, data$V1))
##   (Intercept) V1A12 V1A13 V1A14  
## 1           1     0     0     0 1
## 2           1     1     0     0 2
## 3           1     0     0     1 4
## 4           1     0     0     0 1
## 5           1     0     0     0 1
## 6           1     0     0     1 4

# 因为model.matrix函数对数值型和2水平类别变量没有影响，所以可以将四个变量一起用该函数生成新的数据集modelData，就可以用该数据集建模了。
modelData <- model.matrix(~V1 + V2 + V5 + V8 + V21, data)
head(modelData)
##   (Intercept) V1A12 V1A13 V1A14 V2   V5 V8 V21
## 1           1     0     0     0  6 1169  4   1
## 2           1     1     0     0 48 5951  2   2
## 3           1     0     0     1 12 2096  2   1
## 4           1     0     0     0 42 7882  2   1
## 5           1     0     0     0 24 4870  3   2
## 6           1     0     0     1 36 9055  2   1

# (2) nnet::class.ind()  #Generates Class Indicator Matrix from a Factor
library("nnet")
dummyV12 <- class.ind(data$V1)
head(cbind(dummyV12, data$V1))
#      A11 A12 A13 A14  
# [1,]   1   0   0   0 1
# [2,]   0   1   0   0 2
# [3,]   0   0   0   1 4
# [4,]   1   0   0   0 1
# [5,]   1   0   0   0 1
# [6,]   0   0   0   1 4
# 注：该结果和model.matrix稍有区别，生成了四个虚拟变量。为避免多重共线性，对于level=n的分类变量只需选取其任意n-1个虚拟变量。

# (3) caret::dummyVars()  #Create A Full Set of Dummy Variables
library("caret")
# 利用dummyVars函数对V1变量进行哑变量处理
a <- dummyVars(~V1,data)
dummyV12 <- predict(a,data)
head(cbind(dummyV12, data$V1))
#   V1.A11 V1.A12 V1.A13 V1.A14  
# 1      1      0      0      0 1
# 2      0      1      0      0 2
# 3      0      0      0      1 4
# 4      1      0      0      0 1
# 5      1      0      0      0 1
# 6      0      0      0      1 4

# 对整个数据（中的分类变量）进行哑变量化
a <- dummyVars(~.,data)
b <- predict(a,data)
b[1:6,1:10]
#   V1.A11 V1.A12 V1.A13 V1.A14 V2 V3.A30 V3.A31 V3.A32 V3.A33 V3.A34
# 1      1      0      0      0  6      0      0      0      0      1
# 2      0      1      0      0 48      0      0      1      0      0
# 3      0      0      0      1 12      0      0      0      0      1
# 4      1      0      0      0 42      0      0      1      0      0
# 5      1      0      0      0 24      0      0      0      1      0
# 6      0      0      0      1 36      0      0      1      0      0


# (4) dummies::dummy()  #Flexible, efficient creation of dummy variables
library(dummies)
dummyV12 <- dummy(data$V1, sep = ".")
head(cbind(dummyV12,data$V1))
#      V1.A11 V1.A12 V1.A13 V1.A14  
# [1,]      1      0      0      0 1
# [2,]      0      1      0      0 2
# [3,]      0      0      0      1 4
# [4,]      1      0      0      0 1
# [5,]      1      0      0      0 1
# [6,]      0      0      0      1 4

2.3 线性回归小实例

如上所述，其实R语言在对包括分类变量(factor)的数据建模时，一般会将其自动处理为虚拟变量或哑变量。下面用一个小例子简单说明一下。

library(tidyverse)
library(car)

# Load the data
data("Salaries", package = "carData")
str(Salaries)
# 'data.frame':	397 obs. of  6 variables:
#   $ rank         : Factor w/ 3 levels "AsstProf","AssocProf",..: 3 3 1 3 3 2 3 3 3 3 ...
# $ discipline   : Factor w/ 2 levels "A","B": 2 2 2 2 2 2 2 2 2 2 ...
# $ yrs.since.phd: int  19 20 4 45 40 6 30 45 21 18 ...
# $ yrs.service  : int  18 16 3 39 41 6 23 45 20 18 ...
# $ sex          : Factor w/ 2 levels "Male","Female": 1 1 1 1 1 1 1 1 1 2 ...
# $ salary       : int  139750 173200 79750 115000 141500 97000 175000 147765 119250 129000 ...
sample_n(Salaries, 3)  # Inspect the data

# R 语言默认以因子的第一水平为参考组
# 二分类变量
contrasts(Salaries$sex)
#           Male        
# Female      0
# Male        1
model <- lm(salary ~ sex, data = Salaries)
summary(model)$coef
#             Estimate Std. Error   t value     Pr(>|t|)
# (Intercept) 101002.41   4809.386 21.001103 2.683482e-66
# sexMale      14088.01   5064.579  2.781674 5.667107e-03

# 修改参考组
Salaries <- Salaries %>% mutate(sex = relevel(sex, ref = "Male"))
contrasts(Salaries$sex)
#        Female
# Male        0
# Female      1

model <- lm(salary ~ sex, data = Salaries)
summary(model)$coef
#              Estimate Std. Error   t value      Pr(>|t|)
# (Intercept) 115090.42   1587.378 72.503463 2.459122e-230
# sexFemale   -14088.01   5064.579 -2.781674  5.667107e-03

# 多分类变量
res <- model.matrix(~rank, data = Salaries)
head(cbind(res, Salaries$rank))
#   (Intercept) rankAssocProf rankProf  
# 1           1             0        1 3
# 2           1             0        1 3
# 3           1             0        0 1
# 4           1             0        1 3
# 5           1             0        1 3
# 6           1             1        0 2

model2 <- lm(salary ~ yrs.service + res[,-1] + discipline + sex, data = Salaries)
summary(model2)

model3 <- lm(salary ~ yrs.service + rank + discipline + sex, data = Salaries)
summary(model3)
Anova(model3)
head(model.matrix(model3))  #查看模型的设计矩阵，rank的取值与head(cbind(res, Salaries$rank))中设置的哑变量一致
#   (Intercept) yrs.service rankAssocProf rankProf disciplineB sexFemale
# 1           1          18             0        1           1         0
# 2           1          16             0        1           1         0
# 3           1           3             0        0           1         0
# 4           1          39             0        1           1         0
# 5           1          41             0        1           1         0
# 6           1           6             1        0           1         0

参考阅读：
医咖会 | SPSS教程：手把手教你设置哑变量及解读结果！
R语言中生成虚拟变量/哑变量
Regression with Categorical Variables: Dummy Coding Essentials in R