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

Question

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. — gazza89, Jun 15 '23 at 11:00
Is it possible to do this without computing analytically the distribution? — Julien, Jun 15 '23 at 11:57
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: https://stats.stackexchange.com/a/438402/) — Sextus Empiricus, Jun 16 '23 at 08:33

COOLSerdash · Answer 1 · 2023-06-16T09:47:29.867

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:

Estimate the model and store the coefficients and variance-covariance matrix.
Draw $n$ coefficient vectors from a multivariate normal distribution with the mean and covariance matrix equal to the estimated ones in step 1.
Calculate the quantity of interest for each of the $n$ draws.
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. — COOLSerdash, Jun 16 '23 at 09:28
With car I get a bigger standard error than with a stratified bootstrap, is this expected ? — Julien, Jun 16 '23 at 09: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? — COOLSerdash, Jun 16 '23 at 09:33
@SextusEmpiricus You're probably right since the glm function shows a warning : Warning message: glm.fit: fitted probabilities numerically 0 or 1 occurred — Julien, Jun 16 '23 at 12:30
@Julien When I use your data from the question, I don't get this error. — COOLSerdash, Jun 16 '23 at 13:00
Great answer, and thank you for illustrating a nice use of clarify :) — Noah, Sep 02 '23 at 19:05

Sextus Empiricus · Answer 2 · 2023-06-16T13:59:42.530

5

Aside from

simulations
bootstrapping
delta method

another way is to

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

edited Jun 16 '23 at 13:59

answered Jun 16 '23 at 12:51

Sextus Empiricus

77,915

+1. It's interesting that your standard error is virtually identical to that from the delta method. – COOLSerdash Jun 16 '23 at 13:27
2

@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 https://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) – Sextus Empiricus Jun 16 '23 at 13:29
What is $p$ in the equation ? Is it the proportion of TRUE in the response variable (score_dicho here)? – Julien Jun 16 '23 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) – Sextus Empiricus Jun 16 '23 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}$ ? – Julien Jun 16 '23 at 14:46
$p$ is eventually a function of $x$ then it would be $$P(Y=y|x) = \begin{cases} p(x) &\qquad \text{if y = 1} \ 1-p(x) &\qquad \text{if y = 0} \ 0 &\qquad \text{else} \end{cases}$$ where $$p(x) = \frac{1}{1+e^{-(a+bx)}}$$ so with $p$ I mean the conditional probability, and $p$ is considered a (logistic) function of the features/regressors $x$. – Sextus Empiricus Jun 16 '23 at 15:05

gazza89 · Answer 3 · 2023-06-16T08:19:21.447

To expand upon my comment/answer your question:

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? — Julien, Jun 15 '23 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 — gazza89, Jun 15 '23 at 16:13
Is the first part $\hat{\beta}_0,\hat{\beta}_0$ an error? I'm not sure if I understood that part. — Shawn Hemelstrand, Jun 16 '23 at 00:38
yes, thanks, should beta $\hat{\beta}{0}, \hat{\beta}{1}$, will edit — gazza89, Jun 16 '23 at 08:19

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

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 ***

3 Answers3

Simulate 100000 coefficients from the model

Define the function that calculates the negative ratio of the coefficients

Apply the function to the simulated coefficients earlier

Plot and summaries the results

approach using car package deltaMethod

solution : 617.024

se : 16.157

approach using manual fit with optim

optim function to get the fit

solution : 617.027

se : 16.16775

Linked