{"id":739,"date":"2016-10-22T10:16:15","date_gmt":"2016-10-22T14:16:15","guid":{"rendered":"http:\/\/www.clayford.net\/statistics\/?p=739"},"modified":"2020-09-13T10:17:56","modified_gmt":"2020-09-13T14:17:56","slug":"profile-likelihood-ratio-confidence-intervals","status":"publish","type":"post","link":"https:\/\/www.clayford.net\/statistics\/profile-likelihood-ratio-confidence-intervals\/","title":{"rendered":"Profile likelihood ratio confidence intervals"},"content":{"rendered":"

When you fit a generalized linear model (GLM) in R and call confint<\/code> on the model object, you get confidence intervals for the model coefficients. But you also get an interesting message:<\/p>\n

Waiting for profiling to be done...<\/code><\/p>\n

What's that all about? What exactly is being profiled? Put simply, it's telling you that it's calculating a profile likelihood ratio confidence interval<\/em>.<\/p>\n

The typical way to calculate a 95% confidence interval is to multiply the standard error of an estimate by some normal quantile such as 1.96 and add\/subtract that product to\/from the estimate to get an interval. In the context of GLMs, we sometimes call that a Wald confidence interval. <\/p>\n

Another way to determine an upper and lower bound of plausible values for a model coefficient is to find the minimum and maximum value of the set of all coefficients that satisfy the following:<\/p>\n

\\[-2\\log\\left(\\frac{L(\\beta_{0}, \\beta_{1}|y_{1},…,y_{n})}{L(\\hat{\\beta_{0}}, \\hat{\\beta_{1}}|y_{1},…,y_{n})}\\right) < \\chi_{1,1-\\alpha}^{2}\\]<\/p>\n

Inside the parentheses is a ratio of likelihoods<\/em>. In the denominator is the likelihood of the model we fit. In the numerator is the likelihood of the same model but with different coefficients<\/em>. (More on that in a moment.) We take the log of the ratio and multiply by -2. This gives us a likelihood ratio test (LRT) statistic. This statistic is typically used to test whether a coefficient is equal to some value, such as 0, with the null likelihood in the numerator (model without coefficient, that is, equal to 0) and the alternative or estimated likelihood in the denominator (model with coefficient). If the LRT statistic is less than \\(\\chi_{1,0.95}^{2} \\approx 3.84\\), we fail to reject the null. The coefficient is statisically not much different from 0. That means the likelihood ratio is close to 1. The likelihood of the model without the coefficient is almost as high the model with it. On the other hand, if the ratio is small, that means the likelihood of the model without the coefficient is much smaller than the likelihood of the model with the coefficient. This leads to a larger LRT statistic since it's being log transformed, which leads to a value larger than 3.84 and thus rejection of the null.<\/p>\n

Now in the formula above, we are seeking all such coefficients in the numerator that would make it a true statement. You might say we're “profiling” many different null values and their respective LRT test statistics. Do they fit the profile of a plausible coefficient value in our model?<\/em> The smallest value we can get without violating the condition becomes our lower bound, and likewise with the largest value. When we're done we'll have a range of plausible values for our model coefficient that gives us some indication of the uncertainly of our estimate. <\/p>\n

Let's load some data and fit a binomial GLM to illustrate these concepts. The following R code comes from the help page for confint.glm<\/code>. This is an example from the classic Modern Applied Statistics with S<\/a>. ldose<\/code> is a dosing level and sex<\/code> is self-explanatory. SF<\/code> is number of successes and failures, where success is number of dead worms. We're interested in learning about the effects of dosing level and sex on number of worms killed. Presumably this worm is a pest of some sort.<\/p>\n

# example from Venables and Ripley (2002, pp. 190-2.)\r\nldose <- rep(0:5, 2)\r\nnumdead <- c(1, 4, 9, 13, 18, 20, 0, 2, 6, 10, 12, 16)\r\nsex <- factor(rep(c("M", "F"), c(6, 6)))\r\nSF <- cbind(numdead, numalive = 20-numdead)\r\nbudworm.lg <- glm(SF ~ sex + ldose, family = binomial)\r\nsummary(budworm.lg)\r\n<\/code><\/pre>\n
## \r\n## Call:\r\n## glm(formula = SF ~ sex + ldose, family = binomial)\r\n## \r\n## Deviance Residuals: \r\n##      Min        1Q    Median        3Q       Max  \r\n## -1.10540  -0.65343  -0.02225   0.48471   1.42944  \r\n## \r\n## Coefficients:\r\n##             Estimate Std. Error z value Pr(>|z|)    \r\n## (Intercept)  -3.4732     0.4685  -7.413 1.23e-13 ***\r\n## sexM          1.1007     0.3558   3.093  0.00198 ** \r\n## ldose         1.0642     0.1311   8.119 4.70e-16 ***\r\n## ---\r\n## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\r\n## \r\n## (Dispersion parameter for binomial family taken to be 1)\r\n## \r\n##     Null deviance: 124.8756  on 11  degrees of freedom\r\n## Residual deviance:   6.7571  on  9  degrees of freedom\r\n## AIC: 42.867\r\n## \r\n## Number of Fisher Scoring iterations: 4\r\n<\/pre>\n
The coefficient for ldose<\/code> looks significant. Let's determine a confidence interval for the coefficient using the confint<\/code> function. We call confint<\/code> on our model object, budworm.lg<\/code> and use the parm<\/code> argument to specify that we only want to do it for ldose<\/code>:<\/p>\n
confint(budworm.lg, parm = "ldose")\r\n<\/code><\/pre>\n
## Waiting for profiling to be done...\r\n<\/pre>\n
##     2.5 %    97.5 % \r\n## 0.8228708 1.3390581\r\n<\/pre>\n
We get our “waiting” message though there really was no wait. If we fit a larger model and request multiple confidence intervals, then there might actually be a waiting period of a few seconds. The lower bound is about 0.8 and the upper bound about 1.32. We might say every increase in dosing level increase the log odds of killing worms by at least 0.8. We could also exponentiate to get a CI for an odds ratio estimate:<\/p>\n
exp(confint(budworm.lg, parm = "ldose"))\r\n<\/code><\/pre>\n
## Waiting for profiling to be done...\r\n<\/pre>\n
##    2.5 %   97.5 % \r\n## 2.277027 3.815448\r\n<\/pre>\n
The odds of “success” (killing worms) is at least 2.3 times higher at one dosing level versus the next lower dosing level.<\/p>\n
To better understand the profile likelihood ratio confidence interval, let's do it “manually”. Recall the denominator in the formula above was the likelihood of our fitted model. We can extract that with the logLik<\/code> function:<\/p>\n
den <- logLik(budworm.lg)\r\nden\r\n<\/code><\/pre>\n
## 'log Lik.' -18.43373 (df=3)\r\n<\/pre>\n
The numerator was the likelihood of a model with a different<\/em> coefficient. Here's the likelihood of a model with a coefficient of 1.05:<\/p>\n
num <- logLik(glm(SF ~ sex + offset(1.05*ldose), family = binomial))\r\nnum\r\n<\/code><\/pre>\n
## 'log Lik.' -18.43965 (df=2)\r\n<\/pre>\n
Notice we used the offset<\/code> function. That allows us to fix the coefficient to 1.05 and not have it estimated.<\/p>\n
Since we already extracted the log<\/em> likelihoods, we need to subtract them. Remember this rule from algebra?<\/p>\n
\\[\\log\\frac{M}{N} = \\log M – \\log N\\]<\/p>\n
So we subtract the denominator from the numerator, multiply by -2, and check if it's less than 3.84, which we calculate with qchisq(p = 0.95, df = 1)<\/code><\/p>\n
-2*(num - den)\r\n<\/code><\/pre>\n
## 'log Lik.' 0.01184421 (df=2)\r\n<\/pre>\n
-2*(num - den) < qchisq(p = 0.95, df = 1)\r\n<\/code><\/pre>\n
## [1] TRUE\r\n<\/pre>\n
It is. 1.05 seems like a plausible value for the ldose<\/code> coefficient. That makes sense since the estimated value was 1.0642. Let's try it with a larger value, like 1.5:<\/p>\n
num <- logLik(glm(SF ~ sex + offset(1.5*ldose), family = binomial))\r\n-2*(num - den) < qchisq(p = 0.95, df = 1)\r\n<\/code><\/pre>\n
## [1] FALSE\r\n<\/pre>\n
FALSE. 1.5 seems too big to be a plausible value for the ldose<\/code> coefficient. <\/p>\n
Now that we have the general idea, we can program a while<\/code> loop to check different values until we exceed our threshold of 3.84. <\/p>\n
cf <- budworm.lg$coefficients[3]  # fitted coefficient 1.0642\r\ncut <- qchisq(p = 0.95, df = 1)   # about 3.84\r\ne <- 0.001                        # increment to add to coefficient\r\nLR <- 0                           # to kick start our while loop \r\nwhile(LR < cut){\r\n  cf <- cf + e\r\n  num <- logLik(glm(SF ~ sex + offset(cf*ldose), family = binomial))\r\n  LR <- -2*(num - den)\r\n}\r\n(upper <- cf)\r\n<\/code><\/pre>\n
##    ldose \r\n## 1.339214\r\n<\/pre>\n
To begin we save the original coefficient to cf<\/code>, store the cutoff value to cut<\/code>, define our increment of 0.001 as e<\/code>, and set LR<\/code> to an initial value of 0. In the loop we increment our coefficient estimate which is used in the offset<\/code> function in the estimation step. There we extract the log likelihood and then calculate LR<\/code>. If LR<\/code> is less than cut<\/code> (3.84), the loop starts again with a new coefficient that is 0.001 higher. We see that our upper bound of 1.339214 is very close to what we got above using confint<\/code> (1.3390581). If we set e<\/code> to smaller values we'll get closer.<\/p>\n
We can find the LR profile lower bound in a similar way. Instead of adding the increment we subtract it:<\/p>\n
cf <- budworm.lg$coefficients[3]  # reset cf\r\nLR <- 0                           # reset LR \r\nwhile(LR < cut){\r\n  cf <- cf - e\r\n  num <- logLik(glm(SF ~ sex + offset(cf*ldose), family = binomial))\r\n  LR <- -2*(num - den)\r\n}\r\n(lower <- cf)\r\n<\/code><\/pre>\n
##    ldose \r\n## 0.822214\r\n<\/pre>\n
The result, 0.822214, is very close to the lower bound we got from confint<\/code> (0.8228708).<\/p>\n
This is a very<\/em> basic implementation of calculating a likelihood ratio confidence interval. It is only meant to give a general sense of what's happening when you see that message Waiting for profiling to be done...<\/code>. I hope you found it helpful. To see how R does it, enter getAnywhere(profile.glm)<\/code> in the console and inspect the code. It's not for the faint of heart.<\/p>\n
I have to mention the book Analysis of Categorical Data with R<\/a>, from which I gained a better understanding of the material in this post. The authors have kindly shared their R code at the following web site if you want to have a look: http:\/\/www.chrisbilder.com\/categorical\/<\/a><\/p>\n
To see how they “manually” calculate likelihood ratio confidence intervals, go to the following R script and see the section “Examples of how to find profile likelihood ratio intervals without confint()”: http:\/\/www.chrisbilder.com\/categorical\/Chapter2\/Placekick.R<\/a><\/p>\n","protected":false},"excerpt":{"rendered":"
When you fit a generalized linear model (GLM) in R and call confint on the model object, you get confidence… Read more<\/a><\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[22,7,12,13],"tags":[54,64,33],"class_list":["post-739","post","type-post","status-publish","format-standard","hentry","category-hypothesis-testing","category-probability","category-regression","category-using-r","tag-glm","tag-likelihood-profile","tag-logistic-regression"],"_links":{"self":[{"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/posts\/739","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/comments?post=739"}],"version-history":[{"count":2,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/posts\/739\/revisions"}],"predecessor-version":[{"id":742,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/posts\/739\/revisions\/742"}],"wp:attachment":[{"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/media?parent=739"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/categories?post=739"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/tags?post=739"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}