人大经济论坛

标题: 带置信区间的画图问题 [打印本页]
作者: ywh19860616    时间: 2011-1-10 08:54:22     标题: 带置信区间的画图问题

我想使用R画出带置信区间的图像，类似于了分位数回归图，见图例。
    文字描述：
    假如我有两个向量X和Y，每个X对应一个Y
   X=[0.1, 0.2, 0.3, 0.4, 0.6, 0.7];
   Y=[0.2,0.3,0.4,0.5,0.6,0.7];
   如果只画X和Y的对应图，那很简单，直接plot(X,Y)就能解决。
  但是现在问题是Y中每个数值都对应一个置信区间，
如[0.15,0.22],[0.25,0.32],[0.35,0.42],[0.45,0.52],[0.55,0.62],[0.65,0.7]
  怎么画如下图形？请大家指点，谢谢

图片附件: 未命名.jpg (2011-1-10 08:56:46, 15.64 KB)
https://bbs.pinggu.org/forum.php?mod=attachment&aid=ODMyMDg5fDQxNTgzOGUxfDE3NjY0NzQ0MTR8MHww

作者: ltx5151    时间: 2011-1-10 10:03:42

三个选择
1. lines
2. abline
如果真的想要途中的那种填充效果，用contour。
这些函数help里面都介绍的比较详细，lz可以自己查阅一下。
作者: shenbaiseshatan    时间: 2011-1-10 15:02:30

简单点可以实现类似功能的程序
  1. n<-100
  2. set.seed(123);x<-rnorm(n,5,1);y<-x+rnorm(n,3,.5)
  3. lm.xy<-lm(y~x)
  4. lm.for.plot<-as.data.frame(predict.lm(lm.xy,interval = "confidence"))
  5. attach(lm.for.plot)
  6. index<-order(fit)
  7. plot(sort(fit),pch=20,col="blue")
  8. points(1:n,lwr[index],pch=20,col="red",type='l')
  9. points(1:n,upr[index],pch=20,col="red",type='l')
  10. segments(1:n,lwr[index],1:n,upr[index],col="brown")
复制代码

  1. n<-100
  2. set.seed(123);x<-rnorm(n,5,1);y<-x+rnorm(n,3,.5)
  3. lm.xy<-lm(y~x)
  4. lm.for.plot<-as.data.frame(predict.lm(lm.xy,interval = "confidence"))
  5. attach(lm.for.plot)
  6. index<-order(fit)
  7. index.fill<-seq(1,n,by=2)
  8. points.fill<-c(upr[index],rev(lwr[index]))
  9. plot(c(1:n,n:1),points.fill,type='n')
  10. polygon(c(1:n,n:1),points.fill,border=NA,col="gray")
  11. points(sort(fit),pch=20,col="blue",add=T)
  12. points(1:n,lwr[index],pch=20,col="red",type='l')
  13. points(1:n,upr[index],pch=20,col="red",type='l')
  14. #segments(1:n,lwr[index],1:n,upr[index],col="brown")
复制代码
作者: 楚韵荆风    时间: 2011-1-11 13:35:49

把置信带下面的数据和上面的数据分别用两个数组y1,y2存储，再用lines命令就可以了，
plot(x,y,type="l")
lines(x,y1)
lines(x,y2)
作者: trier2006    时间: 2011-1-13 08:58:47

嗯，大家的建议本质好像都是计算出区间值，然后画上去。似乎没有现成的函数用。
作者: shenbaiseshatan    时间: 2011-1-14 14:12:08

5# trier2006
现成的函数有类似的，比如Hmisc包的Dotplot
作者: epoh    时间: 2011-3-26 11:18:23

function rq.fit.panel()
is now superseded by the use of the function rqss() in the R quantreg package,
or the use of method = "lasso" for the function rq() in that package.

Package‘quantreg’Version 4.57  Date 2011-2-26  

#############
library(quantreg)
data(engel)
attach(engel)
summary(rq(foodexp~income,tau = 1:49/50,data=engel))

Call: rq(formula = foodexp ~ income, tau = 1:49/50, data = engel)
tau: [1] 0.02
Coefficients:
            coefficients lower bd  upper bd
(Intercept) 108.60499    103.39624 141.26345
income        0.34664      0.28526   0.34664

Call: rq(formula = foodexp ~ income, tau = 1:49/50, data = engel)
tau: [1] 0.04
............

#Example of plotting of coefficients and their confidence bands
plot(summary(rq(foodexp~income,tau = 1:49/50,data=engel)))
作者: ywh19860616    时间: 2011-3-26 13:02:39

##epoh ,哈哈，谢谢啊
我看过rqss（）函数帮助，但是找不到识别面板结构的程序，还不知道怎么使用。
作者: epoh    时间: 2011-3-28 19:58:36

panel data:
  package "Ecdat"
  Produc Us States Production
  a panel of 48 observations from 1970 to 1986
  number of observations : 816

##############
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))
w=1

####case 1
####use function rq.fit.panel()
require(SparseM)
require(quantreg)

X <- as.matrix(X)
p <- ncol(X)
n <- length(levels(as.factor(s)))
N <- length(y)
if(N != length(s) || N != nrow(X))
    stop("dimensions of y,X,s must match")
Z <- as.matrix.csr(model.matrix(~as.factor(s)-1))
xx<- cbind(as(w,"matrix.diag.csr") %x% X,w %x% Z)
rq.fit.sfn(xx,y,tau=0.5)
>$coef
[1] -0.001856897  0.227955745  0.806905899 -0.003254027  2.348518722
[6]  2.531511125  2.416402785  2.646310826  2.554288408  2.598468250
[11]  2.540988532  2.473361510  2.401692284  2.497764134  2.537963266
[16]  2.425204261  2.482897448  2.510841856  2.553142722  2.703008416
[21]  2.404225489  2.537962766  2.486090983  2.567535469  2.455404207
[26]  2.403301315  2.463135907  2.545757570  2.474922409  2.515871136
[31]  2.443668078  2.580372182  2.643673666  2.605560174  2.384878447
[36]  2.541337214  2.468143678  2.587964518  2.509402455  2.432777812
[41]  2.495204924  2.258353874  2.465798076  2.372996183  2.600392620
[46]  2.471721790  2.457974490  2.519673773  2.590197854  2.472063213
[51]  2.470425047  2.862672642


####case 2
####use function rqss()
tt=as.matrix(xx)
fit=rqss(y~tt-1)
summary(fit)
>
Formula:
y ~ tt - 1

Parametric coefficients:
      Estimate    Std. Error   t value    Pr(>|t|)   
