了解为什么线性回归没有按预期处理我的分类变量？

Question

我正在阅读有关公式和线性回归的信息，但我无法理解如何解释 lm 的输出以用于具有多个参数和分类变量的线性回归。

我想我了解如何解释简单 y = 的输出 a + bx 公式（如果我在下面说的有误，请纠正我）。

#library(tidyverse)
#library(modelr)
require(ggplot2)
require(dplyr)

diamonds2 <- diamonds %>%
  mutate(lprice = log2(price), lcarat = log2(carat))

mod <- lm(
  lprice ~ lcarat,
  data = diamonds2
)

#diamonds2 %>% modelr::add_predictions(mod, "pred_price")
diamonds2$pred_price <- predict(mod, diamonds2) # if you don't have modelr

型号（mod）是

Call:
lm(formula = lprice ~ lcarat, data = diamonds2)

Coefficients:
(Intercept)       lcarat  
     12.189        1.676  

据我了解，这意味着当我添加预测时，我生成预测的公式是

pred_price = 12.189 + (1.676 * lcarat)

当我将分类变量添加到我的公式时，我感到困惑

diamonds2 <- diamonds %>%
  mutate(lprice = log2(price), lcarat = log2(carat))

mod <- lm(
  lprice ~ lcarat + cut,   # I added a categorical variable here
  data = diamonds2
)

diamonds2 %>%
  add_predictions(mod, "pred_price")

现在模型是

Call:
lm(formula = lprice ~ lcarat + cut, data = diamonds2)

Coefficients:
(Intercept)       lcarat        cut.L        cut.Q        cut.C        cut^4  
   12.10711      1.69577      0.32364     -0.09583      0.07631      0.02688  

我对一些事情感到困惑。

1) diamonds$cut 有五个可能的值（一般、良好、非常好、优质、理想），那么为什么模型只显示四个切割值？

2) 据我了解，R 在线性回归方程中将分类变量视为 1 或 0，因此在评估数据行时，每个“切割”系数将乘以 1 或 0。对吗？

3) 我如何根据上面给出的系数写出y = a_0 + (a_1 * x_1) + (a_2 * x_2)...？在这种情况下可能吗？

Answer 1

1) lm 做了一些不需要的事情，它将 diamonds$cut 变量视为有序分类而不是分类（即不使用通常的 1/0 虚拟变量对比处理）。

最初，我认为您只需要通过编写 lm(formula = lprice ~ lcarat + factor(cut)) 或修复数据框 diamonds2$cut <- factor(diamonds2$cut) 来确保 lm 得到分类。

您希望看到分类cut 的对比度级别。 （5 级分类会给出 4 个对比（和一个截距）；请参阅帮助文档（对比）。但是您没有从 lm 获得对比级别，而是得到多项式系数。

深入研究为什么会发生这种情况，我们注意到str(cut) 告诉我们cut 是一个有序分类（这是罪魁祸首）：

> str(diamonds$cut)
 Ord.factor w/ 5 levels "Fair"<"Good"<..: 5 4 2 4 2 3 3 3 1 3 ...

进一步挖掘 lm 的原因，我查看了 help(contrast), help(lm) and help(model.matrix.default) 的页面，这导致我找到了 options('contrasts')：

> getOption('contrasts')
        unordered           ordered 
"contr.treatment"      "contr.poly" 

这意味着 lm 用于在 有序 分类（如 diamonds$cut）上生成对比的默认行为是不需要的 contr.poly()。所以要么更改默认值，要么将cut 转换为无序分类，或者创建一个新变量diamonds$cut <- factor(diamonds$cut, ordered=F) 两种解决方案的代码都在底部。

2) 不，如果您将 diamonds$cut 作为分类传递，就会发生这种情况。但是，您得到了cut.L, .Q, .C, ^4...，它们是 cut 值中（不需要的）多项式的线性、二次、三次、四次系数，它试图拟合 lcarat 的观察值（请参阅contrasts 的帮助，它正在调用@987654336 @ 并使用 n=4 级）。这不是你想要的。

