# Using a broken stick model with random forests

I'm using random forests in R for classification problems (via package randomForest). I have two objectives: to get the classification accuracy, and to find out which variables contribute most to the discrimination. My question is about the latter: how to find a consistent way to judge which variables are important.

The importance scores can be extracted either raw or divided by the SD, but in either case they are correlated with the numbers of trees in the forest and with the number of observations in the dataset. It occurred to me that the problem was analogous to choosing the number of components from a PCA, so I have been experimenting with using a broken stick model. Here is my example:

library(randomForest)
library(vegan)
data("imports85")

#Run a random forest
set.seed(47)
rf<-randomForest(imports85[,10:18], imports85\$bodyStyle,
importance = T, ntree=5000, na.action = na.exclude)
rf

imp<-importance(rf)[,"MeanDecreaseAccuracy"]

#Broken stick model where 'total variance' = sum of importance values > 0
stick<-bstick(length(imp), sum(imp[imp>0]))
stick<-stick[order(stick, decreasing = T)]

#Plot the predicted importance and overlay the real data
scatter.smooth(seq(9,1)~stick, type = "n",
xlim = c(min(imp), max(imp)), xlab = "Importance", ylab = "Variable")
points(seq(9,1)~sort(imp, decreasing = T))

#Count how many variables had greater than predicted importance
sum(sort(imp, decreasing = T) > stick)


My question is whether this is a valid approach. It makes sense to me, as the broken stick model represents a null model to compare which variables are more important than would be expected by chance. Negative values are excluded when calculating 'total importance', as otherwise the broken stick model would underestimate the value of the most important variables. I would be very grateful for insights/expertise (please be gentle - I'm an ecologist rather than a mathematician)!