# How to compute the standard error of the variable $- \frac{\hat{\beta}_0}{\hat{\beta}_1}$ of a logistic regression?

I've got the following data with a logistic regression model

# data
score_dicho score
<lgl>       <dbl>
1 FALSE         175
2 FALSE         175
3 FALSE         175
4 FALSE         175
5 FALSE         189
6 FALSE         189
7 FALSE         189
8 FALSE         189
9 FALSE         210
10 FALSE         210

data <- structure(list(score_dicho = c(FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, TRUE, TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, TRUE, TRUE, FALSE, TRUE, TRUE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE), score = c(175, 175, 175, 175, 189, 189, 189, 189, 210, 210, 210, 210, 221, 221, 221, 221, 235, 235, 235, 235, 247, 247, 247, 247, 275, 275, 275, 275, 278, 278, 278, 278, 288, 288, 288, 288, 317, 317, 317, 317, 329, 329, 329, 329, 348, 348, 348, 348, 362, 362, 362, 362, 375, 375, 375, 375, 387, 387, 387, 387, 403, 403, 403, 403, 412, 412, 412, 412, 431, 431, 431, 431, 445, 445, 445, 445, 451, 451, 451, 451, 462, 462, 462, 462, 472, 472, 472, 472, 485, 485, 485, 485, 497, 497, 497, 497, 503, 503, 503, 503, 512, 512, 512, 512, 525, 525, 525, 525, 535, 535, 535, 535, 547, 547, 547, 547, 556, 556, 556, 556, 578, 578, 578, 578, 593, 593, 593, 593, 601, 601, 601, 601, 614, 614, 614, 614, 628, 628, 628, 628, 650, 650, 650, 650)), row.names = c(NA, -144L), class = c("tbl_df", "tbl", "data.frame"))

mod_log <- glm(score_dicho ~ score, data = data, family = "binomial")

summary(mod_log)
# Call:
# glm(formula = score_dicho ~ score, family = "binomial", data = data)
#
# Coefficients:
#               Estimate Std. Error z value Pr(>|z|)
# (Intercept) -16.354969   4.024865  -4.063 4.83e-05 ***
# score         0.026506   0.006896   3.844 0.000121 ***


How to compute the standard error of the estimates of the variable $$- coef(Intercept)/coef(score)$$ ?

Or the confidence intervals.

NB : I'd like another method than the bootstrap.

• Unlike linear regression, there's not a bunch of known analytical results for the parameters in logistic regression, so it's much harder to compute the standard error. The most obvious is indeed to use bootstrapping, which you said you did not want to do. I think your next best option is to try to calculate the joint posterior for your two regression coefficients, which will involve a pretty arduous calculation of the evidence term. This will then give you a bivariate distribution you can sample from and you can simply take many samples, calculate the ratio and average over the samples. Jun 15, 2023 at 11:00
• Is it possible to do this without computing analytically the distribution? Jun 15, 2023 at 11:57
• If I have to compute the joint distribution, how to do so ? Jun 15, 2023 at 12:00
• The variable $\hat\beta_0/\hat\beta_1$ follows a ratio distribution with an infinite mean. So computing a standard error is not so trivial. You could however compute the standard error for the limiting/approximate distribution that follows a normal distribution (in a similar way as here: stats.stackexchange.com/a/438402) Jun 16, 2023 at 8:33

One possibility is the delta method (see also this post). Another is simulation-based inference as described in the paper by King et al. (2000) with a recent exploration of simulation-based inference by Rainey (2023). There is also the bootstrap, although you mention that you don't want to use it.

But note that the ratio estimator will be biased.

Delta method

The car package for R makes applying the delta method very easy. For your example:

library(car)

deltaMethod(mod_log, "-b0/b1", parameterNames = paste0("b", 0:1))
Estimate      SE   2.5 % 97.5 %
-b0/b1  617.024  16.157 585.358 648.69


The point estimate is $$617.02$$ with an approximate 95% confidence interval of $$(585.36;\;648.69)$$.

Simulation-based inference (CLARIFY)

The simulation-based inference involves the following steps:

1. Estimate the model and store the coefficients and variance-covariance matrix.
2. Draw $$n$$ coefficient vectors from a multivariate normal distribution with the mean and covariance matrix equal to the estimated ones in step 1.
3. Calculate the quantity of interest for each of the $$n$$ draws.
4. Calculate the point estimate and 95% confidence intervals based on the $$n$$ calculated quantities of interest by taking the average and $$0.025/0.975$$ quantiles, respectively.

The clarify package for R makes this easy. Again for your data and model:

# Reproducibility
set.seed(142857)

