[R-sig-eco] Anova in quantile regressions

[Tomas]

Dear list members

I'm using the "quantreq" package to perform quantile regressions on leaf 
traits (photosynthesis and leaf nutrients). The rq function is great for 
this because its clear that there are other things that limits 
photosynthesis)
I'm interested in determining the 90% quantile regression for the 
Photosynthesis vs. leaf nitrogen for two independent data sets. I was 
able to do that as bellow:

 > head(data1)
   TRIP      SLA Vcmax.mass  N.mass P.mass
1  CAM 42.10000 0.24216835 11.2530 0.6800
2  CAM 42.40000 0.09373972  4.7430 2.4800
3  CAM 44.80000 0.19999955  6.6360 3.1500
4  CAM 46.60000 0.11043004  5.1780 0.7800
5  CAM 47.00000 0.24968015 11.4540 0.7200
6   WA 47.15013 0.41529054 11.0038 0.7901

# Here I divide my data set into 2 subsets (samples from Cameroon and
# from West Africa). This was a tip I got from Nikhil Kaza on I previous
# post I placed on the R-Help list.

 > data_sub1 <- data1[data1$TRIP=="CAM",]
 > data_sub2 <- data1[data1$TRIP=="WA",]

 > ab1 <- rq(Vcmax.mass~N.mass, data=data_sub1, tau=0.9)
 > ab2 <- rq(Vcmax.mass~N.mass, data=data_sub2, tau=0.9)

 > summary(ab1)

Call: rq(formula = Vcmax.mass ~ N.mass, tau = 0.9, data = data_sub1)

tau: [1] 0.9

Coefficients:
             coefficients lower bd upper bd
(Intercept) -0.01013     -0.02926  0.07401
N.mass       0.03181      0.02808  0.03444
Warning message:
In rq.fit.br(x, y, tau = tau, ci = TRUE, ...) : Solution may be nonunique

 > summary(ab2)

Call: rq(formula = Vcmax.mass ~ N.mass, tau = 0.9, data = data_sub2)

tau: [1] 0.9

Coefficients:
             coefficients lower bd upper bd
(Intercept) 0.28323      0.11667  0.37078
N.mass      0.02496      0.02226  0.03583


So far so good, but when I try to do an ANOVA I get an error message:

 > anova(ab1, ab2)
Error in drop(.Call("La_dgesv", a, as.matrix(b), tol, PACKAGE = "base")) 
: 'a' is 0-diml
 > traceback()
8: .Call("La_dgesv", a, as.matrix(b), tol, PACKAGE = "base")
7: drop(.Call("La_dgesv", a, as.matrix(b), tol, PACKAGE = "base"))
6: solve.default((V[[1]])[nullH, nullH], (coef[[1]])[nullH])
5: solve((V[[1]])[nullH, nullH], (coef[[1]])[nullH])
4: solve((V[[1]])[nullH, nullH], (coef[[1]])[nullH])
3: anova.rqlist(object, ...)
2: anova.rq(ab1, ab2)
1: anova(ab1, ab2)

I have no problem comparing slopes between different quantiles for the 
complete data set (or any of the sub data sets):

 > ab.5<-rq(Vcmax.mass~N.mass, data=data1, tau=.5)
 > ab.9<-rq(Vcmax.mass~N.mass, data=data1, tau=.9)
 > anova(ab.5,ab.9)
trace: anova(ab.5, ab.9)
Quantile Regression Analysis of Deviance Table

Model: Vcmax.mass ~ N.mass
Joint Test of Equality of Slopes: tau in {  0.5 0.9  }

     Df Resid Df F value    Pr(>F)
1    1     1061  12.454 0.0004349 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Any ideas would be most appreciated.
Tomas Ferreira Domingues

-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.