"anova.rq" <-
function (object, ..., test = "Wald", joint = TRUE,
       score = "tau", se = "nid", R = 200, trim = NULL)

{
    if (length(list(object, ...)) > 1) {
	objects <- list(object, ...)
        return(anova.rqlist(objects, ..., test = test, joint = joint,
             score = score, se = se, R = R, trim = trim))
    }
    stop("Anova is only defined (yet) for lists of rq objects")
}
anova.rqs <- function(object, ..., se = "nid", joint = TRUE){
    class(object) <- "rq"
    m <- length(object$tau)
    z <- rep(list(object), m)
    for(i in 1:m){
        z[[i]]$coefficients <- object$coefficients[,i]
        z[[i]]$tau <- object$tau[i]
        z[[i]]$rho <- object$rho[i]
    }
    return(anova.rqlist(z, ..., se = se, joint = joint))
}

"anova.rqlist" <-
function (object, ..., test = "Wald", joint = TRUE, 
		score = "tau", se = "nid", R = 200, trim = NULL) 
{
    objects <- object
    responses <- as.character(lapply(objects, function(x) formula(x)[[2]]))
    sameresp <- responses == responses[1]
    if (!all(sameresp)) 
        stop("Models don't all have the same response variable")
    n <- length(objects[[1]]$y)
    models <- as.character(lapply(objects, function(x) formula(x)))
    nobjects <- length(objects)
    dimp <- lapply(objects, function(x) length(coef(x)))
    objects <- objects[order(-unlist(dimp))]
    mf <- model.frame(objects[[1]])
    models <- as.character(lapply(objects, function(x) formula(x)))
    taus <- unlist(lapply(objects, function(x) x$tau))
    if(is.matrix(coef(objects[[1]])))
    	names <- lapply(objects, function(x) dimnames(coef(x))[[1]])
    else
    	names <- lapply(objects, function(x) names(coef(x)))
    if (test == "Wald") 
        objects <- lapply(objects, function(x) summary(x,se=se,covariance = TRUE))
    sametaus <- taus == taus[[1]]
    if (all(sametaus)) {
        Tn <- rep(0, nobjects - 1)
        ndf <- Tn
        ddf <- Tn
        pvalue <- Tn
        topnote <- paste("Model ", format(1:nobjects), ": ", 
            models, sep = "", collapse = "\n")
        if (test == "anowar") {
            X1 <- model.matrix(objects[[1]],mf,contrasts=objects[[1]]$contrasts)
            y <- model.response(mf)
	    weights <- as.vector(model.weights(mf))
            tau <- taus[[1]]
            for (i in 2:nobjects) {
                if (!all(names[[i]] %in% names[[1]])) 
                  stop("Models aren't nested")
                mf <- model.frame(objects[[i]])
                X0 <- model.matrix(objects[[i]], mf,contrasts=objects[[i]]$contrasts)
                Htest <- rq.test.anowar(X0, X1, y, tau = tau, R = R)
                ndf[i - 1] <- Htest$ndf
                Tn[i - 1] <- Htest$Tn
                ddf[i - 1] <- Htest$ddf
                pvalue[i - 1] <- Htest$pvalue
            }
    	table <- data.frame(ndf, ddf, Tn, pvalue)
        }
        else if (test == "rank") {
            x1 <- model.matrix(objects[[1]],mf,contrasts=objects[[1]]$contrasts)
            y <- model.response(mf)
	    weights <- as.vector(model.weights(mf))
            for (i in 2:nobjects) {
                if (!all(names[[i]] %in% names[[1]])) 
                  stop("Models aren't nested")
                nullH <- is.na(match(names[[1]], names[[i]]))
		X1 <- as.matrix(x1[, nullH])
                mf <- model.frame(objects[[i]])
                X0 <- model.matrix(objects[[i]], mf,contrasts=objects[[i]]$contrasts)
		if(score == "tau") tau <- taus[[1]]
                Htest <- rq.test.rank(X0, X1, y, score = score, weights = weights, 
			tau = tau, trim = trim)
                ndf[i - 1] <- Htest$ndf
                Tn[i - 1] <- Htest$Tn
                ddf[i - 1] <- Htest$ddf
                pvalue[i - 1] <- Htest$pvalue
            }
    	table <- data.frame(ndf, ddf, Tn, pvalue)
        }
        else if (test == "Wald") {
            V <- lapply(objects, function(x) x$cov)
            coef <- lapply(objects, function(x) coef(x)[,1])
            for (i in 2:nobjects) {
                if (!all(names[[i]] %in% names[[1]])) 
                  stop("Models aren't nested")
                nullH <- is.na(match(names[[1]], names[[i]]))
                ndf[i - 1] <- sum(nullH)
                Tn[i - 1] <- t((coef[[1]])[nullH]) %*% solve((V[[1]])[nullH, 
                  nullH], (coef[[1]])[nullH])/ndf[i - 1]
                ddf[i - 1] <- n - length(names[[1]])
                pvalue[i - 1] <- 1 - pf(Tn[i - 1], ndf[i - 1], 
                  ddf[i - 1])
            }
    	table <- data.frame(ndf, ddf, Tn, pvalue)
        }
        else stop("test only defined for anowar, Wald and rank")
    }
    else {
        m <- length(taus)
	n <- NROW(objects[[1]]$residuals)
        for (i in 2:m) {
            if (!setequal(names[[i]], names[[1]])) 
                stop("Models with common tau don't have same X")
        }
        if (names[[1]][1] != "(Intercept)") 
            stop("Intercept required in common tau testing")
        Omega <- outer(taus, taus, pmin) - outer(taus, taus)
        J <- objects[[1]]$J
        p <- dim(J)[1]
        H <- array(unlist(lapply(objects, function(x) x$Hinv)), 
            c(p, p, m))
        H <- matrix(aperm(H, c(1, 3, 2)), p * m, p) %*% t(chol(J))
        W <- (H %*% t(H)) * (kronecker(Omega, outer(rep(1, p), 
            rep(1, p))))
        coef <- unlist(lapply(objects, function(x) coef(x)[,1]))
	if(joint){
        	D <- kronecker(diff(diag(m)), cbind(0, diag(p - 1)))
        	ndf <- (p - 1) * (m - 1)
        	Tn <- t(D %*% coef) %*% solve(D %*% W %*% t(D), D %*% coef)/ndf
        	ddf <- n * m - (p-1) * (m - 1)
        	pvalue <- 1 - pf(Tn, ndf, ddf)
        	nobjects <- 1
        	tnote1 <- paste("Model: ", models[[1]], "\n", sep = "")
        	tnote2 <- paste("Joint Test of Equality of Slopes: tau in { ", 
        	    	paste(taus, collapse = " "), " }\n")
        	topnote <- paste(tnote1, tnote2, sep = "")
    		table <- data.frame(ndf, ddf, Tn, pvalue)
		}
	else{
	   Tn <- pvalue <- rep(0,p-1)
	   ndf <- m-1
	   ddf <- n*m - (m-1)
	   for(i in 2:p){
		E <- matrix(0, 1, p)
		E[1,i] <- 1
		D <- kronecker(diff(diag(m)),E)
		Tn[i-1] <-  t(D %*% coef) %*% solve(D %*% W %*% t(D), D %*% coef)/ndf
        	pvalue[i-1] <- 1 - pf(Tn[i-1], ndf, ddf)
		}
      	   tnote1 <- paste("Model: ", models[[1]], "\n", sep = "")
       	   tnote2 <- paste("Tests of Equality of Distinct Slopes: tau in { ", 
            	paste(taus, collapse = " "), " }\n")
       	   topnote <- paste(tnote1, tnote2, sep = "")
    	   table <- data.frame(ndf, ddf, Tn, pvalue)
	   dimnames(table)[[1]] <- names[[1]][2:p]
          }
    }
    x <- list(table=table,topnote=topnote)
    class(x) <- "anova.rq"
    return(x)
}
"print.anova.rq" <-
function(x,...){
    table <- x$table
    topnote <- x$topnote
    dimnames(table)[[2]] <- c("Df", "Resid Df", "F value", "Pr(>F)")
    title <- "Quantile Regression Analysis of Deviance Table\n"
    a <- structure(table, heading = c(title, topnote), class = c("anova", 
        "data.frame"))
    print(a)
}
rq.test.anowar <- function(x0,x1,y,tau,R){
        if(!requireNamespace("logspline", quietly = TRUE))
	    stop("anowar test requires logspline package")
        n <- length(y)
        f0 <- rq(y ~ x0 - 1, tau = tau)
        f1 <- rq(y ~ x1 - 1, tau = tau)
        Rho <- function(u,tau) u * (tau - (u < 0))
        Mn <- sum(Rho(f0$resid,tau)) - sum(Rho(f1$resid,tau))
	if(!hasArg(R)) R <- 200
        W  <- matrix(rexp(n*R,1),n,R)
        B1 <- t(boot.rq.wxy(x1,y,W,tau))
        B0 <- t(boot.rq.wxy(x0,y,W,tau))
        R0 <- Rho(y - x0 %*% B0,tau)
        R1 <- Rho(y - x1 %*% B1,tau)
        W0 <- apply(R0*W, 2, sum) - apply(W * Rho(f0$resid,tau), 2, sum)
        W1 <- apply(R1*W, 2, sum) - apply(W * Rho(f1$resid,tau), 2, sum)
        RefDistn <- logspline::logspline(W0 - W1)
        pvalue <- 1 - logspline::plogspline(Mn,RefDistn)
        ndf <- ncol(x1) - ncol(x0)
        ddf <- n - ncol(x1)
        list(Tn = Mn, ndf=ndf, ddf=ddf, pvalue=pvalue)
        }
