# illustration of using boot() in the boot library to calculate
#  a bootstrap confidence interval for a ROS estimate of the mean

library(boot)    # need to install the boot library
library(NADA)    # and NADA to use ros()
# illustration data is the heron data, which seem to be logN

heron <- read.table('c:/philip/stat 505/data/heron.txt',as.is=T,header=T)

# to use boot() need to write a function that takes two arguments
#   original data, and a vector of indices for bootstrap sample
#   and returns the desired statistic
# since the value and censoring indicator 'go together', it is
#  easiest to pass the data as a matrix with 2 columns
# the first column is the observation; the second is the censoring flag

ros.mean <- function(xc,i) {
# get the values
  x <- xc[,1]
# get the censoring indicator
  cens <- xc[,2]

# then use the indices to construct the bootstrap data set
  xb <- x[i]
  cb <- cens[i]
  
# and use ros to get the mean
# since this is the value that is 'printed', 
#  it is the value of the function

  mean(ros(xb,cb))
  }
  
# test our function:
#  use the original data with indices from 1 to last obs
n <- length(heron$Pb)
print(ros.mean(cbind(heron$Pb,heron$Cens),1:n))

#  print() was used to make sure R prints the results, even if
#   you source the file

# now do the bootstrap
#  minimum arguments to boot() are:
#    the data (a matrix with col 1 = values, col2 = censoring flag)
#    the name of the function that computes the statistic
#    the number of bootstrap replicates

heron.boot <- boot(cbind(heron$Pb,heron$Cens), ros.mean, 1000)

# get confidence intervals from boot.ci()
print(boot.ci(heron.boot))