在R中使用分类变量的线性模型

Question

在R中使用分类变量的线性模型

4

我正试图用一些分类变量去拟合一个线性模型

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,..."，似乎出了些问题，但我不知道具体是哪里出了问题。我已经尝试使用函数对每个变量进行因子化。

- Rashid Dergal Rufeil

检查分类变量的类别。它必须是因子(factor)。因子中的每个级别都将有一个系数。 - vagabond

4

这是一个序数词。?ordered 可能表示“有序的”。 - Neal Fultz

我刚刚解决了我的问题，谢谢。 - Rashid Dergal Rufeil

3个回答

5

由于@Roland没有发布他认为更好的方法（我有点同意他的观点），所以我需要自己学习如何在R中使用真正的统计学方法。最终，我在SO上找到了正确的编码建议，这是由@SvenHohenstein在一个帖子中发布的：如何在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

因此，我们来看一下在每个克拉重量的五等分内削减的影响：

> 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

持有克拉恒定的情况下，公平与理想之间的估计价格差异平均为(-0.6324555 -0.6324555)*815.8，即价格差异为负1031.91（美元或欧元...无论价格变量的单位是什么）。

我已经决定删除我原本要放在这里的其他内容，因为我认为这足以证明我的基本观点：在正确解释和传达“效应”的大小时，需要了解底层编码。单独的系数并没有意义。来自poly的线性对比会创建一个“全”范围差异的效应系数。如果使用R中的poly()，则需要使用对比矩阵值和估计系数进行比较。对比范围通常在1左右，线性对比以0为中心。

- IRTFM

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

- tim

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- IRTFM · Accepted Answer

您遇到了如何处理“有序”（ordinal）因子变量的回归函数和默认对比集是正交多项式对比，最高次数为n-1，其中n是该因子的级别数。如果没有自然顺序，那么解释这个结果将不会很容易...即使有，而且在这种情况下可能确实存在，您可能不想要默认排序（按因子水平字母顺序），并且您可能不想在多项式对比中有更多的度数。

在ggplot2的钻石数据集的情况下，因子水平被正确设置，但大多数新手在遇到有序因子时会得到像“优秀”＆lt;“一般”<“良好”＜“差”这样的有序级别。（失败）

> 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