18
$\begingroup$

This is similar to Bootstrap: estimate is outside of confidence interval

I have some data that represents counts of genotypes in a population. I want to estimate genetic diversity using Shannon's index and also generate a confidence interval using bootstrapping. I've noticed, however, that the estimate via bootstrapping tends to be extremely biased and results in a confidence interval that lies outside of my observed statistic.

Below is an example.

# Shannon's index
H <- function(x){
  x <- x/sum(x)
  x <- -x * log(x, exp(1))
  return(sum(x, na.rm = TRUE))
}
# The version for bootstrapping
H.boot <- function(x, i){
  H(tabulate(x[i]))
}

Data generation

set.seed(5000)
X <- rmultinom(1, 100, prob = rep(1, 50))[, 1]

Calculation

H(X)

## [1] 3.67948

xi <- rep(1:length(X), X)
H.boot(xi)

## [1] 3.67948

library("boot")
types <- c("norm", "perc", "basic")
(boot.out <- boot::boot(xi, statistic = H.boot, R = 1000L))

## 
## CASE RESAMPLING BOOTSTRAP FOR CENSORED DATA
## 
## 
## Call:
## boot::boot(data = xi, statistic = H.boot, R = 1000)
## 
## 
## Bootstrap Statistics :
##     original     bias    std. error
## t1*  3.67948 -0.2456241  0.06363903

Generating the CIs with bias-correction

boot.ci(boot.out, type = types)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boot.out, type = types)
## 
## Intervals : 
## Level      Normal              Basic              Percentile     
## 95%   ( 3.800,  4.050 )   ( 3.810,  4.051 )   ( 3.308,  3.549 )  
## Calculations and Intervals on Original Scale

Assuming that the variance of t can be used for the variance of t0.

norm.ci(t0 = boot.out$t0, var.t0 = var(boot.out$t[, 1]))[-1]

## [1] 3.55475 3.80421

Would it be correct to report the CI centered around t0? Is there a better way to generate the bootstrap?

Improve this question
$\endgroup$

2 Answers 2

Reset to default
16
+25
$\begingroup$

In the setup given by the OP the parameter of interest is the Shannon entropy $$\theta(\mathbf{p}) = - \sum_{i = 1}^{50} p_i \log p_i,$$ which is a function of the probability vector $\mathbf{p} \in \mathbb{R}^{50}$. The estimator based on $n$ samples ($n = 100$ in the simulation) is the plug-in estimator $$\hat{\theta}_n = \theta(\hat{\mathbf{p}}_n) = - \sum_{i=1}^{50} \hat{p}_{n,i} \log \hat{p}_{n,i}.$$ The samples were generated using the uniform distribution for which the Shannon entropy is $\log(50) = 3.912.$ Since the Shannon entropy is maximized in the uniform distribution, the plug-in estimator must be downward biased. A simulation shows that $\mathrm{bias}(\hat{\theta}_{100}) \simeq -0.28$ whereas $\mathrm{bias}(\hat{\theta}_{500}) \simeq -0.05$. The plug-in estimator is consistent, but the $\Delta$-method does not apply for $\mathbf{p}$ being the uniform distribution, because the derivative of the Shannon entropy is 0. Thus for this particular choice of $\mathbf{p}$, confidence intervals based on asymptotic arguments are not obvious.

The percentile interval is based on the distribution of $\theta(\mathbf{p}_n^*)$ where $\mathbf{p}_n^*$ is the estimator obtained from sampling $n$ observations from $\hat{\mathbf{p}}_n$. Specifically, it is the interval from the 2.5% quantile to the 97.5% quantile for the distribution of $\theta(\mathbf{p}_n^*)$. As the OP's bootstrap simulation shows, $\theta(\mathbf{p}_n^*)$ is clearly also downward biased as an estimator of $\theta(\hat{\mathbf{p}}_n)$, which results in the percentile interval being completely wrong.

For the basic (and normal) interval, the roles of the quantiles are interchanged. This implies that the interval does seem to be reasonable (it covers 3.912), though intervals extending beyond 3.912 are not logically meaningful. Moreover, I don't know if the basic interval will have the correct coverage. Its justification is based on the following approximate distributional identity:

$$\theta(\mathbf{p}_n^*) - \theta(\hat{\mathbf{p}}_n) \overset{\mathcal{D}}{\simeq} \theta(\hat{\mathbf{p}}_n) - \theta(\mathbf{p}),$$ which might be questionable for (relatively) small $n$ like $n = 100$.

The OP's last suggestion of a standard error based interval $\theta(\hat{\mathbf{p}}_n) \pm 1.96\hat{\mathrm{se}}_n$ will not work either because of the large bias. It might work for a bias-corrected estimator, but then you first of all need correct standard errors for the bias-corrected estimator.

