Bootstrap sesgado: ¿está bien centrar el IC alrededor de la estadística observada?

Tengo algunos datos que representan recuentos de genotipos en una población. Quiero estimar la diversidad genética usando el índice de Shannon y también generar un intervalo de confianza usando bootstrapping. Sin embargo, he notado que la estimación mediante bootstrapping tiende a ser extremadamente sesgada y da como resultado un intervalo de confianza que se encuentra fuera de mi estadística observada.

A continuación se muestra 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]))
}

Generacion 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

Generando los ICs con corrección de sesgos

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 usarse 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 informar el IC centrado en t0 ? ¿Hay una mejor manera de generar el bootstrap?

r confidence-interval bootstrap bias diversity ZNK
fuente

Respuestas:

En la configuración dada por el OP, el parámetro de interés es la entropía de Shannon que es una función del vector de probabilidad . El estimador basado en muestras ( en la simulación) es el estimador de complemento Las muestras se generaron usando la distribución uniforme para la cual la entropía de Shannon esDado que la entropía de Shannon se maximiza en la distribución uniforme, el estimador del complemento debe estar sesgado hacia abajo . Una simulación muestra que

θ (p) = - \sum_{i = 1}^{50} p_{i} \log p_{i},

$\theta(\mathbf{p}) = - \sum_{i = 1}^{50} p_i \log p_i,$

p \in R^{50}

$\mathbf{p} \in \mathbb{R}^{50}$

n

$n$

n = 100

$n = 100$

{\hat{θ}}_{n} = θ ({\hat{p}}_{n}) = - \sum_{i = 1}^{50} {\hat{p}}_{n, i} \log {\hat{p}}_{n, i} .

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

\log (50) = 3.912.

$\log(50) = 3.912.$

b i a s ({\hat{θ}}_{100}) ≃ - 0.28

$\mathrm{bias}(\hat{\theta}_{100}) \simeq -0.28$ mientras que . El estimador del complemento es coherente, pero el método no se aplica a siendo la distribución uniforme, porque la derivada de la entropía de Shannon es 0. Por lo tanto, para esta elección particular de , Los intervalos de confianza basados en argumentos asintóticos no son obvios.

b i a s ({\hat{θ}}_{500}) ≃ - 0.05

$\mathrm{bias}(\hat{\theta}_{500}) \simeq -0.05$

Δ

$\Delta$

p

$\mathbf{p}$

p

$\mathbf{p}$

El intervalo de percentil se basa en la distribución de donde es el estimador obtenido del muestreo de observaciones de . Específicamente, es el intervalo desde el cuantil del 2.5% hasta el cuantil del 97.5% para la distribución de . Como muestra la simulación de arranque del OP, también está claramente sesgada hacia abajo como un estimador de , lo que da como resultado que el intervalo de percentil sea Completamente mal. $\theta(\mathbf{p}_n^*)$ $\mathbf{p}_n^*$ $n$ $\hat{\mathbf{p}}_n$ $\theta(\mathbf{p}_n^*)$ $\theta(\mathbf{p}_n^*)$ $\theta(\hat{\mathbf{p}}_n)$

Para el intervalo básico (y normal), los roles de los cuantiles se intercambian. Esto implica que el intervalo parece ser razonable (cubre 3.912), aunque los intervalos que se extienden más allá de 3.912 no son lógicamente significativos. Además, no sé si el intervalo básico tendrá la cobertura correcta. Su justificación se basa en la siguiente identidad distribucional aproximada:

θ (p_{n}^{*}) - θ ({\hat{p}}_{n}) \overset{D}{≃} θ ({\hat{p}}_{n}) - θ (p),

$\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ños como .

n

$n$

n = 100

$n = 100$

Última propuesta de la OP de un error estándar basa intervalo de será no trabajar ya sea por el sesgo grande. Es posible que funcione para un estimador con corrección de sesgo, pero antes que nada necesita errores estándar correctos para el estimador con corrección de sesgo. $\theta(\hat{\mathbf{p}}_n) \pm 1.96\hat{\mathrm{se}}_n$

Consideraría un intervalo de probabilidad basado en el perfil de probabilidad de registro para . Me temo que no conozco ninguna forma simple de calcular la probabilidad de registro del perfil para este ejemplo, excepto que necesita maximizar la probabilidad de registro sobre para diferentes valores fijos de . $\theta(\mathbf{p})$ $\mathbf{p}$ $\theta(\mathbf{p})$

NRH
fuente

