Example: Simulation Study of elm vs. Bootstrap

## Define elm function and its service functions
f <- scan('http://www-stat.stanford.edu/~owen/empirical/el.S',
          what=list(x=''), sep='\n')$x
## Author used _ as assignment operator, not a good idea
## Otherwise could have used source('http:...') directly
f <- gsub('_', '<-', f)
source(textConnection(f))

library(boot)

elcl <- function(x, pl=FALSE) {
  guess <- quantile(x, c(.1,.9))
  mu <- seq(guess[1],guess[2],length=100)
  n <- length(mu)
  w <- sapply(mu, function(X) elm(x, X)$logelr)
  i <- which(w == max(w))
  cls <- c(approx(x=w[1:i], y=mu[1:i], xout=max(w)-3.84/2, rule=2)$y,
           approx(x=w[i:n], y=mu[i:n], xout=max(w)-3.84/2, rule=2)$y)
  if(pl) {
    plot(mu, w)
    abline(v=mean(x))
    abline(v=cls)
    abline(h=max(w)-3.84/2, lty=2)
  }
  cls
}
m <- 1000
bnames <- c('normal','basic','percent','bca')
cl <- array(NA, dim=c(m,1+length(bnames),2),
            dimnames=list(NULL,c('el',bnames),c('lower','upper')))
set.seed(1)
for(i in 1:m) {
  cat(i,'')
  x <- exp(rnorm(25))
  cl[i,'el',] <- elcl(x)
  ## clb[i,] <- smean.cl.boot(x, na.rm=FALSE)[-1]
  b <- boot(x, function(x,i) mean(x[i]), R=999)
  cb <- boot.ci(b)
  ## For each type of bootstrap confidence interval, fetch last
  ## two values in vector - these are the CLs; put into matrix
  cb <- t(sapply(cb[bnames], 
                 function(w) {
                  l <- length(w)
                  c(w[l-1],w[l])
                 }))
  cl[i,rownames(cb),] <- cb
}

mu <- exp(.5)

cat('\nMean CL lengths:\n\n')
apply(apply(cl,1:2,function(x)x[2]-x[1]), 2, mean)

cat('\nLeft Tail Areas:\n\n')
apply(cl[,,'lower'] > mu, 2, mean)

cat('\nRight Tail Areas:\n\n')
apply(cl[,,'upper'] < mu, 2, mean)

Mean CL lengths:

> apply(apply(cl,1:2,function(x)x[2]-x[1]), 2, mean)
      el   normal    basic  percent      bca 
1.385881 1.440112 1.416760 1.416760 1.640999 
> 
> cat('\nLeft Tail Areas:\n\n')

Left Tail Areas:

> apply(cl[,,'lower'] > mu, 2, mean)
     el  normal   basic percent     bca 
  0.024   0.003   0.002   0.008   0.026 
> 
> cat('\nRight Tail Areas:\n\n')

Right Tail Areas:

> apply(cl[,,'upper'] < mu, 2, mean)
     el  normal   basic percent     bca 
  0.097   0.148   0.172   0.131   0.090 

This topic: Main > WebHome > CE > EmpLike > EmpLikeEx1
Topic revision: 24 Oct 2004, FrankHarrell
 
This site is powered by FoswikiCopyright © 2013-2022 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding Vanderbilt Biostatistics Wiki? Send feedback