带置信区间的画图问题

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

我想使用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]
  怎么画如下图形？请大家指点，谢谢

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏10 回帖

关键词：画图问题 置信区间 分位数回归 plot 分位数 画图 置信区间

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

简单点可以实现类似功能的程序

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

复制代码

把置信带下面的数据和上面的数据分别用两个数组y1,y2存储，再用lines命令就可以了，
plot(x,y,type="l")
lines(x,y1)
lines(x,y2)

嗯，大家的建议本质好像都是计算出区间值，然后画上去。似乎没有现成的函数用。

5# trier2006
现成的函数有类似的，比如Hmisc包的Dotplot

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

##epoh ,哈哈，谢谢啊
我看过rqss（）函数帮助，但是找不到识别面板结构的程序，还不知道怎么使用。

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

#（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 )