6
$\begingroup$

I'm a fish biologist and we often use a logistic regression to estimate what we refer to as the L50, i.e. the length at which you expect one fish out of two (50%) to have developed gonads.

How to assess the uncertainty around the L50 estimate is not trivial. Here is an example based on 95 sampled Walleye females that vary in length (LT) from 165 to 680 mm, of which 20 had developed gonads (MATURITE coded as 1) and thus, 75 had undeveloped gonads (MATURITE coded as 0).

The data are referred to as saviMEGI14, and the binary response variable MATURITE is analyzed according to the continuous independent variable LT. The logistic regression model is coded in R like this:


summary(m.saviMEGI14.LT <- glm(MATURITE ~ LT, 
   family = binomial, data = saviMEGI14))

Using the dose.p() function of the MASS package, I can easily get the estimation of the L50 and its associated SE:


library(MASS)

dose.p(m.saviMEGI14.LT)

p = 0.5

Dose = 491.9017

SE = 20.76949


And from these estimates, calculate the Wald-based CIs, simply multipling the SE by 1.96 and adding/substracting this value from the "Dose" estimate obtained, i.e. the L50 of 492 mm. The Wald CI would thus be [451, 532].

However, in such a relatively small sample size, the Wald CIs are not ideal because they are based on normal theory, so the default method in R, i.e. the profile likelihood function, should be used instead (see for instance Royston 2007 The Stata Journal)

I am aware that I can get the parameter estimates (central value and CI) for the Intercept and the LT by using the confint() function = profile likelihood or the confint.default() = Wald; but this does not allow me to compute the SE/CIs around the L50 from the logistic regression model.

To get the profile likelihood CIs around the L50 of 492 mm, I use the predict() function applied to the range of LT values by first creating a new data frame called below nd14_LT:


nd14_LT <- data.frame(LT=seq(from=165, to=680, 
              by=1))

and then get the predicted values for the nd14_LT data, asking also to obtain the associated predicted SE on the logit scale (type="link"):


pred14_LT <- predict(m.saviMEGI14.LT, nd14_LT, 
               type="link", se.fit=TRUE)

With the predicted values and their associated predicted SE, one can then convert them from the logit to the response scale and then find at which LT a probability of 0.5 is found for both the lower and upper CI bounds.

Doing so provide a profile likelihood-based CI of [457, 551], which is quite different than the Wald one [452, 532], especially for its upper portion.

Here is the plot showing the regression curve, as well as the central value and profile likelihood CIs (which is called a maturity ogive in fisheries science):

MATURITY OGIVE

All this to come up with this question:

How can someone obtain the lower and upper bounds of the profile likelihood function CI from a logistic regression conducted in R in a rapid manner?

With colleagues, we are using simulations to test different methods for the estimation of the uncertainty around the L50 (i.e., parametric and non-parametric bootstrapping, Fieller analytical method, credible intervals from the Bayesian approach, etc…) and we’d like to find a way to estimate the profile likelihood CI of a given dataset in a timely, valid manner.

$\endgroup$
14
  • 2
    $\begingroup$ A warning: profile-likelihood CI aren't necessarily the default in R. IF the MASS package is installed then calling confint() on a glm or nls object will end up calling the corresponding MASS functions, which do use profile likelihoods. Otherwise, the Wald CI based on the covariance matrix of the coefficient estimates will be used. $\endgroup$
    – EdM
    Oct 21, 2021 at 21:58
  • 1
    $\begingroup$ In addition to what @EdM said: Do you have a reference that confirms that confidence intervals obtained by predict are profile likelihood intervals? I've never heard that term in that context. Thank you. $\endgroup$ Oct 22, 2021 at 6:54
  • 1
    $\begingroup$ @Edm I've ran the same model on the same data, MASS being installed, but not loaded into active R session, and use first the confint() and obtain the message "Waiting for profiling to be done..." indicating that profile likelihood CIs were computed. Using the confint.default() provided me with narrower CIs for the parameter estimates. I couldn't use dose.p(), as MASS was not loaded into R. After doing so, I've got the same results. But thanks for the warning note: someone who does not have MASS installed could get different result, but I don't know for this. $\endgroup$ Oct 22, 2021 at 10:56
  • 1
    $\begingroup$ I think these are simply backtransformed Wald-based confidence intervals on the logit scale (see Hosmer DW, Lemeshow S, Sturdivant RX (2013): Applied Logistic Regression. 3rd ed. Wiley, page 17). $\endgroup$ Oct 22, 2021 at 12:02
  • 1
    $\begingroup$ The standard error of dose.p is based on the delta method. The resulting confidence interval could also be described as a Wald-based interval, but the standard error is estimated differently. $\endgroup$ Oct 22, 2021 at 12:32

