原文鏈接：http://tecdat.cn/?p=14569?

每次我們?cè)趹?yīng)用計(jì)量經(jīng)濟(jì)學(xué)課程中遇到實(shí)際應(yīng)用時(shí)，我們都要處理類別變量。學(xué)生也提出了同樣的問題：我們?nèi)绾巫詣?dòng)組合因子水平？有簡單的R函數(shù)嗎？

因此我想編寫一個(gè)R函數(shù)。為了說明這一點(diǎn)，請(qǐng)考慮以下內(nèi)容

1. 


2. 


3. 'data.frame': 200 obs. of ?3 variables:

4. $ y : num ?1.345 1.863 1.946 2.481 0.765 ...

5. $ x1: num ?0.266 0.372 0.573 0.908 0.202 ...

6. $ x2: Factor w/ 10 levels "I","A","H","F",..: 4 4 6 4 3 6 7 3 4 8 ...

7. table(b$x2)[LETTERS[1:10]]

8. 


9. A ?B ?C ?D ?E ?F ?G ?H ?I ?J

10. 11 12 23 34 23 36 12 32 ?3 14

沒有定義一個(gè)（連續(xù)的）因變量，沒有定義一個(gè)連續(xù)的協(xié)變量，也沒有定義一個(gè)分類變量，此處有十個(gè)級(jí)別。我們可以使用

1. plot(b$x1,y,col="white",xlim=c(0,1.1))

2. text(b$x1,y,as.character(b$x2),cex=.5)

?

線性回歸的輸出得出以下預(yù)測

1. for(i in 1:10){

2. 


3. 


4. 


5. lines(u,v)}

?

?

斜率是相同的，我們只需為每個(gè)級(jí)別添加一個(gè)不同的常數(shù)。如我們所見，某些級(jí)別非常接近，因此將它們組合為一個(gè)類別。這是線性回歸的輸出，

1. 


2. Coefficients:

3. Estimate Std. Error t value Pr(>|t|)

4. (Intercept) ?0.843802 ? 0.119655 ? 7.052 3.23e-11 ***

5. x1 ? ? ? ? ? 1.992878 ? 0.053838 ?37.016 ?< 2e-16 ***

6. x2A ? ? ? ? ?0.055500 ? 0.131173 ? 0.423 ? 0.6727

7. x2H ? ? ? ? ?0.009293 ? 0.121626 ? 0.076 ? 0.9392

8. x2F ? ? ? ? -0.177002 ? 0.121020 ?-1.463 ? 0.1452

9. x2B ? ? ? ? -0.218152 ? 0.130192 ?-1.676 ? 0.0955 .

10. x2D ? ? ? ? -0.206970 ? 0.121294 ?-1.706 ? 0.0896 .

11. x2G ? ? ? ? -0.407417 ? 0.129999 ?-3.134 ? 0.0020 **

12. x2C ? ? ? ? -0.526708 ? 0.123690 ?-4.258 3.24e-05 ***

13. x2J ? ? ? ? -0.664281 ? 0.128126 ?-5.185 5.54e-07 ***

14. x2E ? ? ? ? -0.816454 ? 0.123625 ?-6.604 3.94e-10 ***

15. ---

16. Signif. codes: ?0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

17. 


18. Residual standard error: 0.2014 on 189 degrees of freedom

19. Multiple R-squared: ?0.8995, Adjusted R-squared: ?0.8942

20. F-statistic: 169.1 on 10 and 189 DF, ?p-value: < 2.2e-16

21. 


22. AIC

23. 


24. [1] -60.74443

25. 


26. BIC

27. 


28. [1] -21.16463

這里的參考類別是“ I”。看起來我們實(shí)際上可以將該類別與其他幾個(gè)類別結(jié)合起來。這里的一種策略是選擇似乎沒有顯著差異的所有類別，然后運(yùn)行（多個(gè)）測試

1. 


2. 


3. Hypothesis:

4. x2A = 0

5. x2H = 0

6. x2F = 0

7. 


8. Model 1: restricted model

9. Model 2: y ~ x1 + x2

10. 


11. Res.Df ? ?RSS Df Sum of Sq ? ? ?F Pr(>F)

12. 1 ? ?192 8.4651

13. 2 ? ?189 7.6654 ?3 ? 0.79971 6.5726 ?3e-04 ***

14. ---

15. Signif. codes: ?0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

我們可以將這四個(gè)類別結(jié)合在一起。

我們看到更改參考類別時(shí)的情況（在所有類別上循環(huán)）

1. 


2. 


3. plot(1:nlevels(b$x2),1:nlevels(b$x2),col="white",xlab="",ylab="",axes=F,xlim=c(0,10.5),

4. ylim=c(0,10.5))

5. 


6. 


7. 


8. 


9. 


