¿Bootstrap parcial: está bien que el centro de la CI en la estadística observada?

Question

Esto es similar a la de Bootstrap: estimación está fuera del intervalo de confianza

Tengo algunos datos que representa los recuentos de los genotipos en una población. Yo queremos estimar la diversidad genética utilizando de Shannon índice y también generar un intervalo de confianza utilizando el arranque. Me he dado cuenta, sin embargo, que la estimación a través de bootstrapping tiende a ser extremadamente y sesgada de los resultados en un intervalo de confianza que se encuentra fuera de mi observado de la estadística.

A continuación es un ejemplo.

# 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]))
}

La generación de datos

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

Cálculo

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

La generación de la Cei con el sesgo de corrección

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

Suponiendo que la varianza de t puede ser utilizado para la varianza de t0.

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

## [1] 3.55475 3.80421

¿Sería correcto informe de la CI centrado alrededor de t0? Hay un mejor manera de generar la secuencia de arranque?

Answer 1

En la configuración dada por el OP el parámetro de interés es la entropía de Shannon $$\theta(\mathbf{p}) = - \sum_{i = 1}^{50} p_i \log p_i,$$ que es una función de la probabilidad de vectores $\mathbf{p} \in \mathbb{R}^{50}$. El estimador basado en $n$ de las muestras ($n = 100$ en la simulación) es el plug-in estimador de $$\hat{\theta}_n = \theta(\hat{\mathbf{p}}_n) = - \sum_{i=1}^{50} \hat{p}_{n,i} \log \hat{p}_{n,i}.$$ Las muestras fueron generados utilizando la distribución uniforme para que la entropía de Shannon es $\log(50) = 3.912.$ Ya que la entropía de Shannon es maximizada en la distribución uniforme, el plug-in estimador debe ser sesgada a la baja. Una simulación muestra que $\mathrm{bias}(\hat{\theta}_{100}) \simeq -0.28$, mientras que $\mathrm{bias}(\hat{\theta}_{500}) \simeq -0.05$. El plug-in estimador es consistente, pero el $\Delta$-método no se aplica para $\mathbf{p}$ la distribución uniforme, debido a que la derivada de la entropía de Shannon es de 0. Por tanto, para esta elección particular de $\mathbf{p}$, intervalo de confianza asintótico basado en argumentos que no son evidentes.

El percentil intervalo se basa en la distribución de $\theta(\mathbf{p}_n^*)$ donde $\mathbf{p}_n^*$ es el estimador obtenido a partir de un muestreo $n$ observaciones $\hat{\mathbf{p}}_n$. Específicamente, es el intervalo desde el 2,5% de los cuantiles para el 97.5% de los cuantiles de la distribución de $\theta(\mathbf{p}_n^*)$. Como el OP bootstrap de simulación muestra, $\theta(\mathbf{p}_n^*)$ es claramente a la baja sesgada como un estimador de $\theta(\hat{\mathbf{p}}_n)$, lo que resulta en el percentil intervalo de estar completamente equivocado.

Para lo básico (y normal) el intervalo, la función de los cuantiles son intercambiados. Esto implica que el intervalo no parece ser razonable (que cubre 3.912), a pesar de que los intervalos se extiende más allá de 3.912 no son lógicamente significativo. Por otra parte, no sé si el basic intervalo tendrá la correcta cobertura. Su justificación se basa en los siguientes aproximada de la distribución de identidad:

$$\theta(\mathbf{p}_n^*) - \theta(\hat{\mathbf{p}}_n) \overset{\mathcal{D}}{\simeq} \theta(\hat{\mathbf{p}}_n) - \theta(\mathbf{p}),$$ que podría ser cuestionable para (relativamente) pequeño $n$ como $n = 100$.

El OP de la última sugerencia de un error estándar intervalo basado en $\theta(\hat{\mathbf{p}}_n) \pm 1.96\hat{\mathrm{se}}_n$ no funciona bien debido a la gran sesgo. Se podría trabajar para una corrección de sesgo del estimador, pero entonces lo primero de todo necesitas corregir los errores estándar para la corrección de sesgo del estimador.

Me volvería a considerar un intervalo de probabilidad basada en el perfil de la log-verosimilitud para $\theta(\mathbf{p})$. Me temo que no sé de alguna forma más sencilla para calcular el perfil de la log-verosimilitud para este ejemplo, excepto que usted necesita para maximizar el logaritmo de la probabilidad de más de $\mathbf{p}$ para diferentes valores fijos de $\theta(\mathbf{p})$.

Answer 2

Como la respuesta @NRH señala, el problema no es que el bootstrap dio una visión sesgada resultado. Es que el simple "plug in" estimación de la entropía de Shannon, basado en datos de una muestra, está sesgada a la baja desde el verdadero valor poblacional.

Este problema fue reconocido en la década de 1950, dentro de un par de años de la definición de este índice. Este artículo describe los problemas de fondo, con referencias a la literatura relacionada.

El problema surge de la no lineal de la relación individual de probabilidades a esta medida de la entropía. En este caso, la observación de genotipo fracción de genes i en la muestra n, $\hat{p}_{n,i}$, es un estimador imparcial de la verdadera probabilidad, $p_{n,i}$. Pero cuando el valor observado se aplica a la "enchufe" la fórmula de la entropía sobre 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$. De orden superior correcciones son evaluados en el documento enlazado más arriba.

Hay paquetes en R que lidiar con este problema. El simboot paquete en particular, tiene una función estShannonf que hace que estos sesgos correcciones, y una función sbdiv para el cálculo de intervalos de confianza. Será mejor utilizar tales establecido herramientas de código abierto para su análisis, en lugar de tratar de empezar de cero.