Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add confidence intervals to model outcomes #1

Open
tavareshugo opened this issue Mar 27, 2024 · 1 comment
Open

add confidence intervals to model outcomes #1

tavareshugo opened this issue Mar 27, 2024 · 1 comment

Comments

@tavareshugo
Copy link

tavareshugo commented Mar 27, 2024

currently, there is code to plot the fitted line (using broom::augment()), but it would be nice to also plot the confidence interval of the estimate.
I'm not sure broom::augment() can do this, but here is some "manual" code for the logistic case (edit: ignore the code below, see my next comment):

library(ggplot2)

#### simulate data ####

# sample size
n_samples <- 30

# x values uniformly sampled from -1 to 1
sim <- data.frame(x = runif(n_samples, -1, 1))

# successes sampled from binomial
# model is:
# Y ~ Binom(n, p)
# logit(p) = beta0 + beta1*x
# number of trials 'n' is simulated as 10
# using beta0 = 0 and beta1 = 2
# the probability 'p' is calculated using the inverse logit (plogis)
sim$success <- rbinom(n_samples, 
                      size = 10, 
                      prob = plogis(2*sim$x))
sim$fail <- 10 - sim$success


#### fit the model ####

fit <- glm(cbind(success, fail) ~ x, data = sim, family = "binomial")

# the confidence interval of the estimate
confint(fit)

# create a table of predicted values - requires knowing about the link function (and its inverse)
fit_pred <- data.frame(x = seq(min(sim$x), max(sim$x), length.out = 500)) |> 
  transform(pred = plogis(fit$coefficients[1] + fit$coefficients[2]*x),
            lo = plogis(confint(fit)[1,1] + confint(fit)[2,1]*x),
            hi = plogis(confint(fit)[1,2] + confint(fit)[2,2]*x))

# visualise
sim |> 
  ggplot(aes(x)) + 
  geom_ribbon(data = fit_pred, 
              aes(ymin = lo, ymax = hi), alpha = 0.2) +
  geom_line(data = fit_pred, aes(y = pred)) +
  geom_point(aes(y = success/(success+fail)))

A similar thing should work for the Poisson model, using exp() as the inverse of the link function.

@tavareshugo
Copy link
Author

I guess that manual approach will be harder if the model is more complex (multiple predictors, interactions, etc.).

Maybe insight::get_predicted_ci() (here) is an easier alternative (haven't tested it).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant