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
```
Topic revision: r2 - 24 Oct 2004, FrankHarrell

• Biostatistics Webs

Copyright © 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