10. points(((1:10))[idx],rep(i,length(idx)),pch=1,cex=2)

11. 


12. points(((1:10))[idx],rep(i,length(idx)),pch=19,cex=2)}

?

?

?

我們看到它是對(duì)稱的：如果將“ H”與“ I”組合，則“ I”也應(yīng)與“ H”組合。

我們可以手動(dòng)預(yù)定義一些順序

1. 


2. 


3. for(i in 1:nlevels(b$x2)){

4. 


5. 


6. 


7. 


8. 


9. points(((1:10))[idx],rep(i,length(idx)),pch=1,cex=2)

10. idx=which(p>.1)

11. points(((1:10))[idx],rep(i,length(idx)),pch=19,cex=2)

12. }

我們得到

?

我們已經(jīng)合并了類別。

實(shí)際上，可以使用其他策略。我們從某個(gè)級(jí)別開始，說“ A”。然后，我們將其與所有不顯著不同的級(jí)別合并。如果“ B”不是其中之一，我們將其用作新參考。

1. for(i in 1:nlevels(b$x2)){

2. 


3. 


4. 


5. 


6. 


7. 


8. b$x2=recode(b$x2, paste("c('",paste(mix,collapse = "','"),"')='",paste(mix,collapse = "+"),"'",sep=""))

9. }}

最后的類別是

1. table(b$x2)

2. 


3. A+I+H B+D+F ? C+G ? ? E ? ? J

4. 46 ? ?82 ? ?35 ? ?23 ? ?14

有以下回歸輸出

1. 


2. Coefficients:

3. Estimate Std. Error t value Pr(>|t|)

4. (Intercept) ?0.86407 ? ?0.03950 ?21.877 ?< 2e-16 ***

5. x1 ? ? ? ? ? 1.99180 ? ?0.05323 ?37.417 ?< 2e-16 ***

6. x2B+D+F ? ? -0.21517 ? ?0.03699 ?-5.817 2.44e-08 ***

7. x2C+G ? ? ? -0.50545 ? ?0.04528 -11.164 ?< 2e-16 ***

8. x2E ? ? ? ? -0.83617 ? ?0.05128 -16.305 ?< 2e-16 ***

9. x2J ? ? ? ? -0.68398 ? ?0.06131 -11.156 ?< 2e-16 ***

10. ---

11. Signif. codes: ?0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

12. 


13. Residual standard error: 0.2008 on 194 degrees of freedom

14. Multiple R-squared: ?0.8975, Adjusted R-squared: ?0.8948

15. F-statistic: 339.6 on 5 and 194 DF, ?p-value: < 2.2e-16

16. 


17. AIC

18. [1] -66.76939

19. BIC

20. [1] -43.68117

這與我們之前得到的組一致。但是，如果我們更改順序，我們可以得到不同的組合。例如，如果我們從“ J”到“ A”，而不是從“ A”到“ J”，我們得到

1. for(i in nlevels(b$x2):1){

2. 


3. 


4. 


5. 


6. 


7. mix=c(LETTERS[i],names(p)[idx])

8. b$x2=recode(b$x2, paste("c('",paste(mix,collapse = "','"),"')='",paste(mix,collapse = "+"),"'",sep=""))

9. }}

10. table(b$x2)

11. 


12. E ? ? ? ? G+C I+A+B+D+F+H ? ? ? ? ? J

13. 23 ? ? ? ? ?35 ? ? ? ? 128 ? ? ? ? ?14

這里有不同的信息標(biāo)準(zhǔn)

1. AIC(lm(y~x1+x2,data=b))

2. [1] -36.61665

3. BIC(lm(y~x1+x2,data=b))

4. [1] -16.82675

最后但重要的一點(diǎn)是，可以使用回歸樹。問題是還有另一個(gè)可能相互干擾的解釋變量。所以我建議（1）擬合線性模型，計(jì)算殘差（2）運(yùn)行回歸樹，解釋未定義分類變量

?

?

觀察葉子與我們得到的葉子具有相同的組。

1. arbre

2. n= 200

3. 


4. node), split, n, deviance, yval

5. * denotes terminal node

6. 


7. 1) root 200 22.563500 ?7.771561e-18

8. 2) x2=G,C,J,E 72 ?4.441495 -3.232525e-01

9. 4) x2=J,E 37 ?1.553520 -4.578492e-01 *

10. 5) x2=G,C 35 ?1.509068 -1.809646e-01 *

11. 3) x2=I,A,H,F,B,D 128 ?6.366628 ?1.818295e-01

12. 6) x2=F,B,D 82 ?2.983381 ?1.048246e-01 *

13. 7) x2=I,A,H 46 ?2.030229 ?3.190993e-01 *

我想有可能改善回歸的水平組合。請(qǐng)選中你要保存的內(nèi)容，粘貼到此文本框