R中具有分类变量的线性模型

Question

我正在尝试用一些分类变量拟合线性模型

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, ...”，有问题但我不知道有什么问题，我已经尝试过每个变量的函数因子。

Answer 1

您会遇到回归函数如何处理“有序”（有序）因子变量，并且默认的对比集是直到 n-1 次的正交多项式对比，其中 n 是该因子的水平数。解释这个结果并不容易……尤其是在没有自然顺序的情况下。即使有，而且在这种情况下很可能有，您可能不想要默认排序（按因子级别按字母顺序排列），并且您可能不希望多项式对比中的度数超过几个。

在 ggplot2 的 diamonds 数据集的情况下，因子水平设置正确，但大多数新手在偶然发现有序因子时会得到有序水平，例如“优秀”

> 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

Answer 2

由于@Roland 没有发布他认为更好的方法（我有点同意他的观点），我需要自学真正的统计学家如何在 R 中做到这一点。我最终找到了正确的方法在@SvenHohenstein 的帖子中关于 SO 的编码建议：How to properly set contrasts in R 我喜欢使用as.numeric 的原因是我知道如何解释系数。该系数是 LHS 结果或因变量水平上的一级差异的“效果”（请记住，工作“效果”并不意味着因果关系）。看看我目前高于这个答案的第一个答案，您会看到 cut2-5 的系数值约为 1000，而 cut1 没有值。 cut==1 的“值”的贡献隐藏在 '(Intercept)' 中。估计值如下：

> 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

Fair 与 Ideal 的估计价格保持 carat 常数的平均差异为 (-0.6324555 -0.6324555)*815.8 或负 1031.91 的价格差异（美元或欧元....无论价格变量的单位是多少)

我决定删除我要放在这里的一堆其他东西，因为我认为这充分证明了我的基本观点，即人们需要理解底层编码才能正确解释和传达“影响”。单独的系数没有意义。来自poly 的线性对比创建了一个效果系数，该系数基本上用于“全”范围差异。如果使用 R poly()，则需要使用对比矩阵值和估计系数进行比较。对比度范围通常在 1 左右，线性对比度以 0 为中心。

Answer 3

此处节省功率的合理先验方法是评估线性、四次和三次对比。这允许大多数合理的模型，并避免测试大量级别所允许的那些高阶多项式，但如果在理论上依赖这将使奥卡姆的威廉感到不适:-)

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