I would consider a likelihood interval based of the profile log-likelihood for $\theta(\mathbf{p})$. I'm afraid that I don't know any simple way to compute the profile log-likelihood for this example except that you need to maximize the log-likelihood over $\mathbf{p}$ for different fixed values of $\theta(\mathbf{p})$.

Improve this answer
$\endgroup$
7
  • 7
    $\begingroup$ The bias problem with using the "plug-in" estimator for entropy has been appreciated for decades. This paper analyzes less-biased estimates. A bias correction up to order $1/n$, which dates to 1955 (see eq. 4 of the linked paper), can be applied to the case presented by the OP. The correction is 0.245, almost identical to the bias identified by the bootstrap. Perhaps the bootstrap should be used here for estimating the entropy itself, not just its confidence limits. $\endgroup$
    – EdM
    Commented Jun 24, 2015 at 20:56
  • $\begingroup$ @EdM this is very useful information. I did not know the literature on this particular bias problem. It could be really useful if you could turn the comment into an answer that explains the bias correction and how it could be used with bootstrapping, say, to obtain confidence intervals. $\endgroup$
    – NRH
    Commented Jun 25, 2015 at 7:17
  • $\begingroup$ I didn't know this literature either, until this question and your answer came up. Which is somewhat embarrassing, since Shannon entropy is often used as a measure in my area of biomedical science. I'll see what I can put together as an additional answer. $\endgroup$
    – EdM
    Commented Jun 25, 2015 at 12:45
  • 1
    $\begingroup$ Increasing the number of bootstrap samples will not really help. It has to be large enough so that you can reliable estimate the quantities of interest for the distribution of $\theta(\mathbf{p}_n^*)$, say, but otherwise increasing the number of bootstrap samples will not remove the bias or make the confidence any more appropriate. $\endgroup$
    – NRH
    Commented Jun 25, 2015 at 21:23
  • 1
    $\begingroup$ Sorry ZNK, I misunderstood your question. If you increase the sample size $n$, the bias will be smaller, yes! The estimator is consistent. Precisely for the uniform distribution I would be somewhat sceptical about the actual coverage of the confidence intervals even for large $n$ for the reasons I described in the answer. For all other distributions the CLT applies, and the different methods will produce asymptotically correct coverage for $n \to \infty$. $\endgroup$
    – NRH
    Commented Jun 26, 2015 at 7:14
13
$\begingroup$

As the answer by @NRH points out, the problem is not that the bootstrapping gave a biased result. It's that the simple "plug in" estimate of the Shannon entropy, based on data from a sample, is biased downward from the true population value.

This problem was recognized in the 1950s, within a few years of the definition of this index. This paper discusses the underlying issues, with references to associated literature.

The problem arises from the nonlinear relation of the individual probabilities to this entropy measure. In this case, the observed genotype fraction for gene i in sample n, $\hat{p}_{n,i}$, is an unbiased estimator of the true probability, $p_{n,i}$. But when that observed value is applied to the "plug in" formula for entropy over M genes:

$$\hat{\theta}_n = \theta(\hat{\mathbf{p}}_n) = - \sum_{i=1}^{M} \hat{p}_{n,i} \log \hat{p}_{n,i}.$$

the non-linear relation means that the resulting value is a biased under-estimate of the true genetic diversity.

The bias depends on the number of genes, $M$ and the number of observations, $N$. To first order, the plug-in estimate will be lower than the true entropy by an amount $(M -1)/2N$. Higher order corrections are evaluated in the paper linked above.

There are packages in R that deal with this issue. The simboot package in particular has a function estShannonf that makes these bias corrections, and a function sbdiv for calculating confidence intervals. It will be better to use such established open-source tools for your analysis rather than try to start over from scratch.

Improve this answer
$\endgroup$
3
  • $\begingroup$ So the estimator in and of itself is erroneous due to sample size? The simboot package looks promising, but doesn't seem suited for my purposes as it need a control sample to estimate confidence intervals. $\endgroup$
    – ZNK
    Commented Jun 25, 2015 at 23:58
  • 1
    $\begingroup$ "Erroneous" isn't quite right; the estimator is "biased" in that its expected value is not the same as the actual population value. That doesn't mean it's "erroneous"; biased estimators can be useful, as illustrated by the bias-variance tradeoff in selecting estimators. If simboot doesn't meet your needs, Google "shannon entropy bias r" for links to other R packages like entropy, entropart, and EntropyEstimation. $\endgroup$
    – EdM
    Commented Jun 26, 2015 at 13:30
  • 1
    $\begingroup$ There are additional issues arising from the fact that some genotypes present in the population are likely to be missed in any particular sample. Some of the population- and ecology-based R packages appear to have ways to deal with this problem. $\endgroup$
    – EdM
    Commented Jun 26, 2015 at 13:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.