Solved – reliable nonparametric confidence interval for the mean of a skewed distribution

bootstrapconfidence intervalmeanmediannonparametric

Very skewed distributions such as the log-normal do not result in accurate bootstrap confidence intervals. Here is an example showing that the left and right tail areas are far from the ideal 0.025 no matter which bootstrap method you try in R:

require(boot)
n    <- 25
B    <- 1000
nsim <- 1000
set.seed(1)
which <- c('basic', 'perc', 'norm', 'bca', 'stud')
mul <- 0; sdl <- 1.65   # on log scale
dist <- c('normal', 'lognormal')[2]
switch(dist, normal    = {g <- function(x) x; mu <- mul},
             lognormal = {g <- exp; mu <- exp(mul + sdl * sdl / 2)})
count <- matrix(0, nrow=length(which), ncol=2,
                dimnames=list(which, c('lower', 'upper')))
stat <- function(x, j) {
## See http://www.psychology.mcmaster.ca/bennett/boot09/percentileT.pdf
  x <- x[j]
  m <- mean(x)
  s <- sd(x)
  n <- length(x)
  sem <- s / sqrt(n)
  m.var <- sem ^ 2
  c(m, m.var)
}
for(i in 1 : nsim) {
  if(i %% 100 == 0) cat(i, '')
  x <- g(rnorm(n, mul, sdl))
  b  <- boot(x, stat, R=B)
  ci <- boot.ci(b, type=which)
  for(w in which) {
    nam <- switch(w, perc='percent', norm='normal', basic='basic',
                  stud='student', bca='bca')
    z <- rev(rev(ci[[nam]])[1:2])
    count[w, 'lower'] <- count[w, 'lower'] + (z[1] > mu)
    count[w, 'upper'] <- count[w, 'upper'] + (z[2] < mu)
  }
}
cat('\n')
count / nsim

The result is below:

      lower upper
basic 0.000 0.329
perc  0.003 0.257
norm  0.000 0.287
bca   0.015 0.185
stud  0.005 0.129

For $n=400$ single bootstraps still do not provide adequately accurate coverage:

      lower upper
basic 0.001 0.114
perc  0.005 0.093
norm  0.002 0.102
bca   0.017 0.067
stud  0.011 0.058

Empirical likelihood also fails to provide accurate confidence intervals when sampling from the lognormal distribution.

Is there a general-purpose approach out there that does not depend on knowing the distribution in advance? Has anyone tried to get confidence intervals for the mean by fitting the data to the Tukey generalized $\lambda$ distribution (this distribution is highly flexible)? What about using Kolmogorov-Smirnov confidence bands for the CDF? Would computing the mean on the upper and lower bounds on the CDF be horribly conservative? I would settle for some conservatism if a method has wide applicability.

To restate the goals, I'm seeking a generally applicable approach for getting a confidence interval for a population mean such that

the interval is asymmetric if the raw data distribution is asymmetric
the interval has correct coverage in both tails (e.g., 0.025 error probability in both)
the procedure does not require the analyst to specify anything about the underlying distribution or the transformation needed to make the distribution symmetric

Note that the central limit theorem is irrelevant here; I have a fixed small sample size and the confidence interval must be asymmetric to be accurate in both tails. The parametric $t$-based confidence interval under a lognormal model with $\mu=0, \sigma=1.65$ and $n=20000$ still has bad coverage (left tail error 0.012, right 0.047 when both should be 0.025).

In continuing to think about this there are two broad ways of conceptualizing the problem that I'd like to discuss.

The mean is not a quantity that lends itself to nonparametric inference, at least when exactness of inference is required. The sample median is meaningful for any continuous distribution and we have a simple exact confidence interval for the median. In a sample of size $n=20$ from a normal distribution the confidence interval for the median is $1.28 \times$ longer than the exact $t$-based confidence interval for the mean (see code below). Perhaps this factor of 1.28 is a reasonable price to pay for robustness and complete distributional freedom.
Even though no single bootstrap will give adequately accurate confidence limits for samples from extremely skewed distributions, the double bootstrap can significantly improve the confidence coverage in both tails. Nankervis has some nice results and provides an excellent computational algorithm. But no software I could find implements this.

