[R] Confidence interval from resampling

Adriana Bejarano
2011-06-23T13:44:20Z

Dear R gurus,

I have the following code, but I still not know how to estimate and extract
confidence intervals (95%CI) from resampling.

Thanks!

~Adriana

#data
penta<-c(770,729,640,486,450,410,400,340,306,283,278,260,253,242,240,229,201,198,190,186,180,170,168,151,150,148,147,125,117,110,107,104,85,83,80,74,70,66,54,46,45,43,40,38,10)
x<-log(penta+1)
plot(ecdf(x), ylab="Probability", xlab="Concentration
(Ln+1)")

x.wei<-fitdistr(x,"weibull")
x.wei
     shape        scale
  6.7291685   5.3769965
 (0.7807718) (0.1254696)
xwei.shape <- x.wei$estimate[[1]]
xwei.scale <-  x.wei$estimate[[2]]

x.wei<-fitdistr(x,"weibull")
x.wei
xwei.shape <- x.wei$estimate[[1]]
xwei.scale <-  x.wei$estimate[[2]]
curve(pweibull(x, shape=xwei.shape, scale = xwei.scale,lower.tail=TRUE,
log.p=FALSE), add=TRUE,col='green',lwd=3)

#draw random numbers from a weibull distribution 100 times with
shape=xwei.shape, scale = xwei.scale
draw <- lapply(1:100, function(.x){
out<-rweibull(x, shape=xwei.shape, scale = xwei.scale)
})
newx<- data.frame(draw)

colnames(newx)<-paste("x", 1:100, sep = "")
newmat<-data.matrix(newx)

# matrix of coefficients
rownum=2
colnum=100
ResultMat<-matrix(NA, ncol=colnum, nrow=rownum)
rownum2=45
colnum2=100
ResultMat2<-matrix(NA, ncol=colnum2, nrow=rownum2)

#loop through each column in the source matrix
for (i in 1:100)
                {
        sel_col<-newmat[col(newmat)==i]
  {ResultMat[,i]<-coef(fitdistr(sel_col,"weibull"))}
                 xwei.shape<- ResultMat[1,i]
   xwei.scale<- ResultMat[2,i]
 curve(pweibull(x, shape=xwei.shape, scale=xwei.scale, lower.tail=TRUE,
log.p = FALSE), add=TRUE, col='grey',lwd=0.5)
 ResultMat2[,i]<-pweibull(x, shape=xwei.shape, scale
xwei.scale,lower.tail=TRUE, log.p=FALSE)
}

## convert dataframe to numeric
MatOut<- as.matrix(ResultMat2)
mode(MatOut) <- "numeric"

# initiate variables
mm<-ml<-mu<-rep(0,length(MatOut[,1]))

# mean and upper/lower quantiles
for(i in 1:length(MatOut[,1])){
 mm[i]<- mean(MatOut[i,])
 ml[i]<- quantile(MatOut[i,], 0.025, na.rm=TRUE)
 mu[i]<- quantile(MatOut[i,], 0.975, na.rm=TRUE)
}
#lines(x, mm, col="black")
lines(x, ml, col="blue", lwd=2)
lines(x, mu, col="blue", lwd=2)

	[[alternative HTML version deleted]]