tt1  -0.001857   0.040474  -0.046   0.9634   
tt2   0.227956   0.033949   6.715 3.68e-11 ***
tt3   0.806906   0.041542  19.424  < 2e-16 ***
tt4  -0.003254   0.001402  -2.321   0.0205 *  
tt5   2.348519   0.227838  10.308  < 2e-16 ***
tt6   2.531511   0.226443  11.179  < 2e-16 ***
tt7   2.416403   0.216806  11.145  < 2e-16 ***
tt8   2.646311   0.254619  10.393  < 2e-16 ***
tt9   2.554288   0.222018  11.505  < 2e-16 ***
tt10  2.598468   0.220106  11.806  < 2e-16 ***
tt11  2.540989   0.210428  12.075  < 2e-16 ***
tt12  2.473362   0.234957  10.527  < 2e-16 ***
tt13  2.401692   0.226512  10.603  < 2e-16 ***
tt14  2.497764   0.208589  11.975  < 2e-16 ***
tt15  2.537963   0.241893  10.492  < 2e-16 ***
tt16  2.425204   0.227976  10.638  < 2e-16 ***
tt17  2.482897   0.232136  10.696  < 2e-16 ***
tt18  2.510842   0.231957  10.825  < 2e-16 ***
tt19  2.553143   0.231057  11.050  < 2e-16 ***
tt20  2.703008   0.244407  11.059  < 2e-16 ***
tt21  2.404225   0.201550  11.929  < 2e-16 ***
tt22  2.537963   0.228998  11.083  < 2e-16 ***
tt23  2.486091   0.220021  11.299  < 2e-16 ***
tt24  2.567535   0.239121  10.737  < 2e-16 ***
tt25  2.455404   0.235207  10.439  < 2e-16 ***
tt26  2.403301   0.223782  10.739  < 2e-16 ***
tt27  2.463136   0.226028  10.897  < 2e-16 ***
tt28  2.545758   0.227966  11.167  < 2e-16 ***
tt29  2.474922   0.232918  10.626  < 2e-16 ***
tt30  2.515871   0.211395  11.901  < 2e-16 ***
tt31  2.443668   0.201243  12.143  < 2e-16 ***
tt32  2.580372   0.224698  11.484  < 2e-16 ***
tt33  2.643674   0.222489  11.882  < 2e-16 ***
tt34  2.605560   0.253096  10.295  < 2e-16 ***
tt35  2.384878   0.221027  10.790  < 2e-16 ***
tt36  2.541337   0.225034  11.293  < 2e-16 ***
tt37  2.468144   0.238798  10.336  < 2e-16 ***
tt38  2.587965   0.227848  11.358  < 2e-16 ***
tt39  2.509402   0.221705  11.319  < 2e-16 ***
tt40  2.432778   0.239579  10.154  < 2e-16 ***
tt41  2.495205   0.191899  13.003  < 2e-16 ***
tt42  2.258354   0.211827  10.661  < 2e-16 ***
tt43  2.465798   0.222536  11.080  < 2e-16 ***
tt44  2.372996   0.232450  10.209  < 2e-16 ***
tt45  2.600393   0.251702  10.331  < 2e-16 ***
tt46  2.471722   0.214506  11.523  < 2e-16 ***
tt47  2.457974   0.199837  12.300  < 2e-16 ***
tt48  2.519674   0.226890  11.105  < 2e-16 ***
tt49  2.590198   0.241976  10.704  < 2e-16 ***
tt50  2.472063   0.224931  10.990  < 2e-16 ***
tt51  2.470425   0.229340  10.772  < 2e-16 ***
tt52  2.862673   0.230734  12.407  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


  Quantile Fidelity at tau = 0.5  is      10.8697
  Effective Degrees of Freedom = 52        Sample Size = 816
作者: king64    时间: 2011-3-28 21:36:19

#（1）先回归
   lmfit <- lm( dat$Sunday ~ dat$Daily, data=dat )
# （2）  计算置信区间（95% ）---------------------------
     Ypred <- predict( lmfit, interval = "confidence", newdata= dat, level = 0.95 )
#（3） 画图
        plot(dat$Daily, dat$Sunday,
              xlab="Daily", ylab="Sunday",
             main="Figure:  Yhat Predict confidence")
        lines ( Ypred[ ,1 ] ~ dat$Daily, col="red", lwd=1 )
        lines ( Ypred[ ,2 ] ~ dat$Daily, col="gray", lwd=1 )
        lines ( Ypred[ ,3 ] ~ dat$Daily, col="gray", lwd=1 )
作者: ywh19860616    时间: 2011-3-28 22:49:55

epoh，谢谢您了
呵呵，感谢楼上兄弟
作者: casperyc    时间: 2011-3-30 00:36:47

收藏了。。。。。。。。。。。。。。
作者: ywh19860616    时间: 2011-3-31 13:59:02

谢谢各位了，谢谢epoh老师
作者: ywh19860616    时间: 2011-4-1 12:40:54

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个分位同时估计吗？
作者: epoh    时间: 2011-4-1 13:37:58

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
作者: ywh19860616    时间: 2011-4-1 13:53:02

##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个分位点，而就直接可以？
作者: ywh19860616    时间: 2011-4-1 13:54:34

我试加了，程序也是可以运行的，只是结果不一样
> 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

您觉得我哪里理解错了？
作者: epoh    时间: 2011-4-1 21:48:23

糊涂了,封掉a没发觉.
这问题,我再仔细看一下function rqss().
作者: ywh19860616    时间: 2011-4-1 22:12:25

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

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

谢谢您
作者: epoh    时间: 2011-4-2 15:07:15

仔细把底下三个程序看完
function quantreg()
  function rqss()
  function sfn()

  1.function rq()  method = "lasso"
     清楚的告诉你
     if(object$method == "lasso")
        stop("no inference for lasso'd rq fitting: try rqss (if brave)")

  2.看来panel data还是只能适用function rq.fit.sfn()
    a = structure of the design matrix A' stored in csr format
    毕竟是一个  "design matrix"--------s,taus,w

    不过只估计出coefficients & resid
    楼主需要的统计量,可能要用bootstrap
.
作者: ywh19860616    时间: 2011-4-2 16:30:58

##epoh前辈，根据作者说明是rq和rqss都可以实现rq.fit.sfn功能的了
我都试过了
1、rq函数选择lasso时根本不可行，结果和您一样
2、rqss（）函数我觉得应该可以的，从您第一次用一个分位点试验的结果可以看出，他的结果和rq.fit.sfn结果是一样的，所以应该是可以的，我只是不知道怎么修改成同时估计多个分位点。
3、因为rq.fit.sfn函数要经过bootstrap抽样，感觉不是很稳定

前辈，若使用rq.fit.sfn函数，再bootstrap得出标准误和t值，可以计算其置信区间得到 帖子中的图像吗？
作者: ywh19860616    时间: 2011-4-2 16:32:26

呵呵，那个bootstrap实现统计量的程序已经让您编写了，我只是觉得这个程序不稳定了
作者: epoh    时间: 2011-4-2 17:06:52

rqss（）函数只适用一个分位点
无法估计同时估计多个分位点
最主要是在程序中summary.rqss()
已经写明
if(tau + h > 1)
    stop("tau + h > 1:  error in summary.rq")

ex:tau=0.5
     h=0.1
      if(tau + h > 1)
      stop("tau + h > 1:  error in summary.rq")
      #OK
   
    tau=c(0.5,0.75)
      if(tau + h > 1)
      stop("tau + h > 1:  error in summary.rq")
      #PROBLEM
   
  function rqss()呼叫
  rqss.fit(X, Y, tau = tau, rhs = rhs, method = "sfnc")
  再呼叫
  fit <- sfn = rq.fit.sfn(x, y, tau = tau,  nsubmax = nsubmax, nnzlmax = nnzlmax, tmpmax = tmpmax, ...)
  再由fortran srqfnb.o <- .Fortran("srqfn",....)

  所产生的object fit 就是提供 coef &resid

  使用rq.fit.sfn函数，再bootstrap得出标准误和t值，
  有了标准误,就可以计算其置信区间,得到帖子中的图像.
作者: ywh19860616    时间: 2011-4-2 20:22:04

epoh前辈，谢谢您了
画图程序我自己看看，再向您请教

rqss不能实现多个分位点，那就没有什么意义了，因为quantile regression要的就是对几个分位进行加权处理
作者: epoh    时间: 2011-4-2 22:00:27

Stata:
command qreg:
  qreg depvar indepvars , quantile(choice of quantile)

  The standard errors estimated by qreg are valid for the
  location model only. To obtain better standard errors
  estimates (given with the bootstrap), you should use
bsqreg instead.

  To obtain simultaneously several quantile regressions, use
  sqreg:
  sqreg depvar indepvars , quantiles(choice of quantiles)

  N.B: The standard errors estimates provided by sqreg
         Are computed with the bootstrap.
作者: ywh19860616    时间: 2011-4-2 22:17:56

epoh前辈，谢谢您
呵呵，可能我刚才说的用bootstrap结果不稳定，是不正确的
那个stata命令我看过了，如果是用截面数据，那很完美，图形也非常漂亮

可惜stata目前还未有估计面板分位的命令


呵呵，前辈，您能否把求置信区间的命令也写出来？
作者: epoh    时间: 2011-4-3 16:56:23

Grunfeld   Grunfeld Investment Data
  a panel of 10 observations from 1935 to 1954
  number of observations : 200
