了解 GLM 中的虚拟(手动或自动)变量创建

机器算法验证 r 广义线性模型 分类数据 分类编码
2022-02-12 09:02:30

如果在 glm 公式中使用了因子变量(例如,具有 M 和 F 水平的性别),则会创建虚拟变量,并且可以在 glm 模型摘要中找到它们及其相关系数(例如,genderM)

如果不是依靠 R 以这种方式拆分因子,而是将因子编码为一系列数字 0/1 变量(例如,genderM(1 表示 M,0 表示 F),genderF(1 表示 F,0 表示M)然后这些变量在glm公式中用作数值变量,系数结果会有所不同吗?

基本上问题是:在使用因子变量和数值变量时,R 是否使用不同的系数计算?

后续问题(可能由上述回答):除了让 R 创建虚拟变量的效率之外,将因子重新编码为一系列数字 0,1 变量并使用模型中的变量是否有任何问题?

2个回答

分类变量(在 R 中称为“因子”)在多元回归模型中需要用数字代码表示。有很多可能的方法可以适当地构造数字代码(请参阅UCLA 统计帮助网站上的这个伟大列表)。默认情况下,R 使用参考级别编码(R 称之为“contr.treatment”),这几乎是默认的统计范围。可以使用?options更改整个 R 会话的所有对比,或者使用?contrasts?C(注意大写)更改特定分析/变量。如果您需要有关参考级别编码的更多信息,我会在此处进行解释:Regression based for example on days of week.

有些人发现参考级别编码令人困惑,您不必使用它。如果你愿意,你可以有男性和女性的两个变量;这称为级别意味着编码。但是,如果你这样做,你将需要抑制截距,否则模型矩阵将是奇异的,并且回归不能像上面的 @Affine 注释那样拟合,正如我在这里解释的那样:定性变量编码导致奇异性要抑制截距,您可以通过添加-1or来修改公式+0y~... -1y~... +0.

使用水平意味着编码而不是参考水平编码将改变估计的系数以及与输出一起打印的假设检验的含义。当您有一个两级因子(例如,男性与女性)并使用参考级编码时,您将看到调用的截距(constant)并且在输出中仅列出一个变量(可能是sexM)。截距是参考组(可能是女性)sexM的平均值,是男性平均值与女性平均值之间的差异。与截距相关的 p 值是一个样本-测试参考水平是否具有平均值0与相关的 p 值会sexM告诉您性别是否在您的回答中有所不同。但是,如果您改用水平均值编码,您将列出两个变量,每个 p 值将对应一个样本-测试该水平的平均值是否为0. 也就是说,没有一个 p 值可以测试性别是否不同。

set.seed(1)
y    = c(    rnorm(30), rnorm(30, mean=1)         )
sex  = rep(c("Female",  "Male"          ), each=30)
fem  = ifelse(sex=="Female", 1, 0)
male = ifelse(sex=="Male", 1, 0)

ref.level.coding.model   = lm(y~sex)
level.means.coding.model = lm(y~fem+male+0)

summary(ref.level.coding.model)
# ...
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  0.08246    0.15740   0.524    0.602    
# sexMale      1.05032    0.22260   4.718 1.54e-05 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# ...
summary(level.means.coding.model)
# ...
# Coefficients:
#      Estimate Std. Error t value Pr(>|t|)    
# fem   0.08246    0.15740   0.524    0.602    
# male  1.13277    0.15740   7.197 1.37e-09 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# ...

在您创建与 R 一致的虚拟变量(即数字变量)的条件下,估计的系数将是相同的。例如:让我们创建一个假数据并使用因子拟合 Poisson glm。请注意,gl函数会创建一个因子变量。

> counts <- c(18,17,15,20,10,20,25,13,12)
> outcome <- gl(3,1,9)
> outcome
[1] 1 2 3 1 2 3 1 2 3
Levels: 1 2 3
> class(outcome)
[1] "factor"
> glm.1<- glm(counts ~ outcome, family = poisson())
> summary(glm.1)

Call:
glm(formula = counts ~ outcome, family = poisson())

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.9666  -0.6713  -0.1696   0.8471   1.0494  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   3.0445     0.1260  24.165   <2e-16 ***
outcome2     -0.4543     0.2022  -2.247   0.0246 *  
outcome3     -0.2930     0.1927  -1.520   0.1285    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 10.5814  on 8  degrees of freedom
Residual deviance:  5.1291  on 6  degrees of freedom
AIC: 52.761

Number of Fisher Scoring iterations: 4

由于结果具有三个级别,因此我创建了两个虚拟变量(如果结果=2,则为 dummy.1=0,如果结果=3,则为 dummy.2=1)并使用这些数值进行重新拟合:

> dummy.1=rep(0,9)
> dummy.2=rep(0,9)
> dummy.1[outcome==2]=1
> dummy.2[outcome==3]=1
> glm.2<- glm(counts ~ dummy.1+dummy.2, family = poisson())
> summary(glm.2)

Call:
glm(formula = counts ~ dummy.1 + dummy.2, family = poisson())

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.9666  -0.6713  -0.1696   0.8471   1.0494  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   3.0445     0.1260  24.165   <2e-16 ***
dummy.1      -0.4543     0.2022  -2.247   0.0246 *  
dummy.2      -0.2930     0.1927  -1.520   0.1285    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 10.5814  on 8  degrees of freedom
Residual deviance:  5.1291  on 6  degrees of freedom
AIC: 52.761

Number of Fisher Scoring iterations: 4

如您所见,估计的系数是相同的。但是,如果您想获得相同的结果,则在创建虚拟变量时需要小心。例如,如果我创建两个虚拟变量 (dummy.1=0 if output=1 and dummy.2=1 if output=2) 那么估计结果不同如下:

> dummy.1=rep(0,9)
> dummy.2=rep(0,9)
> dummy.1[outcome==1]=1
> dummy.2[outcome==2]=1
> glm.3<- glm(counts ~ dummy.1+dummy.2, family = poisson())
> summary(glm.3)

Call:
glm(formula = counts ~ dummy.1 + dummy.2, family = poisson())

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.9666  -0.6713  -0.1696   0.8471   1.0494  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   2.7515     0.1459   18.86   <2e-16 ***
dummy.1       0.2930     0.1927    1.52    0.128    
dummy.2      -0.1613     0.2151   -0.75    0.453    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 10.5814  on 8  degrees of freedom
Residual deviance:  5.1291  on 6  degrees of freedom
AIC: 52.761

Number of Fisher Scoring iterations: 4

这是因为当您outcome在 glm.1 中添加变量时,R 默认情况下会创建两个虚拟变量,即outcome2outcome3定义它们dummy.1dummy.2glm.2 中的和类似,即第一级结果是所有其他虚拟变量(outcome2outcome3)设置为零。