带置信区间的画图问题

epoh，谢谢您了
呵呵，感谢楼上兄弟

收藏了。。。。。。。。。。。。。。

谢谢各位了，谢谢epoh老师

epoh老师，您好，能否帮我看看这个该如何实现？
您编写的case1，我修改了一下可以实现多个分位同时估计
case1：
x1=log(Produc$pcap)
x2=log(Produc$pc)
x3=log(Produc$emp)
x4=Produc$unemp
X=cbind(x1,x2,x3,x4)
y=log(Produc$gsp)
N<-length(y)
s= rep(1:48,rep(17,48))
X <- as.matrix(X)
p <- ncol(X)
n <- length(levels(as.factor(s)))
Z <- as.matrix.csr(model.matrix(~as.factor(s)-1))
lambda=1
w=c(0.1,0.25,0.5,0.25,0.1)
taus=c(0.1,0.25,0.5,0.75,0.9)
K <- length(w)
xx1<- cbind(as(w,"matrix.diag.csr") %x% X,w %x% Z)
Penalty <- cbind(as.matrix.csr(0,n,K*p),lambda*as(n,"matrix.diag.csr"))
xx<-rbind(xx1,Penalty)
y <- 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/2,n))
##rhs = (1 - tau) * c(t(a) %*% rep(1, length(y)))
rq.fit.sfn(xx,y,a)
case2： 我修改了一点
library(SparseM)
library(quantreg)
library(Ecdat)
data("Produc")  
#n=48, T=17, N=816
x1=log(Produc$pcap)
x2=log(Produc$pc)
x3=log(Produc$emp)
x4=Produc$unemp
X=cbind(x1,x2,x3,x4)
y=log(Produc$gsp)
s= rep(1:48,rep(17,48))
X <- as.matrix(X)
n <- length(levels(as.factor(s)))
Z <- as.matrix.csr(model.matrix(~as.factor(s)-1))
w<-1
xx<- cbind(as(w,"matrix.diag.csr") %x% X,w %x% Z)  tt<-as.matrix(xx)
taus<-c(0.1,0.25,0.75,0.5,0.9)
for(j in 1:length(taus)){
fit<-rqss(y~tt-1,tau=taus[j])
}
但是case不能实现多个分位同时运行
就假如我一次一次运行，如现在运行5个分位点，就5次，每一次都估计出x1-x4系数，而每一次的固定效应结果不同。
在rq.fit.sfn中，我同时估计5个分位点，得出的只有一组固定效应。那两者肯定不一样
在rqss中能实现5个分位同时估计吗？

library(SparseM)
library(quantreg)
library(Ecdat)
data("Produc")  
x1=log(Produc$pcap)
x2=log(Produc$pc)
x3=log(Produc$emp)
x4=Produc$unemp
X=cbind(x1,x2,x3,x4)
y=log(Produc$gsp)
N<-length(y)
s= rep(1:48,rep(17,48))
X <- as.matrix(X)
p <- ncol(X)
n <- length(levels(as.factor(s)))
Z <- as.matrix.csr(model.matrix(~as.factor(s)-1))
lambda=1
w=c(0.1,0.25,0.5,0.25,0.1)
taus=c(0.1,0.25,0.5,0.75,0.9)
K <- length(w)
xx1<- cbind(as(w,"matrix.diag.csr") %x% X,w %x% Z)
Penalty <- cbind(as.matrix.csr(0,n,K*p),lambda*as(n,"matrix.diag.csr"))
xx<-rbind(xx1,Penalty)
y <- 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/2,n))
##############case1
rq.fit.sfn(xx,y)
################
$coef
[1]  2.836199e-01  4.602578e-01  4.218402e-01 -1.293751e-02  2.836199e-01
[6]  4.602578e-01  4.218402e-01 -1.293751e-02  2.836199e-01  4.602578e-01
[11]  4.218402e-01 -1.293751e-02  2.836199e-01  4.602578e-01  4.218402e-01
[16] -1.293751e-02  2.836199e-01  4.602578e-01  4.218402e-01 -1.293751e-02
[21] -1.412773e-01  6.149486e-02  2.992886e-02 -4.000146e-02  1.343490e-01
[26]  2.414043e-01  2.334016e-01 -4.028156e-02 -4.768962e-02  1.829905e-01
[31] -5.174577e-02 -6.450664e-02 -5.472290e-02 -4.615042e-02  3.985560e-02
[36]  2.038454e-10  1.901025e-01  7.812054e-02  1.347637e-01  2.923434e-11
[41] -8.490129e-02 -3.631191e-02  2.705989e-02  1.312208e-03 -1.009822e-01
[46]  1.433116e-01  2.792632e-01  1.822704e-01  1.587613e-01 -5.398773e-02
[51]  3.816106e-10  1.061240e-02 -9.002397e-02  7.444380e-02  9.299111e-02
[56] -1.201060e-01  4.259159e-01 -5.967792e-02 -4.599368e-03 -1.326447e-01
[61] -1.550668e-01  1.112986e-01  2.706360e-01  9.057896e-02 -1.589800e-02
[66] -2.147877e-02 -3.657946e-03  2.225002e-01


