Short history:
As I understood can the normalized (profile) likelihood be used to get the CI of an estimate. "Normalized" means that the function will have a unit integral.
This works beautifully for the binomial distribution (the normalized likelihood is a beta-density). For the expectation of the normal distribution the dispersion parameter needs to be considered -> profile likelihood, and the normalized profile likelihood is the t-density.
My problem:
Now I thought I should be able go a similar way with gamma-distributed data. The ready solution in R seems to be fitting a gamma-glm and use confint(). Here is a working example:
y = c(269, 299, 437, 658, 1326)
mod = glm(y~1, family=Gamma(link=identity))
(coef(mod)-mean(y) < 1E-4) # TRUE
confint(mod, level=0.9)
Ok, now my attempt to solve this problem step-by-step by hand:
The likelihood function, reparametrized to mu = scale*shape:
L = Vectorize(function(mu,s) prod(dgamma(y,shape=mu/s,scale=s)))
L.opt = function(p) 1/L(p[1],p[2])
optim(par=c(100,100), fn=L.opt, control=c(reltol=1E-20))
The MLE for mu is equal to the coefficient of the model. The next lines produce a contour plot of the likelihood function with a red line indicating the maximum depending on mu:
m = seq(200,1500,len=100)
s = seq( 50,1500,len=100)
Lmat = outer(m,s,L)
Rmat = Lmat/max(Lmat)
contour(m,s,Rmat,levels=seq(0.1,1,len=10))
s.max = function(mu) optimize(function(s) L(mu,s), lower=1, upper=2000, maximum=TRUE)$maximum
lines(m,Vectorize(s.max)(m),col=2,lty=3)
The next step is to get the profile likelihood (and look at a plot):
profileL = Vectorize(function(mu) L(mu,s.max(mu)))
plot(profileL, xlim=c(200,2000))
Then this profile likelihood is normalized to get a unit integral:
AUC = integrate(profileL,0,Inf,rel.tol=1E-20)$value
normProfileL = function(mu) profileL(mu)/AUC
(integrate(normProfileL,0,Inf)$value-1 < 1E-4) # TRUE
plot(normProfileL, xlim=c(200,2000))
Now obtain the 90% CI as the quantiles for which the integral of this normalized profile likelihood is 0.05 and 0.95:
uniroot(function(mu) integrate(normProfileL,0,mu)$value-0.05, lower=0, upper=5000)$root
uniroot(function(mu) integrate(normProfileL,0,mu)$value-0.95, lower=0, upper=5000)$root
This is 413.9 ... 1557.
confint(mod,level=0.9) # gives 364.8 ... 1076
This seems to be not just a rounding error.
Where am I wrong?
Hi,
Your question: How to get a confidence interval (CI) for gamma-distributed data?
http://www.ncbi.nlm.nih.gov/pubmed/9131766
http://vmtz-won-app1.cdc.gov/wonder/help/cancer/FayFeuerConfidenceIntervals.pdf
http://www.stat.umn.edu/geyer/old03/5102/notes/ci.pdf
http://stats.stackexchange.com/questions/89230/confidence-interval-for-a-random-sample-selected-from-gamma-distribution
http://www4.stat.ncsu.edu/~boos/slab/slabe/lab7.html
http://dl.sjp.ac.lk/dspace/bitstream/123456789/1007/1/Confidence%20intervals%20for%20the%20median.pdf
Thank you Senthilvel.
The links were not very helpful, unfortunately.
But I think I solved the problem with help of the book "In all Likelihood" (Yudi Pawitan, Oxfort University Press, 2013). On page 94 there is in fact an example of getting the profile likelihood of a gamma sample (example 4.9). I tried this in R and get the same results.
I think one of my problems was a wrong parametrization of the gamma density.
Then I noticed that the confidence interval of the glm in R is based on the inversion of the likelihood-ratio test and, for mu, simply finds the values of the likelihood (at phi = phi.hat [= the MLE for phi]), not of the profile likelihood. This gives areas under the normalized (profile)likelihood that are lower than the given confidence level (i tried levels of for instance 0.9 and got areas of 0.83). This (the difference) seems to be the result of the chi-quare approximation (the larger the sample size the better is the match between desired conf.level and area under the norm. likelihood). I would be greatful if an expert could comment on this (am I right or wrong? how wrong? wrong in what way?...)
Here is the R code reproducing the results from example 4.9 in the book, showing the likelihood contour, the likelihood function for mu at phi.hat and the profile likelihood, plus the confidence limits obtained by glm:
y = c(23.91, 27.33, 0.15, 3.65, 5.99, 0.88, 0.93,0.53,0.17,14.17,6.18,0.05,3.89,0.24,0.08)
n = length(y)
# re-parametrized gamma desnity with mu and phi
# scale (s) = phi*mu
# shape (a) = 1/phi
# f(x, shape=a, scale=s) = 1/(s^a Gamma(a)) x^(a-1) e^-(x/s)
# --> mu = scale/phi = scale*shape
dgamma2 = function(x,mu,phi,...) dgamma(x, shape=1/phi, scale=phi*mu, ...)
# fit the parameters of the distribution to get MLEs
fitted = MASS::fitdistr(y,dgamma2,start=list(mu=mean(y),phi=sd(y)/sqrt(n)))
mu.hat = fitted$estimate["mu"]
phi.hat = fitted$estimate["phi"]
# logLikelihood function
logL = function(mu,phi) sum(dgamma2(y, mu, phi, log=TRUE))
logL = Vectorize(logL)
# Relative logLikelihood function
maxLogL = logL(mu.hat, phi.hat)
logR = function(mu,phi) logL(mu,phi)-maxLogL
logR = Vectorize(logR)
# Find the ranges for mu and phi to plot
mu.range = c(
uniroot(function(mu) logR(mu, phi.hat)+4, c(min(y)/10, mu.hat))$root,
uniroot(function(mu) logR(mu, phi.hat)+4, c(mu.hat,10*max(y)))$root)
phi.range = c(
uniroot(function(phi) logR(mu.hat, phi)+4, c(1E-8, phi.hat))$root,
uniroot(function(phi) logR(mu.hat, phi)+4, c(phi.hat,sd(y)))$root)
# relative Likelihood contour plot
mu = seq(mu.range[1], mu.range[2], len=100)
phi = seq(phi.range[1], phi.range[2], len=100)
mat = outer(mu,phi,logR)
contour(mu,phi, exp(mat), xlab=bquote(mu), ylab=bquote(phi))
abline(h=phi.hat, col=2, lty=3)
abline(v=mu.hat, col=2, lty=3)
# profile logL
profileLogL = Vectorize(function(mu) {
phi = optimize(function(phi) logL(mu,phi), phi.range, maximum=TRUE)$maximum
logL(mu, phi)
})
# compare logL for phi=phi.hat and profileLogL
plot(function(mu) logL(mu,phi.hat), xlim=mu.range, xlab=bquote(mu), ylab="log relative Likelihood")
plot(function(mu) profileLogL(mu), xlim=mu.range, col=2, add=TRUE)
# GLM
conf = 0.9
model = glm(y~1, family=Gamma("identity"))
glm.ci = confint(model, level=conf)
abline(v=exp(coef(model)),lty=3)
abline(v=glm.ci)
abline(h= maxLogL-0.5*qchisq(conf,1), lty=2) # cuts logL, not profileLogL
Dear Dr.,
There is a relation between Normal distribution and Gamma distribution, so I suggest using the transformation technique into the Normal distribution then, estimate the CI for Normal distribution. Finally, returne the CI into Gamma distribution.
Best Regards.
This URL will be helpful
http://stats.stackexchange.com/questions/89230/confidence-interval-for-a-random-sample-selected-from-gamma-distribution
There is very nice book fo R written by Michael J. Crawley.
https://www.cs.upc.edu/~robert/teaching/estadistica/TheRBook.pdf
- Badges
- Science topic