What is the relationship between profile likelihood and confidence intervals?

Question

To make this chart I generated random samples of different size from a normal distribution with mean=0 and sd=1. Confidence intervals were then calculated using alpha cutoffs ranging from .001 to .999 (red line) with the t.test() function, the profile likelihood was calculated using the code below which I found in lecture notes put on line (I can't find the link at the moment Edit:Found it), this is shown by the blue lines. Green lines show the normalized density using the R density() function and the data is shown by the boxplots at the bottom of each chart. On the right is a caterpillar plot of the 95% confidence intervals (red) and 1/20th of max likelihood intervals (blue).

R Code used for profile likelihood:

  #mn=mean(dat)
  muVals <- seq(low,high, length = 1000)
  likVals <- sapply(muVals,
                    function(mu){
                      (sum((dat - mu)^2) /
                         sum((dat - mn)^2)) ^ (-n/2)
                    }
  )

My specific question is whether there is a known relationship between these two types of intervals and why the confidence interval appears to be more conservative for all cases except when n=3. Comments/answers about whether my calculations are valid (and a better way to do this) and the general relationship between these two types of intervals are also desired.

R code:

samp.size=c(3,4,5,10,20,1000)
cnt2<-1
ints=matrix(nrow=length(samp.size),ncol=4)
layout(matrix(c(1,2,7,3,4,7,5,6,7),nrow=3,ncol=3, byrow=T))
par(mar=c(5.1,4.1,4.1,4.1))
for(j in samp.size){


  #set.seed(200)
  dat<-rnorm(j,0,1)
  vals<-seq(.001,.999, by=.001)
  cis<-matrix(nrow=length(vals),ncol=3)
  cnt<-1
  for(ci in vals){
    x<-t.test(dat,conf.level=ci)$conf.int[1:2]
    cis[cnt,]<-cbind(ci,x[1],x[2])
    cnt<-cnt+1
  }


  mn=mean(dat)
  n=length(dat)
  high<-max(c(dat,cis[970,3]), na.rm=T)
  low<-min(c(dat,cis[970,2]), na.rm=T)
  #high<-max(abs(c(dat,cis[970,2],cis[970,3])), na.rm=T)
  #low<--high


  muVals <- seq(low,high, length = 1000)
  likVals <- sapply(muVals,
                    function(mu){
                      (sum((dat - mu)^2) /
                         sum((dat - mn)^2)) ^ (-n/2)
                    }
  )


  plot(muVals, likVals, type = "l", lwd=3, col="Blue", xlim=c(low,high),
       ylim=c(-.1,1), ylab="Likelihood/Alpha", xlab="Values",
       main=c(paste("n=",n), 
              "True Mean=0 True sd=1", 
              paste("Sample Mean=", round(mn,2), "Sample sd=", round(sd(dat),2)))
  )
  axis(side=4,at=seq(0,1,length=6),
       labels=round(seq(0,max(density(dat)$y),length=6),2))
  mtext(4, text="Density", line=2.2,cex=.8)

  lines(density(dat)$x,density(dat)$y/max(density(dat)$y), lwd=2, col="Green")
  lines(range(muVals[likVals>1/20]), c(1/20,1/20), col="Blue", lwd=4)
  lines(cis[,2],1-cis[,1], lwd=3, col="Red")
  lines(cis[,3],1-cis[,1], lwd=3, col="Red")
  lines(cis[which(round(cis[,1],3)==.95),2:3],rep(.05,2), 
        lty=3, lwd=4, col="Red")
  abline(v=mn, lty=2, lwd=2)
  #abline(h=.05, lty=3, lwd=4, col="Red")
  abline(h=0, lty=1, lwd=3)
  abline(v=0, lty=3, lwd=1)

  boxplot(dat,at=-.1,add=T, horizontal=T, boxwex=.1, col="Green")
  stripchart(dat,at=-.1,add=T, pch=16, cex=1.1)

  legend("topleft", legend=c("Likelihood"," Confidence Interval", "Sample Density"),
         col=c("Blue","Red", "Green"), lwd=3,bty="n")

  ints[cnt2,]<-cbind(range(muVals[likVals>1/20])[1],range(muVals[likVals>1/20])[2],
                     cis[which(round(cis[,1],3)==.95),2],cis[which(round(cis[,1],3)==.95),3])
  cnt2<-cnt2+1
}
par(mar=c(5.1,4.1,4.1,2.1))


plot(0,0, type="n", ylim=c(1,nrow(ints)+.5), xlim=c(min(ints),max(ints)), 
     yaxt="n", ylab="Sample Size", xlab="Values")
for(i in 1:nrow(ints)){
  segments(ints[i,1],i+.2,ints[i,2],i+.2, lwd=3, col="Blue")
  segments(ints[i,3],i+.3,ints[i,4],i+.3, lwd=3, col="Red")
}
axis(side=2, at=seq(1.25,nrow(ints)+.25,by=1), samp.size)

Answer 1

I will not give a complete answer (I have a hard time trying to understand what you are doing exactly), but I will try to clarify how profile likelihood is built. I may complete my answer later.

The full likelihood for a normal sample of size $n$ is $$L(\mu, \sigma^2) = \left( \sigma^2 \right)^{-n/2} \exp\left( - \sum_i (x_i-\mu)^2/2\sigma^2 \right).$$

If $\mu$ is your parameter of interest, and $\sigma^2$ is a nuisance parameter, a solution to make inference only on $\mu$ is to define the profile likelihood $$L_P(\mu) = L\left(\mu, \widehat{\sigma^2}(\mu) \right)$$ where $\widehat{\sigma^2}(\mu)$ is the MLE for $\mu$ fixed: $$\widehat{\sigma^2}(\mu) = \text{argmax}_{\sigma^2} L(\mu, \sigma^2).$$

One checks that $$\widehat{\sigma^2}(\mu) = {1\over n} \sum_k (x_k - \mu)^2.$$

Hence the profile likelihood is $$L_P(\mu) = \left( {1\over n} \sum_k (x_k - \mu)^2 \right)^{-n/2} \exp( -n/2 ).$$

Here is some R code to compute and plot the profile likelihood (I removed the constant term $\exp(-n/2)$):

> data(sleep)
> difference <- sleep$extra[11:20]-sleep$extra[1:10]
> Lp <- function(mu, x) {n <- length(x); mean( (x-mu)**2 )**(-n/2) }
> mu <- seq(0,3, length=501)
> plot(mu, sapply(mu, Lp, x = difference), type="l")

Link with the likelihood I’ll try to highlight the link with the likelihood with the following graph.

First define the likelihood:

L <- function(mu,s2,x) {n <- length(x); s2**(-n/2)*exp( -sum((x-mu)**2)/2/s2 )}

Then do a contour plot:

sigma <- seq(0.5,4, length=501)
mu <- seq(0,3, length=501)

z <- matrix( nrow=length(mu), ncol=length(sigma))
for(i in 1:length(mu))
  for(j in 1:length(sigma))
    z[i,j] <- L(mu[i], sigma[j], difference)

# shorter version
# z <- outer(mu, sigma, Vectorize(function(a,b) L(a,b,difference)))

contour(mu, sigma, z, levels=c(1e-10,1e-6,2e-5,1e-4,2e-4,4e-4,6e-4,8e-4,1e-3,1.2e-3,1.4e-3))

And then superpose the graph of $\widehat{\sigma^2}(\mu)$:

hats2mu <- sapply(mu, function(mu0) mean( (difference-mu0)**2 ))
lines(mu, hats2mu, col="red", lwd=2)

The values of the profile likelihood are the values taken by the likelihood along the red parabola.

You can use the profile likelihood just as a univariate classical likelihood (cf @Prokofiev’s answer). For example, the MLE $\hat\mu$ is the same.

For your confidence interval, the results will differ a little because of the curvature of the function $\widehat{\sigma^2}(\mu)$, but as long that you deal only with a short segment of it, it’s almost linear, and the difference will be very small.

You can also use the profile likelihood to build score tests, for example.

Answer 2

In a general framework, profile likelihood intervals are approximate confidence intervals. The proof of this result is essentially the same as proving that the likelihood ratio statistic is (asymptotically) approximately distributed as a $\chi^2_k$ distribution. The idea consists of inverting the hypothesis test obtained from a likelihood ratio statistic.

For example, a $0.147$-level profile likelihood interval has an approximate confidence of $95\%$.

These are classical results and therefore I will simply provide some references on this:

http://www.jstor.org/stable/2347496

http://www.stata-journal.com/sjpdf.html?articlenum=st0132

http://www.unc.edu/courses/2010fall/ecol/563/001/docs/lectures/lecture11.htm

http://en.wikipedia.org/wiki/Likelihood-ratio_test

http://en.wikipedia.org/wiki/Likelihood_function#Profile_likelihood

The following R code shows that, even for small samples, the intervals obtained with both approaches are similar (I am re-using Elvis example):

Note that you have to use the normalised profile likelihood.

data(sleep)
x <- sleep$extra[11:20]-sleep$extra[1:10]
n <- length(x)
Rp <- function(mu) {mean( (x-mean(x))^2 )^(n/2)/mean( (x-mu)^2 )^(n/2) }
Rp(mean(x))

mu <- seq(0,3, length=501)
plot(mu, sapply(mu, Rp), type="l")


Rpt<- function(mu) Rp(mu)-0.147 # Just an instrumental function

# Likelihood-confidence interval of 95% level

c(uniroot(Rpt,c(0.5,1.5))$root,uniroot(Rpt,c(1.51,3))$root)

# t confidence interval

t.test(x,conf.level=0.95)$conf.int

If we use a larger sample size, the confidence intervals are even closer:

set.seed(123)
x <- rnorm(100)
n <- length(x)
Rp <- function(mu) {mean( (x-mean(x))^2 )^(n/2)/mean( (x-mu)^2 )^(n/2) }
Rp(mean(x))

mu <- seq(-0.5,0.5, length=501)
plot(mu, sapply(mu, Rp), type="l")


Rpt<- function(mu) Rp(mu)-0.147 # Just an instrumental function

# Likelihood-confidence interval of 95% level

c(uniroot(Rpt,c(-0.4,0))$root,uniroot(Rpt,c(0,0.4))$root)

# t confidence interval

t.test(x,conf.level=0.95)$conf.int

AN IMPORTANT POINT:

Note that for specific samples different kinds of confidence intervals may differ in terms of their length or location, what really matters is their coverage. In the long run, all of them should provide the same coverage, independently on how much they differ for specific samples.

Answer 3

I will not give an overly mathematical answer, but I would like to address your central question about the relationship between CI's and profile likelihood intervals. As the other respondents have pointed out, CI's can be constructed from a profile likelihood by using the $\chi^2$ approximation to the $normalized$ likelihood ratio. The accuracty of this approach depends on one of two things being approximately true:

1. The profile log-likelihood is approximatley quadratic
2. There exists a parameter transform that makes the profile log-likelihood approximately quadratic.

The quadratic is important because it defines a normal distribution in log-scale. The more quadratic it is, the better the approximation and the resulting CIs'. Your choice of 1/20th cutoff for the likelihood intervals is equivalent to more than a 95% CI in the asymptotic limit, whcih is why the blue intervals are generally longer than the red ones.

Now, there is another issue with profile likelihood that needs some attention. If you have a lot of variables that you are profiling over, then if the number of data points per dimension is low, the profile likelihood can be very biased and optimistic. Marginal, conditional, and modified profile likelihoods are then used to reduce this bias.

So, the answer to your question is YES...the connection is the asymptotic normality of most maximum likelihood estimators, as manifested in the chi-squared distribution of the likelihood ratio.