Bagaimana cara menghitung interval kepercayaan untuk rata-rata kumpulan data log-normal?

Saya pernah mendengar / melihat di beberapa tempat bahwa Anda dapat mengubah set data menjadi sesuatu yang terdistribusi normal dengan mengambil logaritma masing-masing sampel, menghitung interval kepercayaan untuk data yang diubah, dan mengubah interval kepercayaan kembali menggunakan operasi terbalik (mis. naikkan 10 pangkat masing-masing untuk batas bawah dan atas, untuk ).$\log_{10}$

Namun, saya agak curiga terhadap metode ini, hanya karena metode ini tidak berfungsi untuk rata-rata itu sendiri:$10^{\operatorname{mean}(\log_{10}(X))} \ne \operatorname{mean}(X)$

Apa cara yang benar untuk melakukan ini? Jika itu tidak bekerja untuk mean itu sendiri, bagaimana mungkin itu bekerja untuk interval kepercayaan untuk mean?

Ada beberapa cara untuk menghitung interval kepercayaan untuk rata-rata distribusi lognormal. Saya akan menyajikan dua metode: Bootstrap dan kemungkinan profil. Saya juga akan mempresentasikan diskusi tentang Jeffrey sebelumnya.

Bootstrap

Untuk MLE

Dalam hal ini, MLE dari $(\mu,\sigma)$ untuk sampel $(x_1,...,x_n)$ adalah

Kemudian, MLE dari mean adalah δ = exp ( μ + σ 2 / 2 ) . Dengan resampling kita dapat memperoleh sampel bootstrap dari δ dan, dengan menggunakan ini, kita dapat menghitung beberapa bootstrap interval keyakinan. Kode berikut menunjukkan cara mendapatkan ini.$\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")

Untuk mean sampel

Sekarang, pertimbangkan estimator $\tilde{\delta}=\bar{x}$ daripada MLE. Jenis penduga lain mungkin dipertimbangkan juga.

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")

Kemungkinan profil

Untuk definisi fungsi kemungkinan dan kemungkinan profil, lihat . Menggunakan properti invarian dari kemungkinan kita bisa reparameterise sebagai berikut $(\mu,\sigma)\rightarrow(\delta,\sigma)$ , di mana $\delta=\exp(\mu+\sigma^2/2)$ dan kemudian menghitung numerik kemungkinan profil dari $\delta$ .

Fungsi ini mengambil nilai dalam $(0,1]$ ; interval level $0.147$ memiliki perkiraan kepercayaan $95\%$ . Kita akan menggunakan properti ini untuk membangun interval kepercayaan untuk $\delta$ . RKode berikut menunjukkan cara mendapatkan interval ini.

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

Pada bagian ini, disajikan algoritma alternatif, berdasarkan pada pengambilan sampel Metropolis-Hastings dan penggunaan Jeffrey, untuk menghitung interval kredibilitas untuk $\delta$ .

Ingat bahwa Jeffreys sebelum untuk $(\mu,\sigma)$ dalam model lognormal adalah

$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))

Perhatikan bahwa mereka sangat mirip.

Anda dapat mencoba pendekatan Bayesian dengan Jeffreys sebelumnya. Ini harus menghasilkan interval kredibilitas dengan properti pencocokan frequentist yang benar: tingkat kepercayaan interval kredibilitas dekat dengan tingkat kredibilitasnya.

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

Namun, saya agak curiga dengan metode ini, hanya karena tidak bekerja untuk rata-rata itu sendiri: 10 berarti (log10 (X)) ≠ berarti (X)

Anda benar - itu rumus untuk mean geometrik, bukan mean aritmatika. Mean aritmatika adalah parameter dari distribusi normal, dan seringkali tidak terlalu berarti untuk data lognormal. Mean geometrik adalah parameter yang sesuai dari distribusi lognormal jika Anda ingin berbicara lebih bermakna tentang kecenderungan sentral untuk data Anda.

Dan Anda memang akan menghitung CI tentang mean geometrik dengan mengambil logaritma data, menghitung mean dan CI seperti biasa, dan mengubah-balik. Anda benar bahwa Anda benar-benar tidak ingin mencampur distribusi Anda dengan meletakkan CI untuk rata-rata geometris di sekitar rata-rata aritmatika .... yeowch!