理解为什么线性回归没有按预期处理我的分类变量？
Understanding why linear regression isn't treating my categorical variable as expected?
我正在阅读有关公式和线性回归的内容，但我无法理解如何解释 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 有五个可能的值（fair、good、very good、premium、ideal），那么为什么模型只显示 cut 的四个值？
2) 据我了解，R 在线性回归方程中将分类变量视为 1 或 0，因此在评估数据行时每个 "cut" 系数将乘以 1 或 0。对吗？
3) 如何根据上面给出的系数写出 y = a_0 + (a_1 * x_1) + (a_2 * x_2)...？在这种情况下可以吗？
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 转换为无序分类，或者创建一个新的 var diamonds$cut <- factor(diamonds$cut, ordered=F)
两种解决方案的代码都在底部。
2) 不，如果您将 diamonds$cut 作为分类传递，就会发生这种情况。但是，您得到了 cut.L, .Q, .C, ^4...，它们是 cut 值中（不需要的）多项式的线性、二次、三次、四次系数，它试图拟合 lcarat 的观察值（请参阅 contrasts 上的帮助正在调用 contr.poly 并使用 n=4 个级别）。这不是你想要的。
另一个线索是那些级别应该被命名
cut.Fair, cut.Good, cut.Very_Good, cut.Premium, cut.Ideal（你会从这 5 个中得到 4 个；另一个将被删除）。
3) 修复它以将 cut 视为（无序）因子后，您应该得到：
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
TL;DR
diamonds$cut has five possible values (fair, good, very good, premium, ideal), so why does the model only show four values for cut?
因子通常用比因子水平少一个的系数表示。那是因为您还在估计模型的截距。您本应在因子的另一个级别获得的信息在截距中表示。
From my understanding, R treats a categorical variable as being either 1 or 0 in the linear regression equation, so each "cut" coefficient will either be multiplied by 1 or 0 when evaluating a data row. Is that correct?
情况并非总是如此。这对于传统的虚拟编码 (contr.treatment) 来说是正确的，但是还有很多其他方法可以将因子输入到模型中。在您提供的模型中，您有 orthogonal polynomial contrast codes.
3) How do I write a y = a_0 + (a_1 * x_1) + (a_2 * x_2)... from that coefficients given above? Is that possible in this case?
不是不可能，而是比较难（详见下文）。多项式对比变量不能总是在单组比较中整齐地表示，因为它们代表了水平上的总体趋势，因此它们更难根据回归方程来考虑。一个合适的近似值是：
lprice = 12.10711 + 1.69577*lcarat + 0.32364*Lin_cut + -0.09583*Qua_cut + 0.07631*Cub_cut + 0.02688*4_cut + error
其中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 ...
由于有序因子的分析和解释通常与其他因子略有不同，因此在 lm() 模型中输入时，R 的默认值以不同方式对待它们：
> 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"
在此示例中，Petal.Width 的平均值在参考组 (setosa) 中为 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 形，中间有一个尖峰。 contrast codes 可以解释为那些趋势中的每一个，所以 L 代码被解释为数据中的线性趋势，Q 代码是数据中的二次趋势等。
数据中的每个个案都获取这四个对比代码中每一个的值，这些对比代码变量用于估计模型。例如，对于 cut="Fair" 的变量，线性对比代码变量的值为 -0.632，二次值为 0.534，三次值为 -.316，^4 值为 0.119。
对于您的模型，您最终得到了正线性切割趋势（cut.L 的系数为正，并且显着不同于零，您可以通过 运行 summary(mod)).这意味着切割越好，控制对数（克拉），对数（价格）越高：理想的价格高于优质，价格高于非常好，等等。你也看到负二次趋势，虽然，这表明倒U形。这表明中等质量的切割比线性趋势预期的 log(price) 更高——正线性加上负二次。正立方表明 log(price) 从 Good 到 Premium 有一些下降，或者至少比线性和二次趋势预期的增长要少。 ^4 趋势表明 Very Good 切工的对数（价格）相对于 Good 和 Premium 而言高于预期。
进一步阅读
有关多项式趋势对比的更深入的解释，请参阅 this excellent answer on Cross Validated。
我正在阅读有关公式和线性回归的内容，但我无法理解如何解释 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 有五个可能的值（fair、good、very good、premium、ideal），那么为什么模型只显示 cut 的四个值？
2) 据我了解，R 在线性回归方程中将分类变量视为 1 或 0，因此在评估数据行时每个 "cut" 系数将乘以 1 或 0。对吗？
3) 如何根据上面给出的系数写出 y = a_0 + (a_1 * x_1) + (a_2 * x_2)...？在这种情况下可以吗？
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 转换为无序分类，或者创建一个新的 var diamonds$cut <- factor(diamonds$cut, ordered=F)
两种解决方案的代码都在底部。
2) 不，如果您将 diamonds$cut 作为分类传递，就会发生这种情况。但是，您得到了 cut.L, .Q, .C, ^4...，它们是 cut 值中（不需要的）多项式的线性、二次、三次、四次系数，它试图拟合 lcarat 的观察值（请参阅 contrasts 上的帮助正在调用 contr.poly 并使用 n=4 个级别）。这不是你想要的。
另一个线索是那些级别应该被命名
cut.Fair, cut.Good, cut.Very_Good, cut.Premium, cut.Ideal（你会从这 5 个中得到 4 个；另一个将被删除）。
3) 修复它以将 cut 视为（无序）因子后，您应该得到：
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
TL;DR
diamonds$cut has five possible values (fair, good, very good, premium, ideal), so why does the model only show four values for cut?
因子通常用比因子水平少一个的系数表示。那是因为您还在估计模型的截距。您本应在因子的另一个级别获得的信息在截距中表示。
From my understanding, R treats a categorical variable as being either 1 or 0 in the linear regression equation, so each "cut" coefficient will either be multiplied by 1 or 0 when evaluating a data row. Is that correct?
情况并非总是如此。这对于传统的虚拟编码 (contr.treatment) 来说是正确的，但是还有很多其他方法可以将因子输入到模型中。在您提供的模型中，您有 orthogonal polynomial contrast codes.
3) How do I write a y = a_0 + (a_1 * x_1) + (a_2 * x_2)... from that coefficients given above? Is that possible in this case?
不是不可能，而是比较难（详见下文）。多项式对比变量不能总是在单组比较中整齐地表示，因为它们代表了水平上的总体趋势，因此它们更难根据回归方程来考虑。一个合适的近似值是：
lprice = 12.10711 + 1.69577*lcarat + 0.32364*Lin_cut + -0.09583*Qua_cut + 0.07631*Cub_cut + 0.02688*4_cut + error
其中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 ...
由于有序因子的分析和解释通常与其他因子略有不同，因此在 lm() 模型中输入时，R 的默认值以不同方式对待它们：
> 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"
在此示例中，Petal.Width 的平均值在参考组 (setosa) 中为 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 形，中间有一个尖峰。 contrast codes 可以解释为那些趋势中的每一个，所以 L 代码被解释为数据中的线性趋势，Q 代码是数据中的二次趋势等。
数据中的每个个案都获取这四个对比代码中每一个的值，这些对比代码变量用于估计模型。例如，对于 cut="Fair" 的变量，线性对比代码变量的值为 -0.632，二次值为 0.534，三次值为 -.316，^4 值为 0.119。
对于您的模型，您最终得到了正线性切割趋势（cut.L 的系数为正，并且显着不同于零，您可以通过 运行 summary(mod)).这意味着切割越好，控制对数（克拉），对数（价格）越高：理想的价格高于优质，价格高于非常好，等等。你也看到负二次趋势，虽然，这表明倒U形。这表明中等质量的切割比线性趋势预期的 log(price) 更高——正线性加上负二次。正立方表明 log(price) 从 Good 到 Premium 有一些下降，或者至少比线性和二次趋势预期的增长要少。 ^4 趋势表明 Very Good 切工的对数（价格）相对于 Good 和 Premium 而言高于预期。
进一步阅读
有关多项式趋势对比的更深入的解释，请参阅 this excellent answer on Cross Validated。