Bias in bagged decision tree estimates of probability (classification) As far as I understand it, random forests are as biased as any tree in the forest.  Bagging a non-random forest (using all available variables at any given split) is unbiased -- bias in random forests comes from sampling the independent variables.  
However, is 
$$
\widehat{pr}(y)_{bagging} = N^{-1}\displaystyle\sum_{trees}vote(y=1)
$$
a unbiased measure of the probability of $y=1$?
The following suggests not:
par(mfrow=c(1,1))
    N=3000
    x = rnorm(N)
    yt = pnorm(x+rnorm(N))
    y = as.factor(yt>.5)
    m = randomForest(y=y,x = cbind(x),mtry=1)
    yhat = predict(m,type='prob')[,2]
plot(yt,yhat)
abline(lm(yhat~yt),col='red',lwd=2)


Am I missing something?  I ask because I've recently heard some loose talk about bagging replacing nonparametric regression in contexts where causality is important.
 A: [Edited] Oh very important detail. I guess you intend to verify that the a RF estimate produce similar probabilistic predictions as an unbiased normal distribution estimate. But you define $y_t$ = pnorm(x + rnorm(N)), instead of $y_t$ = pnorm(x). Then the unbiased normal distribtion estimate is 'over fitted'.
Your RF model is also overfitted. The single trees of your ensemble are fully grown and overfitted (of course). To ensure the ensemble is not also overfitted, RF algorithm make use of both bootstrapping and random variable subspace to grow decorrelated trees. In this rather special case (1 feature + 50% noise component), it is difficult to grow adequately decorrelated trees. Then, default random forest hyper parameters are not favorable. You can either grow more robust trees by increasing terminal nodesizes or by decorrelating trees further by bootstrapping fewer observations for each tree.
Bias (understood as lack-of-fit) in random forest certainly also comes the trees themselves, if the data structure contains a curvature which is difficult for trees to fit (e.g. saddle points or spherical class borders in feature space).
library(randomForest)
par(mfrow=c(1,1))
#make data
N=3000
x = rnorm(N)
y.value = x + rnorm(N) #training and raw variable
yt = pnorm(x)
y = as.factor(y.value>0)

#tune sampsize
sampsizes = c(2,2*c(2:25)^2,1500,2000,2500)
oob.prob.cor.sampsize = sapply(sampsizes,function(i) {
  m = randomForest(y=y,x = cbind(x),mtry=1,sampsize = i)
  plot(x=yt,y=predict(m,type="prob")[,2],main=paste("sampsize=",i),xlab="yt",ylab="yhat") #movie time
  cor(y.value,predict(m,type="prob")[,2])
})
plot(sampsizes,oob.prob.cor.sampsize,log="x",type="p",main="sampsize")





#...or tune nodesize
nodesizes = c(1:5,6,8,10,15,20,25,50,100,250,500)
oob.prob.cor = sapply(nodesizes, function(i) {
  m =randomForest(y=y,x = cbind(x),nodesize=i)
  plot(x=yt,y=predict(m,type="prob")[,2],main=paste("sampsize=",i),xlab="yt",ylab="yhat") #movie time
  cor(y.value,predict(m,type="prob")[,2])
})
plot(nodesizes,oob.prob.cor,log="x",type="p",main="optimal nodesize")

.