另一个迹象是这些关卡应该被命名 cut.Fair, cut.Good, cut.Very_Good, cut.Premium, cut.Ideal（你会从这 5 个中得到 4 个；另一个将被丢弃）。

3）一旦你修复它以将切割视为（无序）因素，你应该得到： lprice = coeff.price * lcarat + coeff.Fair * cut.Fair + coeff.Good * cut.Good + ... + coeff.Ideal * cut.Ideal

FIX/WORKAROUND 1:

# Save the old ordered-categorical, then clobber it with the unordered one so that lm() Does The Right Thing (tm)
diamonds$cut.ordered <- diamonds$cut
diamonds$cut <- factor(diamonds$cut, ordered=F)

OR ELSE FIX 2:
# Make lm treat all categoricals as unordered categoricals, even ordered ones
#options('contrasts' = c('contr.treatment','contr.treatment') )

lm(log(price) ~ log(carat) + cut, data=diamonds)

Coefficients:
 (Intercept)    log(carat)       cutGood  cutVery Good    cutPremium  
      8.2001        1.6958        0.1632        0.2408        0.2382  
    cutIdeal  
      0.3172  

Answer 2

TL;DR

diamonds$cut 有五个可能的值（一般、好、非常好、优质、理想），那么为什么模型只显示四个切割值？

因子通常用比因子中的水平少一个的系数来表示。那是因为您还在估计模型的截距。相反，您在因子的其他级别上获得的信息会显示在截距中。

据我了解，R 在线性回归方程中将分类变量视为 1 或 0，因此在评估数据行时，每个“截断”系数将乘以 1 或 0。对吗？

情况并非总是如此。这对于传统的虚拟编码 (contr.treatment) 来说是正确的，但是还有很多其他方法可以将因子输入到模型中。在您展示的模型中，您有orthogonal polynomial contrast codes。

3) 如何根据上面给出的系数写出 y = a_0 + (a_1 * x_1) + (a_2 * x_2)...？在这种情况下可能吗？

并非不可能，但难度更大（详情见下文）。多项式对比变量不能总是在单组比较中整齐地表示，因为它们代表了各个级别的总体趋势，因此很难根据回归方程来考虑它们。一个不错的近似值是：

lprice = 12.10711 + 1.69577*lcarat + 0.32364*Lin_cut + -0.09583*Qua_cut + 0.07631*Cub_cut + 0.02688*4_cut + 误差

其中Lin_cut是cut的线性趋势，qua_cut是cut的二次趋势，Cub_cut是cut的三次趋势，4_cut是cut的4^趋势。

更长的解释

cut 是一个有序因子，这意味着它是分类的，但它代表一些潜在的连续变量，因此级别的顺序很重要。请注意 R 描述 cut 的方式与另一个因素相比的差异：

> str(diamonds$cut)
 Ord.factor w/ 5 levels "Fair"<"Good"<..: 5 4 2 4 2 3 3 3 1 3 ...
> str(iris$Species)
 Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...

由于有序因子的分析和解释通常与其他因子略有不同，因此 R 的默认值在输入 lm() 模型时会以不同方式处理它们：

> options("contrasts")
$contrasts
        unordered           ordered 
"contr.treatment"      "contr.poly" 

要在lm() 中输入具有 k 个级别的因子作为预测变量，您需要将其转换为 k-1 代码。有几种方法可以做到这一点，它们都是不错的选择；不同之处在于它们会改变您从模型中获得的系数的解释，因此根据您想要回答的问题类型，您需要选择一种编码分类变量的策略而不是另一种。

contr.treatment

contr.treatment 创建有时称为"traditional" dummy codes 的内容。一个水平（因子的第一水平，默认情况下）被视为参考组，然后每个代码代表该参考组与其他水平之间的差异。

> lm(Petal.Width ~ Species, data = iris)

Call:
lm(formula = Petal.Width ~ Species, data = iris)

Coefficients:
      (Intercept)  Speciesversicolor   Speciesvirginica  
            0.246              1.080              1.780  
