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

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

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.

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.
