R中因子水平的自动组合
每次我们在应用计量经济学课程中面对实际应用时，我们都必须处理分类变量。同样的问题也发生在学生身上：我们怎样才能自动地结合因素水平呢？有简单的R函数吗？
在过去的几年里，我确实上传了一些博客文章。但到目前为止没有什么令人满意的。让我写下几行关于可以做什么的话。如果有人想写一个很好的R函数，那就太棒了。为了说明这一想法，请考虑以下(模拟数据集)：
n=200
set.seed(1)
x1=runif(n)
x2=runif(n)
y=1+2*x1-x2+rnorm(n,0,.2)
LB=sample(LETTERS[1:10])
b=data.frame(y=y,x1=x1,
x2=cut(x2,breaks=
c(-1,.05,.1,.2,.35,.4,.55,.65,.8,.9,2),
labels=LB))
str(b)
'data.frame':200 obs. of 3 variables:
$ y : num 1.345 1.863 1.946 2.481 0.765 ...
$ x1: num 0.266 0.372 0.573 0.908 0.202 ...
$ x2: Factor w/ 10 levels "I","A","H","F",..: 4 4 6 4 3 6 7 3 4 8 ...
table(b$x2)[LETTERS[1:10]]
A B C D E F G H I J
11 12 23 34 23 36 12 32 3 14
有一个(连续)因变量y，一个连续协变量x_1和一个范畴变量x_2，具有十个水平。我们可以使用以下方法绘制数据：
plot(b$x1,y,col="white",xlim=c(0,1.1))
text(b$x1,y,as.character(b$x2),cex=.5)
线性回归的输出产生以下预测：
for(i in 1:10){
p=function(x) predict(lm(y~x1+x2,data=b),newdata=data.frame(x1=x,x2=LETTERS[i]))
u=seq(-1,1.065,by=.01)
v=Vectorize(p)(u)
lines(u,v)}
x_1的斜率是相同的，我们只需为每个级别添加一个不同的常数。正如我们所看到的，一些级别非常接近，因此将它们合并成一个类别似乎是合理的。以下是线性回归的输出：
summary(lm(y~x1+x2,data=b))
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.843802 0.119655 7.052 3.23e-11 ***
x1 1.992878 0.053838 37.016 < 2e-16 ***
x2A 0.055500 0.131173 0.423 0.6727
x2H 0.009293 0.121626 0.076 0.9392
x2F -0.177002 0.121020 -1.463 0.1452
x2B -0.218152 0.130192 -1.676 0.0955 .
x2D -0.206970 0.121294 -1.706 0.0896 .
x2G -0.407417 0.129999 -3.134 0.0020 **
x2C -0.526708 0.123690 -4.258 3.24e-05 ***
x2J -0.664281 0.128126 -5.185 5.54e-07 ***
x2E -0.816454 0.123625 -6.604 3.94e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.2014 on 189 degrees of freedom
Multiple R-squared: 0.8995,Adjusted R-squared: 0.8942
F-statistic: 169.1 on 10 and 189 DF, p-value: < 2.2e-16
AIC(lm(y~x1+x2,data=b))
[1] -60.74443
BIC(lm(y~x1+x2,data=b))
[1] -21.16463
在这里，参考类别是“I”。看起来我们可以把这个类别和其他几个类别结合起来。这里的一种策略是选择似乎没有明显区别的所有类别，并运行一个(多个)测试：
library(car)
linearHypothesis(lm(y~x1+x2,data=b), c("x2A = 0", "x2H = 0", "x2F = 0"))
Hypothesis:
x2A = 0
x2H = 0
x2F = 0
Model 1: restricted model
Model 2: y ~ x1 + x2
Res.Df RSS Df Sum of Sq F Pr(>F)
1 192 8.4651
2 189 7.6654 3 0.79971 6.5726 3e-04 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
看来我们可以把这四个类别结合起来。
在这里，我们可以看到当我们更改引用类别时发生了什么(实际上，在所有类别上循环)：
P=matrix(NA,nlevels(b$x2),nlevels(b$x2))
colnames(P)=rownames(P)=LETTERS[1:10]
plot(1:nlevels(b$x2),1:nlevels(b$x2),col="white",xlab="",ylab="",axes=F,xlim=c(0,10.5),
ylim=c(0,10.5))
text(1:10,0,LETTERS[1:10])
text(0,1:10,LETTERS[1:10])
for(i in 1:nlevels(b$x2)){
#levels(b$x2)=LETTERS[1:10]
b$x2=relevel(b$x2,LETTERS[i])
p=summary(lm(y~x1+x2,data=b))$coefficients[-(1:2),4]
names(p)=substr(names(p),3,3)
P[LETTERS[i],names(p)]=p
p=P[LETTERS[i],]
idx=which(p>.05)
points(((1:10))[idx],rep(i,length(idx)),pch=1,cex=2)
idx=which(p>.1)
points(((1:10))[idx],rep(i,length(idx)),pch=19,cex=2)}
我们很高兴看到它是对称的：如果“H”应该与“i”结合，那么“i”也应该与“H”结合。
在这里，黑点与10%p值有关，白色点与5%p值有关.这张图其实很难读.。事实上，这让我们想起Bertin(1967年).
在这里，我们可以手动预定义一些排序(下面我们将看到它是如何自动化的)：
LETTERSord=c("I","A","H","F","B","D","G","C","J","E")
P=matrix(NA,nlevels(b$x2),nlevels(b$x2))
colnames(P)=rownames(P)=LETTERSord
plot(1:nlevels(b$x2),1:nlevels(b$x2),col="white",xlab="",ylab="",axes=F,xlim=c(0,10.5),
ylim=c(0,10.5))
ct=c(3,3,2,1,1)
abline(v=.5+c(0,cumsum(ct)),lty=2)
abline(h=.5+c(0,cumsum(ct)),lty=2)
text(1:10,0,LETTERSord)
text(0,1:10,LETTERSord)
for(i in 1:nlevels(b$x2)){
#levels(b$x2)=LETTERS[1:10]
b$x2=relevel(b$x2,LETTERSord[i])
p=summary(lm(y~x1+x2,data=b))$coefficients[-(1:2),4]
names(p)=substr(names(p),3,3)
P[LETTERSord[i],names(p)]=p
p=P[LETTERSord[i],]
idx=which(p>.05)
points(((1:10))[idx],rep(i,length(idx)),pch=1,cex=2)
idx=which(p>.1)
points(((1:10))[idx],rep(i,length(idx)),pch=19,cex=2)
}
在这里，我们得到以下信息：
It looks like we have our combined categories...
实际上，使用另一种策略是可能的。我们从某种程度上说“A”。然后，我们将它与所有非显著的不同级别合并起来。如果“B”不是其中之一，我们使用它作为新的参考。等
for(i in 1:nlevels(b$x2)){
if(LETTERS[i]%in%levels(b$x2)){
b$x2=relevel(b$x2,LETTERS[i])
p=summary(lm(y~x1+x2,data=b))$coefficients[-(1:2),4]
names(p)=substr(names(p),3,nchar(p))
idx=which(p>.05)
mix=c(LETTERS[i],names(p)[idx])
b$x2=recode(b$x2, paste("c('",paste(mix,collapse = "','"),"')='",paste(mix,collapse = "+"),"'",sep=""))
}}
最后的类别是：
table(b$x2)
A+I+H B+D+F C+G E J
46 82 35 23 14
使用以下回归输出：
summary(lm(y~x1+x2,data=b))
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.86407 0.03950 21.877 < 2e-16 ***
x1 1.99180 0.05323 37.417 < 2e-16 ***
x2B+D+F -0.21517 0.03699 -5.817 2.44e-08 ***
x2C+G -0.50545 0.04528 -11.164 < 2e-16 ***
x2E -0.83617 0.05128 -16.305 < 2e-16 ***
x2J -0.68398 0.06131 -11.156 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.2008 on 194 degrees of freedom
Multiple R-squared: 0.8975,Adjusted R-squared: 0.8948
F-statistic: 339.6 on 5 and 194 DF, p-value: < 2.2e-16
AIC(lm(y~x1+x2,data=b))
[1] -66.76939
BIC(lm(y~x1+x2,data=b))
[1] -43.68117
这和我们以前的团队是一致的。但实际上，如果我们改变顺序，我们可以得到不同的组合。例如，如果我们从“J”到“A”，而不是“A”到“J”，我们得到：
for(i in nlevels(b$x2):1){
#levels(b$x2)=LETTERS[1:10]
if(LETTERS[i]%in%levels(b$x2)){
b$x2=relevel(b$x2,LETTERS[i])
p=summary(lm(y~x1+x2,data=b))$coefficients[-(1:2),4]
names(p)=substr(names(p),3,nchar(p))
idx=which(p>.05)
mix=c(LETTERS[i],names(p)[idx])
b$x2=recode(b$x2, paste("c('",paste(mix,collapse = "','"),"')='",paste(mix,collapse = "+"),"'",sep=""))
}}
table(b$x2)
E G+C I+A+B+D+F+H J
23 35 128 14
这里有不同的信息标准：
AIC(lm(y~x1+x2,data=b))
[1] -36.61665
BIC(lm(y~x1+x2,data=b))
[1] -16.82675
我想有必要随机运行我们通过这些级别的顺序。最后，但同样重要的是，我们可以使用回归树。问题是，还有另一个解释变量可能会插入。所以我建议(1)拟合一个线性模型
为了计算残差，
(2)运行回归树，解释
使用范畴变量x_2(我解释了树是如何构建的，而解释变量是以前的职位):
library(rpart)
library(rpart.plot)
b$e=residuals(lm(y~x1,data=b))
arbre=rpart(e~x2,data=b)
prp(arbre,type=2,extra=1)
观察到叶子和我们得到的叶子是一样的。
arbre
n= 200
node), split, n, deviance, yval
* denotes terminal node
1) root 200 22.563500 7.771561e-18
2) x2=G,C,J,E 72 4.441495 -3.232525e-01
4) x2=J,E 37 1.553520 -4.578492e-01 *
5) x2=G,C 35 1.509068 -1.809646e-01 *
3) x2=I,A,H,F,B,D 128 6.366628 1.818295e-01
6) x2=F,B,D 82 2.983381 1.048246e-01 *
7) x2=I,A,H 46 2.030229 3.190993e-01 *
我想应该可以把所有这些放在一个R函数中，建议可能改进回归的水平组合。