# Simulate 100000 coefficients from the model
s <- sim(mod_log, n = 1e5)

# Define the function that calculates the negative ratio of the coefficients
my_fun <- function(fit) {
betas <- unname(coef(fit))
c("-b0/b1" = -betas[1]/betas[2])
}

# Apply the function to the simulated coefficients earlier
est1 <- sim_apply(s, FUN = my_fun)

# Plot and summaries the results
plot(est1) # Not shown
print(summary(est1), digits = 5)
Estimate  2.5 % 97.5 %
-b0/b1   617.02 590.77 667.40

confint(est1)
2.5 %   97.5 %
-b0/b1 590.7698 667.3992


The simulation-based estimate is $$617.02$$ with a 95% confidence interval of $$(590.77;\; 667.40)$$.

Nonparametric bootstrap

To complete the answer, here is an analysis using the (nonparametric, case-resampling) bootstrap, also implemented in the car package:

library(car)
set.seed(142857) # Reproducibility
# Bootstrap
res_boot <- Boot(mod_log, f = function(obj){-coef(obj)[1]/coef(obj)[2]}, R = 999)
summary(res_boot)
R original bootBias bootSE bootMed
(Intercept) 999   617.02   1.3104 18.998  616.14

confint(res_boot)
Bootstrap bca confidence intervals

2.5 %   97.5 %
(Intercept) 588.2415 664.4407


The bootstrap standard error is $$18.998$$ with a 95% confidence interval of $$(588.24;\;664.44)$$. This is very close to the simulation-based confidence intervals.

References

King, G., Tomz, M., & Wittenberg, J. (2000). Making the most of statistical analyses: Improving interpretation and presentation. American journal of political science, 347-361.

Rainey, C. (2023). A careful consideration of CLARIFY: simulation-induced bias in point estimates of quantities of interest. Political Science Research and Methods, 1-10.

• @Julien I've used the car package because it's so easy. I also know of the msm package that also implements the delta method. But it's also fairly easy to code it "by hand" for this case. Jun 16, 2023 at 9:28
• With car I get a bigger standard error than with a stratified bootstrap, is this expected ? Jun 16, 2023 at 9:32
• @Julien It's possible because the delta method relies on a normal approximation whereas the bootstrap uses the empirical sampling distribution which can be far from normal. What exactly did you do? Jun 16, 2023 at 9:33
• Jun 16, 2023 at 9:34
• Great answer, and thank you for illustrating a nice use of clarify :)
– Noah
Sep 2, 2023 at 19:05

Aside from

1. simulations
2. bootstrapping
3. delta method

another way is to

1. re-parameterize the model

In your case this can be done by modeling the Bernoulli distribution probability parameter $$p$$ as

$$p = \frac{1}{1+e^{-(\beta_0+\beta_1x)}} = \frac{1}{1+e^{-\beta_1 (x+\beta_0/\beta_1)}} = \frac{1}{1+e^{-\beta_1 (x-\beta_2)}}$$

where the parameter $$\beta_2 = \beta_0/\beta_1$$ on the right-hand side is what you are looking for.

This parameter value and it's standard error can be estimated by using any optimization algorithm. The Hessian of the log likelihood can be used to estimate the covariance and the standard error.

Example code:

set.seed(1)
y = data$$score_dicho x = data$$score

###
### approach using car package deltaMethod
###

mod_log <- glm(y ~ x, family = binomial())
mod_log
car::deltaMethod(mod_log, "-b0/b1", parameterNames = c("b0","b1"))
### solution : 617.024
### se       :  16.157

###
### approach using manual fit with optim
###
logL2 = function(par) {
ab = par[1]
b = par[2]
p = 1/(1+exp(-b*(x-ab)))
-sum(dbinom(y,1,p, log = TRUE))
}

# optim function to get the fit
opt = optim(c(600,0.02),logL2, hessian = TRUE, control = list(parscale = c(600,0.03), ndeps = c(10^-6,10^-6)))
opt
sqrt(solve(opt$hessian)[1,1]) ### <--- this thing here is the SE estimate ### solution : 617.027 ### se : 16.16775  • +1. It's interesting that your standard error is virtually identical to that from the delta method. Jun 16, 2023 at 13:27 • @COOLSerdash I believe that, being based on the gradient, a first order approximation, it should be exactly identical and the discrepancies are only in the computation. (Here is a case where there is a difference stats.stackexchange.com/a/481446 but that is the case when the entire variable is transformed and not just the model for the parameters of the distribution) Jun 16, 2023 at 13:29 • What is$p$in the equation ? Is it the proportion of TRUE in the response variable (score_dicho here)? Jun 16, 2023 at 13:57 • @Julien I used$p$as the parameter of the assumed Bernoulli distribution $$P(Y=y|p) = \begin{cases}p &\qquad \text{if y = 1} \\ 1-p &\qquad \text{if y = 0} \\ 0 &\qquad \text{else} \end{cases}$$ for the binary observations$Y$(which in your code are, I believe data$score_dicho) Jun 16, 2023 at 14:02
• You mean $P(Y=y|x) = \begin{cases}p &\qquad \text{if y = 1} \\ 1-p &\qquad \text{if y = 0} \\ 0 &\qquad \text{else} \end{cases}$ ? Jun 16, 2023 at 14:46

