I have two variables (continuous $x_1$, control/treatment $x_2$) that I want to use to predict a probability. Domain knowledge suggests that the relationship is roughly linear in the log-odds, so I am going with a logistic regression on the two variables and their interaction.

$$ \log\Bigg( \dfrac{p}{1-p} \Bigg) = \beta_0 + \beta_1x_1 + \beta_2 x_2 + \beta_3 x_1x_2 $$

Unlike for my previous work with GLM inference, we are interested in a particular probability $p_0$ and by how much the two groups differ in their $x_1$ value when $p_0$ is reached.

For a model of each group individually, this is easy. Just fit the model for that group and calculate the $x_1$ for which $p_0$ is achieved.

For the two-variable model, however, I am stumped. Worse, I would like a confidence interval, and my naïve approach of considering the upper and lower confidence limits of the parameters$^{\dagger}$ is giving me confidence intervals that are too wide ($95\%$ confidence interval has $100\%$ coverage in $500,000$ Monte Carlo iterations).

How can I estimate the difference in $x_1$ for which $p_0$ is achieved and give a corresponding confidence interval?

$^{\dagger}$A reference my boss sent my way calculates the confidence interval by taking the lower confidence limits for each parameter, writing a model with those as the parameter estimates, finding $x_1$ that gives $p_0$, and then doing it all again for the upper confidence limits. I had a feeling that this method would be inadequate for a confidence interval of $x_1$ giving $p_0$.

  • $\begingroup$ Are you groups independent, or do the same participants/whatever appear in both groups? $\endgroup$
    – Eoin
    Jan 19, 2021 at 15:36
  • $\begingroup$ @Eoin They are totally independent; subjects can be control or treatment, but not both. $\endgroup$
    – Dave
    Jan 19, 2021 at 15:37

1 Answer 1


I think this can be solved using posterior samples (in a Bayesian model), or, similarly, using Monte Carlo samples from the maximum likelihood parameters and their covariance matrix. I'm 95% sure this is correct, but am always happy to be corrected.

First, confidence intervals for the one group case.

You already know how to calculate your quantity of interest (let's call it $x(p_0)$) given estimates of $\beta_0, \beta_1$. The arm package provides a simple function, arm::sim (source code) that uses the model covariance matrix to draw samples for each parameter such that the sample mean is the ML value, the 2.5th and 97.5th percentiles correspond to the 95% CI, and covariances between parameters are preserved. You can thus estimate the 95% CI for $x(p_0)$ by calculating $x(p_0)$ for each sampled pair of values for of the $\beta$s, and calculating the percentiles.

Now, because your groups are independent, you can do this separately for each group, yielding a collection of samples of $x(p_0)$ for group 1, and similar samples for group 2. Subtracting these sets of samples pairwise yields samples from your estimate of the difference, $\delta$, in $x(p_0)$ between the two groups. The maximum likelihood estimate of $\delta$ is the mean of these samples (assuming they remain Normally distributed -- something I'm not sure of), and the 95% CI can be estimated by taking the percentiles, as before.

  • $\begingroup$ Why split it into one regression per group? This seems like it should be able to be accomplished in one regression fit. $\endgroup$
    – Dave
    Jan 19, 2021 at 16:53
  • $\begingroup$ It probably can be, but that makes the equation for calculating $x(p_0)$ from the $\beta$ more complicated, and doesn't change the result if the groups are independent (since data from the control group doesn't affect parameter estimates from the treatment group, and vice verse). $\endgroup$
    – Eoin
    Jan 19, 2021 at 17:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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