##########
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))
w=1
n <- length(levels(as.factor(s)))
Z <- as.matrix.csr(model.matrix(~as.factor(s)-1))
xx<- cbind(as(w,"matrix.diag.csr") %x% X,w %x% Z)
####################
tt=as.matrix(xx)
taus=seq(0.1,0.9,0.01)
ntau=length(taus)  #81
coef1mat=matrix(NA,ntau,5)
coef2mat=matrix(NA,ntau,5)
rq.fit=rq(y~tt-1,taus)
sum=summary(rq.fit,"ker")
#If n is large(>1000), t=1.96
t=1.96
for (i in 1:ntau){
#coefficient1
coef1mat[i,1]=sum[]$tau
coef1mat[i,2]=sum[]$coefficients[1,1]
coef1mat[i,3]=sum[]$coefficients[1,2]
coef1mat[i,4]=sum[]$coefficients[1,1]+t*sum[]$coefficients[1,2]
coef1mat[i,5]=sum[]$coefficients[1,1]-t*sum[]$coefficients[1,2]
#coefficient2
coef2mat[i,1]=sum[]$tau
coef2mat[i,2]=sum[]$coefficients[2,1]
coef2mat[i,3]=sum[]$coefficients[2,2]
coef2mat[i,4]=sum[]$coefficients[2,1]+t*sum[]$coefficients[2,2]
coef2mat[i,5]=sum[]$coefficients[2,1]-t*sum[]$coefficients[2,2]
}
coef1mat
coef2mat
###plot coef1,coef2 95%CI
par(mfrow = c(1, 2))
####coefficient1
upperbd1=coef1mat[,4]
lowerbd1=coef1mat[,5]
coef1=coef1mat[,2]
plot(taus,upperbd1,type = 'n',  ylim = c(0.00, 0.20),
     ylab = 'Coef1 Value', xlab = 'Quantile')
lines(taus, lowerbd1, col = 'grey')
lines(taus, upperbd1, col = 'grey')
polygon(c(taus, rev(taus)), c(upperbd1, rev(lowerbd1)),
     col = "lightblue", border = 'black')
points(taus, coef1,col='blue',pch=21)
lines(taus, coef1,col='blue')
title("gross Investment vs value of the firm")
####coefficient2
upperbd2=coef2mat[,4]
lowerbd2=coef2mat[,5]
coef2=coef2mat[,2]
plot(taus,upperbd2,type = 'n',  ylim = c(0.00, 0.50),
     ylab = 'Coef2 Capital', xlab = 'Quantile')
lines(taus, lowerbd2, col = 'grey')
lines(taus, upperbd2, col = 'grey')
polygon(c(taus, rev(taus)), c(upperbd2, rev(lowerbd2)),
     col = "lightblue", border = 'black')
points(taus, coef2,col='blue',pch=21)
lines(taus, coef2,col='blue')
title("gross Investment vs stock of plant and equipment")

         graph1.jpeg

图片附件: graph1.jpeg (2011-4-3 16:56:23, 75.74 KB)
https://bbs.pinggu.org/forum.php?mod=attachment&aid=ODgwNzgwfDVhOTlmZDI1fDE3NjY0NzQ0MTR8MHww

作者: ywh19860616    时间: 2011-4-3 22:33:23

epoh前辈，谢谢您了
您真是太厉害了，这应该和R内部的程序一样了
不过您还是使用rq函数，而不是在rq.fit.panel运行基础上的
谢谢
作者: epoh    时间: 2011-4-3 23:11:31

不太懂你的题意
作者: ywh19860616    时间: 2011-4-4 00:09:51

epoh前辈，您说过rq和rqss根本不能实现多个分位数同时估计的，而且使用rq函数，内部就有画图程序，效果和您的是一样的。如plot（summary（...））就可以。
我的意思是在rq.fit.sfn函数基础上（即使用rq.fit.panel函数），使用bootstrap出来的标准误直接计算出置信区间，进而画出图形



n1=length(y)  #30
b=1000     #Number of bootstrap samples
ncoef=16
boot_bhat=matrix(NA,b, ncoef)
block_length = 10  
num_blocks = n1/block_length    #n/block_length  
Indices = seq(1:n1)  # All of the indices from 1 to n
Indices = matrix(Indices,block_length,num_blocks)
for (i in 1:b){          #Number of bootstrap samples
randblock =sample(seq(1:num_blocks),num_blocks,replace = TRUE) # Choose which blocks to use
Ind_sim = Indices[,randblock]       #Find which data are in each block  
Ind_sim = c(Ind_sim)      
Xsim = X[Ind_sim,1:dim(X)[2]]      #Construct the x data
Ysim = y[Ind_sim]          #Construct the y data
boot_bhat[i,] = rq.fit.panel(Xsim,Ysim,s)$coef
}
bhat=colMeans(boot_bhat)
#bhat

cov=cov(boot_bhat)
serr = sqrt(diag(cov))
#serr   

p <- length(bhat)
rdf <- n1 - p
vnames<- dimnames(x)[[2]]
coef <- array(bhat, c(p, 4))
dimnames(coef) <- list(vnames, c("Value", "Std. Error", "t value","Pr(>|t|)"))
coef[, 2] <- serr
coef[, 3] <- coef[, 1]/coef[, 2]
coef[, 4] <- if (rdf > 0) 2 * (1 - pt(abs(coef[, 3]), rdf))
coef
就是在上面计算出系数coef后，直接计算置信区间，进而再画出图形（尽管现在估计的分位点较少，但是我只是想看大致趋势，没有关系）
作者: epoh    时间: 2011-4-4 09:59:57

在function rq.fit.sfn()基础上
底下结果都相同,你要的是哪种图形
fit1=rq.fit.sfn(xx,yy,tau=0.5,a)
fit1
fit2=rq.fit.sfn(xx,yy,tau=0.3,a)
fit2
fit3=rq.fit.sfn(xx,yy,tau=0.7,a)
fit3
fit4=rq.fit.sfn(xx,yy,tau=c(0.1,0.25,0.5,0.75,0.9) ,a)
fit4


另针你说用bootstrap结果不稳定,
应该说bootstrap结果是稳定,
是我block bootstrap程序,忙中出错.

还是回到我提供的BOOTSTRAP FOR PANEL DATA MODELS.pdf
page 6/30 block bootstrap resampling is operating with columns
所以根据page 2/30
panel dataset with N individuals and T time periods is represented
              by a matrix Y of N rows and T columns

就以提供的数据为例:
year=c(1991,1992,1993,1994,1995,1996,1991,1992,1993,1994,1995,1996)
x1=c(10,20,30,40,50,60,70,80,90,100,110,120)
x2=c(101,201,301,401,501,601,701,801,901,1001,1101,1201)
x=cbind(x1,x2)
y=year
#x
> x
      year  x1   x2
[1,] 1991  10  101
[2,] 1992  20  201
[3,] 1993  30  301
[4,] 1994  40  401
[5,] 1995  50  501
[6,] 1996  60  601
[7,] 1991  70  701
[8,] 1992  80  801
[9,] 1993  90  901
[10,] 1994 100 1001
[11,] 1995 110 1101
[12,] 1996 120 1201

数据需先转换成下列pattern, 再以 column 抽样

   T 1991 1992 1993 1994 1995 1996
              10      20     30      40      50      60
      id  701    801   901 1001  1101 1201

然后再还原成paneldata格式,进行回归.
作者: ywh19860616    时间: 2011-4-4 10:50:21

epoh前辈，呵呵，我还是没有表达清楚，勿怪
我要的图形应该是fit4=rq.fit.sfn(xx,yy,tau=c(0.1,0.25,0.5,0.75,0.9) ,a)
意思就是说假如现在我有4个解释变量（包括intercept），根据上面命令，每个解释变量都估计了5个分位点
我想直观观察每个解释变量随着分位点不同，其值的变化规律。
就是以分位点（0.1,0.25,0.5,0.75,0.9）为横轴，而4个解释变量的系数值为纵轴，这样就画出了4幅反映
各解释变量规律的图像了。如
n1=length(y)  #30
b=1000     #Number of bootstrap samples
ncoef=16
boot_bhat=matrix(NA,b, ncoef)
block_length = 10  
num_blocks = n1/block_length    #n/block_length  
Indices = seq(1:n1)  # All of the indices from 1 to n
Indices = matrix(Indices,block_length,num_blocks)
for (i in 1:b){          #Number of bootstrap samples
randblock =sample(seq(1:num_blocks),num_blocks,replace = TRUE) # Choose which blocks to use
Ind_sim = Indices[,randblock]       #Find which data are in each block  
Ind_sim = c(Ind_sim)      
Xsim = X[Ind_sim,1:dim(X)[2]]      #Construct the x data
Ysim = y[Ind_sim]          #Construct the y data
boot_bhat[i,] = rq.fit.panel(Xsim,Ysim,s)$coef
}
bhat=colMeans(boot_bhat)
#bhat

