R中具有分类变量的线性模型
Linear model with categorical variables in R
我正在尝试用一些分类变量拟合线性模型
model <- lm(price ~ carat+cut+color+clarity)
summary(model)
答案是：
Call:
lm(formula = price ~ carat + cut + color + clarity)
Residuals:
Min 1Q Median 3Q Max
-11495.7 -688.5 -204.1 458.2 9305.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3696.818 47.948 -77.100 < 2e-16 ***
carat 8843.877 40.885 216.311 < 2e-16 ***
cut.L 755.474 68.378 11.049 < 2e-16 ***
cut.Q -349.587 60.432 -5.785 7.74e-09 ***
cut.C 200.008 52.260 3.827 0.000131 ***
cut^4 12.748 42.642 0.299 0.764994
color.L 1905.109 61.050 31.206 < 2e-16 ***
color.Q -675.265 56.056 -12.046 < 2e-16 ***
color.C 197.903 51.932 3.811 0.000140 ***
color^4 71.054 46.940 1.514 0.130165
color^5 2.867 44.586 0.064 0.948729
color^6 50.531 40.771 1.239 0.215268
clarity.L 4045.728 108.363 37.335 < 2e-16 ***
clarity.Q -1545.178 102.668 -15.050 < 2e-16 ***
clarity.C 999.911 88.301 11.324 < 2e-16 ***
clarity^4 -665.130 66.212 -10.045 < 2e-16 ***
clarity^5 920.987 55.012 16.742 < 2e-16 ***
clarity^6 -712.168 52.346 -13.605 < 2e-16 ***
clarity^7 1008.604 45.842 22.002 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1167 on 4639 degrees of freedom
Multiple R-squared: 0.9162, Adjusted R-squared: 0.9159
F-statistic: 2817 on 18 and 4639 DF, p-value: < 2.2e-16
但是我不明白为什么答案是“.L,.Q,.C,^4, ...”，有问题但我不知道是什么问题，我已经试过了每个变量的函数因子。
您遇到了回归函数如何处理“有序”（有序）因子变量的问题，默认的对比集是最高 n-1 次的正交多项式对比，其中 n 是该因子的水平数。解释这个结果并不是一件容易的事……尤其是在没有自然规律的情况下。即使存在，并且在这种情况下很可能存在，您可能不希望默认排序（按因子级别按字母顺序排列）并且您可能不希望多项式对比中的度数超过几度。
对于 ggplot2 的钻石数据集，因子级别设置正确，但大多数新手在偶然发现有序因子时会获得有序级别，例如 "Excellent" <"Fair" < "Good" < "Poor"。 （失败）
> levels(diamonds$cut)
[1] "Fair" "Good" "Very Good" "Premium" "Ideal"
> levels(diamonds$clarity)
[1] "I1" "SI2" "SI1" "VS2" "VS1" "VVS2" "VVS1" "IF"
> levels(diamonds$color)
[1] "D" "E" "F" "G" "H" "I" "J"
一种在正确设置后使用有序因子的方法是将它们包裹在 as.numeric 中，这样可以对趋势进行线性测试。
> contrasts(diamonds$cut) <- contr.treatment(5) # Removes ordering
> model <- lm(price ~ carat+cut+as.numeric(color)+as.numeric(clarity), diamonds)
> summary(model)
Call:
lm(formula = price ~ carat + cut + as.numeric(color) + as.numeric(clarity),
data = diamonds)
Residuals:
Min 1Q Median 3Q Max
-19130.3 -696.1 -176.8 556.9 9599.8
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5189.460 36.577 -141.88 <2e-16 ***
carat 8791.452 12.659 694.46 <2e-16 ***
cut2 909.433 35.346 25.73 <2e-16 ***
cut3 1129.518 32.772 34.47 <2e-16 ***
cut4 1156.989 32.427 35.68 <2e-16 ***
cut5 1264.128 32.160 39.31 <2e-16 ***
as.numeric(color) -318.518 3.282 -97.05 <2e-16 ***
as.numeric(clarity) 522.198 3.521 148.31 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1227 on 53932 degrees of freedom
Multiple R-squared: 0.9054, Adjusted R-squared: 0.9054
F-statistic: 7.371e+04 on 7 and 53932 DF, p-value: < 2.2e-16
因为@Roland 没有 post 他认为更好的方法（我有点同意他的看法），所以我需要自学一个真正的统计学家如何在 R 中做到这一点. 我最终在@SvenHohenstein 的 post 中找到了关于 SO 的正确编码建议： How to properly set contrasts in R 我喜欢使用 as.numeric 的原因是我知道如何解释系数。系数是 LHS 结果或因变量水平差异的 'effect'（记住工作 'effect' 并不意味着因果关系）。查看我目前位于该答案上方的第一个答案，您会看到 cut2-5 的系数值约为 1000，而 cut1 的系数没有值。 "value" 对 cut==1 的贡献隐藏在“（拦截）”中。估计值如下：
> cbind( levels(diamonds$cut), c(coef(model.cut)[grep('Intercept|cut', names(coef(model.cut)))] ))
[,1] [,2]
(Intercept) "Fair" "-5189.46034442502"
cut2 "Good" "909.432743872746"
cut3 "Very Good" "1129.51839934007"
cut4 "Premium" "1156.98898349819"
cut5 "Ideal" "1264.12800574865"
您可以绘制未调整的均值，但未调整的值并没有真正意义，（因此强调回归分析的必要性）：
> with(diamonds, tapply(price, cut, mean))
Fair Good Very Good Premium Ideal
4358.758 3928.864 3981.760 4584.258 3457.542
所以看看在carat的五分位数内切割的效果：
> with(diamonds, tapply(price, list(cut, cut2(carat, g=5) ), mean))
[0.20,0.36) [0.36,0.54) [0.54,0.91) [0.91,1.14) [1.14,5.01]
Fair 802.4528 1193.162 2336.543 4001.972 8682.351
Good 574.7482 1101.406 2701.412 4872.072 9788.294
Very Good 597.9258 1151.537 2727.251 5464.223 10158.057
Premium 717.1096 1149.550 2537.446 5214.787 10131.999
Ideal 739.8972 1254.229 2624.180 6050.358 10317.725
那么影响...什么？对于双向分析，'cut' 的整个值范围内的平均值可能为 800？
contrasts(diamonds$cut, how.many=1) <- poly(1:5)
> model.cut2 <- lm(price ~ carat+cut, diamonds)
> model.cut2
Call:
lm(formula = price ~ carat + cut, data = diamonds)
Coefficients:
(Intercept) carat cut1
-2555.1 7838.5 815.8
> contrasts(diamonds$cut)
1
Fair -0.6324555
Good -0.3162278
Very Good 0.0000000
Premium 0.3162278
Ideal 0.6324555
公平与理想的估计价格持有 carat 常数的平均差异为 (-0.6324555 -0.6324555)*815.8 或负 1031.91 的价格差异（美元或欧元......无论单位如何价格变量）
我决定删除我要放在这里的一堆其他东西，因为我认为这充分证明了我的要点，即人们需要理解底层编码才能正确解释和传达量级"effects"。单独的系数没有意义。 poly 的线性对比创建了一个效果系数，本质上是 "full" 范围差异。如果使用 R poly()，则需要使用对比矩阵值 和估计系数 进行比较。对比度的范围通常在 1 左右，线性对比度以 0 为中心。
此处节省功率的合理先验方法是评估线性、四次和三次对比。这允许最合理的模型，并避免测试大量级别允许的那些高阶多项式，但如果在理论上依赖，奥卡姆的威廉会不适:-)
library(ggplot2)
df = diamonds[1:1000, ] # a chunk of data
contrasts(df$cut , how.many=3) = contr.poly(nlevels(df$cut))
contrasts(df$color , how.many=3) = contr.poly(nlevels(df$color))
contrasts(df$clarity, how.many=3) = contr.poly(nlevels(df$clarity))
model <- lm(price ~ carat+cut+color+clarity, data = df)
summary(model)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -692.74 30.99 -22.353 < 2e-16 ***
carat 4444.79 41.37 107.431 < 2e-16 ***
cut.L 286.35 22.31 12.835 < 2e-16 ***
cut.Q -88.61 20.26 -4.374 1.35e-05 ***
cut.C 120.91 18.51 6.532 1.03e-10 ***
color.L -660.17 24.93 -26.476 < 2e-16 ***
color.Q -119.34 23.90 -4.993 7.03e-07 ***
color.C 37.18 20.90 1.779 0.0756 .
clarity.L 1356.12 43.22 31.375 < 2e-16 ***
clarity.Q -220.86 33.48 -6.596 6.87e-11 ***
clarity.C 375.47 31.10 12.073 < 2e-16 ***
Multiple R-squared: 0.929, Adjusted R-squared: 0.9283
F-statistic: 1293 on 10 and 989 DF, p-value: < 2.2e-16
我正在尝试用一些分类变量拟合线性模型
model <- lm(price ~ carat+cut+color+clarity)
summary(model)
答案是：
Call:
lm(formula = price ~ carat + cut + color + clarity)
Residuals:
Min 1Q Median 3Q Max
-11495.7 -688.5 -204.1 458.2 9305.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3696.818 47.948 -77.100 < 2e-16 ***
carat 8843.877 40.885 216.311 < 2e-16 ***
cut.L 755.474 68.378 11.049 < 2e-16 ***
cut.Q -349.587 60.432 -5.785 7.74e-09 ***
cut.C 200.008 52.260 3.827 0.000131 ***
cut^4 12.748 42.642 0.299 0.764994
color.L 1905.109 61.050 31.206 < 2e-16 ***
color.Q -675.265 56.056 -12.046 < 2e-16 ***
color.C 197.903 51.932 3.811 0.000140 ***
color^4 71.054 46.940 1.514 0.130165
color^5 2.867 44.586 0.064 0.948729
color^6 50.531 40.771 1.239 0.215268
clarity.L 4045.728 108.363 37.335 < 2e-16 ***
clarity.Q -1545.178 102.668 -15.050 < 2e-16 ***
clarity.C 999.911 88.301 11.324 < 2e-16 ***
clarity^4 -665.130 66.212 -10.045 < 2e-16 ***
clarity^5 920.987 55.012 16.742 < 2e-16 ***
clarity^6 -712.168 52.346 -13.605 < 2e-16 ***
clarity^7 1008.604 45.842 22.002 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1167 on 4639 degrees of freedom
Multiple R-squared: 0.9162, Adjusted R-squared: 0.9159
F-statistic: 2817 on 18 and 4639 DF, p-value: < 2.2e-16
但是我不明白为什么答案是“.L,.Q,.C,^4, ...”，有问题但我不知道是什么问题，我已经试过了每个变量的函数因子。
您遇到了回归函数如何处理“有序”（有序）因子变量的问题，默认的对比集是最高 n-1 次的正交多项式对比，其中 n 是该因子的水平数。解释这个结果并不是一件容易的事……尤其是在没有自然规律的情况下。即使存在，并且在这种情况下很可能存在，您可能不希望默认排序（按因子级别按字母顺序排列）并且您可能不希望多项式对比中的度数超过几度。
对于 ggplot2 的钻石数据集，因子级别设置正确，但大多数新手在偶然发现有序因子时会获得有序级别，例如 "Excellent" <"Fair" < "Good" < "Poor"。 （失败）
> levels(diamonds$cut)
[1] "Fair" "Good" "Very Good" "Premium" "Ideal"
> levels(diamonds$clarity)
[1] "I1" "SI2" "SI1" "VS2" "VS1" "VVS2" "VVS1" "IF"
> levels(diamonds$color)
[1] "D" "E" "F" "G" "H" "I" "J"
一种在正确设置后使用有序因子的方法是将它们包裹在 as.numeric 中，这样可以对趋势进行线性测试。
> contrasts(diamonds$cut) <- contr.treatment(5) # Removes ordering
> model <- lm(price ~ carat+cut+as.numeric(color)+as.numeric(clarity), diamonds)
> summary(model)
Call:
lm(formula = price ~ carat + cut + as.numeric(color) + as.numeric(clarity),
data = diamonds)
Residuals:
Min 1Q Median 3Q Max
-19130.3 -696.1 -176.8 556.9 9599.8
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5189.460 36.577 -141.88 <2e-16 ***
carat 8791.452 12.659 694.46 <2e-16 ***
cut2 909.433 35.346 25.73 <2e-16 ***
cut3 1129.518 32.772 34.47 <2e-16 ***
cut4 1156.989 32.427 35.68 <2e-16 ***
cut5 1264.128 32.160 39.31 <2e-16 ***
as.numeric(color) -318.518 3.282 -97.05 <2e-16 ***
as.numeric(clarity) 522.198 3.521 148.31 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1227 on 53932 degrees of freedom
Multiple R-squared: 0.9054, Adjusted R-squared: 0.9054
F-statistic: 7.371e+04 on 7 and 53932 DF, p-value: < 2.2e-16
因为@Roland 没有 post 他认为更好的方法（我有点同意他的看法），所以我需要自学一个真正的统计学家如何在 R 中做到这一点. 我最终在@SvenHohenstein 的 post 中找到了关于 SO 的正确编码建议： How to properly set contrasts in R 我喜欢使用 as.numeric 的原因是我知道如何解释系数。系数是 LHS 结果或因变量水平差异的 'effect'（记住工作 'effect' 并不意味着因果关系）。查看我目前位于该答案上方的第一个答案，您会看到 cut2-5 的系数值约为 1000，而 cut1 的系数没有值。 "value" 对 cut==1 的贡献隐藏在“（拦截）”中。估计值如下：
> cbind( levels(diamonds$cut), c(coef(model.cut)[grep('Intercept|cut', names(coef(model.cut)))] ))
[,1] [,2]
(Intercept) "Fair" "-5189.46034442502"
cut2 "Good" "909.432743872746"
cut3 "Very Good" "1129.51839934007"
cut4 "Premium" "1156.98898349819"
cut5 "Ideal" "1264.12800574865"
您可以绘制未调整的均值，但未调整的值并没有真正意义，（因此强调回归分析的必要性）：
> with(diamonds, tapply(price, cut, mean))
Fair Good Very Good Premium Ideal
4358.758 3928.864 3981.760 4584.258 3457.542
所以看看在carat的五分位数内切割的效果：
> with(diamonds, tapply(price, list(cut, cut2(carat, g=5) ), mean))
[0.20,0.36) [0.36,0.54) [0.54,0.91) [0.91,1.14) [1.14,5.01]
Fair 802.4528 1193.162 2336.543 4001.972 8682.351
Good 574.7482 1101.406 2701.412 4872.072 9788.294
Very Good 597.9258 1151.537 2727.251 5464.223 10158.057
Premium 717.1096 1149.550 2537.446 5214.787 10131.999
Ideal 739.8972 1254.229 2624.180 6050.358 10317.725
那么影响...什么？对于双向分析，'cut' 的整个值范围内的平均值可能为 800？
contrasts(diamonds$cut, how.many=1) <- poly(1:5)
> model.cut2 <- lm(price ~ carat+cut, diamonds)
> model.cut2
Call:
lm(formula = price ~ carat + cut, data = diamonds)
Coefficients:
(Intercept) carat cut1
-2555.1 7838.5 815.8
> contrasts(diamonds$cut)
1
Fair -0.6324555
Good -0.3162278
Very Good 0.0000000
Premium 0.3162278
Ideal 0.6324555
公平与理想的估计价格持有 carat 常数的平均差异为 (-0.6324555 -0.6324555)*815.8 或负 1031.91 的价格差异（美元或欧元......无论单位如何价格变量）
我决定删除我要放在这里的一堆其他东西，因为我认为这充分证明了我的要点，即人们需要理解底层编码才能正确解释和传达量级"effects"。单独的系数没有意义。 poly 的线性对比创建了一个效果系数，本质上是 "full" 范围差异。如果使用 R poly()，则需要使用对比矩阵值 和估计系数 进行比较。对比度的范围通常在 1 左右，线性对比度以 0 为中心。
此处节省功率的合理先验方法是评估线性、四次和三次对比。这允许最合理的模型，并避免测试大量级别允许的那些高阶多项式，但如果在理论上依赖，奥卡姆的威廉会不适:-)
library(ggplot2)
df = diamonds[1:1000, ] # a chunk of data
contrasts(df$cut , how.many=3) = contr.poly(nlevels(df$cut))
contrasts(df$color , how.many=3) = contr.poly(nlevels(df$color))
contrasts(df$clarity, how.many=3) = contr.poly(nlevels(df$clarity))
model <- lm(price ~ carat+cut+color+clarity, data = df)
summary(model)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -692.74 30.99 -22.353 < 2e-16 ***
carat 4444.79 41.37 107.431 < 2e-16 ***
cut.L 286.35 22.31 12.835 < 2e-16 ***
cut.Q -88.61 20.26 -4.374 1.35e-05 ***
cut.C 120.91 18.51 6.532 1.03e-10 ***
color.L -660.17 24.93 -26.476 < 2e-16 ***
color.Q -119.34 23.90 -4.993 7.03e-07 ***
color.C 37.18 20.90 1.779 0.0756 .
clarity.L 1356.12 43.22 31.375 < 2e-16 ***
clarity.Q -220.86 33.48 -6.596 6.87e-11 ***
clarity.C 375.47 31.10 12.073 < 2e-16 ***
Multiple R-squared: 0.929, Adjusted R-squared: 0.9283
F-statistic: 1293 on 10 and 989 DF, p-value: < 2.2e-16