## install ebci package
## if (!requireNamespace("remotes")) {
##   install.packages("remotes")
## }
## remotes::install_github("kolesarm/ebci")


rm(list=ls())

library(ebci)
library(lme4)

attach(sleepstudy)

S = unique(Subject)
k = length(S)
n = nrow(sleepstudy)

left = (0:8)*13 - 1
right = (0:8)*13 + 11


beta = rep(0,18)
se = rep(0,18)
for(i in 1:18){
    I = (1:n)[Subject == S[i]]
    x = Days[I]
    y = Reaction[I]
    out = lm(y ~ x)
    beta[i] = out$coef[2]
    se[i] = sqrt(vcov(out)[2,2])
    }

tmp = ebci(beta ~ 1,data.frame(beta=beta),se,alpha=.05)
hat = tmp$df$th_eb
left = hat - tmp$df$len_eb
right = hat + tmp$df$len_eb
print(hat)
print(cbind(left,right))

pdf("ConfidenceSleep.pdf")
plot(1:18,hat,pch=21,bg="lightpink",cex=1.6,ylim=c(min(left),max(right)),
                xlab="",ylab="",cex.axis=1.4)
for(i in 1:18){
     segments(i,left[i],i,right[i],col="purple",lwd=3)
     }
title(ylab=expression(beta[1]),line=2.5,cex=2.5)
title(xlab="Subjects",line=3,cex=2.5)    
dev.off()