cov=cov(boot_bhat)
serr = sqrt(diag(cov))
#serr   

p <- length(bhat)
rdf <- n1 - p
vnames<- dimnames(x)[[2]]
coef <- array(bhat, c(p, 4))
dimnames(coef) <- list(vnames, c("Value", "Std. Error", "t value","Pr(>|t|)"))
coef[, 2] <- serr
coef[, 3] <- coef[, 1]/coef[, 2]
coef[, 4] <- if (rdf > 0) 2 * (1 - pt(abs(coef[, 3]), rdf))
coef
如上面您已经bootstrap得出了各解释变量在不同分位点的系数，那就想把这些系数趋势画在图形中反映出来


我前面已经说了，呵呵，不是bootstrap不稳定问题，而是我自己参数设置问题
啊，bootstrap抽样有错误吗？那前辈能修改吗？
作者: epoh    时间: 2011-4-4 13:48:34

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)  
Penalty <- cbind(as.matrix.csr(0,n,K*p),lambda*as(n,"matrix.diag.csr"))  
D <- rbind(Fidelity,Penalty)   
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)
#####block bootstrap
#n=10, T=20, N=200
N=length(y)  #200
T=20
n=10
blockmat=matrix(NA,n,T)
for(i in 1:n){
   for(j in 1:T){
  #    blockmat[i,j]=(j-1)*T+i
        blockmat[i,j]=(i-1)*T+j
      }
}
b=5000   #Number of bootstrap samples
ncoef=28
boot_bhat=matrix(NA,b, ncoef)
indices = seq(1:T)  # All of the indices from 1 to 20
####

for (j in 1:b){          #Number of bootstrap samples
randblock =sample(indices,replace = TRUE) # Choose which blocks to use
ind_sim=0*seq(1:n)
   for(i in 1:n){
    ind_sim[((i-1)*T+1):(i*T)]=(i-1)*T+randblock
          }   
xsim=cbind(x1[ind_sim],x2[ind_sim]) #Construct the X data
xx1=cbind(as(w,"matrix.diag.csr") %x% xsim,w %x% Z)
xx<-rbind(xx1,Penalty)
ysim = y[ind_sim]          #Construct the y data
yy <- c(w %x% ysim,rep(0,n))
a <- c((w*(1-taus)) %x% (t(xsim)%*%rep(1,N)),
                sum(w*(1-taus)) * (t(Z) %*% rep(1,N)) + lambda * rep(1/2,n))
boot_bhat[j,] = rq.fit.sfn(xx,yy,taus,a)$coef
}

bhat=colMeans(boot_bhat)
#bhat
cov=cov(boot_bhat)
serr = sqrt(diag(cov))
#serr   
p <- length(bhat)
rdf <- N-ncoef
coef <- array(bhat, c(p, 4))
dimnames(coef) <- list(seq(1:ncoef), c("Value", "Std. Error", "t value","Pr(>|t|)"))
coef[, 2] <- serr
coef[, 3] <- coef[, 1]/coef[, 2]
coef[, 4] <- if (rdf > 0) 2 * (1 - pt(abs(coef[, 3]), rdf))
coef
> coef
           Value  Std. Error    t value     Pr(>|t|)
1     0.06107249  0.01386247  4.4055984 1.850521e-05
2     0.15487193  0.02205991  7.0205165 4.903056e-11
3     0.07234995  0.01482287  4.8809663 2.396156e-06
4     0.16567351  0.02101250  7.8845210 3.477219e-13
5     0.08053246  0.01565537  5.1440794 7.267875e-07
6     0.17087047  0.02093900  8.1603941 6.794565e-14
7     0.08633521  0.01577335  5.4734852 1.539569e-07
8     0.17589873  0.02234190  7.8730426 3.721468e-13
9     0.08989036  0.01543300  5.8245559 2.752130e-08
10    0.18514027  0.02556903  7.2408017 1.425859e-11
11    0.09275225  0.01509710  6.1437128 5.434659e-09
12    0.19914111  0.03041511  6.5474402 6.494829e-10
13    0.09563246  0.01524897  6.2714042 2.799298e-09
14    0.21819694  0.03619383  6.0285673 9.816594e-09
15    0.09971964  0.01643330  6.0681427 8.017352e-09
16    0.24902179  0.05148496  4.8367874 2.915302e-06
17    0.10538962  0.02081095  5.0641423 1.048965e-06
18    0.31456336  0.07023116  4.4789711 1.362888e-05
19   70.87586275 53.31286209  1.3294327 1.854658e-01
20  178.11811868 39.54191280  4.5045398 1.224064e-05
21 -145.88272883 33.85621760 -4.3088903 2.754051e-05
22   -1.70397795  8.63275211 -0.1973853 8.437591e-01
23  -51.72378988 17.42847433 -2.9677750 3.427315e-03
24   -1.91490600  6.15461311 -0.3111334 7.560759e-01
25  -24.71403762  8.68786452 -2.8446619 4.985616e-03
26  -31.82235648 11.51383878 -2.7638355 6.335172e-03
27  -42.08427967 10.15349925 -4.1448055 5.328559e-05
28   -4.48096796  0.96710385 -4.6333886 7.076866e-06
###############
coef1mat=matrix(NA,ntaus,5)
coef2mat=matrix(NA,ntaus,5)
#If n is large(>1000), t=1.96
t=1.96
for (i in 1:ntaus){
#coefficient1
coef1mat[i,1]=taus
coef1mat[i,2]=coef[(2*(i-1)+1),1]
coef1mat[i,3]=coef[(2*(i-1)+1),2]
coef1mat[i,4]=coef[(2*(i-1)+1),1]+t*coef[(2*(i-1)+1),2]
coef1mat[i,5]=coef[(2*(i-1)+1),1]-t*coef[(2*(i-1)+1),2]
#coefficient2
coef2mat[i,1]=taus
coef2mat[i,2]=coef[2*i,1]
coef2mat[i,3]=coef[2*i,2]
coef2mat[i,4]=coef[2*i,1]+t*coef[2*i,2]
coef2mat[i,5]=coef[2*i,1]-t*coef[2*i,2]
}
coef1mat
coef2mat
> coef1mat
      [,1]       [,2]       [,3]       [,4]       [,5]
[1,]  0.1 0.06107249 0.01386247 0.08824294 0.03390204
[2,]  0.2 0.07234995 0.01482287 0.10140278 0.04329712
[3,]  0.3 0.08053246 0.01565537 0.11121698 0.04984793
[4,]  0.4 0.08633521 0.01577335 0.11725098 0.05541944
[5,]  0.5 0.08989036 0.01543300 0.12013903 0.05964168
[6,]  0.6 0.09275225 0.01509710 0.12234257 0.06316193
[7,]  0.7 0.09563246 0.01524897 0.12552045 0.06574448
[8,]  0.8 0.09971964 0.01643330 0.13192891 0.06751036
[9,]  0.9 0.10538962 0.02081095 0.14617909 0.06460016
> coef2mat
      [,1]      [,2]       [,3]      [,4]      [,5]
