Jak obliczyć przedział ufności dla średniej logarytmicznego zestawu danych?

Słyszałem / widziałem w kilku miejscach, że możesz przekształcić zestaw danych w coś o rozkładzie normalnym, biorąc logarytm każdej próbki, oblicz przedział ufności dla przekształconych danych i przekształcaj przedział ufności z powrotem za pomocą operacji odwrotnej (np. podnieś 10 do potęgi dolnej i górnej granicy odpowiednio dla ). $\log_{10}$

Jednak jestem nieco podejrzliwy wobec tej metody, po prostu dlatego, że nie działa ona dla samego środka: $10^{\operatorname{mean}(\log_{10}(X))} \ne \operatorname{mean}(X)$

Jaki jest właściwy sposób to zrobić? Jeśli to nie działa dla samego środka, to jak może działać dla przedziału ufności dla środka?

confidence-interval mean lognormal Vegard
źródło

Masz rację. Takie podejście na ogół nie działa i często daje przedziały ufności, które nie obejmują średniej populacji, a nawet średniej próby. Oto kilka dyskusji na ten temat: amstat.org/publications/jse/v13n1/olsson.html To nie jest odpowiedź, ponieważ nie przyjrzałem się tej sprawie na tyle, aby szczegółowo skomentować link.

Erik

Ten problem ma klasyczne rozwiązanie: projecteuclid.org/… . Niektóre inne rozwiązania, w tym kod, są dostępne na stronie epa.gov/oswer/riskassessment/pdf/ucl.pdf - ale przeczytaj to z ciężkim ziarnem soli, ponieważ przynajmniej jedna z opisanych tam metod („Metoda nierówności Czebyszewa”) jest po prostu źle.

whuber

Odpowiedzi:

Istnieje kilka sposobów obliczania przedziałów ufności dla średniej rozkładu logarytmicznego. Przedstawię dwie metody: Bootstrap i prawdopodobieństwo profilu. Przedstawię również dyskusję na temat przeora Jeffreysa.

Bootstrap

Dla MLE

W tym przypadku MLE dla próbki wynosi $(\mu,\sigma)$ $(x_1,...,x_n)$

\hat{μ} = \frac{1}{n} \sum_{jot = 1}^{n} \log (x_{jot}); {\hat{σ}}^{2)} = \frac{1}{n} \sum_{jot = 1}^{n} (\log (x_{jot}) - \hat{μ})^{2)} .

$\hat\mu= \dfrac{1}{n}\sum_{j=1}^n\log(x_j);\,\,\,\hat\sigma^2=\dfrac{1}{n}\sum_{j=1}^n(\log(x_j)-\hat\mu)^2.$

Zatem MLE średniej to . Przez resampling możemy otrzymać próbkę bootstrap z i za pomocą tego, możemy obliczyć kilka bootstrap przedziały ufności. Poniższe kody pokazują, jak je uzyskać. $\hat\delta=\exp(\hat\mu+\hat\sigma^2/2)$ $\hat\delta$ R

rm(list=ls())
library(boot)

set.seed(1)

# Simulated data
data0 = exp(rnorm(100))

# Statistic (MLE)

mle = function(dat){
m = mean(log(dat))
s = mean((log(dat)-m)^2)
return(exp(m+s/2))
}

# Bootstrap
boots.out = boot(data=data0, statistic=function(d, ind){mle(d[ind])}, R = 10000)
plot(density(boots.out$t))

# 4 types of Bootstrap confidence intervals
boot.ci(boots.out, conf = 0.95, type = "all")

Dla próbki średniej

Teraz, biorąc pod uwagę estymator zamiast MLE. Można również rozważyć inny rodzaj estymatorów. $\tilde{\delta}=\bar{x}$

rm(list=ls())
library(boot)

set.seed(1)

# Simulated data
data0 = exp(rnorm(100))

# Statistic (MLE)

samp.mean = function(dat) return(mean(dat))

# Bootstrap
boots.out = boot(data=data0, statistic=function(d, ind){samp.mean(d[ind])}, R = 10000)
plot(density(boots.out$t))

# 4 types of Bootstrap confidence intervals
boot.ci(boots.out, conf = 0.95, type = "all")

Prawdopodobieństwo profilu

Aby zapoznać się z definicją funkcji wiarygodności i funkcji wiarygodności profilu, patrz . Korzystając z właściwości niezmienniczości prawdopodobieństwa, możemy ponownie sparametryzować w następujący sposób , gdzie a następnie obliczyć numerycznie prawdopodobieństwo profilu . $(\mu,\sigma)\rightarrow(\delta,\sigma)$ $\delta=\exp(\mu+\sigma^2/2)$ $\delta$

