[R] Confidence interval from resampling

[David Winsemius]

On Jun 23, 2011, at 9:44 AM, Adriana Bejarano wrote:

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

If you have a distribution of values, say "resamp.stat", of a  
statistic from a properly performed resampling operation you can  
extract and display easily the 5th and 95th percentiles.

CI.stat <- quantile(resamp.stat, c(0.05, 0.95) )
CI.stat

Note: I do not think that 100 replications would generally be  
sufficient for final work, although its probably acceptable for testing.

Your code as posted initially threw a bunch of errors since you did  
not include library(MASS), but fixing that fairly obvious problem  
shows that you can draw a nice plot. However, it remains unclear what  
statistic of what distribution you desire to assess. Mean, median, ...  
of what?

I do not think the error that arose on my machine from the wrapped  
text here:

#draw random numbers from a weibull distribution 100 times with
  ... shape=xwei.shape, scale = xwei.scale   -> error

..... was causing any problem, but there were a bunch of warnings that  
ought to be investigated:

Right after the loop I see ten of these:
Warning messages:
1: In dweibull(x, shape, scale, log) : NaNs produced

-- 
David

> 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]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

David Winsemius, MD
West Hartford, CT