«График важности переменной» и выбор переменной

Классификация деревьев хороша. Они предоставляют интересную альтернативу логистической регрессии. Я начал включать их в свои курсы, может быть, семь или восемь лет назад. Хороший вопрос (как получить оптимальное разбиение), хорошая алгоритмическая процедура (хитрость разделения по одной переменной и только по одной в каждом узле, а затем двигаться вперед, а не назад) и визуальный вывод просто идеально (с этой древовидной структурой). Но прогноз может быть довольно плохим. Производительность этого алгоритма едва ли может конкурировать с (точно определенной) логистической регрессией.

Затем я обнаружил леса (см  . Страницу Лео Бреймана  для подробной презентации). Будучи большим поклонником процедур Boostrap, мне понравилась идея. В регрессионных моделях я обычно упоминаю начальную загрузку, чтобы избежать асимптотических приближений: мы загружаем строки (наблюдения). В случае случайного леса я должен признать, что идея случайного выбора набора возможных переменных в каждом узле очень умна. Производительность намного лучше, но интерпретация обычно сложнее. И то, что я люблю, когда есть много ковариаций: график переменной важности. Это то, что мы вряд ли можем получить с эконометрическими моделями (пожалуйста, дайте мне знать, если я ошибаюсь).

Чтобы проиллюстрировать это, давайте сгенерируем большой набор данных. Не обязательно огромный, но большой, так что нам действительно нужно выбирать переменные. Поскольку более интересно, если у нас есть, возможно, коррелированные переменные, нам нужна ковариационная матрица. В R есть хороший пакет для случайной генерации ковариационных матриц.

> set.seed(1)
> n=500
> library(clusterGeneration)
> library(mnormt)
> S=genPositiveDefMat("eigen",dim=15)
> S=genPositiveDefMat("unifcorrmat",dim=15)
> X=rmnorm(n,varcov=S$Sigma)
> library(corrplot)
> corrplot(cor(X), order = "hclust")

> Score=-1+X[,1]/sd(X[,1])+X[,2]/sd(X[,2])+
+ 0.5*X[,3]/sd(X[,3])-0.25*X[,4]/sd(X[,4])
> P=exp(Score)/(1+exp(Score))
> Y=rbinom(n,size=1,prob=P)
> df=data.frame(Y,X)
> allX=paste("X",1:ncol(X),sep="")
> names(df)=c("Y",allX)

> require(randomForest)
> fit=randomForest(factor(Y)~., data=df)

> (VI_F=importance(fit))
    MeanDecreaseGini
X1          31.14309
X2          31.78810
X3          20.95285
X4          13.52398
X5          13.54137
X6          10.53621
X7          10.96553
X8          15.79248
X9          14.19013
X10         10.02330
X11         11.46241
X12         11.36008
X13         10.82368
X14         10.17462
X15         10.45530

> library(caret)
> varImp(fit)
     Overall
X1  31.14309
X2  31.78810
X3  20.95285
X4  13.52398
X5  13.54137
X6  10.53621
X7  10.96553
X8  15.79248
X9  14.19013
X10 10.02330
X11 11.46241
X12 11.36008
X13 10.82368
X14 10.17462
X15 10.45530

> varImpPlot(fit,type=2)

> importanceOrder=order(-fit$importance)
> names=rownames(fit$importance)[importanceOrder][1:15]
> par(mfrow=c(5, 3), xpd=NA)
> for (name in names)
+   partialPlot(fit, df, eval(name), main=name, xlab=name,ylim=c(-.2,.9))

> library(rpart)
> fit=rpart(factor(Y)~., df)
> plot(fit)
> text(fit)

> tmp=rownames(fit$splits)
> allVars=colnames(attributes(fit$terms)$factors)  
> rownames(fit$splits)=1:nrow(fit$splits)
> splits=data.frame(fit$splits)
> splits$var=tmp
> splits$type=""
> frame=as.data.frame(fit$frame)
> index=0
> for(i in 1:nrow(frame)){
+   if(frame$var[i] != "<leaf>"){
+   index=index + 1
+   splits$type[index]="primary"
+   if(frame$ncompete[i] > 0){
+   for(j in 1:frame$ncompete[i]){
+   index=index + 1
+   splits$type[index]="competing"}}
+   if(frame$nsurrogate[i] > 0){
+   for(j in 1:frame$nsurrogate[i]){
+      index=index + 1
+      splits$type[index]="surrogate"}}}}
> splits$var=factor(as.character(splits$var))
> splits=subset(splits, type != "surrogate")
> out=aggregate(splits$improve,
+     list(Variable = splits$var),
+     sum, na.rm = TRUE)
> allVars=colnames(attributes(fit$terms)$factors)
>  if(!all(allVars %in% out$Variable)){
+   missingVars=allVars[!(allVars %in% out$Variable)]
+   zeros=data.frame(x = rep(0, length(missingVars)), Variable = missingVars)
+   out=rbind(out, zeros)}
> out2=data.frame(Overall = out$x)
> rownames(out2)=out$Variable
> out2
      Overall
X1  51.024692
X10  4.328443
X11 19.087255
X12 10.399549
X13 15.248933
X15  9.989834
X2  68.758329
X3  41.986055
X4  15.211913
X5  18.247668
X7  18.857998
X8  43.318540
X9  30.299429
X6   0.000000
X14  0.000000

> VI_T=out2
> barplot(unlist(VI_T/sum(VI_T)),names.arg=1:15)

> barplot(t(VI_F/sum(VI_F)))