Draw confidence interval in R

Question

Can you explain me how, with R code, we draw coverage graph for the CI? For example, we have 40 samples (each one of size 20) and we have calculated the CI for each one of the 40. Our mission right now is to draw something like that:

Where each line is for vector of 20 samples, and the Y axis is for the actual interval.

Note: I have the "pivot" of each interval (the middle point).

Thanks!

Answer 1

You say you have 40 samples and only show 20. From the information given, it is not possible to know from which population distribution you are sampling [I used $\mathsf{Norm}(\mu = 0, \sigma=15)$] or what confidence level you are using [I used 95%]. Also, I used $m =40$ samples of size $n = 20.$

I used standard t confidence intervals of the form $$ \bar X \pm t^* S/\sqrt{n},$$ where $t^* = 2.093$ cuts 2.5% from the upper tail of Student's t distribution with $\nu = 20 - 1 = 19$ degrees of freedom. So the center of each interval is the sample mean. (You may be mis-using the term 'pivot'.)

On average, there will be 2 confidence intervals out of 40 that do not cover. For the seed chosen, there happen to be 38 CIs that cover and two that don't. (If you want a variety of plots, delete the set.seed statement before you make subsequent runs of the program.)

You can adjust constants on the second line of the R code shown below the figure to get exactly what you want. (The R code is entirely from the base of R, with no special libraries required.)

set.seed(1776)
m = 40;  n = 20;  mu = 0;  sigma = 15;  conf.level = .95
x = rnorm(m*n, mu, sigma)
MAT = matrix(x, nrow=m)  # m x n matrix: each row a sample of size n
x.bar = rowMeans(MAT);  s = apply(MAT, 1, sd)

t.crit = qt(1-(1-conf.level)/2, n-1)
LCL = x.bar - t.crit*s/sqrt(n);  UCL = x.bar + t.crit*s/sqrt(n)
cover = LCL < mu & UCL > mu
HI = max(UCL); LO = min(LCL)  # to set dimensions of the plot

plot(c(0,m+1), c(HI, LO), col="white", ylab="", xlab="Sample", main="", xaxs="i")
  abline(h=mu, col="green2")
  for(i in 1:m) {bar="blue"; if (cover[i]==F) {bar="red"}
    lines(c(i,i), c(UCL[i], LCL[i]), col=bar, lwd=2) }

.

i = 1:40; cbind(i, LCL, UCL, cover)
       i         LCL        UCL cover
 [1,]  1 -10.7907965  5.3863254     1
 [2,]  2  -5.5686896  6.8660868     1
 [3,]  3  -9.0946722  3.9352009     1
 [4,]  4  -3.0539188  9.8247816     1
...
[10,] 10  -9.0793201  3.1892928     1
[11,] 11  -5.5498471  9.8125370     1
[12,] 12  -4.3934771  4.6644040     1
[13,] 13   0.7405062 12.8408320     0   # Doesn't cover 0
[14,] 14  -8.6025527  6.0888690     1
[15,] 15  -9.2913948  8.3837346     1
[16,] 16 -12.5804007  3.4733748     1
[17,] 17  -9.9732150  2.4025221     1
[18,] 18  -9.0158737  6.1480452     1
[19,] 19  -9.4274926  8.9666786     1
...
[30,] 30  -8.8345204  8.4811912     1
[31,] 31  -8.5916045  3.1271369     1
[32,] 32  -6.7154571  6.0728731     1
[33,] 33  -8.5203468  5.4368537     1
[34,] 34 -12.0104057 -0.6229541     0   # Doesn't cover 0
[35,] 35  -0.6030757 13.2032218     1
[36,] 36  -7.0561399  5.3908935     1
[37,] 37  -5.5087385  8.1079198     1
[38,] 38  -1.6004350  9.5415809     1
[39,] 39  -6.2662541  7.9475938     1
[40,] 40  -0.4373327 14.7683446     1