[1,]  0.1 0.1548719 0.02205991 0.1981093 0.1116345
[2,]  0.2 0.1656735 0.02101250 0.2068580 0.1244890
[3,]  0.3 0.1708705 0.02093900 0.2119109 0.1298300
[4,]  0.4 0.1758987 0.02234190 0.2196889 0.1321086
[5,]  0.5 0.1851403 0.02556903 0.2352556 0.1350250
[6,]  0.6 0.1991411 0.03041511 0.2587547 0.1395275
[7,]  0.7 0.2181969 0.03619383 0.2891369 0.1472570
[8,]  0.8 0.2490218 0.05148496 0.3499323 0.1481113
[9,]  0.9 0.3145634 0.07023116 0.4522164 0.1769103
###plot coef1,coef2 95%CI
par(mfrow = c(1, 2))
####coefficient1
upperbd1=coef1mat[,4]
lowerbd1=coef1mat[,5]
coef1=coef1mat[,2]
plot(taus,upperbd1,type = 'n',  ylim = c(0.00, 0.20),
     ylab = 'Coef1 Value', xlab = 'Quantile')
lines(taus, lowerbd1, col = 'grey')
lines(taus, upperbd1, col = 'grey')
polygon(c(taus, rev(taus)), c(upperbd1, rev(lowerbd1)),
     col = "lightblue", border = 'black')
points(taus, coef1,col='blue',pch=21)
lines(taus, coef1,col='blue')
title("gross Investment vs value of the firm")
####coefficient2
upperbd2=coef2mat[,4]
lowerbd2=coef2mat[,5]
coef2=coef2mat[,2]
plot(taus,upperbd2,type = 'n',  ylim = c(0.00, 0.50),
     ylab = 'Coef2 Capital', xlab = 'Quantile')
lines(taus, lowerbd2, col = 'grey')
lines(taus, upperbd2, col = 'grey')
polygon(c(taus, rev(taus)), c(upperbd2, rev(lowerbd2)),
     col = "lightblue", border = 'black')
points(taus, coef2,col='blue',pch=21)
lines(taus, coef2,col='blue')
title("gross Investment vs stock of plant and equipment")
          graph.jpeg

图片附件: graph.jpeg (2011-4-4 13:48:34, 87.81 KB)
https://bbs.pinggu.org/forum.php?mod=attachment&aid=ODgxMDY4fDBkYzI1M2I5fDE3NjY0NzQ0MTR8MHww

作者: ywh19860616    时间: 2011-4-4 17:25:52

epoh前辈，谢谢您
可以看出，您很认真负责
作者: ywh19860616    时间: 2011-4-4 23:08:46

epoh前辈，我运行这一段时，会提示错误
t=1.96
for (i in 1:ntaus){
#coefficient1
coef1mat[i,1]=taus
coef1mat[i,2]=coef[(2*(i-1)+1),1]
coef1mat[i,3]=coef[(2*(i-1)+1),2]
coef1mat[i,4]=coef[(2*(i-1)+1),1]+t*coef[(2*(i-1)+1),2]
coef1mat[i,5]=coef[(2*(i-1)+1),1]-t*coef[(2*(i-1)+1),2]
#coefficient2
coef2mat[i,1]=taus
coef2mat[i,2]=coef[2*i,1]
coef2mat[i,3]=coef[2*i,2]
coef2mat[i,4]=coef[2*i,1]+t*coef[2*i,2]
coef2mat[i,5]=coef[2*i,1]-t*coef[2*i,2]
}
coef1mat
coef2mat

错误于coef1mat[i, 1] = taus : 被替换的项目不是替换值长度的倍数

我哪里出错？
作者: epoh    时间: 2011-4-4 23:15:08

coef1mat[i,1]=taus[i]
coef2mat[i,1]=taus[i]
作者: ywh19860616    时间: 2011-4-4 23:48:35

额，epoh老师，已经可以了
谢谢您的指导
作者: ywh19860616    时间: 2011-4-8 18:32:28

epoh老师，又碰到疑问了
我看到一些文章运行出各分位点的系数后，如在0.1,0.2,0.3,0.4,0.5的系数值分别为coeff1，coeff2，coeff3，coeff4，coeff5，文章提到了通过统计量检验这些系数间是否存在显著差异。有些文章提到了系数的联合显著性检验（ie，有些文章提到：对FDI在不同分位点的各系数差异性的联合F检验），这些在R有相关命令吗？
附件中有一篇文章，其中我用红色标志了2个地方，您看在R中有相关命令可以实现吗？

附件: Quantile regression.pdf (2011-4-8 18:32:28, 477.97 KB)
https://bbs.pinggu.org/forum.php?mod=attachment&aid=ODgzNzkxfDY4MDE0M2FlfDE3NjY0NzQ0MTR8MHww
作者: epoh    时间: 2011-4-9 08:56:00

function anova.rq()

There are two (as yet) distinct forms of the test.
In the first the fitted objects all have the same specified quantile (tau) and
the intent is to test the hypothesis that smaller models are adequate
relative to the largest specified model.

In the second form of the test the linear predictor of the fits
are all the same, but the specified quantiles (taus) are different
作者: ywh19860616    时间: 2011-4-9 10:41:29

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（）函数
作者: epoh    时间: 2011-4-9 11:01:21

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

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

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

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

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的矩阵，为什么会提示说行不同呢？？请您指点一下
作者: epoh    时间: 2012-4-9 22:17:34

sun1018 发表于 2012-4-9 17:43
前辈，您好，在论坛里看到您是R语言的高手，我问一下，为什么欲行面板分位数回归时总是出问题：
> libra ...
嘿嘿,一年了,有点忘了
不知ywh兄,还记不记得?
作者: sun1018    时间: 2012-4-9 22:53:50

epoh 发表于 2012-4-9 22:17
嘿嘿,一年了,有点忘了
不知ywh兄,还记不记得?
呵呵，刚才仔细看了一下， 是W的定义有问题
作者: ywh19860616    时间: 2012-4-10 08:39:12

epoh 发表于 2012-4-9 22:17
嘿嘿,一年了,有点忘了
不知ywh兄,还记不记得?
呵呵，epoh老师，楼主自己解决了。
作者: ywh19860616    时间: 2012-4-10 16:46:27

sun1018 发表于 2012-4-9 22:53
呵呵，刚才仔细看了一下， 是W的定义有问题
DB file 是MicroTSP 的数据文件，在TSP中用fetch得到。
epoh老师，DB结尾的数据文件在TSP中如何打开呢？
好像是一个网页数据。
作者: epoh    时间: 2012-4-10 21:07:16

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
作者: ywh19860616    时间: 2012-4-10 21:20:18

epoh 发表于 2012-4-10 21:07
MicroTSP Databank Format
  http://econpy.googlecode.com/svn/trunk/pytrix/opendatabank.txt
呵呵，谢谢epoh老师，可以显示
我开始用了dbprint命令，一直显示错误。
作者: xiaomuyoucun    时间: 2013-3-15 12:06:28

epoh 发表于 2011-4-1 21:48
糊涂了,封掉a没发觉.
这问题,我再仔细看一下function rqss().
epoh 老师，我运行上边程序出现错误：

> xx1<- 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)
作者: ywh19860616    时间: 2013-3-15 12:18:45

xiaomuyoucun 发表于 2013-3-15 12:06
epoh 老师，我运行上边程序出现错误：

> xx1
哈哈，epoh老师太忙了
你可以直接利用rqpd 包了。
作者: xiaomuyoucun    时间: 2013-3-15 12:45:03

ywh19860616 发表于 2013-3-15 12:18
哈哈，epoh老师太忙了
你可以直接利用rqpd 包了。
这个程序之前你也用过，那请问怎么回事？是哪里出问题了，你能修改过来吗？
作者: xiaomuyoucun    时间: 2013-3-15 15:11:29

ywh19860616 发表于 2011-4-1 13:54
我试加了，程序也是可以运行的，只是结果不一样
> library(SparseM)
> library(quantreg)
你好，我发现这一段程序总是运行不过去，你之前运行过没有错误吗？为什么我的总出现以下错误？
急于求解：
> xx1<- 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)
作者: ywh19860616    时间: 2013-3-15 15:15:11

xiaomuyoucun 发表于 2013-3-15 15:11
你好，我发现这一段程序总是运行不过去，你之前运行过没有错误吗？为什么我的总出现以下错误？
急于求解 ...
你的最后目的是什么？运行面板分位回归模型吗？
如果答案是，那你直接利用rqpd包就可以了，作者
已经把这些程序打包为package了，可以直接用
作者: ywh19860616    时间: 2013-3-15 15:15:53