R_{p} (δ) = \frac{\underset{σ}{łyk} L. (δ, σ)}{\underset{δ, σ}{łyk} L. (δ, σ)} .

$R_p(\delta)=\dfrac{\sup_{\sigma}{\mathcal L}(\delta,\sigma)}{\sup_{\delta,\sigma}{\mathcal L}(\delta,\sigma)}.$

Ta funkcja przyjmuje wartości w ; przedział poziomu ma przybliżoną pewność . Będziemy używać tej właściwości do konstruowania przedziału ufności dla . Poniższe kody pokazują, jak uzyskać ten przedział . $(0,1]$ $0.147$ $95\%$ $\delta$ R

set.seed(1)

# Simulated data
data0 = exp(rnorm(100))

# Log likelihood
ll = function(mu,sigma) return( sum(log(dlnorm(data0,mu,sigma))))

# Profile likelihood
Rp = function(delta){
temp = function(sigma) return( sum(log(dlnorm(data0,log(delta)-0.5*sigma^2,sigma)) ))
max=exp(optimize(temp,c(0.25,1.5),maximum=TRUE)$objective     -ll(mean(log(data0)),sqrt(mean((log(data0)-mean(log(data0)))^2))))
return(max)
}

vec = seq(1.2,2.5,0.001)
rvec = lapply(vec,Rp)
plot(vec,rvec,type="l")

# Profile confidence intervals
tr = function(delta) return(Rp(delta)-0.147)
c(uniroot(tr,c(1.2,1.6))$root,uniroot(tr,c(2,2.3))$root)

$\star$ Bayesian

W tej sekcji przedstawiono alternatywny algorytm oparty na próbkowaniu Metropolis-Hastings i wykorzystaniu wcześniejszych Jeffreysów do obliczania przedziału wiarygodności dla . $\delta$

Przypomnijmy, że w Jeffreys przed o w modelu logarytmiczno-normalnego jest $(\mu,\sigma)$

π (μ, σ) \propto σ^{- 2)},

$\pi(\mu,\sigma)\propto \sigma^{-2},$

i że ten przełożony jest niezmienny podczas reparametryzacji. To wcześniejsze jest niewłaściwe, ale tylna część parametrów jest właściwa, jeśli wielkość próbki . Poniższy kod pokazuje, jak uzyskać 95% przedział wiarygodności za pomocą tego modelu Bayesa. $n\geq 2$ R

library(mcmc)

set.seed(1)

# Simulated data
data0 = exp(rnorm(100))

# Log posterior
lp = function(par){
if(par[2]>0) return( sum(log(dlnorm(data0,par[1],par[2]))) - 2*log(par[2]))
else return(-Inf)
}

# Metropolis-Hastings
NMH = 260000
out = metrop(lp, scale = 0.175, initial = c(0.1,0.8), nbatch = NMH)

#Acceptance rate
out$acc

deltap = exp(  out$batch[,1][seq(10000,NMH,25)] + 0.5*(out$batch[,2][seq(10000,NMH,25)])^2  )

plot(density(deltap))

# 95% credibility interval
c(quantile(deltap,0.025),quantile(deltap,0.975))

Zauważ, że są bardzo podobne.

kjetil b halvorsen
źródło

(+1) Myślę, że można również uzyskać przedziały ufności w oparciu o teorię największego prawdopodobieństwa z pakietem distrMod R

Stéphane Laurent

@ StéphaneLaurent Dzięki za informację. Chciałbym zobaczyć wynik twojego kodu z nowym przeorem. Nie wiedziałem o komendach i pakiecie, którego używasz.

Piękna odpowiedź @ Procrastinator. Innym podejściem jest estymator rozmazywania, który wykorzystuje wszystkie pozostałości poza średnią na skali logarytmicznej, aby uzyskać przewidywanych wartości na oryginalnej skali i po prostu je uśrednić. Nie jestem jednak na bieżąco z przedziałami ufności stosującymi to podejście, z wyjątkiem standardowej metody percentyla bootstrap.

n

$n$

Frank Harrell,

Znakomita odpowiedź! Proponowane tutaj podejścia zakładają błędy modelu homoscedastycznego - pracowałem nad projektami, w których to założenie nie było możliwe do utrzymania. Sugerowałbym również zastosowanie regresji gamma jako alternatywy, która pomijałaby potrzebę korekcji odchylenia.