3 Answers 3

6
$\begingroup$

My collaborator Rafael de Andrade Moral has written a R function that allows to estimate the uncertainty for L50 estimates according to 13 different methods, including the profile-likelihood method. We used this R function in our recently-published paper in Fisheries Research:

Monitoring reproduction in fish: Assessing the adequacy of ogives and the predicted uncertainty of their L50 estimates for more reliable biological inferences

to evaluate the coverage probability of different approaches, such as the Delta method, non-parametric bootstrapping and Fieller (1944)'s analytical approach. The profile-likehood method did pretty good but the Monte Carlo approach with BCa interval led to the best nominal coverage for the simulated dataset that we have analyzed.

The R scripts for the profile-likelihood method (and others) to assess uncertainty in the L50 (or else such as the L95) is available at:

https://github.com/rafamoral/L50/blob/main/confint_L.R

The supplementary material of our MS contains information on how to use the confint_L() function.

Cheers,

$\endgroup$
5
$\begingroup$

EDIT: I see you mentioned the Fieller method in your original post. Perhaps you were referring to the solution I provided below.

Here is a great paper on the topic. Using a logistic regression with a logit link function you can model the proportion of fish as a function of length, with $\lambda:=$LD50. Based on the asymptotic normality of

$$ \frac{(\hat{\beta}_0 + \lambda\hat{\beta}_1)-\text{ln}\big(\frac{0.5}{1-0.5}\big)}{\sqrt{\hat{\text{se}}_{0}^2 + \lambda^2\hat{\text{se}}_1^2 + 2\lambda\hat{\text{cov}}_{01}}}$$

a $100(1-\alpha)\%$ confidence interval is found by identifying the set of $\lambda$ that satisfy

$$ \frac{\left[(\hat{\beta}_0 + \lambda\hat{\beta}_1)-\text{ln}\big(\frac{0.5}{1-0.5}\big)\right]^2}{\hat{\text{se}}_{0}^2 + \lambda^2\hat{\text{se}}_1^2 + 2\lambda\hat{\text{cov}}_{01}}<z_{\alpha}^2$$

where $\hat{\text{se}}_{0}$ is the estimated standard error of $\hat{\beta}_0$, $\hat{\text{se}}_{1}$ is the estimated standard error of $\hat{\beta}_1$, and $\hat{\text{cov}}_{01}$ is the estimated covariance between $\hat{\beta}_0$ and $\hat{\beta}_1$. This works well even in small sample sizes and is a much better normal approximation than a Wald interval for $\lambda$ based on an identity link using the dose.p() output. The confidence interval above can be calculated using standard output from the logistic regression without calling dose.p(), and should perform similarly to the likelihood ratio CI you are interested in. The only part that would require some work is numerically inverting the quantity above. You can create a sequence of values for $\lambda$, evaluate the quantity above for each value of $\lambda$, and identify those values that satisfy the inequality.

A great way to visualize this is to define and plot the following functions $$H(\lambda)=1-\Phi\Bigg(\frac{(\hat{\beta}_0 + \lambda\hat{\beta}_1)-\text{ln}\big(\frac{0.5}{1-0.5}\big)}{\sqrt{\hat{\text{se}}_{0}^2 + \lambda^2\hat{\text{se}}_1^2 + 2\lambda\hat{\text{cov}}_{01}}}\Bigg)$$

