0
点赞
收藏
分享

微信扫一扫

R中因子水平的自动组合

耶也夜 2023-06-15 阅读 73


每次我们在应用计量经济学课程中面对实际应用时,我们都必须处理分类变量。同样的问题也发生在学生身上:我们怎样才能自动地结合因素水平呢?有简单的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)


R中因子水平的自动组合_5e

线性回归的输出产生以下预测:

 

R中因子水平的自动组合_线性回归_02

 

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)}


R中因子水平的自动组合_5e_03

我们很高兴看到它是对称的:如果“H”应该与“i”结合,那么“i”也应该与“H”结合。

在这里,黑点与10%p值有关,白色点与5%p值有关.这张图其实很难读.。事实上,这让我们想起Bertin(1967年).

R中因子水平的自动组合_线性回归_04

在这里,我们可以手动预定义一些排序(下面我们将看到它是如何自动化的):

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)
 
}


在这里,我们得到以下信息:

R中因子水平的自动组合_5e_05

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)拟合一个线性模型

R中因子水平的自动组合_5e_06

为了计算残差,

R中因子水平的自动组合_5e_07

(2)运行回归树,解释

R中因子水平的自动组合_5e_08

使用范畴变量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)


R中因子水平的自动组合_线性回归_09

观察到叶子和我们得到的叶子是一样的。

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函数中,建议可能改进回归的水平组合。

举报

相关推荐

11-R因子

R语言:因子分析 factor analysis

0 条评论