How to easily obtain the profile likelihood 95% confidence interval for a predicted value in a logistic regression model in R?

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"):

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): 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.

Thanks, merci.

Julien

• Can you do Quantile (logistic) regression on your logistic regression model? Then you would just choose the 2.5 and 97.5 percentile curves to be drawn in addition to the 50 percentile curve. See journals.sagepub.com/doi/pdf/10.1177/1536867X1101100301 Oct 21 '21 at 19:43
• Thanks for the suggestion @Ggjj11. I was aware about the quantile regression, but not the logistic quantile regression. It looks like the lrq package in R can allow such an analysis. I will definitely have a look at it. But regarding the profile likelhood CIs, do you have a suggestion to rapidly estimate the uncertainty without going through all the steps I've described? Thanks. Oct 21 '21 at 20:15
• 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.
– EdM
Oct 21 '21 at 21:58
• 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. Oct 22 '21 at 6:54
• @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. Oct 22 '21 at 10:56

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{[(\hat{\beta}_0 + \lambda\hat{\beta}_1)-\text{ln}\big(\frac{0.5}{1-0.5}\big)]^2}{\hat{\text{se}}_{0}^2 + \lambda^2\hat{\text{se}}_1^2 + 2\lambda\hat{\text{cov}}_{01}}

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.

• 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. Oct 22 '21 at 11:20
• 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?". Oct 22 '21 at 11:31
• 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. Oct 22 '21 at 23:29
• Here is a related thread on stackoverflow. One of the answers offers a drc package, but I think this produces the standard Wald interval. Oct 23 '21 at 0:03
• @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). Oct 23 '21 at 0:14

After some research on the profile likelihood function, it seems that R packages exist to estimate the confidence intervals of parameter (Beta) estimates of a logistic regression model with this approach (e.g., ProfileLikelihood), on top of the confint() function that is already available for this. However, for the predicted probability I haven't been able to find any.

For our needs, my statistician collaborator will write our own R scripts for this and if everything works fine, we'll make them available.

Two last things:

1. the dose.p() function of the MASS package seems to provide an SE at the response scale that corresponds to the Delta method, not the Wald method. We have also encountered some problems with the dose.p() function with small sample sizes. The deltavar() function of the emdbook package seems a better option. For more information, visit for instance the following B. Bolker's webpage: https://bbolker.github.io/stat4c03/HW/hw3_sol.html

2. The plot that I have included with my initial question shows the predicted probability to observe developed gonads according to fish length but contrary to what I've originally thought, the CI shown were obtained with the Wald method, i.e. the predicted SE on the logit scale is multiplied by + or - 1.96 and then these values are back-transformed on the response scale (e.g., Xu and Long 2005), which seems to be the most commonly-used approach, but under some circumstances it may provide less reliable results (Brown et al. 2003).

Our main objective will be to compare the perfomance of alternative methods to estimate the uncertainty of logistic regression model predictions for the L50 in fish.