如何画出这样的可信区间图

http://www.sr.bham.ac.uk/~ajrs/R/gallery/plot_cluster_z_vs_kT.txt
中的source("http://www.sr.bham.ac.uk/~ajrs/R/scripts/errorbar.R")
#--Add errorbars (using function in sourced errorbar.R script):
with(A, errorbar(z, 0, 0, Tx, -Tx.le, Tx.ue, bnd=F, col=col))        # bnd=F means input data are errors not confidence bounds

自己研究下

1. # Make up some dummy data
   
2. set.seed(1705190905)
   
3. 
   
4. dfDmmy <-
   
5.   merge(
   
6.     data.frame( quarter = c(1:4) ),
   
7.     data.frame( year = c(2006:2009) )
   
8.        )
   
9. 
   
10. dfDmmy$mid <-
    
11.   1.5*exp( -0.2*( seq.int( nrow(dfDmmy) ) - 1 )*
    
12.              exp( rnorm( nrow(dfDmmy), sd = 0.05 ) ) )
    
13. 
    
14. dfDmmy$upper <- with( dfDmmy, mid + runif( length(mid) )*mid/2 )
    
15. 
    
16. dfDmmy$lower <- with( dfDmmy, mid - runif( length(mid) )*mid )
    
17. 
    
18. # Set up plot
    
19. par(family = "serif")
    
20. 
    
21. plot(
    
22.   NA, type = "n",
    
23.   xlim = range( seq.int( nrow(dfDmmy) ) ), ylim = c(-0.5, 2.5),
    
24.   xlab = NA_character_, ylab = NA_character_, axes = FALSE
    
25.     )
    
26. 
    
27. # Axes
    
28. with(
    
29.   dfDmmy,
    
30.   {
    
31.     axis( side = 1, at = seq_along(quarter), labels = quarter,
    
32.           lwd = 0, lwd.ticks = 1 )
    
33.    
    
34.     lapply(
    
35.       unique(year),
    
36.       function(i)
    
37.         mtext( text = as.character(i), side = 1, line = 2,
    
38.                at = mean( which( year == i ) ), adj = 0.5, padj = 1 )
    
39.           )
    
40.   }
    
41.     )
    
42. 
    
43. abline( h = par("usr")[3] )
    
44. 
    
45. axis( side = 2, at = c(-1:5/2), labels = sprintf("%5.2f", -1:5/2),
    
46.       las = 1, lwd = 0, lwd.ticks = 1 )
    
47. 
    
48. abline( v = par("usr")[1] )
    
49. 
    
50. # Title
    
51. # mtext(
    
52. #   text = expression( paste( "Statistics with BCPNN methods, ", italic(IC) ) ),
    
53. #   side = 3, adj = 0, xpd = FALSE
    
54. #      )
    
55. # text(
    
56. #   x = par("usr")[1] - par("mar")[2]*par("cxy")[1],
    
57. #   y = par("usr")[4] + par("cxy")[2]/2,
    
58. #   labels = expression( paste( "Statistics with BCPNN methods, ", italic(IC) ) ),
    
59. #   adj = c(0, 0), xpd = NA
    
60. #     )
    
61. 
    
62. windowsFonts( SONG = windowsFont("TT SimSun") )
    
63. 
    
64. text(
    
65.   x = par("usr")[1] - par("mar")[2]*par("cxy")[1],
    
66.   y = par("usr")[4] + par("cxy")[2]/2,
    
67.   labels = "BCPNN",
    
68.   family = "serif", adj = c(0, 0), xpd = NA
    
69.     )
    
70. 
    
71. text(
    
72.   x = par("usr")[1] - par("mar")[2]*par("cxy")[1] + strwidth("BCPNN"),
    
73.   y = par("usr")[4] + par("cxy")[2]/2,
    
74.   labels = "\u{6CD5}\u{7EDF}\u{8BA1}\u{91CF}\u{FF0C}",
    
75.   family = "SONG", adj = c(0, 0), xpd = NA
    
76.     )
    
77. 
    
78. text(
    
79.   x =
    
80.     par("usr")[1] - par("mar")[2]*par("cxy")[1] +
    
81.     strwidth("BCPNN\u{6CD5}\u{7EDF}\u{8BA1}\u{91CF}\u{FF0C}") + par("cxy")[1],
    
82.   y = par("usr")[4] + par("cxy")[2]/2,
    
83.   labels = "IC",
    
84.   family = "serif", font = 3, adj = c(0, 0), xpd = NA
    
85.     )
    
86. 
    
87. # Horizontal reference line
    
88. abline( h = 0, lty = 3, col = "blue" )
    
89. 
    
90. # Plot the data
    
91. with(
    
92.   dfDmmy,
    
93.   {
    
94.     # Intervals
    
95.     dvX <- seq_along(mid)
    
96.    
    
97.     # mapply(
    
98.     #   lines,
    
99.     #   x =
    
100.     #     apply( cbind( seq_along(lower) - 0.1, seq_along(lower) + 0.1 ), 1, I ),
     
101.     #   y = apply( cbind( lower, lower ), 1, I )
     
102.     #       )
     
103.     mapply(
     
104.       segments,
     
105.       x0 = lapply(dvX - 0.1, I), x1 = lapply(dvX + 0.1, I),
     
106.       y0 = lapply(lower, I)
     
107.           )
     
108.    
     
109.    
     
110.     # mapply(
     
111.     #   lines,
     
112.     #   x =
     
113.     #     apply( cbind( seq_along(upper) - 0.1, seq_along(upper) + 0.1 ), 1, I ),
     
114.     #   y = apply( cbind( upper, upper ), 1, I )
     
115.     #       )
     
116.     mapply(
     
117.       segments,
     
118.       x0 = lapply(dvX - 0.1, I), x1 = lapply(dvX + 0.1, I),
     
119.       y0 = lapply(upper, I)
     
120.           )
     
121.    
     
122.    
     
123.     # mapply(
     
124.     #   lines,
     
125.     #   x = apply( cbind( seq_along(lower), seq_along(upper) ), 1, I ),
     
126.     #   y = apply( cbind( lower, upper ), 1, I )
     
127.     #       )
     
128.     mapply(
     
129.       segments,
     
130.       x0 = lapply(dvX, I), y0 = lapply(lower, I), y1 = lapply(upper, I)
     
131.           )
     
132. 
     
133.     # Points and lines
     
134.     points( mid, type = "o" )
     
135.   }
     
136.     )

复制代码

Output: