Chapter 7 Example 2.3

The likelihood for a set of data, a model, and a set of parameters for the model is simply the probability (or probability density) for those data in that model with those parameters. For a variety of reasons, we generally work with the (natural) log of the likelihood rather than the likelihood itself. We can use log likelihood to compare models and to create confidence regions for parameters.

The R function logLik() (yes, they left off the “e”) is the basic tool for computing the log likelihood of a model fitted to data. We will use it explicitly and implicitly when we do likelihood procedures on data.

In this example, we are doing inference on the mean of the differences (Standard - Ergonomic) for the runstiching data. The statistical model is simply that the data follow a normal distribution with some mean \(\mu\) and some variance \(\sigma^2\). We use the lm() function in R (lm() for linear model) to do the estimation, and we use logLik() of the resulting fitted model to get the likelihood.

The primary argument of lm() is a “formula” that takes the form response ~ linmod, where response is the data to be fit, and linmod is a specification of how the means of response are structured. By changing linmod, we can specify that everything has mean 0, or everything has a common mean \(\mu\), or the means depend linearly on a predictor variable x, or there are two groups that have their own means, among countless possibilities.

For the runstitching differences, we need to be able to compute the likelihood for a model where we assume a particular value for \(\mu\) but let \(\sigma^2\) vary. The formula that does that for zero mean is response ~ 0. If you have some other mean, say .456, then you need to modify the response and then fit a mean of zero via response - .456 ~ 0. We also need to be able to compute the likelihood when the mean is unconstrained and estimated from the data; this model is response ~ 1. In this setting, we are estimating a coefficient or multiplier for a column of all ones, giving us the same mean for every response.

Load the runstitch data from the package and create differences between Standard and Ergonomic.

> library(cfcdae)
> data(RunStitch)
> differences <- RunStitch[,"Standard"] - RunStitch[,"Ergonomic"]

7.1 Likelihood and Likelihood Ratio Test

The likelihood ratio test (LRT) compares a statistic to a chi-squared distribution. The statistic is twice the difference in log likelihoods between the model at the null value of the parameters and the model with the parameters set to the values that maximize the likelihood. The chi-squared distribution has degrees of freedom equal to the difference in number of parameters between the two models (that is, how many parameters does the null hypothesis restrict).

The p-values from this test are approximate, with the approximation getting better as the sample size increases. In modest sized samples and simple models such as this, you would typically use a t-test.

We have two null hypotheses. Mr. Skeptical thinks that the mean should be 0, and Mr. Enthusiastic thinks that the mean should be .5; we should fit both of those models.
> null0model <- lm(differences~0)
> nullp5model <- lm(differences-.5~0)
We also need to fit the model where the mean is unconstrained and estimated from the data.
> fullmodel <- lm(differences~1)
The likelihood for the model assuming a mean of 0 is:
> null0Like <- logLik(null0model); null0Like
'log Lik.' -29.98793 (df=1)
The value -29.99 is the log likelihood, and the df=1 says that this model involves one estimated parameter (\(\sigma^2\)). The log likelihood for the model assuming a mean of .5 is lower:
> nullp5Like <- logLik(nullp5model); nullp5Like
'log Lik.' -32.37791 (df=1)
Finally, the log likelihood for the model where we can estimate the mean is highest, but we can see that it uses an additional parameter (the mean we just estimated) to achieve that higher likelihood.
> fullLike <- logLik(fullmodel); fullLike
'log Lik.' -28.88137 (df=2)
The likelihood ratio test looks at twice the difference of the log likelihoods (full model minus null model) and compares that the a chi-squared distribution with degrees of freedom equal to the additional degrees of freedom in the full model. For our null hypotheses we find:
> LRT0 <- as.numeric(2*(fullLike-null0Like));LRT0
[1] 2.213124
> LRTp5 <- as.numeric(2*(fullLike-nullp5Like));LRTp5
[1] 6.993092

Note: the as.numeric() is not technically required, but using it strips off the labeling associated with fullLike, making the printed output cleaner.

The models differ by one parameter (set the mean to zero or let it be fit), so we should compare this to a chi-squared with one degree of freedom, using the upper tail area to get the p-value.
> pchisq(LRT0, 1, lower.tail=FALSE)
[1] 0.1368413
> pchisq(LRTp5, 1, lower.tail=FALSE)
[1] 0.008182487

The p-value is about .14 for the null of 0, so there is no evidence against the null that the mean difference is 0. However, the p-value for the null of .5 is just less than .01, so there is some, but not terribly strong, evidence against that null.

7.2 Likelihood Confidence Interval

The values of a parameter in a likelihood confidence interval with coverage 1-\(\alpha\) are those null values that would not be rejected by the LRT at level \(\alpha\). This process is not automated in R for this simple model, so we need to do a little work.

Create a function that computes the log likelihood for a given null value:
> myLogLike <- function(null) {
+   logLik(lm(differences-null~0))
+ }

Next we create a dense grid of potential null mean values and evaluate the log likelihood at each potential null. The grid used here is probably more dense than necessary, and the denser the grid, the longer it takes to do the work.

> x <- seq(-.2, .6, .001) # steps of .001
> loglike <- x
> for(i in 1:length(x)) {
+   loglike[i] <- myLogLike(x[i])
+ }
Those null values that do not get rejected by the LRT are in the interval. A null value does not get rejected if its log likelihood is bigger than a cutoff. The cutoff is the maximized log likelihood minus one half of the 99th percentile (for 99% coverage). Note that the LRT, which is compared to chi-squared, uses twice the log likelihood, so we use .5 times the chi-squared when using just the likelihood. To get our interval, we compute our log likelihood cutoff and then find the range of those parameter values that correspond to log likelihoods above the cutoff. Our 99% likelihood confidence interval is -.139 to .490.
> cutoff <- logLik(fullmodel) - qchisq(.99,1)/2
> range(x[loglike>cutoff])
[1] -0.139  0.490
It’s also nice to see a plot of the log likelihood with a horizontal line at the cutoff.
> plot(x, loglike, xlab="Mean", ylab="log likelihood", type='l')
> abline(h=cutoff)
Log likelihood curve for run stitch differences with 99% cutoff

Figure 7.1: Log likelihood curve for run stitch differences with 99% cutoff