R code illustrating 1. above:

## Exact CI for median from DescTools package SignTest.default
## See also ttp://www.stat.umn.edu/geyer/old03/5102/notes/rank.pdf,
## http://de.scribd.com/doc/75941305/Confidence-Interval-for-Median-Based-on-Sign-Test
cimed <- function(x, alpha=0.05, na.rm=FALSE) {
  if(na.rm) x <- x[! is.na(x)]
  n <- length(x)
  k <- qbinom(p=alpha / 2, size=n, prob=0.5, lower.tail=TRUE)
  ## Actual CL: 1 - 2 * pbinom(k - 1, size=n, prob=0.5) >= 1 - alpha
  sort(x)[c(k, n - k + 1)]
}

n <- 20
m <- 20000
cil <- cilt <- 0
z <- qt(0.975, n - 1)

for(i in 1 : m) {
  x <- rnorm(n)
  cil  <- cil + diff(cimed(x))
  cilt <- cilt + 2 * z * sqrt(var(x) / n)
}
cil  <- cil / m
cilt <- cilt / m

c(cil, cilt, cilt / cil, cil / cilt)

Best Answer

I am somewhat pessimistic about a such non-parametric method, at least without the introduction of some sort of constraints on the underlying distribution.

My reasoning for this is that there will always be a distribution that breaks the true coverage probability for any finite $n$ (although as $n \rightarrow \infty$, this distribution will become more and more pathological), or the confidence interval will have to be arbitrarily large.

To illustrate, you could imagine a distribution that looks like a normal up to some value $\alpha$, but after $\alpha$ becomes extremely right skewed. This can have unbounded influence on the distribution's mean and as you push $\alpha$ out as far as possible, this can have arbitrarily small probability of making it into your sample. So you can imagine that for any $n$, you could pick an $\alpha$ to be so large that all points in your sample have extremely high probability of looking like it comes from a normal distribution with mean = 0, sd = 1, but you can also have any true mean.

So if you're looking for proper asymptotic coverage, of course this can be achieved by the CLT. However, your question implies that you are (quite reasonably) interested in the finite coverage. As my example shows, there will always be a pathological case that ruins any finite length CI.

Now, you still could have a non-parametric CI that achieves good finite coverage by adding constraints to your distribution. For example, the log-concave constraint is a non-parametric constraint. However, it seems inadequate for your problem, as log-normal is not log-concave.

Perhaps to help illustrate how difficult your problem could be, I've done unpublished work on a different constraint: inverse convex (if you click on my profile, I have a link to a personal page that has a preprint). This constraint includes most, but not all log-normals. You can also see that for this constraint, the tails can be "arbitrarily heavy", i.e. for any inverse convex distribution up to some $\alpha$, you can have heavy enough tails that the mean will be as large as you like.

Related Solutions

Solved – Coverage probabilities of the basic bootstrap confidence Interval

The terminology is probably not used consistently, so the following is only how I understand the original question. From my understanding, the normal CIs you computed are not what was asked for. Each set of bootstrap replicates gives you one confidence interval, not many. The way to compute different CI-types from the results of a set of bootstrap replicates is as follows:

B    <- 999                  # number of replicates
muH0 <- 100                  # for generating data: true mean
sdH0 <- 40                   # for generating data: true sd
N    <- 200                  # sample size
DV   <- rnorm(N, muH0, sdH0) # simulated data: original sample

Since I want to compare the calculations against the results from package boot, I first define a function that will be called for each replicate. Its arguments are the original sample, and an index vector specifying the cases for a single replicate. It returns $M^{\star}$, the plug-in estimate for $\mu$, as well as $S_{M}^{2\star}$, the plug-in estimate for the variance of the mean $\sigma_{M}^{2}$. The latter will be required only for the bootstrap $t$-CI.