"rq.test.rank" <-
function (x0, x1, y, v = NULL, score = "wilcoxon", weights = NULL, 
    tau = 0.5, iid = TRUE, delta0 = rep(0, NCOL(x1)), omega = 1, trim = NULL, 
    pvalue = "F") 
{
    if (length(weights) > 0) {
        y <- weights * y
        x0 <- weights * x0
        x1 <- weights * x1
    }
    n <- length(y)
    if (!length(v) > 0) 
        if(score == "tau"){
           r <- rq.fit.br(x0,y,tau = tau)$dual - (1 - tau)
           r <- list(ranks = r, A2 =  tau * (1 - tau))
           }
        else{
           v <- rq(y ~ x0 - 1, tau = -1)
           r <- ranks(v, score, tau, trim)
           }
    else 
        r <- ranks(v, score, tau, trim)
    if(iid == FALSE && score == "tau"){
         h <- bandwidth.rq(tau, n, hs = TRUE)
         bhi <- rq.fit.br(x0, y, tau + h, ci = FALSE)
         bhi <- coefficients(bhi)
         blo <- rq.fit.br(x0, y, tau - h, ci = FALSE)
         blo <- coefficients(blo)
         dyhat <- x0 %*% (bhi - blo)
         if(any(dyhat <= 0)) {
              pfis <- (100 * sum(dyhat <= 0))/n
              warning(paste(pfis, "percent fis <=0"))
              }
         eps <- .Machine$double.eps^(2/3)
         f <- pmax(eps, (2 * h)/(dyhat - eps))
         x1hat <- resid(lm.wfit(x0, x1, w = f))
        }
    else
        x1hat <- as.matrix(qr.resid(qr(x0), x1))
    Tn <- as.matrix(t(x1hat) %*% r$ranks)
    Tn <- t(Tn) %*% solve(crossprod(x1hat)) %*% Tn/r$A2
    ncp <- (omega^2) * t(delta0) %*% (crossprod(x1hat)/n) %*% 
        delta0/r$A2
    ndf <- NCOL(x1)
    Tn <- Tn/ndf
    ddf <- length(y) - NCOL(x0) - NCOL(x1)
    if (pvalue == "F") 
        pvalue <- 1 - pf(Tn, ndf, ddf, ncp)
    else 
        pvalue <- 1 - pchisq(Tn*ndf, ndf, ncp)
    list(Tn = Tn, ndf = ndf, ddf = ddf, pvalue = pvalue)
}
"ranks" <- 
function (v, score = "wilcoxon", tau = 0.5, trim = NULL) 
{
    A2 <- 1
    if(length(trim) & !(score == "wilcoxon"))
        stop("trimming is only permitted with the wilcoxon score functions") 
    if (score == "wilcoxon") {
        phibar <- 1/2
        A2 <- 1/12
        u <- v$sol[1,]
        J <- length(u)
        d <- v$dsol
        if(length(trim)){
             if(length(trim) == 1) trim <- c(trim, 1 - trim)
             else if(length(trim) > 2) stop("Only 2 trimming proportions allowed")
             else {
                if(any(trim < 0) || any(trim > 1)) stop("trim must lie in [0,1]")
                trim <- sort(trim)
                klo <- findInterval(trim[1],u)
                khi <- findInterval(trim[2],u)
                wlo <- (trim[1] - u[klo])/(u[klo+1] - u[klo])
                dlo <- wlo * d[,klo] + (1-wlo) * d[,klo + 1]
                s <-  (klo+1):khi
                s <-  (klo+1):khi
                if(khi < J) {
                   whi <- (trim[2] - u[khi])/(u[khi+1] - u[khi])
                   dhi <- whi * d[,khi] + (1-whi) * d[,khi + 1]
                   u <- c(trim[1], u[s], trim[2])
                   d <- cbind(dlo, d[,s], dhi)
                   }
                else {
                   u <- c(trim[1],u[s])
                   d <- cbind(dlo,  d[,s])
                   }
                a <- trim[1]; b <- trim[2]
                phibar <- c <- ( a^2 + 2*b - b^2)/2
                A2 <- (2/3) *(a^3 - b^3) - c * a^2 - 2 * b * c + c^2 + (1+c) * b^2 
                }
            }
        D <- v$dsol
        J <- length(u)
        ranks <- c((0.5 * (d[, 2:J] + d[, 1:(J - 1)]) %*% diff(u)) - phibar)
        return(list(ranks = ranks, A2 = A2))
    }
    else if (score == "normal") {
        J <- ncol(v$sol)
        dt <- v$sol[1, 2:J] - v$sol[1, 1:(J - 1)]
        dphi <- c(0, dnorm(qnorm(v$sol[1, 2:(J - 1)])), 0)
        dphi <- diff(dphi)
        ranks <- as.vector((((v$dsol[, 2:J] - v$dsol[, 1:(J - 
            1)]))) %*% (dphi/dt))
        return(list(ranks = ranks, A2 = A2))
    }
    else if (score == "sign") {
        j.5 <- sum(v$sol[1, ] < 0.5)
        w <- (0.5 - v$sol[1, j.5])/(v$sol[1, j.5 + 1] - v$sol[1, 
            j.5])
        r <- w * v$dsol[, j.5 + 1] + (1 - w) * v$dsol[, j.5]
        ranks <- 2 * r - 1
        return(list(ranks = ranks, A2 = A2))
    }
    else if (score == "tau") {
        j.tau <- sum(v$sol[1, ] < tau)
        w <- (tau - v$sol[1, j.tau])/(v$sol[1, j.tau + 1] - v$sol[1, 
            j.tau])
        r <- w * v$dsol[, j.tau + 1] + (1 - w) * v$dsol[, j.tau]
        ranks <- r - (1 - tau)
        A2 <- tau * (1 - tau)
        return(list(ranks = ranks, A2 = A2))
    }
    else if (score == "normalscale") {
	J <- ncol(v$sol)
	taus <- v$sol[1,]
        dt <- taus[2:J] - taus[1:(J - 1)]
	Qt <- qnorm(taus[2:(J-1)])
        phi <- c(0, -dnorm(Qt)*Qt, 0)
        dphi <- diff(phi)
        ranks <- as.vector((((v$dsol[, 2:J] - v$dsol[, 1:(J -
            1)]))) %*% (dphi/dt))
        return(list(ranks = ranks, A2 = 2))
	}
    else if (score == "halfnormalscale") {
	J <- ncol(v$sol)
	taus <- v$sol[1,]
        dt <- taus[2:J] - taus[1:(J - 1)]
	Qt <- qnorm(taus[2:(J-1)])
        phi <- c(0, (taus[2:(J-1)] - dnorm(Qt)*Qt)*(taus[2:(J-1)] > .5), 1)
        dphi <- diff(phi)
	dphi[dphi > .5] <- 0 #Kludge to zap jump point.
        ranks <- (((v$dsol[, 2:J] - v$dsol[, 1:(J-1)]))) 
	ranks <- as.vector(ranks %*% (dphi/dt))
        return(list(ranks = ranks, A2 = 1.25))
	}
    else if(score == "lehmann"){
	J <- ncol(v$sol)
        taus <- v$sol[1,]
        dt <- taus[2:J] - taus[1:(J - 1)]
	taus <- taus[2:(J-1)]
	phi <- c(0, -(taus -1)*log(1-taus),0)
	dphi <- diff(phi)
        ranks <- (((v$dsol[, 2:J] - v$dsol[, 1:(J-1)])))
        ranks <- as.vector(ranks %*% (dphi/dt))
        return(list(ranks = ranks, A2 = 1))
	}
    else if (score == "interquartile") {
        j.25 <- sum(v$sol[1, ] < 0.25)
        w <- (0.25 - v$sol[1, j.25])/(v$sol[1, j.25 + 1] - v$sol[1, 
            j.25])
        r.25 <- w * v$dsol[, j.25 + 1] + (1 - w) * v$dsol[, j.25]
        j.75 <- sum(v$sol[1, ] < 0.75)
        w <- (0.75 - v$sol[1, j.75])/(v$sol[1, j.75 + 1] - v$sol[1, 
            j.75])
        r.75 <- w * v$dsol[, j.75 + 1] + (1 - w) * v$dsol[, j.75]
        ranks <- 0.5 + r.75 - r.25
        A2 <- 1/4
        return(list(ranks = ranks, A2 = A2))
    }
    else stop("invalid score function")
}