When you do standard out-of-the-box logistic regression, you are calculating the maximum likelihood estimators $$\hat{\beta}_{0}, \hat{\beta}_{1}$$

which are given by

$$argmin_{\beta_{0}, \beta_{1}} \sum_{i=1}^{N}y_{i}\ln\sigma(\beta_{0} + \beta_{1}x_{1}) + (1-y_{i})\ln (1 - \sigma(\beta_{0} + \beta_{1}x_{1}))$$

where I'm using $$\sigma(x)$$ to denote $$\frac{1}{1+e^{-x}}$$

For what's to come, it will be instructive to note that this is equivalent to writing:

$$argmin_{\beta_{0}, \beta_{1}} \sum_{i=1}^{N}y_{i}\ln p(y_{i}|x_{i}, \beta_{0}, \beta_{1})=argmin_{\beta_{0}, \beta_{1}} \ln P(Data|\beta_{0}, \beta_{1})$$

If instead of having a maximum-likelihood estimator but wanted to come up with a posterior distribution for $$(\beta_{0}, \beta_{1})$$, you would do this using Bayes' rule

$$p(\beta_{0}, \beta_{1}| D) = \frac{p(D|\beta_{0}, \beta_{1})p(\beta_{0}, \beta_{1})}{\int p(D|\beta_{0}, \beta_{1})p(\beta_{0}, \beta_{1})d\beta_{0}\beta_{1}}$$

You already have an expression for the likelihood term, albeit it's a bit more complicated if you don't take the log, it's given by

$$\prod_{i=1}^{N}\sigma(\beta_{0}+ \beta_{1}x_{i})^{y_{i}}(1 - \sigma( \beta_{0}+\beta_{1}x_{i}))^{1-y_{i}}$$

and you will have to choose some suitable prior $$p(\beta_{0}, \beta_{1})$$ that corresponds to your intuition.

So in theory, I've written down an analytic expression for your posterior function, but this isn't entirely trivial to evaluate (which is why Bayesian inference and MCMC are a discipline unto themselves). The denominator is not something you can calculate analytically, you'll have to solve it numerically.

Note however that you don't technically need to know the posterior, you just need to be able to sample from it, as what you care about is getting a posterior for the quantity $$\frac{\beta_{0}}{\beta_{1}}$$. Because $$(\beta_{0},\beta_{1})$$ will be correlated, even if you had access to the joint posterior it would be non-trivial to calculate the distribution of their ratio. A more pragmatic way to proceed would be to sample $$(\beta_{0},\beta_{1})$$ pairs from the distribution, calculate the ratios every time and take an average and std over many such ratios.

This only requires you to be able to sample from the posterior, and there is a whole literature on how to sample from un-normalised distributions. In this case, that means treating your posterior as

$$p(\beta_{0}, \beta_{1}) = \alpha \cdot \prod_{i=1}^{N}\sigma(\beta_{0}+ \beta_{1}x_{i})^{y_{i}}(1 - \sigma( \beta_{0}+\beta_{1}x_{i}))^{1-y_{i}}$$

In conclusion, you can see why people tend to solve this problem by bootstrapping

• Once you have the distribution, do you use the inverse CDF technique to produce samples? Jun 15, 2023 at 15:47
• no, because it's a 2 dimensional pdf. If you have the full distribution, you can use a sampling technique that works well in low dimensions like Metropolis-Hastings or Gibbs sampling. But as I said, it's hard to calculate the full pdf as you have to evaluate that integral, it's much easier to get the un-normalised pdf (or even easier to get the log thereof), and there are tweaks one can make to the aforementioned sampling techniques that allow you to sample from unnormalised distributions Jun 15, 2023 at 16:13
• Is the first part $\hat{\beta}_0,\hat{\beta}_0$ an error? I'm not sure if I understood that part. Jun 16, 2023 at 0:38
• yes, thanks, should beta $\hat{\beta}_{0}, \hat{\beta}_{1}$, will edit Jun 16, 2023 at 8:19