###############case2
tt<-as.matrix(xx)
fit<-rqss(y~tt-1)
summary(fit)
Parametric coefficients:
       Estimate Std. Error t value Pr(>|t|)   
tt1   2.836e-01  1.726e-02  16.433  < 2e-16 ***
tt2   4.603e-01  1.473e-02  31.244  < 2e-16 ***
tt3   4.218e-01  1.098e-02  38.425  < 2e-16 ***
tt4  -1.294e-02  1.002e-03 -12.908  < 2e-16 ***
tt5   2.836e-01  1.686e-02  16.827  < 2e-16 ***
tt6   4.603e-01  1.441e-02  31.933  < 2e-16 ***
tt7   4.218e-01  1.077e-02  39.172  < 2e-16 ***
tt8  -1.294e-02  9.440e-04 -13.705  < 2e-16 ***
tt9   2.836e-01  1.629e-02  17.411  < 2e-16 ***
tt10  4.603e-01  1.413e-02  32.565  < 2e-16 ***
tt11  4.218e-01  1.048e-02  40.260  < 2e-16 ***
tt12 -1.294e-02  1.007e-03 -12.850  < 2e-16 ***
tt13  2.836e-01  1.686e-02  16.827  < 2e-16 ***
tt14  4.603e-01  1.441e-02  31.933  < 2e-16 ***
tt15  4.218e-01  1.077e-02  39.172  < 2e-16 ***
tt16 -1.294e-02  9.440e-04 -13.705  < 2e-16 ***
tt17  2.836e-01  1.726e-02  16.433  < 2e-16 ***
tt18  4.603e-01  1.473e-02  31.244  < 2e-16 ***
tt19  4.218e-01  1.098e-02  38.425  < 2e-16 ***
tt20 -1.294e-02  1.002e-03 -12.908  < 2e-16 ***
tt21 -1.413e-01  8.672e-03 -16.291  < 2e-16 ***
tt22  6.149e-02  9.322e-03   6.596 4.78e-11 ***
tt23  2.993e-02  8.580e-03   3.488 0.000492 ***
tt24 -4.000e-02  1.679e-02  -2.382 0.017265 *  
tt25  1.343e-01  9.265e-03  14.501  < 2e-16 ***
tt26  2.414e-01  1.648e-02  14.649  < 2e-16 ***
tt27  2.334e-01  1.483e-02  15.738  < 2e-16 ***
tt28 -4.028e-02  1.087e-02  -3.707 0.000213 ***
tt29 -4.769e-02  9.713e-03  -4.910 9.49e-07 ***
tt30  1.830e-01  9.089e-03  20.133  < 2e-16 ***
tt31 -5.175e-02  1.210e-02  -4.278 1.94e-05 ***
tt32 -6.451e-02  8.223e-03  -7.845 5.55e-15 ***
tt33 -5.472e-02  1.170e-02  -4.679 2.98e-06 ***
tt34 -4.615e-02  1.009e-02  -4.576 4.90e-06 ***
tt35  3.986e-02  1.140e-02   3.497 0.000477 ***
tt36  2.038e-10  4.322e-03   0.000 1.000000   
tt37  1.901e-01  9.868e-03  19.264  < 2e-16 ***
tt38  7.812e-02  1.591e-02   4.910 9.50e-07 ***
tt39  1.348e-01  1.544e-02   8.729  < 2e-16 ***
tt40  2.923e-11  6.683e-03   0.000 1.000000   
tt41 -8.490e-02  1.161e-02  -7.316 3.09e-13 ***
tt42 -3.631e-02  1.556e-02  -2.333 0.019698 *  
tt43  2.706e-02  1.056e-02   2.563 0.010407 *  
tt44  1.312e-03  7.382e-03   0.178 0.858914   
tt45 -1.010e-01  1.695e-02  -5.956 2.82e-09 ***
tt46  1.433e-01  1.238e-02  11.579  < 2e-16 ***
tt47  2.793e-01  1.128e-02  24.764  < 2e-16 ***
tt48  1.823e-01  1.111e-02  16.407  < 2e-16 ***
tt49  1.588e-01  1.160e-02  13.682  < 2e-16 ***
tt50 -5.399e-02  1.946e-02  -2.774 0.005568 **
tt51  3.816e-10  5.112e-03   0.000 1.000000   
tt52  1.061e-02  1.576e-02   0.674 0.500640   
tt53 -9.002e-02  1.151e-02  -7.820 6.66e-15 ***
tt54  7.444e-02  9.354e-03   7.958 2.22e-15 ***
tt55  9.299e-02  1.008e-02   9.227  < 2e-16 ***
tt56 -1.201e-01  1.349e-02  -8.906  < 2e-16 ***
tt57  4.259e-01  1.282e-02  33.218  < 2e-16 ***
tt58 -5.968e-02  9.418e-03  -6.337 2.61e-10 ***
tt59 -4.599e-03  1.166e-02  -0.395 0.693160   
tt60 -1.326e-01  1.287e-02 -10.304  < 2e-16 ***
tt61 -1.551e-01  9.147e-03 -16.952  < 2e-16 ***
tt62  1.113e-01  1.079e-02  10.312  < 2e-16 ***
tt63  2.706e-01  1.195e-02  22.653  < 2e-16 ***
tt64  9.058e-02  1.167e-02   7.759 1.09e-14 ***
tt65 -1.590e-02  1.465e-02  -1.085 0.277769   
tt66 -2.148e-02  1.096e-02  -1.959 0.050195 .  
tt67 -3.658e-03  1.035e-02  -0.354 0.723707   
tt68  2.225e-01  1.932e-02  11.514  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


  Quantile Fidelity at tau = 0.5  is      17.5209
  Effective Degrees of Freedom = 248       Sample Size = 4128

