Return to Answer

The delta method states

$$ \operatorname{Var}(g(X)) = [g'(X)]^2 \operatorname{Var}(X)$$

Because this problem involves two parameters, we can extend this to the multivariate delta rule

$$ =\nabla g(X)^T \, \Sigma \, \nabla g $$

Here,

$$ g = \left[e^{\beta_{0}+\beta_{1}}-e^{\beta_{0}}\right] /\left[\left(1+e^{\beta_{0}+\beta_{1}}\right)\left(1+e^{\beta_{0}}\right)\right] $$

and $\Sigma$ is the variance covariance matrix from your model. $\nabla g$ is...gross. I'm not going to do that by hand, and computer algebra while fast yields a mess of symbols. You can however use autodifferentiation compute the gradient. Once you calculate the variance, then its simply your estimate of the difference in probs plus/minus 1.96 times the standard deviation (root of the variance). Caution, this approach will yield answers below 0 or above 1.

We can do this in R in the following way (note you need to install the autodiffr package).

library(autodiffr)

g = function(b)  (exp(b[1] + b[2]) - exp(b[1])) / ((1+ exp(b[1] + b[2]))*(1+exp(b[1])))

x = rbinom(100, 1, 0.5)
eta = -0.8 + 0.2*x
p = plogis(eta)
y = rbinom(100, 1, p)

model = glm(y~x, family=binomial())
Bigma = vcov(model)

grad_g = makeGradFunc(g)
nabla_g = grad_g(coef(model))


se = as.numeric(sqrt(nabla_g %*% Bigma %*% nabla_g))


estimate = diff(predict(model, newdata=list(x=c(0, 1)), type='response'))

estimate + c(-1, 1)*1.96*se

Repeating this procedure for this modest example shows that the resulting confidence interval has near nominal coverage, which is a good thing, but I imagine things would become worse as the probabilities approach 0 or 1.

The delta method states

$$ \operatorname{Var}(g(X)) = [g'(X)]^2 \operatorname{Var}(X)$$

Because this problem involves two parameters, we can extend this to the multivariate delta method

$$ =\nabla g^T \, \Sigma \, \nabla g $$

Here,

$$ g = \left[e^{\beta_{0}+\beta_{1}}-e^{\beta_{0}}\right] /\left[\left(1+e^{\beta_{0}+\beta_{1}}\right)\left(1+e^{\beta_{0}}\right)\right] $$

and $\Sigma$ is the variance covariance matrix from your model. $\nabla g$ is...gross. I'm not going to do that by hand, and computer algebra while fast yields a mess of symbols. You can however use autodifferentiation compute the gradient. Once you calculate the variance, then its simply your estimate of the difference in probs plus/minus 1.96 times the standard deviation (root of the variance). Caution, this approach will yield answers below 0 or above 1.

We can do this in R in the following way (note you need to install the autodiffr package).

library(autodiffr)

g = function(b)  (exp(b[1] + b[2]) - exp(b[1])) / ((1+ exp(b[1] + b[2]))*(1+exp(b[1])))

x = rbinom(100, 1, 0.5)
eta = -0.8 + 0.2*x
p = plogis(eta)
y = rbinom(100, 1, p)

model = glm(y~x, family=binomial())
Bigma = vcov(model)

grad_g = makeGradFunc(g)
nabla_g = grad_g(coef(model))


se = as.numeric(sqrt(nabla_g %*% Bigma %*% nabla_g))


estimate = diff(predict(model, newdata=list(x=c(0, 1)), type='response'))

estimate + c(-1, 1)*1.96*se

Repeating this procedure for this modest example shows that the resulting confidence interval has near nominal coverage, which is a good thing, but I imagine things would become worse as the probabilities approach 0 or 1.

The delta method states

$$ \operatorname{Var}(g(X)) = [g'(X)]^2 \operatorname{Var}(X)$$

Because this problem involves two parameters, we can extend this to the multivariate delta rule

$$ =\nabla g(X)^T \, \Sigma \, \nabla g $$

Here,

$$ g = \left[e^{\beta_{0}+\beta_{1}}-e^{\beta_{0}}\right] /\left[\left(1+e^{\beta_{0}+\beta_{1}}\right)\left(1+e^{\beta_{0}}\right)\right] $$

and $\Sigma$ is the variance covariance matrix from your model. $\nabla g$ is...gross. I'm not going to do that by hand, and computer algebra while fast yields a mess of symbols. You can however use autodifferentiation compute the gradient. Once you calculate the variance, then its simply your estimate of the difference in probs plus/minus 1.96 times the standard deviation (root of the variance). Caution, this approach will yield answers below 0 or above 1.

We can do this in R in the following way (note you need to install the autodiffr package).

library(autodiffr)

g = function(b)  (exp(b[1] + b[2]) - exp(b[1])) / ((1+ exp(b[1] + b[2]))*(1+exp(b[1])))

x = rbinom(100, 1, 0.5)
eta = -0.8 + 0.2*x
p = plogis(eta)
y = rbinom(100, 1, p)

model = glm(y~x, family=binomial())
Bigma = vcov(model)

grad_g = makeGradFunc(g)
nabla_g = grad_g(coef(model))


se = as.numeric(sqrt(nabla_g %*% Bigma %*% nabla_g))


estimate = diff(predict(model, newdata=list(x=c(0, 1)), type='response'))

estimate + c(-1, 1)*1.96*se

```

The delta method states

$$ \operatorname{Var}(g(X)) = [g'(X)]^2 \operatorname{Var}(X)$$

Because this problem involves two parameters, we can extend this to the multivariate delta rule

$$ =\nabla g(X)^T \, \Sigma \, \nabla g $$

Here,

$$ g = \left[e^{\beta_{0}+\beta_{1}}-e^{\beta_{0}}\right] /\left[\left(1+e^{\beta_{0}+\beta_{1}}\right)\left(1+e^{\beta_{0}}\right)\right] $$

and $\Sigma$ is the variance covariance matrix from your model. $\nabla g$ is...gross. I'm not going to do that by hand, and computer algebra while fast yields a mess of symbols. You can however use autodifferentiation compute the gradient. Once you calculate the variance, then its simply your estimate of the difference in probs plus/minus 1.96 times the standard deviation (root of the variance). Caution, this approach will yield answers below 0 or above 1.

We can do this in R in the following way (note you need to install the autodiffr package).

library(autodiffr)

g = function(b)  (exp(b[1] + b[2]) - exp(b[1])) / ((1+ exp(b[1] + b[2]))*(1+exp(b[1])))

x = rbinom(100, 1, 0.5)
eta = -0.8 + 0.2*x
p = plogis(eta)
y = rbinom(100, 1, p)

model = glm(y~x, family=binomial())
Bigma = vcov(model)

grad_g = makeGradFunc(g)
nabla_g = grad_g(coef(model))


se = as.numeric(sqrt(nabla_g %*% Bigma %*% nabla_g))


estimate = diff(predict(model, newdata=list(x=c(0, 1)), type='response'))

estimate + c(-1, 1)*1.96*se

Repeating this procedure for this modest example shows that the resulting confidence interval has near nominal coverage, which is a good thing, but I imagine things would become worse as the probabilities approach 0 or 1.

The delta method states

$$ \operatorname{Var}(g(X)) = [g'(X)]^2 \operatorname{Var}(X)$$

Because this problem involves two parameters, we can extend this to the multivariate delta rule

$$ =\nabla g(X)^T \, \Sigma \, \nabla g $$

Here,

$$ g = \left[e^{\beta_{0}+\beta_{1}}-e^{\beta_{0}}\right] /\left[\left(1+e^{\beta_{0}+\beta_{1}}\right)\left(1+e^{\beta_{0}}\right)\right] $$

and $\Sigma$ is the variance covariance matrix from your model. $\nabla g$ is...gross. I'm not going to do that by hand, and computer algebra while fast yields a mess of symbols. You can however use autodifferentiation compute the gradient. Once you calculate the variance, then its simply your estimate of the difference in probs plus/minus 1.96 times the standard deviation (root of the variance). Caution, this approach will yield answers below 0 or above 1.

The delta method states

$$ \operatorname{Var}(g(X)) = [g'(X)]^2 \operatorname{Var}(X)$$

Because this problem involves two parameters, we can extend this to the multivariate delta rule

$$ =\nabla g(X)^T \, \Sigma \, \nabla g $$

Here,

$$ g = \left[e^{\beta_{0}+\beta_{1}}-e^{\beta_{0}}\right] /\left[\left(1+e^{\beta_{0}+\beta_{1}}\right)\left(1+e^{\beta_{0}}\right)\right] $$

and $\Sigma$ is the variance covariance matrix from your model. $\nabla g$ is...gross. I'm not going to do that by hand, and computer algebra while fast yields a mess of symbols. You can however use autodifferentiation compute the gradient. Once you calculate the variance, then its simply your estimate of the difference in probs plus/minus 1.96 times the standard deviation (root of the variance). Caution, this approach will yield answers below 0 or above 1.

We can do this in R in the following way (note you need to install the autodiffr package).

library(autodiffr)

g = function(b)  (exp(b[1] + b[2]) - exp(b[1])) / ((1+ exp(b[1] + b[2]))*(1+exp(b[1])))

x = rbinom(100, 1, 0.5)
eta = -0.8 + 0.2*x
p = plogis(eta)
y = rbinom(100, 1, p)

model = glm(y~x, family=binomial())
Bigma = vcov(model)

grad_g = makeGradFunc(g)
nabla_g = grad_g(coef(model))


se = as.numeric(sqrt(nabla_g %*% Bigma %*% nabla_g))


estimate = diff(predict(model, newdata=list(x=c(0, 1)), type='response'))

estimate + c(-1, 1)*1.96*se

```