Biased bootstrap: czy można wyśrodkować CI wokół obserwowanej statystyki?

Jest to podobne do Bootstrap: oszacowanie jest poza przedziałem ufności

Mam pewne dane, które reprezentują liczbę genotypów w populacji. Chcę oszacować różnorodność genetyczną za pomocą indeksu Shannona, a także wygenerować przedział ufności za pomocą ładowania początkowego. Zauważyłem jednak, że oszacowanie za pomocą ładowania początkowego jest zwykle bardzo stronnicze i skutkuje przedziałem ufności, który leży poza moją obserwowaną statystyką.

Poniżej znajduje się przykład.

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

Generowanie danych

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

Obliczenie

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

Generowanie elementów CI z korekcją błędów

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

Zakładając, że wariancja t może być wykorzystana do wariancji t0 .

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

## [1] 3.55475 3.80421

Czy poprawne byłoby zgłoszenie CI wyśrodkowanego wokół t0 ? Czy jest lepszy sposób na wygenerowanie bootstrapu?

r confidence-interval bootstrap bias diversity ZNK
źródło

Odpowiedzi:

W konfiguracji podanej przez OP parametrem będącym przedmiotem zainteresowania jest entropia Shannona która jest funkcją wektora prawdopodobieństwa . Estymator oparty na próbkach ( w symulacji) to estymator wtyczki Próbki zostały wygenerowane przy użyciu rozkładu jednolitego, dla którego entropia Shannona wynosiPonieważ entropia Shannona jest zmaksymalizowana w rozkładzie równomiernym, estymator wtyczki musi być tendencyjny w dół . Symulacja pokazuje to

θ (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$ podczas gdy . Estymator wtyczek jest spójny, ale metoda nie ma zastosowania dla będącego rozkładem równomiernym, ponieważ pochodną entropii Shannona jest 0. Zatem dla tego konkretnego wyboru , przedziały ufności oparte na asymptotycznych argumentach nie są oczywiste.

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

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

Δ

$\Delta$

p

$\mathbf{p}$

p

$\mathbf{p}$

Interwał percentyla jest oparty na rozkładzie gdzie jest estymatorem uzyskanym z próbkowania obserwacji z . W szczególności jest to przedział od kwantyla 2,5% do kwantyla 97,5% dla rozkładu . Jak pokazuje symulacja bootstrapu PO, jest również wyraźnie tendencyjna w dół jako estymator , co powoduje, że interwał percentyla jest całkowicie źle. $\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)$

W podstawowym (i normalnym) przedziale role kwantyli są zamienione. Oznacza to, że przedział wydaje się rozsądny (obejmuje 3,912), chociaż przedziały przekraczające 3,912 nie mają logicznego znaczenia. Co więcej, nie wiem, czy podstawowy interwał będzie miał odpowiedni zasięg. Jego uzasadnienie opiera się na następującej przybliżonej tożsamości dystrybucyjnej:

θ (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}),$ co może budzić wątpliwości (względnie) małe np. .

n

$n$

n = 100

$n = 100$

Ostatnia sugestia OP dotycząca standardowego interwału opartego na błędach również nie zadziała z powodu dużego błędu. Może to działać w przypadku estymatora z korekcją błędów, ale najpierw potrzebujesz poprawnych standardowych błędów dla estymatora z korekcją błędów. $\theta(\hat{\mathbf{p}}_n) \pm 1.96\hat{\mathrm{se}}_n$

Rozważę przedział prawdopodobieństwa oparty na prawdopodobieństwie dziennika profilu dla . Obawiam się, że nie znam żadnego prostego sposobu obliczenia prawdopodobieństwa logarytmu profilu dla tego przykładu, z wyjątkiem tego, że musisz zmaksymalizować prawdopodobieństwo logarytmu nad dla różnych stałych wartości . $\theta(\mathbf{p})$ $\mathbf{p}$ $\theta(\mathbf{p})$

NRH
źródło

Problem odchylenia przy użyciu estymatora „plug-in” do entropii jest doceniany od dziesięcioleci. Ten artykuł analizuje mniej tendencyjne szacunki. W przypadku przedstawionym przez PO można zastosować korektę błędu systematycznego do rzędu , która pochodzi z 1955 r. (Patrz równanie 4 w powiązanym dokumencie). Korekta wynosi 0,245, prawie identycznie z odchyleniem zidentyfikowanym przez pasek startowy. Być może należy zastosować tutaj bootstrap do oszacowania samej entropii, a nie tylko jej granic ufności.

1 / n

$1/n$

EdM

