[R] odd behaviour in quantreg::rq

[roger koenker]

Admittedly this seemed quite peculiar....  but if you look at the  
entrails
of the following code you will see that with the weights the first and
second levels of your x$method variable have the same (weighted) median
so the contrast that you are estimating SHOULD be zero.  Perhaps
there is something fishy about the data construction that would have
allowed us to anticipate this?  Regarding the "fn" option, and the
non-uniqueness warning,  this is covered in the (admittedly obscure)
faq on quantile regression available at:

	http://www.econ.uiuc.edu/~roger/research/rq/FAQ

# example:
library(quantreg)

# load data
x <- read.csv(url('http://169.237.35.250/~dylan/temp/test.csv'))

# with weights
summary(rq(sand ~ method, data=x, weights=area_fraction, tau=0.5),  
se='ker')

#Reproduction with more convenient notation:

X0 <- model.matrix(~method, data = x)
y <- x$sand
w <- x$area_fraction
f0 <- summary(rq(y ~ X0 - 1, weights = w),se = "ker")

#Second reproduction with orthogonal design:

X1 <- model.matrix(~method - 1, data = x)
f1 <- summary(rq(y ~ X1 - 1, weights = w),se = "ker")

#Comparing f0 and f1 we see that they are consistent!!  How can that  
be??
#Since the columns of X1 are orthogonal estimation of the 3 parameters  
are separable
#so we can check to see whether the univariate weighted medians are  
reproducible.

s1 <- X1[,1] == 1
s2 <- X1[,2] == 1
g1 <- rq(y[s1] ~ X1[s1,1] - 1, weights = w[s1])
g2 <- rq(y[s2] ~ X1[s2,2] - 1, weights = w[s2])

#Now looking at the g1 and g2 objects we see that they are equal and  
agree with f1.


url:    www.econ.uiuc.edu/~roger            Roger Koenker
email    rkoenker at uiuc.edu            Department of Economics
vox:     217-333-4558                University of Illinois
fax:       217-244-6678                Urbana, IL 61801



On Jun 30, 2009, at 3:54 PM, Dylan Beaudette wrote:

> Hi,
>
> I am trying to use quantile regression to perform weighted- 
> comparisons of the
> median across groups. This works most of the time, however I am  
> seeing some
> odd output in summary(rq()):
>
> Call: rq(formula = sand ~ method, tau = 0.5, data = x, weights =
> area_fraction)
> Coefficients:
>                   Value    Std. Error t value  Pr(>|t|)
> (Intercept)        45.44262  3.64706   12.46007  0.00000
> methodmukey-HRU     0.00000  4.67115    0.00000  1.00000
> 				  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
>
> When I do not include the weights, I get something a little closer  
> to a
> weighted comparison of means, along with an error message:
>
> Call: rq(formula = sand ~ method, tau = 0.5, data = x)
> Coefficients:
>                   Value    Std. Error t value  Pr(>|t|)
> (Intercept)        44.91579  2.46341   18.23318  0.00000
> methodmukey-HRU     9.57601  9.29348    1.03040  0.30380
> Warning message:
> In rq.fit.br(x, y, tau = tau, ...) : Solution may be nonunique
>
>
> I have noticed that the error message goes away when specifying  
> method='fn' to
> rq(). An example is below. Could this have something to do with  
> replication
> in the data?
>
>
> # example:
> library(quantreg)
>
> # load data
> x <- read.csv(url('http://169.237.35.250/~dylan/temp/test.csv'))
>
> # with weights
> summary(rq(sand ~ method, data=x, weights=area_fraction, tau=0.5),  
> se='ker')
>
> # without weights
> # note error message
> summary(rq(sand ~ method, data=x, tau=0.5), se='ker')
>
> # without weights, no error message
> summary(rq(sand ~ method, data=x, tau=0.5, method='fn'), se='ker')
>
> -- 
> Dylan Beaudette
> Soil Resource Laboratory
> http://casoilresource.lawr.ucdavis.edu/
> University of California at Davis
> 530.754.7341
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.