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)
> fullmodel <- lm(differences~1)
> null0Like <- logLik(null0model); null0Like
'log Lik.' -29.98793 (df=1)
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)
> fullLike <- logLik(fullmodel); fullLike
'log Lik.' -28.88137 (df=2)
> 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.
> 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])
+ }
> cutoff <- logLik(fullmodel) - qchisq(.99,1)/2
> range(x[loglike>cutoff])
[1] -0.139 0.490
> plot(x, loglike, xlab="Mean", ylab="log likelihood", type='l')
> abline(h=cutoff)

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