Como calculo um intervalo de confiança para a média de um conjunto de dados log-normal?

19

Ouvi / vi em vários lugares que você pode transformar o conjunto de dados em algo que é distribuído normalmente, pegando o logaritmo de cada amostra, calculando o intervalo de confiança para os dados transformados e transformando o intervalo de confiança novamente usando a operação inversa (por exemplo, eleve 10 à potência dos limites inferior e superior, respectivamente, para o log10 ).

No entanto, desconfio um pouco desse método, simplesmente porque ele não funciona para a média em si: 10mean(log10(X))mean(X)

Qual é a maneira correta de fazer isso? Se não funciona para a média em si, como pode funcionar para o intervalo de confiança da média?

Vegard
fonte
3
Você está certo. Essa abordagem geralmente não funciona e geralmente gera intervalos de confiança que não incluem a média da população ou mesmo a média da amostra. Aqui está uma discussão sobre isso: amstat.org/publications/jse/v13n1/olsson.html Esta não é uma resposta, já que não examinei o assunto o suficiente para comentar o link em detalhes.
Erik
3
Esse problema tem uma solução clássica: projecteuclid.org/… . Algumas outras soluções, incluindo código, são fornecidas em epa.gov/oswer/riskassessment/pdf/ucl.pdf-- mas leia isso com um grão pesado de sal, porque pelo menos um método descrito lá (o "Método de desigualdade de Chebyshev") é simplesmente errado.
whuber

Respostas:

11

Existem várias maneiras de calcular intervalos de confiança para a média de uma distribuição lognormal. Vou apresentar dois métodos: Bootstrap e Probabilidade de perfil. Também apresentarei uma discussão sobre os Jeffreys anteriores.

Bootstrap

Para o MLE

Neste caso, a MLE de (μ,σ) para uma amostra (x1,...,xn) são

μ^=1nj=1nlog(xj);σ^2=1nj=1n(log(xj)μ^)2.

Em seguida, a MLE da média é δ = exp ( μ + σ 2 / 2 ) . Redefinindo podemos obter uma amostra de bootstrap de δ e, usando este, podemos calcular várias inicialização intervalos de confiança. Os códigos a seguir mostram como obtê-los.δ^=exp(μ^+σ^2/2)δ^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")

Para a amostra média

Agora, considerando o estimador δ~=x¯ vez do MLE. Outro tipo de estimadores também pode ser considerado.

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

Probabilidade do perfil

Para a definição de funções de verossimilhança e verossimilhança de perfil, consulte . Usando a propriedade de invariância da probabilidade podemos reparameterise como se segue (μ,σ)(δ,σ) , onde δ=exp(μ+σ2/2) e, em seguida, calcular a probabilidade numericamente perfil de δ .

Rp(δ)=supσeu(δ,σ)supδ,σeu(δ,σ).

Esta função aceita valores em (0 0,1] ; um intervalo de nível 0,147 possui uma confiança aproximada de 95% . Vamos usar essa propriedade para construir um intervalo de confiança para δ . Os Rcódigos a seguir mostram como obter esse intervalo.

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)

Bayesiano

δ

(μ,σ)

π(μ,σ)σ-2,

n2R

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

Observe que eles são muito semelhantes.

kjetil b halvorsen
fonte
1
(+1) Eu acho que você também pode obter intervalos de confiança baseado na teoria de probabilidade máxima com o pacote distrMod R
Stéphane Laurent
@ StéphaneLaurent Obrigado pela informação. Gostaria de ver o resultado do seu código com o novo prior. Eu não estava ciente dos comandos e do pacote que você está usando.
4
n
Excelente resposta! As abordagens sugeridas aqui assumem erros no modelo homoscedástico - trabalhei em projetos nos quais essa suposição não era sustentável. Eu também sugeriria o uso da regressão gama como alternativa, o que ignoraria a necessidade de uma correção de viés.
Isabella Ghement
4

Você pode tentar a abordagem bayesiana com o de Jeffreys. Deveria gerar intervalos de credibilidade com uma propriedade de correspondência frequente correta: o nível de confiança do intervalo de credibilidade é próximo ao seu nível de credibilidade.

 # 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
fonte
Isso parece muito interessante e, como eu costumo gostar de métodos bayesianos, votei de maneira positiva. Ainda pode ser melhorado adicionando algumas referências ou, de preferência, até uma explicação compreensível sobre o porquê de funcionar.
Erik
μσ2μσ2μσ2f(μ,σ2)μσ2. Não sei se existem algumas referências, mas caso contrário, você pode verificar com simulações.
Stéphane Laurent
Muito obrigado pela discussão. Eu apaguei todos os meus comentários para maior clareza e para evitar qualquer confusão. (+1)
1
@ Procrastinator Obrigado também. Eu também apaguei meus comentários e adicionei o ponto sobre o Jeffreys antes no meu código.
Stéphane Laurent
Alguém poderia me explicar como o boot.out = boot (dados = dados0, estatística = função (d, ind) {mle (d [ind])}, R = 10000) funciona. Vejo que "ind" é um índice, mas não entendo como encontrar "ind". Onde está o segundo argumento referenciando? Eu tentei com funções alternativas e não funcionou. Olhando para a inicialização real da função, também não vejo uma referência a Ind.
precisa saber é o seguinte
0

No entanto, desconfio um pouco desse método, simplesmente porque ele não funciona para a média em si: 10mean (log10 (X)) ≠ mean (X)

Você está certo - essa é a fórmula para a média geométrica, não a média aritmética. A média aritmética é um parâmetro da distribuição normal e geralmente não é muito significativa para dados lognormal. A média geométrica é o parâmetro correspondente da distribuição lognormal, se você quiser falar mais significativamente sobre uma tendência central para seus dados.

E você realmente calcularia os ICs sobre a média geométrica, tomando os logaritmos dos dados, calculando a média e os ICs como de costume, e transformando-os novamente. Você está certo que realmente não deseja misturar suas distribuições colocando os ICs da média geométrica em torno da média aritmética ... uau!

dnidz
fonte