Solved – Wald test and Likelihood ratio test, where do the confidence intervals on the regression coefficients come from

confidence intervalhypothesis testinglikelihood-ratiologisticregression coefficients

So I'm trying to build my own Wald test and likelihood ratio test code within a machine learning pipeline. I can get the final fitted logistic regression coefficients from liblinear. I'm coding in MATLAB.

How would I get the variance of the regression coefficients and also the confidence intervals of the regression coefficients? Clearly for a variance and confidence interval, you need a sample of multiple sets of coefficients. But I thought you only get a single set of coefficients at the end of the log likelihood optimization.

Basically trying to replicate the results in the following link
http://www.ats.ucla.edu/stat/mult_pkg/faq/general/nested_tests.htm

Best Answer

When you fit a logistic regression model, there is no closed form solution for the parameter estimates unlike in linear regression. So instead, you search over the parameter space for a set of parameter estimates that maximize the log likelihood or minimize the deviance. (I usually prefer to think in terms of the deviance, but in this case it might be better to think of maximizing the log likelihood.) The most common search procedure is the Newton-Raphson algorithm.

As you search, you could map out the shape of the log likelihood, but this isn't really done. In the process of running the Newton-Raphson algorithm, you calculate the Hessian matrix (and it shouldn't be difficult to get if you run a search algorithm that doesn't use the Hessian instead). That provides a picture of the shape in the region of the parameter space where you are currently. The Wald test of the parameter is based on the assumption that the log likelihood has the shape of a normal distribution (which it would with infinite data but may not with small samples). An estimate of the standard deviation is calculated, and that is used as the standard error. This is typically what is used to form confidence intervals for betas.

The likelihood ratio test works differently. The ratio of the likelihoods is the difference of the log likelihoods. That is, it is the difference between the likelihoods of the model when a parameter is set at two different values. Essentially always, these two values are the maximum likelihood estimate and the null value (0). The difference between these two values should be distributed as chi-squared. It is less common to use this approach to determine the confidence interval for a beta, but it can be done. You need to search over possible beta estimates and work backward to find values that constitute the limits of the interval.

There is a very useful figure (actually taken from John Fox) at the bottom of the linked page that is extremely helpful for understanding this topic.

Related Solutions

Solved – Understanding confidence intervals in Firth penalized logistic regression

The fact that firth=FALSE doesn't give similar results to glm is puzzling to me -- hopefully someone else can answer. As far as pl goes, though, you're almost always better off with profile confidence intervals. The Wald confidence intervals assume that the (implicit) log-likelihood surface is locally quadratic, which is often a bad approximation. Except that they're more computationally intensive, profile confidence intervals are always (? I would welcome counterexamples ?) more accurate. The "improved" p values you get from the Wald estimates are likely overoptimistic.

Generate data:

dd <- data.frame(X=rep(c("yes","no"),c(22,363)),
             Y=rep(c("no","yes","no"),c(22,7,356)))
with(dd,table(X,Y))

Replicate:

m_glm <-glm(Y~X,family=binomial,data=dd)
library("logistf")
m_fp <-logistf(Y~X,data=dd,pl=TRUE,firth=TRUE)
m_mp <- logistf(Y~X,data=dd,pl=TRUE,firth=FALSE)
m_fw <-logistf(Y~X,data=dd,pl=FALSE,firth=TRUE)
m_mw <-logistf(Y~X,data=dd,pl=FALSE,firth=FALSE)

Compare Wald (confint.default) with profile CIs for glm (in this case the profile intervals are actually narrower).

confint.default(m_glm)  ## {-2740, 2710}
confint(m_glm)          ## {NA, 118}

Comparing with the glm2 package (just to make sure that glm isn't doing something wonky).

library("glm2")
glm2(Y~X,family=binomial,data=dd)
## similar results to glm(...)

Solved – Confidence intervals for GLMM: bootstrap vs likelihood profile

tl;dr parametric bootstrap intervals are slightly more reliable, but much slower to compute. I would guess that either would be adequate in your case.

Likelihood profile: Likelihood profile confidence intervals are limited by the accuracy of the asymptotic approximation that differences in deviance (-2 * log-likelihood, possibly with an offset based on the saturated model) are $\chi^2$-distributed. In general they should be pretty good if the minimum number of groups in any random-effects grouping variable is large. $n=26$ is "fairly large" (various sources quote $n=40$ or $n=50$ as "large enough not to worry about it all"; Angrist and Pischke's Mostly Harmless Econometrics, which has a Hitchhiker's Guide to the Galaxy theme, gives $n=42$). So for your case these should be reasonably good, although maybe not so good if you want to compute high-precision confidence intervals (e.g. 99.9% CIs).
Bootstrapping: parametric bootstrapping is the usual approach taken for GLMMs (non-parametric bootstrap, i.e. resampling the data with replacement, has to be done in careful/non-standard ways in order to preserve the grouping structure). It doesn't make any asymptotic assumptions, but [like likelihood profile CIs] it does assume the model assumptions are adequate (e.g. conditional distribution, Normality of random effects, etc.). The main problem is that it's very slow - you have to refit the model once for every bootstrap sample. Also, there's no built-in method for PB in glmmTMB (unlike lme4), although there is a simulate() method for fitted models, so it shouldn't be too hard to put one together if you know what you're doing.