@EdM to bardzo przydatna informacja. Nie znałem literatury na temat tego konkretnego problemu uprzedzeń. Może to być bardzo przydatne, jeśli możesz zamienić komentarz w odpowiedź wyjaśniającą korektę błędu i jak można go użyć, na przykład, w celu uzyskania przedziałów ufności.

NRH,

Nie znałem też tej literatury, dopóki nie pojawiło się to pytanie i odpowiedź. Co jest nieco zawstydzające, ponieważ entropię Shannona często stosuje się jako miarę w mojej dziedzinie nauk biomedycznych. Zobaczę, co mogę złożyć jako dodatkową odpowiedź.

EdM

Zwiększenie liczby próbek bootstrap tak naprawdę nie pomoże. Musi być wystarczająco duży, aby można było w wiarygodny sposób oszacować interesujące ilości dla dystrybucji , powiedzmy, ale w przeciwnym razie zwiększenie liczby próbek bootstrap nie usunie błędu lub nie spowoduje zaufanie bardziej odpowiednie.

θ (p_{n}^{*})

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

NRH,

Przepraszam ZNK, źle zrozumiałem twoje pytanie. Jeśli zwiększysz wielkość próby , obciążenie będzie mniejsze, tak! Estymator jest spójny. Właśnie w odniesieniu do rozkładu równomiernego byłbym nieco sceptycznie nastawiony do faktycznego pokrycia przedziałów ufności nawet dla dużej liczby z powodów opisanych w odpowiedzi. W przypadku wszystkich innych dystrybucji obowiązuje CLT, a różne metody zapewnią asymptotycznie poprawne pokrycie dla .

n

$n$

n

$n$

n \to \infty

$n \to \infty$

NRH

Jak wskazuje odpowiedź @NRH, problemem nie jest to, że ładowanie początkowe dało stronniczy wynik. Chodzi o to, że proste oszacowanie entropii Shannona metodą „plug-in”, oparte na danych z próbki, jest tendencyjne w dół w stosunku do prawdziwej wartości populacji.

Problem ten został rozpoznany w latach 50. XX wieku, w ciągu kilku lat od zdefiniowania tego wskaźnika. W tym artykule omówiono podstawowe problemy z odniesieniami do powiązanej literatury.

Problem wynika z nieliniowej relacji indywidualnych prawdopodobieństw do tej miary entropii. W tym przypadku zaobserwowana frakcja genotypowa dla genu i w próbce n , , jest bezstronnym estymatorem prawdziwego prawdopodobieństwa . Ale kiedy ta zaobserwowana wartość zostanie zastosowana do formuły „plug in” dla entropii ponad genami 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}.$

relacja nieliniowa oznacza, że wynikowa wartość jest tendencyjnym niedoszacowaniem prawdziwej różnorodności genetycznej.

Odchylenie zależy od liczby genów i liczbę obserwacji, . Przy pierwszym zamówieniu oszacowanie wtyczki będzie niższe niż prawdziwa entropia o kwotę . Korekty wyższych rzędów są oceniane w dokumencie powyżej. $M$ $N$ $(M -1)/2N$

Istnieją pakiety w R, które rozwiązują ten problem. W simbootszczególności pakiet ma funkcję, estShannonfktóra dokonuje tych poprawek odchylenia, oraz funkcję sbdivobliczania przedziałów ufności. Lepiej będzie używać do analizy takich uznanych narzędzi typu open source, niż próbować zaczynać od zera.

EdM
źródło

Czyli estymator sam w sobie jest błędny ze względu na wielkość próby? Na simbootwygląd opakowania obiecujące, ale nie wydaje się odpowiednie dla moich celów, gdyż potrzebował próbkę kontrolną do oceny przedziałów ufności.

ZNK,

„Błędny” nie jest w porządku; estymator jest „stronniczy”, ponieważ jego oczekiwana wartość nie jest taka sama jak rzeczywista wartość populacji. To nie znaczy, że jest „błędne”; tendencyjne estymatory mogą być przydatne, co ilustruje kompromis wariancji odchylenia przy wyborze estymatorów. Jeśli simbootnie spełnia Twoich potrzeb, Google „entropii Shannona bias r” dla linków do innych pakietów, takich jak R entropy, entroparti EntropyEstimation.

EdM

Istnieją dodatkowe problemy wynikające z faktu, że niektóre genotypy obecne w populacji prawdopodobnie zostaną pominięte w jakiejkolwiek konkretnej próbce. Wydaje się, że niektóre pakiety R oparte na populacji i ekologii mają sposoby poradzenia sobie z tym problemem.

EdM