带置信区间的画图问题

epoh前辈，anova.rq函数是和rq函数一起使用的
ie：
data(barro)
fit0 <- rq(y.net ~  lgdp2 + fse2 + gedy2 , data = barro)
fit1 <- rq(y.net ~  lgdp2 + fse2 + gedy2 + Iy2 + gcony2, data = barro)
fit2 <- rq(y.net ~  lgdp2 + fse2 + gedy2 + Iy2 + gcony2, data = barro,tau=.75)
fit3 <- rq(y.net ~  lgdp2 + fse2 + gedy2 + Iy2 + gcony2, data = barro,tau=.25)
anova(fit1,fit0)
anova(fit1,fit2,fit3)
anova(fit1,fit2,fit3,joint=FALSE)

我要的是anova(fit1,fit2,fit3,joint=FALSE)这个实现结果。我看了anova.rq内部函数，呵呵，好复杂啊
内部程序：
> anova.rq
function (object, ...)
{
    if (length(list(object, ...)) > 1) {
        return(anova.rqlist(object, ...))
    }
    stop("Anova is only defined (yet) for lists of rq objects")
}
> anova.rqlist
function (object, ..., test = "Wald", joint = TRUE, score = "wilcoxon",
    R = 200, trim = NULL)
{
    objects <- list(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 = "nid",
            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 <- 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 <- 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)
        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)
}

呵呵，看来rq.sfn.fit这个函数还是存在很大缺陷，主要还是用rq（）函数

哈哈!老兄,
应该说package 更新的速度,
不及你研究的"深度".
我稍瞄一下:
The standard errors (in parenthesis) are obtained after 1000 panel-bootstrap replications
我抽空,看看能否有解.
你有没这组DATA?

epoh前辈，见笑了，呵呵。我只是对方法很感兴趣。
其实这些估计分位的函数，包括quantreg包都是roger教授编写的。他专门写了一个rq.fit.panel（不属于R内部函数），其实就是调用rq.fit.sfn函数来估计面板数据分位，而据他主页上介绍rqss（）函数和rq（）函数也是可以实现这一功能了。不过前段时间您已经帮忙看了，他只能一次估计一个分位点，而不能同时估计多个分位点。

这篇文章是roger教授的学生写的一篇，这位老师一般都不共享数据的，所以我也没有数据。
额，标准误是通过panel-bootstrap产生的，这个您已经实现了。

哈哈，可惜epoh老师不是专门学习quantile regression的，要不我就可以向您请教疑问了。
我现在还堆积着很多问题没有得到解答，看的深度不够，基础不扎实啊。

epoh 发表于 2011-4-4 13:48
library(SparseM)
library(quantreg)
library(Ecdat)

前辈，您好，在论坛里看到您是R语言的高手，我问一下，为什么欲行面板分位数回归时总是出问题：
> library(SparseM)
> library(quantreg)
> library(Ecdat)
> data(Grunfeld)
> #panel data Grunfeld
> #n=10, T=20, N=200
> x1=Grunfeld$value
> x2=Grunfeld$capital
> X=cbind(x1,x2)
> y=Grunfeld$inv
> s= rep(1:10,rep(20,10))
> lambda=1
> w=c(0.1,0.2,0.3,0.4,0.5,0.4,0.3,0.2,0.1)
> taus=c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9)
> ntaus=length(taus)
> K <- length(w)
> X <- as.matrix(X)
> p <- ncol(X)
> n <- length(levels(as.factor(s)))
> N <- length(y)
> Z <- as.matrix.csr(model.matrix(~as.factor(s)-1))
> Fidelity <- cbind(as(w,"matrix.diag.csr") %x% X,w %x% Z)  
错误于cbind(deparse.level, ...) : args have differing numbers of rows
此外: 警告信息：
In as.matrix.csr(x, nrow = nrow, ncol = ncol, eps = eps) :
  ncol*nrow indivisable by length(x)
> Penalty <- cbind(as.matrix.csr(0,n,K*p),lambda*as(n,"matrix.diag.csr"))  
> D <- rbind(Fidelity,Penalty)   
错误于rbind(Fidelity, Penalty) : 找不到对象'Fidelity'
> yy <- c(w %x% y,rep(0,n))
> a <- c((w*(1-taus)) %x% (t(X)%*%rep(1,N)),
+        sum(w*(1-taus)) * (t(Z) %*% rep(1,N)) + lambda * rep(1,n))  
> rq.fit.sfn(D,yy,tau=taus,rhs=a)
请您指教
我计算过，一个是200*18的矩阵，一个是200*90的矩阵，为什么会提示说行不同呢？？请您指点一下

sun1018 发表于 2012-4-9 17:43
前辈，您好，在论坛里看到您是R语言的高手，我问一下，为什么欲行面板分位数回归时总是出问题：
> libra ...

嘿嘿,一年了,有点忘了
不知ywh兄,还记不记得?

epoh 发表于 2012-4-9 22:17
嘿嘿,一年了,有点忘了
不知ywh兄,还记不记得?

呵呵，刚才仔细看了一下， 是W的定义有问题

epoh 发表于 2012-4-9 22:17
嘿嘿,一年了,有点忘了
不知ywh兄,还记不记得?

呵呵，epoh老师，楼主自己解决了。

sun1018 发表于 2012-4-9 22:53
呵呵，刚才仔细看了一下， 是W的定义有问题

DB file 是MicroTSP 的数据文件，在TSP中用fetch得到。
epoh老师，DB结尾的数据文件在TSP中如何打开呢？
好像是一个网页数据。

ywh19860616 发表于 2012-4-10 16:46
DB file 是MicroTSP 的数据文件，在TSP中用fetch得到。
epoh老师，DB结尾的数据文件在TSP中如何打开呢？ ...

MicroTSP Databank Format
  http://econpy.googlecode.com/svn/trunk/pytrix/opendatabank.txt

?A simple example
?x.db  (x series:1,3,5,7,9)
?Undated Series
?Line 1 : starting index
?Line 2 : ending index
1
5
1
3
5
7
9
? save as x.db
??????
??????
?open x.db
fetch x;
print x;
? save as x.tsp

epoh 发表于 2012-4-10 21:07
MicroTSP Databank Format
  http://econpy.googlecode.com/svn/trunk/pytrix/opendatabank.txt

呵呵，谢谢epoh老师，可以显示
我开始用了dbprint命令，一直显示错误。