xiaomuyoucun 发表于 2013-3-15 15:11
你好，我发现这一段程序总是运行不过去，你之前运行过没有错误吗？为什么我的总出现以下错误？
急于求解 ...
你是为了运行面板分位回归吗？
如果是，那用rqpd包就可以了
作者: xiaomuyoucun    时间: 2013-3-15 15:18:53

ywh19860616 发表于 2013-3-15 15:15
你是为了运行面板分位回归吗？
如果是，那用rqpd包就可以了
我不是很清楚这个包rqpd怎么用。。主要是我还实现惩罚项和工具变量
作者: xiaomuyoucun    时间: 2013-3-15 15:20:00

ywh19860616 发表于 2013-3-15 15:15
你是为了运行面板分位回归吗？
如果是，那用rqpd包就可以了
还有我不清楚为什么我运行那一段程序时总汇报那个错误，你当时没遇到吗？怎么解决的？
作者: ywh19860616    时间: 2013-3-15 15:30:36

xiaomuyoucun 发表于 2013-3-15 15:18
我不是很清楚这个包rqpd怎么用。。主要是我还实现惩罚项和工具变量
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,as.matrix(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)
作者: bourn    时间: 2013-6-14 00:30:30

  1. 请问给出 x,y 向量 以及y 的置信区间的两个向量，怎样画出echo老师那样漂亮的图呢？
