Voreingenommener Bootstrap: Ist es in Ordnung, das CI um die beobachtete Statistik zu zentrieren?

Dies ähnelt Bootstrap: Die Schätzung liegt außerhalb des Konfidenzintervalls

Ich habe einige Daten, die die Anzahl der Genotypen in einer Population darstellen. Ich möchte die genetische Vielfalt mithilfe des Shannon-Index abschätzen und mithilfe von Bootstrapping ein Konfidenzintervall generieren. Ich habe jedoch festgestellt, dass die Schätzung über Bootstrapping tendenziell extrem voreingenommen ist und zu einem Konfidenzintervall führt, das außerhalb meiner beobachteten Statistik liegt.

Unten ist ein Beispiel.

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

Datengenerierung

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

Berechnung

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

Generieren der CIs mit Bias-Korrektur

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

Angenommen, die Varianz von t kann für die Varianz von t0 verwendet werden .

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

## [1] 3.55475 3.80421

Wäre es richtig, das CI um t0 herum zu melden ? Gibt es eine bessere Möglichkeit, den Bootstrap zu generieren?

r confidence-interval bootstrap bias diversity ZNK
quelle

Antworten:

In dem vom OP gegebenen Aufbau ist der interessierende Parameter die Shannon-Entropie die eine Funktion des Wahrscheinlichkeitsvektors . Die Schätzfunktion basierend auf

θ (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}$

- Proben (

in der Simulation) ist die Plug-inSchätzer

n

$n$

n = 100

$n = 100$

Die Proben wurden unter Verwendung der Gleichverteilung erzeugt, für die die Shannon-Entropie

Da die Shannon-Entropie in der Gleichverteilung maximiert ist, muss der Plug-In-Schätzernach unten vorgespannt sein. Eine Simulation zeigtdass

, während

{\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$

Konfidenzintervalle, die auf asymptotischen Argumenten basieren, nicht offensichtlich . . Der Plug-In-Schätzer ist konsistent, aber die

Methode gilt nicht für

als Gleichverteilung, da die Ableitung der Shannon-Entropie 0 ist

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

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

Δ

$\Delta$

p

$\mathbf{p}$

p

$\mathbf{p}$

Die Perzentil - Intervall wird auf der Grundlage der Verteilung von wobei ist der Schätzer von Abtastung erhalten Beobachtungen aus . Insbesondere ist es das Intervall vom 2,5% -Quantil zum 97,5% -Quantil für die Verteilung von $\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)$ , was dazu führt, dass das Perzentilintervall völlig falsch ist.

Für das Basisintervall (und das normale Intervall) werden die Rollen der Quantile vertauscht. Dies impliziert, dass das Intervall angemessen erscheint (es umfasst 3,912), obwohl Intervalle, die über 3,912 hinausgehen, logisch nicht sinnvoll sind. Außerdem weiß ich nicht, ob das Basisintervall die richtige Abdeckung hat. Ihre Begründung basiert auf der folgenden ungefähren Verteilungsidentität:

θ (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}),$

n

$n$

n = 100

$n = 100$

Der letzte Vorschlag des OP für ein standardmäßiges fehlerbasiertes Intervall $\theta(\hat{\mathbf{p}}_n) \pm 1.96\hat{\mathrm{se}}_n$

$\theta(\mathbf{p})$ $\mathbf{p}$ $\theta(\mathbf{p})$

NRH
quelle

Das Bias-Problem bei der Verwendung des "Plug-in" -Schätzers für die Entropie ist seit Jahrzehnten bekannt.

1 / n

$1/n$

@EdM das sind sehr nützliche Informationen. Ich kannte die Literatur zu diesem speziellen Voreingenommenheitsproblem nicht. Es könnte sehr nützlich sein, wenn Sie den Kommentar in eine Antwort umwandeln, die die Bias-Korrektur erklärt und erklärt, wie sie beispielsweise beim Bootstrapping verwendet werden kann, um Konfidenzintervalle zu erhalten.

NRH

Ich kannte diese Literatur auch nicht, bis diese Frage und Ihre Antwort auftauchten. Das ist etwas peinlich, da die Shannon-Entropie in meinem Bereich der biomedizinischen Wissenschaft häufig als Maß verwendet wird. Ich werde sehen, was ich als zusätzliche Antwort zusammenstellen kann.

EdM

θ (p_{n}^{*})

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

n

$n$

n

$n$

n \to \infty

$n \to \infty$

Wie die Antwort von @NRH zeigt, besteht das Problem nicht darin, dass das Bootstrapping ein voreingenommenes Ergebnis lieferte. Es ist so, dass die einfache "Plug-in" -Schätzung der Shannon-Entropie, basierend auf Daten aus einer Stichprobe, vom wahren Populationswert nach unten verzerrt ist.

Dieses Problem wurde in den 1950er Jahren innerhalb weniger Jahre nach der Definition dieses Index erkannt. In diesem Artikel werden die zugrunde liegenden Probleme unter Bezugnahme auf die zugehörige Literatur erörtert.

$\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}.$

Die nichtlineare Beziehung bedeutet, dass der resultierende Wert eine voreingenommene Unterschätzung der wahren genetischen Vielfalt ist.

$M$ $N$ $(M -1)/2N$ . Korrekturen höherer Ordnung werden in dem oben verlinkten Artikel ausgewertet.

Es gibt Pakete in R, die sich mit diesem Problem befassen. Das simbootPaket hat insbesondere eine Funktion estShannonf, die diese Vorspannungskorrekturen vornimmt, und eine Funktion sbdivzum Berechnen von Konfidenzintervallen. Es ist besser, solche etablierten Open-Source-Tools für Ihre Analyse zu verwenden, als von vorne zu beginnen.

EdM
quelle

Der Schätzer an und für sich ist also aufgrund der Stichprobengröße fehlerhaft? Das simbootPaket sieht vielversprechend aus, scheint aber für meine Zwecke nicht geeignet zu sein, da es eine Kontrollprobe zur Schätzung der Konfidenzintervalle benötigt.

ZNK

"Fehlerhaft" ist nicht ganz richtig; Der Schätzer ist insofern "voreingenommen", als sein erwarteter Wert nicht mit dem tatsächlichen Bevölkerungswert übereinstimmt. Das heißt nicht, dass es "falsch" ist; Voreingenommene Schätzer können nützlich sein, wie der Kompromiss zwischen Bias und Varianz bei der Auswahl von Schätzern zeigt. Wenn simbootnicht erfüllt Ihre Bedürfnisse, Google „Shannon - Entropie Bias r“ für Verbindungen zu anderen R - Paketen wie entropy, entropart, und EntropyEstimation.

EdM

Es gibt zusätzliche Probleme, die sich aus der Tatsache ergeben, dass einige in der Population vorhandene Genotypen in einer bestimmten Stichprobe wahrscheinlich übersehen werden. Einige der bevölkerungs- und ökologiebasierten R-Pakete scheinen Möglichkeiten zu haben, mit diesem Problem umzugehen.

EdM