Today, I was trying to implement a weighted bootstrap method in R and ran into a problem involving the 'weights' argument of the 'rq' function. Essentially, when putting it inside a function, it systematically fails to find the object I assign to that argument except in bizarre cases -- e.g., it works if the object exists in the global environment.
I have a minimal working example of the issue:
# Clear everything
rm(list=ls())
# Load quantile regression library
library(quantreg)
# Set seed
SEED = 10983
set.seed(SEED)
# Simulate data
N = 500
K = 2
# Linear gaussian model: Y = XB + e
B = array(1, dim=c(K,1))
X = array( rnorm(N*K), dim=c(N,K) )
e = array( rnorm(N), dim=c(N,1))
Y = X%*%B + e
# Check everything is fine: Eg.
plot(X[,1], Y)
# Turn into a dataframe
df = cbind.data.frame(Y, X)
names(df) = c( 'Y', paste('X', 1:K, sep=''))
# Estimation
eqn = paste('Y ~ 1 +', paste(names(df)[-1], collapse=' + '), collapse=' ')
mdl = rq(formula=eqn, data=df, tau=-1)
# Estimation with weights
Zs = array(rexp(N), dim=c(N,1))
ws = Zs / sum(Zs)
# Exponential weights
mdl2 = rq(formula=eqn, data=df, tau=-1, weights=N*ws)
# Equal weights
mdl3 = rq(formula=eqn, data=df, tau=-1, weights=array(1, dim=c(N,1)))
print(mdl$sol[,100])
print(mdl2$sol[,100])
print(mdl3$sol[,100])
# Now, let's do it in a function
estimate = function(Y, X, W){
# Turn into a dataframe
df = cbind.data.frame(Y, X)
names(df) = c( 'Y', paste('X', 1:K, sep=''))
# Estimation
eqn = paste('Y ~ 1 +', paste(names(df)[-1], collapse=' + '), collapse=' ')
mdl = rq(formula=eqn, data=df, tau=-1, weights=W)
return(mdl$sol[,100])
}
#### AND NOW IT DOESN'T WORK!
estimate(Y,X,ws*N)
### BUT THIS DOES
W = ws*N
estimate(Y,X,ws*N)
### AND THIS RUNS, BUT ODDLY USES THE SAME WEIGHT 'W' AS IN THE GLOBAL ENV.
estimate = function(Y, X){
# Turn into a dataframe
df = cbind.data.frame(Y, X)
names(df) = c( 'Y', paste('X', 1:K, sep=''))
# Get weights
Zs = array(rexp(N), dim=c(N,1))
W = Zs / sum(Zs)
# Estimation
eqn = paste('Y ~ 1 +', paste(names(df)[-1], collapse=' + '), collapse=' ')
mdl = rq(formula=eqn, data=df, tau=-1, weights=W)
return(mdl$sol[,100])
}
estimate(Y,X)
This is somewhat of a problem. I have an ugly work around that can do the trick, but I am wondering if someone has a better idea. My work around is simple enough: the objective is piece-wise linear, so I can just 'force' the weights by using 'Y ~ 0 + X1 + X2 + X3' with X3 being a column of ones and multiply everything by the random weight vector (i.e., data=df*W). (You need to forcibly include the constant because, otherwise, the constant's score isn't re-weighted in the optimization).
Any idea? I mean, it's nice enough being able to 'make it work,' but that might not be the best solution and I am open to suggestions. It's also a public service to share this seemingly unknown problem.