Isabella Ghement

Możesz spróbować podejścia bayesowskiego z przeorem Jeffreysa. Powinien dawać przedziały wiarygodności z poprawną właściwością dopasowania częstości: poziom ufności przedziału wiarygodności jest zbliżony do jego poziomu wiarygodności.

 # required package
 library(bayesm)

 # simulated data
 mu <- 0
 sdv <- 1
 y <- exp(rnorm(1000, mean=mu, sd=sdv))

 # model matrix
 X <- model.matrix(log(y)~1)
 # prior parameters
 Theta0 <- c(0)
 A0 <- 0.0001*diag(1)
 nu0 <- 0 # Jeffreys prior for the normal model; set nu0 to 1 for the lognormal model
 sigam0sq <- 0
 # number of simulations
 n.sims <- 5000

 # run posterior simulations
 Data <- list(y=log(y),X=X)
 Prior <- list(betabar=Theta0, A=A0, nu=nu0, ssq=sigam0sq)
 Mcmc <- list(R=n.sims)
 bayesian.reg <- runireg(Data, Prior, Mcmc)
 mu.sims <- t(bayesian.reg$betadraw) # transpose of bayesian.reg$betadraw
 sigmasq.sims <- bayesian.reg$sigmasqdraw

 # posterior simulations of the mean of y: exp(mu+sigma²/2)
 lmean.sims <- exp(mu.sims+sigmasq.sims/2)

 # credibility interval about lmean:
 quantile(lmean.sims, probs = c(0.025, 0.975))

Stéphane Laurent
źródło

Brzmi bardzo interesująco, a ponieważ lubię metody bayesowskie, podniosłem głos. Można go jeszcze ulepszyć, dodając kilka referencji lub najlepiej zrozumiałe wyjaśnienie, dlaczego to działa.

Erik

Wiadomo, że „it” (właściwość dopasowania częstych) działa dla i . W przypadku właściwość dopasowania częstości jest idealna: przedział wiarygodności jest dokładnie taki sam jak zwykły przedział ufności. Dla Nie wiem, czy jest dokładny, ale łatwo to sprawdzić, ponieważ rozkład tylny jest odwrotną gamma. Fakt, że działa dla i , niekoniecznie oznacza, że działa dla funkcji z i . Nie wiem, czy istnieją jakieś referencje, ale w przeciwnym razie możesz sprawdzić za pomocą symulacji.

μ

$\mu$

σ^{2}

$\sigma^2$

μ

$\mu$

σ^{2}

$\sigma^2$

μ

$\mu$

σ^{2}

$\sigma^2$

f (μ, σ^{2})

$f(\mu, \sigma^2)$

μ

$\mu$

σ^{2}

$\sigma^2$

Stéphane Laurent,

Wielkie dzięki za dyskusję. Usunąłem wszystkie moje komentarze w celu zachowania przejrzystości i uniknięcia nieporozumień. (+1)

@Procrastinator Też dziękuję. Usunąłem również moje komentarze i dodałem punkt o Jeffreys wcześniej w moim kodzie.

Stéphane Laurent,

Czy ktoś mógłby mi wyjaśnić, jak działa boots.out = boot (data = data0, statystyka = funkcja (d, ind) {mle (d [ind])}, R = 10000). Widzę, że „ind” jest indeksem, ale nie rozumiem, jak znaleźć „ind”. Gdzie jest ten drugi argument? Próbowałem z alternatywnymi funkcjami i nie działało. Patrząc na rzeczywisty rozruch funkcji, nie widzę też odniesienia do Ind.

andor kesselman

Jednak jestem nieco podejrzliwy wobec tej metody, po prostu dlatego, że nie działa ona dla samego środka: 10mean (log10 (X)) ≠ mean (X)

Masz rację - to wzór na średnią geometryczną, a nie średnią arytmetyczną. Średnia arytmetyczna jest parametrem z rozkładu normalnego i często nie ma większego znaczenia dla danych logarytmicznych. Średnia geometryczna jest odpowiednim parametrem z rozkładu logarytmicznego, jeśli chcesz bardziej sensownie mówić o centralnej tendencji dla twoich danych.

I rzeczywiście obliczalibyśmy CI względem średniej geometrycznej, biorąc logarytmy danych, obliczając średnią i CI jak zwykle i przekształcając wstecz. Masz rację, że tak naprawdę nie chcesz mieszać swoich rozkładów, umieszczając CI dla średniej geometrycznej wokół średniej arytmetycznej ... tak!

dnidz
źródło