复制代码
  1. x        y        lower        upper
  2. 1        397.2529977        396.7593464        397.7656951
  3. 2        397.4132485        396.9200467        397.9268448
  4. 3        397.3247914        396.8306472        397.836503
  5. 4        397.2968677        396.8028951        397.8089225
  6. 5        397.3177442        396.8234853        397.8292263
  7. 6        397.2909879        396.7964316        397.8018753
  8. 7        397.2202784        396.7259736        397.7316689
  9. 8        397.235463        396.7412418        397.7470205
  10. 9        397.2530363        396.7594596        397.7658828
  11. 10        397.2721299        396.7783726        397.7846153
  12. 11        397.222693        396.7287409        397.7347888
  13. 12        397.2435759        396.7495669        397.755558
  14. 13        397.29418        396.7997861        397.8053922
  15. 14        397.3203146        396.8257034        397.8310922
  16. 15        397.2451708        396.7521961        397.7592215
  17. 16        397.2563634        396.76206        397.7677566
  18. 17        397.2783285        396.7830434        397.7877583
  19. 18        397.272673        396.7778679        397.7830628
  20. 19        397.2491882        396.7551349        397.7610815
  21. 20        397.2466776        396.7526789        397.7586802
  22. 21        397.2791679        396.7844136        397.7896593
  23. 22        397.2882613        396.7935994        397.7989374
  24. 23        397.2785999        396.7836264        397.788653
  25. 24        397.2863792        396.7910587        397.7957383
  26. 25        397.3181032        396.8227456        397.827388
  27. 26        397.3206856        396.8257823        397.830879
  28. 27        397.3097213        396.8150327        397.8203441
  29. 28        397.3164154        396.8214377        397.8264599
  30. 29        397.3211675        396.8256034        397.8300394
  31. 30        397.2947883        396.7991728        397.8035573
  32. 31        397.2645196        396.7691402        397.7737609
  33. 32        397.2447072        396.7488711        397.753035
  34. 33        397.2548759        396.7589193        397.7629626
  35. 34        397.3896628        396.8927863        397.8959098
  36. 35        397.3084707        396.812229        397.8159872
  37. 36        397.2391652        396.742792        397.7464189
  38. 37        397.2311395        396.7346531        397.7381667
  39. 38        397.2386479        396.7422905        397.7459331
  40. 39        397.2150245        396.7189073        397.72279
  41. 40        397.2281909        396.7319131        397.7356353
  42. 41        397.2544918        396.7581481        397.7618043
  43. 42        397.2654263        396.7688229        397.7722195
  44. 43        397.2322621        396.7357601        397.7392581
  45. 44        397.265934        396.7696922        397.7734505
  46. 45        397.3003763        396.8036283        397.8068803
  47. 46        397.3333101        396.8367572        397.8402043
  48. 47        397.319227        396.8229513        397.8266755
  49. 48        397.3179154        396.8212696        397.8246237
  50. 49        397.3280813        396.8312825        397.8344836
  51. 50        397.2956591        396.7989231        397.8021872
  52. 51        397.246898        396.7503246        397.7537513
  53. 52        397.2410562        396.7444577        397.7478592
  54. 53        397.2526465        396.7558576        397.7590688
  55. 54        397.2260762        396.7291005        397.7321248
  56. 55        397.2513019        396.7541653        397.7570288
  57. 56        397.2610098        396.7639057        397.7668017
  58. 57        397.2741034        396.7766502        397.779197
  59. 58        397.2057649        396.7092354        397.7127059
  60. 59        397.2083586        396.7113139        397.7142691
  61. 60        397.2182534        396.7214999        397.7247464
  62. 61        397.4224004        396.9254085        397.9284166
  63. 62        397.4729842        396.9764236        397.9798629
  64. 63        397.5078162        397.0109405        398.0140648
  65. 64        397.4774982        396.9803206        397.983143
  66. 65        397.4746745        396.9774455        397.9802166
  67. 66        397.4436107        396.9467194        397.9498281
  68. 67        397.4421092        396.945141        397.9481728
  69. 68        397.4114569        396.9145322        397.9176075
  70. 69        397.4490819        396.9515781        397.9540744
  71. 70        397.4590143        396.9615753        397.9641364
  72. 71        397.4454084        396.9481286        397.9508489
  73. 72        397.3980864        396.900859        397.9036316
  74. 73        397.3823058        396.8850411        397.8877765
  75. 74        397.3107738        396.8135396        397.8163054
  76. 75        397.2767014        396.7796459        397.7825905
  77. 76        397.2523147        396.7550692        397.7578237
  78. 77        397.2587444        396.7613372        397.76393
  79. 78        397.2467163        396.7491899        397.7516635
  80. 79        397.2175457        396.7202546        397.7229635
  81. 80        397.25351        396.7561282        397.7587463
  82. 81        397.2781211        396.7806535        397.7831859
  83. 82        397.3212792        396.8237367        397.8261942
  84. 83        397.3426944        396.8452583        397.8478223
  85. 84        397.3641306        396.8665626        397.8689945
  86. 85        397.3454709        396.8478673        397.8502637
  87. 86        397.3105994        396.8130274        397.8154554
  88. 87        397.2833158        396.7857399        397.788164
  89. 88        397.3099919        396.8124487        397.8149055
  90. 89        397.3490211        396.8512496        397.8534781
  91. 90        397.3495788        396.8519033        397.8542278
  92. 91        397.3195795        396.8217894        397.8239992
  93. 92        397.2802401        396.7824254        397.7846107
  94. 93        397.2897467        396.791869        397.7939913
  95. 94        397.2833429        396.7856354        397.7879278
  96. 95        397.2786531        396.7809443        397.7832355
  97. 96        397.2999597        396.8021047        397.8042496
  98. 97        397.3490376        396.8510623        397.853087
  99. 98        397.3478331        396.8496396        397.8514462
  100. 99        397.3179968        396.8198363        397.8216757
  101. 100        397.1811206        396.6835288        397.685937
  102. 101        397.2077541        396.7099375        397.7121208
  103. 102        397.1895242        396.6917288        397.6939334
  104. 103        397.1485445        396.6511544        397.6537643
  105. 104        397.0871233        396.5903143        397.5935053
  106. 105        397.1754044        396.6781595        397.6809145
  107. 106        397.1530638        396.6561618        397.6592599
  108. 107        397.1024263        396.6061229        397.6098196
  109. 108        397.1218832        396.6254609        397.6290385
  110. 109        397.197356        396.7002263        397.7030966
  111. 110        397.1807736        396.6843189        397.6878642
  112. 111        397.129608        396.6332436        397.6368792
  113. 112        397.1150715        396.6184655        397.6218595
  114. 113        397.1560094        396.6589613        397.6619132
  115. 114        397.1070776        396.6103967        397.6137158
  116. 115        396.9879452        396.4916884        397.4954316
  117. 116        397.0148937        396.517861        397.5208284
  118. 117        397.1301963        396.633186        397.6361758
  119. 118        397.1987119        396.7015834        397.704455
  120. 119        397.1686358        396.6712246        397.6738134
  121. 120        397.1553741        396.6579044        397.6604346
  122. 121        397.1698632        396.6723473        397.6748313
  123. 122        397.1512732        396.6536367        397.6560003
  124. 123        397.1824391        396.6846925        397.6869459
  125. 124        397.2089191        396.7112997        397.7136804
  126. 125        397.2534481        396.7559063        397.7583644
  127. 126        397.3367331        396.8395354        397.8423377
  128. 127        397.3094845        396.8119873        397.8144901
  129. 128        397.2581761        396.7603647        397.7625533
  130. 129        397.2520332        396.7542354        397.7564376
  131. 130        397.2119141        396.7142677        397.7166213
  132. 131        397.2067979        396.7092994        397.7118009
  133. 132        397.1784564        396.6810026        397.6835487
  134. 133        397.1958888        396.6983848        397.7008808
  135. 134        397.1194415        396.6220662        397.6246909
  136. 135        397.1583911        396.6610885        397.6637859
  137. 136        397.2219642        396.7248128        397.7276614
  138. 137        397.3493721        396.8524075        397.855443
  139. 138        397.492099        396.9959144        397.9997298
  140. 139        397.4629351        396.9663019        397.9696687
  141. 140        397.432991        396.9360319        397.9390727
  142. 141        397.3931297        396.8957882        397.8984466
  143. 142        397.3475364        396.8502059        397.8528755
  144. 143        397.3335259        396.8362897        397.8390536
  145. 144        397.3699547        396.8730708        397.8761869
  146. 145        397.3905989        396.8937217        397.8968444
  147. 146        397.3990867        396.9019562        397.9048257
  148. 147        397.3419722        396.8443836        397.8467951
  149. 148        397.3696852        396.8722114        397.8747376
  150. 149        397.469935        396.973281        397.976627
  151. 150        397.5179582        397.0215314        398.0251046
  152. 151        397.5139909        397.0171799        398.020369
  153. 152        397.4560409        396.9589295        397.9618182
  154. 153        397.5282006        397.0311914        398.0341822
  155. 154        397.5245302        397.0273726        398.0302149
  156. 155        397.4479811        396.9506993        397.9534174
  157. 156        397.4280954        396.9306465        397.9331977
  158. 157        397.4394616        396.9421023        397.9447431
  159. 158        397.4066816        396.9092261        397.9117707
  160. 159        397.3709748        396.8734705        397.8759661
  161. 160        397.3912829        396.8939729        397.896663
  162. 161        397.4104393        396.9128958        397.9153522
  163. 162        397.3027484        396.8053165        397.8078845
  164. 163        397.3125548        396.8150308        397.8175068
  165. 164        397.3665587        396.8691076        397.8716566
  166. 165        397.4108632        396.9133569        397.9158506
  167. 166        397.3595105        396.8617951        397.8640797
  168. 167        397.4838046        396.9866696        397.9895346
  169. 168        397.4605318        396.9631157        397.9656996
  170. 169        397.4601611        396.9626386        397.9651162
  171. 170        397.4458878        396.9484581        397.9510285
  172. 171        397.4281261        396.9305905        397.933055
  173. 172        397.5054504        397.0085066        398.0115628
  174. 173        397.503132        397.0059013        398.0086706
  175. 174        397.479221        396.9819389        397.9846567
  176. 175        397.4236581        396.9260314        397.9284048
  177. 176        397.3912633        396.8934211        397.8955788
  178. 177        397.4465535        396.9489514        397.9513493
  179. 178        397.4262581        396.9286375        397.931017
  180. 179        397.4029338        396.9053212        397.9077086
  181. 180        397.3956809        396.8981381        397.9005952
  182. 181        397.4132356        396.9156404        397.9180452
  183. 182        397.3444547        396.8465355        397.8486163
  184. 183        397.3184295        396.8204938        397.8225582
  185. 184        397.3398975        396.8425842        397.8452709
  186. 185        397.4724447        396.9755453        397.9786458
  187. 186        397.5229083        397.0258348        398.0287614
  188. 187        397.5275233        397.0303923        398.0332614
  189. 188        397.5404257        397.043108        398.0457904
  190. 189        397.5450314        397.0473685        398.0497056
  191. 190        397.4560584        396.9584698        397.9608812
  192. 191        397.6987633        397.2016411        398.204519
  193. 192        397.737445        397.2401059        398.2427668
  194. 193        397.7678697        397.2703277        398.2727857
  195. 194        397.7189267        397.2212082        398.2234897
  196. 195        397.6266206        397.1288835        398.1311464
  197. 196        397.5406846        397.0427143        398.0447439
  198. 197        397.4876274        396.9896626        397.9916979
  199. 198        397.4432281        396.9452424        397.9472568
  200. 199        397.4111675        396.9131983        397.9152292
  201. 200        397.3853098        396.8874406        397.8895715
  202. 201        397.4889912        396.9909884        397.9929855
  203. 202        397.4480621        396.949932        397.9518019
  204. 203        397.499156        397.0012971        398.0034382
  205. 204        397.489891        396.9919379        397.9939848
  206. 205        397.5452877        397.0475015        398.0497153
  207. 206        397.4986016        397.0007833        398.0029651
  208. 207        397.4161721        396.9184861        397.9208002
  209. 208        397.4061382        396.9080483        397.9099584
  210. 209        397.4088925        396.9106556        397.9124186
  211. 210        397.404952        396.906713        397.9084739
  212. 211        397.4092705        396.9110829        397.9128954
  213. 212        397.4360817        396.9379652        397.9398488
  214. 213        397.4205047        396.9223432        397.9241816
  215. 214        397.4789776        396.9807163        397.982455
  216. 215        397.548326        397.0502104        398.0520947
  217. 216        397.6038013        397.1057365        398.1076716
  218. 217        397.558127        397.0597945        398.061462
  219. 218        397.5074574        397.0091948        398.0109321
  220. 219        397.4732688        396.9748871        397.9765053
  221. 220        397.9415925        397.4435284        398.4454643
  222. 221        398.1893405        397.6921727        398.695005
  223. 222        398.0674215        397.5698176        398.5722137
  224. 223        397.9918473        397.4940857        398.496324
  225. 224        397.9513635        397.453546        398.4557285
  226. 225        397.8062134        397.308559        398.3109045
  227. 226        397.6839158        397.1859854        398.1880551
  228. 227        397.5746487        397.076672        398.0786954
  229. 228        397.5072468        397.0093306        398.0114143
  230. 229        397.4933871        396.995409        397.997431
  231. 230        397.4714165        396.973279        397.9751414
  232. 231        397.4304453        396.9324152        397.9343851
  233. 232        397.4259195        396.9279654        397.9300114
  234. 233        397.4715464        396.9734692        397.9753921
  235. 234        397.486836        396.9890613        397.9912867
  236. 235        397.4558612        396.9579095        397.9599578
  237. 236        397.432416        396.9345531        397.9366901
  238. 237        397.4413069        396.943221        397.9451352
  239. 238        397.4419463        396.9436519        397.9453575
  240. 239        397.4282293        396.9300512        397.9318731
  241. 240        397.3818228        396.883609        397.8853951
  242. 241        397.4668889        396.9686969        397.9705048
  243. 242        397.6081599        397.1106195        398.1130791
  244. 243        397.6218243        397.1240594        398.1262946
  245. 244        397.5699594        397.0718145        398.0736695
  246. 245        397.5108057        397.0126139        398.0144222
  247. 246        397.4521798        396.9540176        397.9558554
  248. 247        397.5037395        397.0057412        398.0077429
  249. 248        397.5831826        397.0852856        398.0873887
  250. 249        397.5731748        397.0750896        398.0770044
  251. 250        397.506805        397.0087782        398.0107514
  252. 251        397.4391656        396.9410331        397.9429005
  253. 252        397.4169983        396.918899        397.9207997
  254. 253        397.4293675        396.9311068        397.932846
  255. 254        397.3915341        396.8932238        397.8949135
  256. 255        397.3572664        396.8590205        397.8607747
  257. 256        397.3500254        396.8518307        397.853636
  258. 257        397.3588718        396.8605696        397.8622674
  259. 258        397.332651        396.8344723        397.8362935
  260. 259        397.2881051        396.7898351        397.791565
  261. 260        397.2440953        396.745938        397.7477808
  262. 261        397.3006066        396.8023837        397.8041608
  263. 262        397.3301263        396.8319513        397.8337762
  264. 263        397.3019143        396.8037238        397.8055333
  265. 264        397.2926148        396.7944405        397.7962661
  266. 265        397.3153468        396.817046        397.8187452
  267. 266        397.3116513        396.8132488        397.8148463
  268. 267        397.2691937        396.7708791        397.7725645
  269. 268        397.2410492        396.742925        397.7448008
  270. 269        397.2668801        396.7685205        397.7701609
  271. 270        397.2894216        396.7912626        397.7931037
  272. 271        397.2903285        396.7921227        397.7939168
  273. 272        397.2967892        396.7986293        397.8004694
  274. 273        397.3155564        396.8171151        397.8186738
  275. 274        397.2874543        396.7891201        397.7907859
  276. 275        397.2782818        396.7799625        397.7816433
  277. 276        397.3195522        396.8210977        397.8226431
  278. 277        397.3033504        396.8049387        397.8065269
  279. 278        397.2754763        396.7772506        397.7790249
  280. 279        397.2250845        396.7269991        397.7289138
  281. 280        397.2149708        396.7169374        397.7189039
  282. 281        397.2167358        396.7186125        397.7204893
  283. 282        397.2171779        396.7190406        397.7209033
  284. 283        397.2140945        396.7160088        397.717923
  285. 284        397.2972437        396.7987866        397.8003296
  286. 285        397.2496108        396.7512286        397.7528464
  287. 286        397.1988664        396.7007223        397.7025782
  288. 287        397.1852783        396.6871194        397.6889605
  289. 288        397.1872237        396.6890263        397.6908289
  290. 289        397.2026106        396.7043841        397.7061576
  291. 290        397.1698592        396.6718343        397.6738094
  292. 291        397.1428419        396.6447838        397.6467256
  293. 292        397.149435        396.6514206        397.6534062
  294. 293        397.1693066        396.6712003        397.6730939
  295. 294        397.1613481        396.6633321        397.6653161
  296. 295        397.0846946        396.5868988        397.5891031
  297. 296        397.0840387        396.5861205        397.5882023
  298. 297        397.1053373        396.6070853        397.6088333
  299. 298        397.0617547        396.5637695        397.5657842
  300. 299        397.0474309        396.5494804        397.5515299
  301. 300        397.0628593        396.5647915        397.5667236
  302. 301        397.0840289        396.5858348        397.5876407
  303. 302        397.0683808        396.5702397        397.5720987
  304. 303        397.0277807        396.5297898        397.531799
  305. 304        397.0572035        396.5592935        397.5613835
  306. 305        397.0884536        396.5904141        397.5923745
  307. 306        397.0789538        396.5809706        397.5829874
  308. 307        397.0556879        396.5575592        397.5594305
  309. 308        397.0574076        396.5592468        397.561086
  310. 309        397.0659106        396.5676502        397.5693899
  311. 310        397.0462162        396.547996        397.5497758
  312. 311        397.0119924        396.5139557        397.5159189
  313. 312        397.0270795        396.5290123        397.5309451
  314. 313        397.0657276        396.5675293        397.569331
  315. 314        397.0696473        396.5714637        397.5732801
  316. 315        397.1077863        396.6097573        397.6117283
  317. 316        397.155256        396.6574585        397.6596611
  318. 317        397.1791815        396.6813293        397.6834771
  319. 318        397.1744758        396.676389        397.6783022
  320. 319        397.1904868        396.6924652        397.6944435
  321. 320        397.1783364        396.68042        397.6825036
  322. 321        397.1844115        396.686387        397.6883625
  323. 322        397.1570135        396.6591109        397.6612084
  324. 323        397.1064275        396.6083975        397.6103676
  325. 324        397.1151357        396.6171195        397.6191033
  326. 325        397.1242998        396.6262582        397.6282166
  327. 326        397.1050804        396.6071167        397.609153
  328. 327        397.133569        396.6357531        397.6379373
  329. 328        397.1691505        396.6712959        397.6734414
  330. 329        397.2001518        396.7021351        397.7041185
  331. 330        397.1800827        396.6820869        397.684091
  332. 331        397.1520945        396.654148        397.6562015
  333. 332        397.1701847        396.6721398        397.674095
  334. 333        397.1801097        396.6819758        397.6838419
  335. 334        397.1652954        396.6673086        397.6693218
  336. 335        397.144548        396.6465505        397.6485529
  337. 336        397.1391778        396.6410287        397.6428796
  338. 337        397.2060993        396.707897        397.7096946
  339. 338        397.2319377        396.7337917        397.7356457
  340. 339        397.2192709        396.7212644        397.723258
  341. 340        397.2620889        396.7639827        397.7658766
  342. 341        397.4161831        396.9182588        397.9203345
  343. 342        397.5255857        397.0277235        398.0298613
  344. 343        397.5382468        397.0402045        398.0421623
  345. 344        397.561016        397.0627237        398.0644315
  346. 345        397.6144327        397.1161196        398.1178065
  347. 346        397.6287099        397.1304042        398.1320986
  348. 347        397.5923295        397.0940746        398.0958198
  349. 348        397.5402467        397.04187        398.0434932
  350. 349        397.4852901        396.9869568        397.9886236
  351. 350        397.4115255        396.9130836        397.9146418
  352. 351        397.3536509        396.8551866        397.8567222
  353. 352        397.3133127        396.8148921        397.8164716
  354. 353        397.3057489        396.8073025        397.808856
  355. 354        397.2925814        396.79414        397.7956986
  356. 355        397.2801407        396.7817555        397.7833702
  357. 356        397.2539467        396.7555739        397.757201
  358. 357        397.2335188        396.7351126        397.7367065
  359. 358        397.23334        396.7348029        397.7362659
  360. 359        397.2856904        396.7869334        397.7881764
  361. 360        397.3179414        396.8191596        397.8203778
  362. 361        397.2461574        396.747544        397.7489305
  363. 362        397.1972163        396.6987316        397.7002469
  364. 363        397.2026201        396.7040633        397.7055065
  365. 364        397.1786508        396.6803047        397.6819586
  366. 365        397.1669719        396.6684855        397.6699992
复制代码

作者: ywh19860616    时间: 2013-6-14 06:52:13

bourn 发表于 2013-6-14 00:30
看看plotCI函数
在plotrix包里面

作者: trier2006    时间: 2013-6-14 08:33:04

shenbaiseshatan 发表于 2011-1-14 14:12
5# trier2006  
现成的函数有类似的，比如Hmisc包的Dotplot
学习了，谢谢
作者: tangaibing    时间: 2013-11-16 22:20:19

现成的函数有类似的，比如Hmisc包的Dotplot
作者: 神释    时间: 2014-8-12 14:03:51

很不错
作者: youyuhaiwan    时间: 2016-4-11 21:11:35

shenbaiseshatan 发表于 2011-1-10 15:02
简单点可以实现类似功能的程序
如果要画的是回归出来后的系数的置信区间阴影图形呢？
作者: Evangeline213    时间: 2018-8-12 13:01:01

epoh 发表于 2011-3-26 11:18
function rq.fit.panel()
is now superseded by the use of the function rqss() in the R quantreg pack ...
请问怎么在ggplot里画出分位数回归线的置信区间呢？




欢迎光临 人大经济论坛 (https://bbs.pinggu.org/) Powered by Discuz! X2