> getM <- function(orgDV, idx) {
+     bsM   <- mean(orgDV[idx])                       # M*
+     bsS2M <- (((N-1) / N) * var(orgDV[idx])) / N    # S^2*(M)
+     c(bsM, bsS2M)
+ }

> library(boot)                                       # for boot(), boot.ci()
> bOut <- boot(DV, statistic=getM, R=B)
> boot.ci(bOut, conf=0.95, type=c("basic", "perc", "norm", "stud"))
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 999 bootstrap replicates
CALL : 
boot.ci(boot.out = bOut, conf = 0.95, type = c("basic", "perc", "norm", "stud"))

Intervals : 
Level      Normal            Basic         Studentized        Percentile    
95%   ( 95.6, 106.0 )   ( 95.7, 106.2 )  ( 95.4, 106.2 )   ( 95.4, 106.0 )  
Calculations and Intervals on Original Scale

Without using package boot you can simply use replicate() to get a set of bootstrap replicates.

boots <- t(replicate(B, getM(DV, sample(seq(along=DV), replace=TRUE))))

But let's stick with the results from boot.ci() to have a reference.

boots   <- bOut$t                     # estimates from all replicates
M       <- mean(DV)                   # M from original sample
S2M     <- (((N-1)/N) * var(DV)) / N  # S^2(M) from original sample
Mstar   <- boots[ , 1]                # M* for each replicate
S2Mstar <- boots[ , 2]                # S^2*(M) for each replicate
biasM   <- mean(Mstar) - M            # bias of estimator M

The basic, percentile, and $t$-CI rely on the empirical distribution of bootstrap estimates. To get the $\alpha/2$ and $1 - \alpha/2$ quantiles, we find the corresponding indices to the sorted vector of bootstrap estimates (note that boot.ci() will do a more complicated interpolation to find the empirical quantiles when the indices are not natural numbers).

(idx <- trunc((B + 1) * c(0.05/2, 1 - 0.05/2)) # indices for sorted vector of estimates
[1] 25 975

> (ciBasic <- 2*M - sort(Mstar)[idx])          # basic CI
[1] 106.21826  95.65911

> (ciPerc <- sort(Mstar)[idx])                 # percentile CI
[1] 95.42188 105.98103

For the $t$-CI, we need the bootstrap $t^{\star}$ estimates to calculate the critical $t$-values. For the standard normal CI, the critical value will just be the $z$-value from the standard normal distribution.

# standard normal CI with bias correction
> zCrit   <- qnorm(c(0.025, 0.975))   # z-quantiles from std-normal distribution
> (ciNorm <- M - biasM + zCrit * sqrt(var(Mstar)))
[1] 95.5566 106.0043

> tStar <- (Mstar-M) / sqrt(S2Mstar)  # t*
> tCrit <- sort(tStar)[idx]           # t-quantiles from empirical t* distribution
> (ciT  <- M - tCrit * sqrt(S2M))     # studentized t-CI
[1] 106.20690  95.44878

In order to estimate the coverage probabilities of these CI-types, you will have to run this simulation many times. Just wrap the code into a function, return a list with the CI-results and run it with replicate() like demonstrated in this gist.

Solved – How to choose bootstrap confidence interval type from boot.ci in R

If your replicates are not normally distributed, do not choose $normal$.

$Basic$ can give you intervals that are out of the range of your replicated data; e.g. your bootstrapped replicates range between 2-200 but your lower confidence interval is -5.

For $student's$ CI, you need to pass a variance alongside whatever statistics (e.g. mean, median) you are dealing with. I would prefer this over $bca$ if you cannot generate a large number of replicates that can satisfy $bca$. If the number of replicates are small, the $bca$ intervals become unstable. One way of checking the stability is to generate many sets of replicates and identify corresponding confidence limits - most likely you will notice that the range of confidence limits based on $bca$ is wider than the rest.

I don't even know why they have the $percentile$ method, the most confusing of all five.

Best Answer

Related Solutions

Solved – Coverage probabilities of the basic bootstrap confidence Interval

Solved – How to choose bootstrap confidence interval type from boot.ci in R

Related Question