> levels(iris$Species)
[1] "setosa"     "versicolor" "virginica" 

在此示例中，参考组 (setosa) 中 Petal.Width 的平均值为 0.246，杂色为 0.246 + 1.080 = 1.36，弗吉尼亚为 0.246 + 1.780 = 2.026。

> library(dplyr)
> iris %>% group_by(Species) %>% summarize(Petal.Width = mean(Petal.Width))
# A tibble: 3 × 2
     Species Petal.Width
      <fctr>       <dbl>
1     setosa       0.246
2 versicolor       1.326
3  virginica       2.026

R 在后台自动为您进行虚拟编码，但您可以随时检查：

> mod$contrasts
$Species
[1] "contr.treatment"

这就是那些虚拟编码变量的样子：

> contr.treatment(levels(iris$Species))
           versicolor virginica
setosa              0         0
versicolor          1         0
virginica           0         1

创建了两个虚拟代码（此处为两列），因为因子中有三个级别。对于数据集中物种为 setosa 的每个案例，两个虚拟代码均为 0。当物种为杂色时，杂色虚拟为 1，弗吉尼亚虚拟为 0。当物种为弗吉尼亚时，虚拟代码为 1，而其他为 0。

contr.poly

虽然您当然可以使用传统的虚拟代码来表示任何分类变量，但这并不总是最能提供信息的方式。生成的系数根据参考水平测试每个水平，这可能对您的数据没有任何特别的兴趣。如果您在 R 中执行 ?contr.treatment，您会看到几个方便的选项，但如果内置代码不能满足您的需求，您也可以从头开始编写自己的代码。

对于有序因子，R 假设多项式趋势对比在大多数情况下最有用，这就是为什么它是默认值。你可以看看它是如何工作的：

> contr.poly(levels(diamonds$cut))
             .L         .Q            .C         ^4
[1,] -0.6324555  0.5345225 -3.162278e-01  0.1195229
[2,] -0.3162278 -0.2672612  6.324555e-01 -0.4780914
[3,]  0.0000000 -0.5345225 -4.095972e-16  0.7171372
[4,]  0.3162278 -0.2672612 -6.324555e-01 -0.4780914
[5,]  0.6324555  0.5345225  3.162278e-01  0.1195229

这些代码不像contr.treatment 代码那样简单易懂，但绘图可能会有所帮助：

library(tidyr)
library(ggplot2)

contr.poly(levels(diamonds$cut)) %>% 
as.data.frame() %>% 
mutate(level=1:nrow(codes)) %>% 
gather("key", "value", -level) %>% 
ggplot(aes(x=level, y=value,color = key)) + 
geom_line()

这使得线性趋势 (L) 的代码形成一条直线，而二次趋势的代码形成 U 形，三次趋势的代码形成一种倾斜的 N-形状，^4 趋势形成一个中间有尖峰的 U。对比码可以解释为这些趋势中的每一个，所以L码解释为数据中的线性趋势，Q码解释为数据中的二次趋势等。

数据中的每个案例都获取这四个对比代码中的每一个的值，而这些对比代码变量用于估计模型。例如，对于 cut="Fair" 的变量，线性对比代码变量的值将为 -0.632，二次的值为 0.534，三次的值为 -.316，^4 的值为 0.119。

对于您的模型，您最终会得到一个正的 cut 线性趋势（cut.L 的系数是正的，并且与零显着不同，您可以通过运行 summary(mod) 看到这一点）。这意味着切割得越好，控制对数（克拉），对数（价格）就越高：理想的价格高于溢价，溢价高于非常好的价格，等等。不过，您也会看到负二次趋势，这表示倒置的 U 形。这表明中等质量削减的对数（价格）高于线性趋势的预期——正线性加负二次。正立方表明 log（价格）从 Good 到 Premium 有一些下降，或者至少比线性和二次趋势预期的增长要小。 ^4 趋势表明，非常好切工的对数（价格）相对于好和溢价高于预期。

进一步阅读

有关多项式趋势对比的更深入解释，请参阅this excellent answer on Cross Validated。