$$H^{\text{-}}(\lambda)=\Phi\Bigg(\frac{(\hat{\beta}_0 + \lambda\hat{\beta}_1)-\text{ln}\big(\frac{0.5}{1-0.5}\big)}{\sqrt{\hat{\text{se}}_{0}^2 + \lambda^2\hat{\text{se}}_1^2 + 2\lambda\hat{\text{cov}}_{01}}}\Bigg)$$

\begin{eqnarray} C(\lambda)= \left\{ \begin{array}{cc} H(\lambda) & \text{if } \lambda\le \hat{\lambda}(\boldsymbol{y}) \\ & \nonumber\\ H^{\text{-}}(\lambda) & \text{if } \lambda\ge \hat{\lambda}(\boldsymbol{y}). \end{array} \right.\nonumber \end{eqnarray}.

where $\hat{\lambda}(\boldsymbol{y})$ is the estimate of LD50 based on the observed data. $C(\lambda)$ is called a confidence curve and depicts p-values and confidence intervals of all levels. In small sample sizes the performance of this interval might be improved by referencing a $t$-distribution with $n-1$ degrees of freedom instead of a standard normal distribution.

If you are still interested in the likelihood ratio test you can create a similar confidence curve:

$$p:=\text{logit}^{-1}({\beta}_0 + \lambda{\beta}_1)$$ $$L(\beta_0,\beta_1)\propto \prod_{i=1}^n \text{logit}^{-1}({\beta}_0 + x_i{\beta}_1)^{y_i}\times[1-\text{logit}^{-1}({\beta}_0 + x_i{\beta}_1)]^{1-y_i}$$ $$\text{LR}=\frac{L(\tilde{\beta}_0,\tilde{\beta}_1)}{L(\hat{\beta}_0,\hat{\beta}_1)}$$ where $\tilde{\beta}_0$ and $\tilde{\beta}_1$ are estimates calculated under the restricted null space for $\lambda$.

\begin{eqnarray} H(\lambda)= \left\{ \begin{array}{cc} \big[1-F_{\chi^2_1}\big(-2\text{log(LR)}\big)\big]/2 & \text{if } \lambda\le \hat{\lambda}(\boldsymbol{y}) \\ & \nonumber\\ \big[1+F_{\chi^2_1}\big(-2\text{log(LR)}\big)\big]/2 & \text{if } \lambda\gt \hat{\lambda}(\boldsymbol{y}). \end{array} \right.\nonumber \end{eqnarray}.

\begin{eqnarray} C(\lambda)= \left\{ \begin{array}{cc} H(\lambda) & \text{if } \lambda\le \hat{\lambda}(\boldsymbol{y}) \\ & \nonumber\\ 1-H(\lambda) & \text{if } \lambda\ge \hat{\lambda}(\boldsymbol{y}). \end{array} \right.\nonumber \end{eqnarray}

where $F_{\chi^2_1}$ is the CDF of a chi-square distribution with 1 degree of freedom. Because the likelihood ratio confidence interval requires profiling nuisance parameters it is almost as computationally intensive as iterative methods such as bootstrap and Monte Carlo approaches.

$\endgroup$
6
  • $\begingroup$ Thanks Geoffrey Johnson for the suggestion to what I understand to be a modification of Fieller (1944) analytical method. Our goal, however, is to get more commonly-used profile likelihood function-based CIs from R, on top of those of Wald, to compare them with those obtained with a Bayesian approach, and iterative procedures, such as parametric and non-parametric bootstrapping. All this through simulations for statistical comparisons of different effect and sample sizes. Computing profile likelihood-based CIs manually is not a option. And this thus needs to be straightforward. Thanks. $\endgroup$ Oct 22, 2021 at 11:20
  • $\begingroup$ Here is a good (and possible sole) example of what we have in mind to contrast different methods (here: Fieller, non-parametric bootstrapping and Monte Carlo algorithm) to assess the uncertainty around the L50 in fish: spo.nmfs.noaa.gov/sites/default/files/13roafis.pdf We'd like to test a broader range of methods, mostly to be able to tell: "do the more simple, non-interative procedures such as the profile-likelihood function, allow us to properly conduct statisticals test about the L50?". $\endgroup$ Oct 22, 2021 at 11:31
  • $\begingroup$ I am more familiar with SAS, but I'll continue looking to find a ready-made R function for profile likelihood ratio confidence limits for LD50. It might be that you will need to write your own R function to produce the confidence interval, and call this function within a simulation to assess its performance over many samples. $\endgroup$ Oct 22, 2021 at 23:29
  • $\begingroup$ Here is a related thread on stackoverflow. One of the answers offers a drc package, but I think this produces the standard Wald interval. $\endgroup$ Oct 23, 2021 at 0:03
  • $\begingroup$ @julienbio99, you might be able to improve the performance of the Fieller method by referencing a t-distribution with n-1 degrees of freedom instead of a standard normal distribution and show that this works as well or better than the iterative methods. The likelihood ratio method is almost as computationally intensive as the iterative methods because of the profiling (re-estimating nuisance parameters under the restricted null space). $\endgroup$ Oct 23, 2021 at 0:14
4
$\begingroup$

TL;DR @kjetilbhalvorsen explains how to compute confidence intervals based on the profile likelihood in Confidence interval for difference between two predicted probabilities in R. The entire thread is worth reading as it describes several methods to compute confidence intervals. I summarize the profile likelihood method briefly and then apply it to calculate a CI for the median lethal dose, LD50, on the tobacco budworms dataset. See Section 7.2 of Modern Applied Statistics with S by Venables and Ripley.

Let $\theta(\boldsymbol{\beta})$ be a scalar function of the model parameters $\boldsymbol{\beta}$. Then the profile likelihood of $\theta$ is: $$ \begin{aligned} L_P(\theta) = \max_{\theta:\theta(\boldsymbol{\beta})=\theta}L(\boldsymbol{\beta}) \end{aligned} $$ where $L(\boldsymbol{\beta})$ is the likelihood. For the dose-response example, $\boldsymbol{\beta} = (\beta_0,\beta_1)$ but the theory — and the R code — are straightforward to extend to more parameters. The LD50, the dose expected to kill 50% of individuals, is the point where the log-odds of survival are equal to zero: $\theta = -\beta_0/\beta_1$.

x <- budworm$x # dose
y <- budworm$y # survival (0 or 1)

# The LD50 (dose that is expected to kill 50% of individuals) is defined as
# the point where the log-odds of survival are equal to zero.
theta <- function(beta0, beta1) {
  -beta0 / beta1
}

# The negative log-likelihood of the logistic regression model
# with an intercept and one predictor x.
negloglik_lrm <- function(beta0, beta1) {
  prob <- plogis(beta0 + beta1 * x)
  -sum(dbinom(y, 1, prob, log = TRUE))
}

# Fit a simple linear regression to find reasonable starting values.
initialize_betas <- function(x, y) {
  betas <- coef(lm(y ~ x))
  list(beta0 = betas[1], beta1 = betas[2])
}

The crux of the method is to profile the likelihood, ie, to maximize the likelihood over a grid of $\theta$ values. For a fixed $\theta$ we find $\max L(\boldsymbol{\beta})$ over the set $\left\{\boldsymbol{\beta}:\theta(\boldsymbol{\beta})=\theta\right\}$.

# Minimize the negative log-likelihood instead of maximizing the likelihood
model <- bbmle::mle2(negloglik_lrm, start = initialize_betas(x, y))

# nlp = negative profile log-likelihood
grid <- profile_negloglik(model)
grid
#> # A tibble: 201 × 2
#>    theta   npl
#>    <dbl> <dbl>
#>  1  1.65  141.
#>  2  1.67  140.
#>  3  1.68  139.
#>  4  1.70  138.
#>  5  1.71  137.
#>  6  1.73  136.
#>  7  1.74  135.
#>  8  1.75  134.
#>  9  1.77  133.
#> 10  1.78  132.
#> # … with 191 more rows

Now that we have calculated the profile likelihood for $\theta$ at a grid of values, it remains to find a subset of those to represent a confidence interval.

Under regularity conditions (which mean the log likelihood is well approximated by a quadratic function), the ratio likelihood test states that under the null hypothesis $H_0:\theta=\theta_0$: $$ \begin{aligned} 2\log\frac{L(\widehat{\theta})}{L(\theta_0)} \sim \chi^2_p \end{aligned} $$ where $p$ is the number of parameters $\theta$, $\widehat{\theta}$ is the maximum likelihood estimate (MLE) and $\chi^2_p$ is the Chi-squared distribution with $p$ degrees of freedom. This is the value at which the (profile) likelihood is minimized.

By inverting the hypothesis test we get a 100(1-$\theta$)% confidence interval (or region) for $\theta$.

$$ \begin{aligned} \left\{\theta:2\log\frac{L(\widehat{\theta})}{L(\theta)}<\chi^2_{p,(1-\alpha)}\right\} \end{aligned} $$

For the budworms example, the LD50 $\theta$ is a scalar, so $p=1$ and we can find the profile likelihood 95% and 99% confidence intervals numerically:

#>   conf.level conf.low conf.high      nll
#> 1       0.95 2.415057  3.084742 113.6548
#> 2       0.99 2.312052  3.195448 115.0515

I also use the MASS:dose.p function to compute the Wald 95% and 99% confidence intervals. (See R code attached.)

data.frame(
  conf.level,
  conf.low = dose.hat + dose.se * qnorm((1 - conf.level) / 2),
  conf.high = dose.hat - dose.se * qnorm((1 - conf.level) / 2)
)
#>   conf.level conf.low conf.high
#> 1       0.95 2.422786  3.071986
#> 2       0.99 2.320789  3.173983

In this example, the profile and the Wald confidence intervals for $\theta$ agree well.

References

In addition to @kjetilbhalvorsen's posts about profile likelihood, I'm also reading "In All Likelihood: Statistical Modelling And Inference Using Likelihood" by Yudi Pawitan. This book explains the theory in great detail and with many examples. (Not about budworms though.)


R code to compute profile likelihood confidence intervals:

library("MASS")
library("tidyverse")

# Compute negative log-likelihood-based intervals for a scalar parameter theta
# at the specified alpha levels.
# This implementation is based on the program `li.r` for computing likelihood
# intervals which accompanies the book "In All Likelihood" by Yudi Pawitan.
# https://www.meb.ki.se/sites/yudpaw/book/
confint_negloglik <- function(theta, nll, conf.level = 0.95) {
  nll_min <- min(nll)
  
  # 2*{ log L(theta_mle) - log L(theta) } < qchisq(df=1,p=1-alpha)
  nll_cut <- nll_min + qchisq(conf.level, 1) / 2
  theta_mle <- mean(theta[nll == nll_min])

  theta_below <- theta[theta < theta_mle]
  if (length(theta_below) < 2) {
    conf.low <- min(theta)
  } else {
    nll.below <- nll[theta < theta_mle]
    conf.low <- approx(nll.below, theta_below, xout = nll_cut)$y
  }

  theta_above <- theta[theta > theta_mle]
  if (length(theta_above) < 2) {
    conf.high <- max(theta)
  } else {
    nll.above <- nll[theta > theta_mle]
    conf.high <- approx(nll.above, theta_above, xout = nll_cut)$y
  }

  data.frame(conf.level, conf.low, conf.high, nll = nll_cut)
}

grid_params <- function(model, pts = 100) {
  seq_range <- function(x, n) {
    seq(min(x), max(x), len = n)
  }
  # Use the confidence intervals to find reasonable ranges for the betas.
  range_betas <- bbmle::confint(bbmle::profile(model))
  expand_grid(
    # Create a grid (= an outer product) for beta0 and beta1
    beta0 = seq_range(range_betas["beta0", ], pts),
    beta1 = seq_range(range_betas["beta1", ], pts)
  ) %>%
    # For each point in the grid, compute the parameter of interest theta
    # and the negative log-likelihood of the model. Both are functions of
    # the regression coefficients beta0 and beta1.
    mutate(
      theta = theta(beta0, beta1),
      nll = map2_dbl(beta0, beta1, negloglik_lrm)
    )
}

profile_negloglik <- function(profile, bins = 201) {
  profile %>%
    grid_params() %>%
    group_by(
      # Subdivide the observed range of theta into many small intervals.
      cut_interval(theta, bins),
      .drop = TRUE
    ) %>%
    summarise(
      # For each interval, find its midpoint as well as the minimum of
      # the negative log-likelihood. This is the negative profile
      # log-likelihood (npl).
      theta = (min(theta) + max(theta)) / 2,
      npl = min(nll)
    ) %>%
    select(
      theta, npl
    )
}

plot_profile <- function(grid, confint) {
  grid %>%
    ggplot(
      aes(theta, npl)
    ) +
    geom_line(
      linewidth = 1
    ) +
    geom_segment(
      aes(
        x = conf.low, xend = conf.high,
        y = nll, yend = nll,
        color = factor(conf.level)
      ),
      inherit.aes = FALSE,
      data = confint,
      linewidth = 1
    ) +
    guides(
      color = guide_legend(title = "confidence level")
    )
}

# Tobacco budworms survival data. See `MASS::dose.p`.
ldose <- rep(0:5, 2)
numdead <- c(1, 4, 9, 13, 18, 20, 0, 2, 6, 10, 12, 16)
numalive <- 20 - numdead

budworm <- list(
  x = c(rep(ldose, numdead), rep(ldose, numalive)),
  y = c(rep(0, sum(numdead)), rep(1, sum(numalive)))
)

conf.level <- c(0.95, 0.99)

x <- budworm$x
y <- budworm$y

# The negative log-likelihood of the logistic regression model
# with an intercept and one predictor x.
negloglik_lrm <- function(beta0, beta1) {
  prob <- plogis(beta0 + beta1 * x)
  -sum(dbinom(y, 1, prob, log = TRUE))
}

# The LD50 (dose that is expected to kill 50% of individuals) is defined as
# the point where the log-odds of survival are equal to zero.
theta <- function(beta0, beta1) {
  -beta0 / beta1
}

# Fit a simple linear regression to find reasonable starting values.
initialize_betas <- function(x, y) {
  betas <- coef(lm(y ~ x))
  list(beta0 = betas[1], beta1 = betas[2])
}

# Minimize the negative log-likelihood instead of maximizing the likelihood
model <- bbmle::mle2(negloglik_lrm, start = initialize_betas(x, y))

profile <- bbmle::profile(model)
grid <- profile_negloglik(model)
grid

# Profile confidence intervals for the betas parameters and for theta
ci.betas <- bbmle::confint(profile)
ci.betas
ci.theta <- confint_negloglik(grid$theta, grid$npl, conf.level)
ci.theta

plt <- plot_profile(grid, ci.theta)
plt +
  labs(
    x = expression(paste(θ, "=", -beta[0] / beta[1])),
    y = "negative profile log-likelihood",
    title = "Negative profile log-likelihood for LD50"
  )

model <- glm(
  y ~ x,
  family = binomial
)

dose.p(model)
dose.hat <- 2.747386
dose.se <- 0.1656153

# Wald confidence interval
data.frame(
  conf.level,
  conf.low = dose.hat + dose.se * qnorm((1 - conf.level) / 2),
  conf.high = dose.hat - dose.se * qnorm((1 - conf.level) / 2)
)
$\endgroup$
1
  • $\begingroup$ Thanks for this very detailed answer and the codes provided. Really appreciated. $\endgroup$ Dec 30, 2022 at 3:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.