[R-sig-eco] Anova in quantile regressions

[Gavin Simpson]

On Tue, 2010-03-02 at 14:54 +0000, Tomas wrote:
> 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",]

It would be more R-ish to use the subset argument if rq's formula method
allows (which it does)

 > ab1 <- rq(Vcmax.mass~N.mass, data=data1, 
             subset = TRIP == "CAM", tau=0.9)
 > ab2 <- rq(Vcmax.mass~N.mass, data=data1, 
             subset = TRIP == "CAM", tau=0.9)

This keeps all your data together etc...

> 
>  > ab1 <- rq(Vcmax.mass~N.mass, data=data_sub1, tau=0.9)
>  > ab2 <- rq(Vcmax.mass~N.mass, data=data_sub2, tau=0.9)
> 
<snip />
> 
> So far so good, but when I try to do an ANOVA I get an error message:

That just doesn't make any sense. Those models aren't in any way shape
or form "nested" so you can't compare them that way. Whatever gave you
the impression that you could do this with *any* anova method in R (at
least that which I have come across)??

Read the Warning section in ?anova.rq which explicitly states that the
models have be be fitted to the *same* data.

The two models need different covariates but be fitted to the same data
to use anova this way. For example:

ab1 <- rq(Vcmax.mass~N.mass, data=data1, tau=0.9)
ab2 <- rq(Vcmax.mass~N.mass + foo, data=data1, tau=0.9)

anova(ab1, ab2)

would be fine, because the model described by ab1 can be reached by
setting the coefficient for covariate 'foo' in ab2 to 0. This is
nesting. AFAIK, this is the only model comparison that can be done with
anova methods in R.

> 
>  > 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):

That is because the package author wrote code to handle this particular
feature and documented it. See ?anova.rq

> 
>  > 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

What is it you are trying to do, compare the slopes between the two
trips? If so, the usual way would be to include TRIP as an extra
covariate *in* the model:

## If I have this right myself,
## this model assumes same slope but different intercepts
m1 <- rq(Vcmax.mass ~ N.mass + TRIP, data=data1, tau=0.9)

## whilst this model allows the slopes to vary
m2 <- rq(Vcmax.mass ~ N.mass * TRIP, data=data1, tau=0.9)
## equivalent to
m2 <- rq(Vcmax.mass ~ N.mass + TRIP + N.mass:TRIP, data=data1, tau=0.9)

## we can then test the full nesting of models using:
m0 <- rq(Vcmax.mass ~ N.mass, data=data1, tau=0.9)
anova(m0, m1, m2)

Of course, I haven't tested any of this, but the above should work given
my reading of the quantreg help (and that I understood it correctly, of
course!).

I strongly recommend you read ?anova.rq and the related references etc
as the author is at pains to stress the need to understand what
hypotheses can be tested using the anova method for rq objects.

HTH

G
-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 Dr. Gavin Simpson             [t] +44 (0)20 7679 0522
 ECRC, UCL Geography,          [f] +44 (0)20 7679 0565
 Pearson Building,             [e] gavin.simpsonATNOSPAMucl.ac.uk
 Gower Street, London          [w] http://www.ucl.ac.uk/~ucfagls/
 UK. WC1E 6BT.                 [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%