El problema del sesgo con el uso del estimador "plug-in" para la entropía se ha apreciado durante décadas. Este artículo analiza estimaciones menos sesgadas. Se puede aplicar una corrección de sesgo hasta el orden , que data de 1955 (véase la ecuación 4 del documento vinculado), al caso presentado por el OP. La corrección es 0.245, casi idéntica al sesgo identificado por el bootstrap. Quizás el bootstrap debería usarse aquí para estimar la entropía en sí, no solo sus límites de confianza.

1 / n

$1/n$

EdM

@EdM esta es una información muy útil. No conocía la literatura sobre este problema de sesgo en particular. Podría ser realmente útil si pudiera convertir el comentario en una respuesta que explique la corrección del sesgo y cómo podría usarse con bootstrapping, por ejemplo, para obtener intervalos de confianza.

NRH

Tampoco conocía esta literatura, hasta que surgió esta pregunta y su respuesta. Lo cual es algo vergonzoso, ya que la entropía de Shannon se usa a menudo como medida en mi área de ciencias biomédicas. Veré lo que puedo armar como una respuesta adicional.

EdM

Aumentar el número de muestras de bootstrap realmente no ayudará. Tiene que ser lo suficientemente grande para que pueda estimar de manera confiable las cantidades de interés para la distribución de , por ejemplo, pero aumentar el número de muestras de bootstrap no eliminará el sesgo ni hará La confianza más apropiada.

θ (p_{n}^{*})

$\theta(\mathbf{p}_n^*)$

NRH

Lo siento ZNK, entendí mal tu pregunta. Si aumenta el tamaño de la muestra , el sesgo será menor, ¡sí! El estimador es consistente. Precisamente para la distribución uniforme, sería un tanto escéptico sobre la cobertura real de los intervalos de confianza, incluso para grandes por las razones que describí en la respuesta. Para todas las demás distribuciones, se aplica el CLT, y los diferentes métodos producirán una cobertura asintóticamente correcta para .

n

$n$

n

$n$

n \to \infty

$n \to \infty$

NRH

Como señala la respuesta de @NRH, el problema no es que el arranque haya dado un resultado sesgado. Es que la simple estimación de "enchufable" de la entropía de Shannon, basada en datos de una muestra, está sesgada hacia abajo del valor real de la población.

Este problema fue reconocido en la década de 1950, a los pocos años de la definición de este índice. Este artículo discute los problemas subyacentes, con referencias a la literatura asociada.

El problema surge de la relación no lineal de las probabilidades individuales con esta medida de entropía. En este caso, la fracción de genotipo observada para el gen i en la muestra n , , es un estimador imparcial de la probabilidad real, . Pero cuando ese valor observado se aplica a la fórmula "plug in" para la entropía sobre los genes M: $\hat{p}_{n,i}$ $p_{n,i}$

{\hat{θ}}_{n} = θ ({\hat{p}}_{n}) = - \sum_{i = 1}^{M} {\hat{p}}_{n, i} \log {\hat{p}}_{n, i} .

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

la relación no lineal significa que el valor resultante es una subestimación sesgada de la verdadera diversidad genética.

El sesgo depende del número de genes, y el número de observaciones, . En primer orden, la estimación del complemento será menor que la entropía verdadera en una cantidad . Las correcciones de orden superior se evalúan en el documento vinculado anteriormente. $M$ $N$ $(M -1)/2N$

Hay paquetes en R que se ocupan de este problema. El simbootpaquete en particular tiene una función estShannonfque realiza estas correcciones de sesgo y una función sbdivpara calcular los intervalos de confianza. Será mejor usar herramientas de código abierto establecidas para su análisis en lugar de intentar comenzar desde cero.

EdM
fuente

Entonces, ¿el estimador en sí mismo es erróneo debido al tamaño de la muestra? El simbootpaquete parece prometedor, pero no parece adecuado para mis propósitos, ya que necesita una muestra de control para estimar los intervalos de confianza.

ZNK

"Erroneous" no está del todo bien; el estimador está "sesgado" porque su valor esperado no es el mismo que el valor real de la población. Eso no significa que sea "erróneo"; los estimadores sesgados pueden ser útiles, como lo ilustra la compensación de la varianza de sesgo en la selección de estimadores. Si simbootno cumple con sus necesidades, Google "entropía de Shannon sesgo r" para los enlaces a otros paquetes como R entropy, entroparty EntropyEstimation.

EdM

Hay problemas adicionales que surgen del hecho de que algunos genotipos presentes en la población probablemente se perderán en cualquier muestra en particular. Algunos de los paquetes R basados en la población y la ecología parecen tener formas de abordar este problema.

EdM