19
$\begingroup$

I've heard/seen in several places that you can transform the data set into something that is normal-distributed by taking the logarithm of each sample, calculate the confidence interval for the transformed data, and transform the confidence interval back using the inverse operation (e.g. raise 10 to the power of the lower and upper bounds, respectively, for $\log_{10}$).

However, I'm a bit suspicious of this method, simply because it doesn't work for the mean itself: $10^{\operatorname{mean}(\log_{10}(X))} \ne \operatorname{mean}(X)$

What is the correct way to do this? If it doesn't work for the mean itself, how can it possibly work for the confidence interval for the mean?

$\endgroup$
  • 3
    $\begingroup$ You are quite right. This approach does not generally work and often yields confidence intervals that do not include the population mean or even the sample mean. Here is some discussion on it: amstat.org/publications/jse/v13n1/olsson.html This is not an answer, since I did not look into the matter enough to actually comment on the link in detail. $\endgroup$ – Erik Jul 31 '12 at 8:53
  • 3
    $\begingroup$ This problem has a classic solution: projecteuclid.org/…. Some other solutions, including code, are provided at epa.gov/oswer/riskassessment/pdf/ucl.pdf--but read this with a heavy grain of salt, because at least one method described there (the "Chebyshev Inequality Method") is just plain wrong. $\endgroup$ – whuber Jul 31 '12 at 15:21
11
$\begingroup$

There are several ways for calculating confidence intervals for the mean of a lognormal distribution. I am going to present two methods: Bootstrap and Profile likelihood. I will also present a discussion on the Jeffreys prior.

Bootstrap

For the MLE

In this case, the MLE of $(\mu,\sigma)$ for a sample $(x_1,...,x_n)$ are

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

Then, the MLE of the mean is $\hat\delta=\exp(\hat\mu+\hat\sigma^2/2)$. By resampling we can obtain a bootstrap sample of $\hat\delta$ and, using this, we can calculate several bootstrap confidence intervals. The following R codes shows how to obtain these.

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

For the sample mean

Now, considering the estimator $\tilde{\delta}=\bar{x}$ instead of the MLE. Other type of estimators might be considered as well.

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

Profile likelihood

For the definition of likelihood and profile likelihood functions, see. Using the invariance property of the likelihood we can reparameterise as follows $(\mu,\sigma)\rightarrow(\delta,\sigma)$, where $\delta=\exp(\mu+\sigma^2/2)$ and then calculate numerically the profile likelihood of $\delta$.

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

This function takes values in $(0,1]$; an interval of level $0.147$ has an approximate confidence of $95\%$. We are going to use this property for constructing a confidence interval for $\delta$. The following R codes shows how to obtain this interval.

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

In this section, an alternative algorithm, based on Metropolis-Hastings sampling and the use of the Jeffreys prior, for calculating a credibility interval for $\delta$ is presented.

Recall that the Jeffreys prior for $(\mu,\sigma)$ in a lognormal model is

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

and that this prior is invariant under reparameterisations. This prior is improper, but the posterior of the parameters is proper if the sample size $n\geq 2$. The following R code shows how to obtain a 95% credibility interval using this Bayesian model.

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

Note that they are very similar.

$\endgroup$
  • 1
    $\begingroup$ (+1) I think you can also get confidence intervals based on maximum-likelihood theory with the distrMod R package $\endgroup$ – Stéphane Laurent Jul 31 '12 at 11:51
  • $\begingroup$ @StéphaneLaurent Thanks for the info. I would like to see the outcome of your code with the new prior. I was not aware of the commands and the package that you are using. $\endgroup$ – user10525 Jul 31 '12 at 11:53
  • 4
    $\begingroup$ Beautiful response @Procrastinator. One other approach is the smearing estimator, which uses all the residuals off of the mean on the log scale to get $n$ predicted values on the original scale and simply averages them. I'm less up to date on confidence intervals using this approach though, except for using the standard bootstrap percentile method. $\endgroup$ – Frank Harrell Jul 31 '12 at 12:15
  • $\begingroup$ Superb response! The approaches suggested here assume homoscedastic model errors - I have worked on projects where this assumption was not tenable. I would also suggest the use of gamma regression as an alternative, which would bypass the need for a bias correction. $\endgroup$ – Isabella Ghement May 18 at 16:53
4
$\begingroup$

You might try the Bayesian approach with Jeffreys' prior. It should yield credibility intervals with a correct frequentist-matching property: the confidence level of the credibility interval is close to its credibility level.

 # 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))
$\endgroup$
  • $\begingroup$ This sounds very interesting and since I tend to like Bayesian methods I upvoted it. It could still be improved by adding an some references or preferably even an understandable explanation on why it works. $\endgroup$ – Erik Jul 31 '12 at 8:58
  • $\begingroup$ It is known that "it" (the frequentist-matching property) works for $\mu$ and $\sigma^2$. For $\mu$ the frequentist-matching property is perfect: the credibility interval is exactly the same as the usual confidence interval. For $\sigma^2$ I don't know whether it is exact but it is easy to check because the posterior distribution is an inverse-Gamma. The fact that it works for $\mu$ and $\sigma^2$ does not necessarily implies that it works for a function $f(\mu, \sigma^2)$ of $\mu$ and $\sigma^2$. I don't know whether there are some references but otherwise you can check with simulations. $\endgroup$ – Stéphane Laurent Jul 31 '12 at 9:17
  • $\begingroup$ Many thanks for the discussion. I deleted all my comments for clearness and to avoid any confusion. (+1) $\endgroup$ – user10525 Jul 31 '12 at 15:08
  • 1
    $\begingroup$ @Procrastinator Thanks too. I have also deleted my comments and added the point about the Jeffreys prior in my code. $\endgroup$ – Stéphane Laurent Jul 31 '12 at 15:15
  • $\begingroup$ Could someone please explain to me how boots.out = boot(data=data0, statistic=function(d, ind){mle(d[ind])}, R = 10000) works. I see that "ind" is an index, but I don't understand how to find "ind". Where is this second argument referencing? I've tried it with alternative functions and it did not work. Looking at the actual function boot, I don't see a reference to Ind either. $\endgroup$ – andor kesselman Sep 30 '15 at 2:32
0
$\begingroup$

However, I'm a bit suspicious of this method, simply because it doesn't work for the mean itself: 10mean(log10(X))≠mean(X)

You're right -- that's the formula for the geometric mean, not the arithmetic mean. The arithmetic mean is a parameter from the normal distribution, and is often not very meaningful for lognormal data. The geometric mean is the corresponding parameter from the lognormal distribution if you want to talk more meaningfully about a central tendency for your data.

And you would indeed calculate the CIs about the geometric mean by taking the logarithms of the data, calculating the mean and CIs as usual, and back-transforming. You're right that you really don't want to mix your distributions by putting the CIs for the geometric mean around the arithmetic mean....yeowch!

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.