我正在研究具有大约 4500 个变量的二元类随机森林。其中许多变量是高度相关的,其中一些只是原始变量的分位数。我不太确定应用 PCA 进行降维是否明智。这会提高模型性能吗?
我希望能够知道哪些变量对我的模型更重要,但如果我使用 PCA,我只能知道哪些 PC 更重要。
提前谢谢了。
我正在研究具有大约 4500 个变量的二元类随机森林。其中许多变量是高度相关的,其中一些只是原始变量的分位数。我不太确定应用 PCA 进行降维是否明智。这会提高模型性能吗?
我希望能够知道哪些变量对我的模型更重要,但如果我使用 PCA,我只能知道哪些 PC 更重要。
提前谢谢了。
我的经验是,RF 之前的 PCA 并不是一个很大的优势(如果有的话)。主成分回归 (PCR) 是,例如,当 PCA 在 OLS 线性回归之前帮助正则化训练特征时,这对于稀疏数据集是非常需要的。由于 RF 本身已经在不假设线性的情况下执行了良好/公平的正则化,因此它不一定是优势。也就是说,我发现自己在两周前为 R 编写了一个 PCA-RF 包装器。该代码包括一个由 100 个特征组成的数据集的模拟数据集,该数据集仅包含 5 个真正的线性分量。在这种情况下,使用 PCA 进行预过滤实际上是一个小优势。代码是无缝实现的,因此每个 RF 参数都可以简单地传递给 RF。加载向量保存在 model_fit 中以在预测期间使用。
@我希望能够知道哪些变量对我的模型更重要,但如果我使用 PCA,我只能知道哪些 PC 更重要。
简单的方法是在没有 PCA 的情况下运行并获得可变重要性,并期望找到类似 PCA-RF 的东西。
繁琐的方法是用你自己的可变重要性代码将 PCA-RF 包装在一个新的 bagging 方案中。可以在 50-100 行左右完成。
PCA-RF 的源代码建议:
#wrap PCA around randomForest, forward any other arguments to randomForest
#define as new S3 model class
train_PCA_RF = function(x,y,ncomp=5,...) {
f.args=as.list(match.call()[-1])
pca_obj = princomp(x)
rf_obj = do.call(randomForest,c(alist(x=pca_obj$scores[,1:ncomp]),f.args[-1]))
out=mget(ls())
class(out) = "PCA_RF"
return(out)
}
#print method
print.PCA_RF = function(object) print(object$rf_obj)
#predict method
predict.PCA_RF = function(object,Xtest=NULL,...) {
print("predicting PCA_RF")
f.args=as.list(match.call()[-1])
if(is.null(f.args$Xtest)) stop("cannot predict without newdata parameter")
sXtest = predict(object$pca_obj,Xtest) #scale Xtest as Xtrain was scaled before
return(do.call(predict,c(alist(object = object$rf_obj, #class(x)="randomForest" invokes method predict.randomForest
newdata = sXtest), #newdata input, see help(predict.randomForest)
f.args[-1:-2]))) #any other parameters are passed to predict.randomForest
}
#testTrain predict #
make.component.data = function(
inter.component.variance = .9,
n.real.components = 5,
nVar.per.component = 20,
nObs=600,
noise.factor=.2,
hidden.function = function(x) apply(x,1,mean),
plot_PCA =T
){
Sigma=matrix(inter.component.variance,
ncol=nVar.per.component,
nrow=nVar.per.component)
diag(Sigma) = 1
x = do.call(cbind,replicate(n = n.real.components,
expr = {mvrnorm(n=nObs,
mu=rep(0,nVar.per.component),
Sigma=Sigma)},
simplify = FALSE)
)
if(plot_PCA) plot(prcomp(x,center=T,.scale=T))
y = hidden.function(x)
ynoised = y + rnorm(nObs,sd=sd(y)) * noise.factor
out = list(x=x,y=ynoised)
pars = ls()[!ls() %in% c("x","y","Sigma")]
attr(out,"pars") = mget(pars) #attach all pars as attributes
return(out)
}
运行代码示例:
#start script------------------------------
#source above from separate script
#test
library(MASS)
library(randomForest)
Data = make.component.data(nObs=600)#plots PC variance
train = list(x=Data$x[ 1:300,],y=Data$y[1:300])
test = list(x=Data$x[301:600,],y=Data$y[301:600])
rf = randomForest (train$x, train$y,ntree =50) #regular RF
rf2 = train_PCA_RF(train$x, train$y,ntree= 50,ncomp=12)
rf
rf2
pred_rf = predict(rf ,test$x)
pred_rf2 = predict(rf2,test$x)
cat("rf, R^2:",cor(test$y,pred_rf )^2,"PCA_RF, R^2", cor(test$y,pred_rf2)^2)
cor(test$y,predict(rf ,test$x))^2
cor(test$y,predict(rf2,test$x))^2
pairs(list(trueY = test$y,
native_rf = pred_rf,
PCA_RF = pred_rf2)
)