##epoh，我对您程序还要疑惑，为什么在case1中，函数rq.fit.sfn(xx,y)不加入
a <- c((w*(1-taus)) %x% (t(X)%*%rep(1,N)), sum(w*(1-taus)) * (t(Z) %*% rep(1,N)) + lambda * rep(1/2,n))？？因为我看了rq.fit.sfn内部是rq.fit.sfn（D，x，a）
在您程序中，为什么没有指明是这5个分位点，而就直接可以？

我试加了，程序也是可以运行的，只是结果不一样
> library(SparseM)
> library(quantreg)
> library(Ecdat)
> data("Produc")  
> x1=log(Produc$pcap)
> x2=log(Produc$pc)
> x3=log(Produc$emp)
> x4=Produc$unemp
> X=cbind(x1,x2,x3,x4)
> y=log(Produc$gsp)
> N<-length(y)
> s= rep(1:48,rep(17,48))
> X <- as.matrix(X)
> p <- ncol(X)
> n <- length(levels(as.factor(s)))
> Z <- as.matrix.csr(model.matrix(~as.factor(s)-1))
> lambda=1
> w=c(0.1,0.25,0.5,0.25,0.1)
> taus=c(0.1,0.25,0.5,0.75,0.9)
> K <- length(w)
> xx1<- cbind(as(w,"matrix.diag.csr") %x% X,w %x% Z)
> Penalty <- cbind(as.matrix.csr(0,n,K*p),lambda*as(n,"matrix.diag.csr"))
> xx<-rbind(xx1,Penalty)
> y <- 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/2,n))
> rq.fit.sfn(xx,y,a)
$coef
[1]  0.3479225632  0.4579939079  0.3419116807 -0.0132019291  0.3479225632
[6]  0.4579939079  0.3419116807 -0.0132019291  0.3479225632  0.4579939079
[11]  0.3419116807 -0.0132019291  0.3479225632  0.4579939079  0.3419116807
[16] -0.0132019291  0.3479225632  0.4579939079  0.3419116807 -0.0132019291
[21] -0.1333094818 -0.0001001004  0.0005137676 -0.0406687377  0.0783307839
[26]  0.1897515687  0.1388376698 -0.0243480092 -0.0477144706  0.1091862206
[31] -0.0495600876 -0.0631116318 -0.0821633242 -0.0649543957 -0.0088250215
[36] -0.0295872475  0.1433398180  0.0385951674  0.1092909477 -0.0157361471
[41] -0.1005352958 -0.0707680805  0.0030014620 -0.0804312417 -0.1401517851
[46]  0.0825749965  0.2307316486  0.1470793361  0.0835593363 -0.0659401370
[51]  0.0033954589 -0.0757982221 -0.0809079723  0.0426541558  0.0304444640
[56] -0.1109061177  0.3313202668 -0.0641237261 -0.0821561387 -0.1284287320
[61] -0.1464990138  0.0452016648  0.1755758956  0.0560798608 -0.0777549839
[66] -0.0559941815 -0.0202300682  0.0801947593

您觉得我哪里理解错了？

糊涂了,封掉a没发觉.
这问题,我再仔细看一下function rqss().

谢谢epoh指导
rqss（）函数也有用到方法lasso，所以也是加了惩罚项，因此才可以实现rq.panel.fit功能（他就是通过引入固定效应惩罚项的）。
我看了很久，没有看懂

epoh老师，如果rqss（）可以实现多个分位同时估计，那也直接给出了标准误和系数，能不能直接画出第一帖附件中的图形？

谢谢您