Get 95% confidence interval with glm(..) in R

Calculating 95% Confidence Intervals for GLM in R

Learn how to obtain and interpret 95% confidence intervals for coefficients from Generalized Linear Models (GLM) in R, covering common distributions and methods.

Generalized Linear Models (GLMs) are a powerful and flexible extension of ordinary least squares regression, allowing for response variables that have error distribution models other than a normal distribution. When working with GLMs in R, it's crucial to not only estimate the model coefficients but also to understand their uncertainty. This is typically done by calculating confidence intervals. This article will guide you through the process of obtaining 95% confidence intervals for GLM coefficients in R, discussing different approaches and considerations for various link functions and distributions.

Understanding GLM Coefficients and Their Standard Errors

In a GLM, the coefficients (often denoted as β) represent the change in the linear predictor for a one-unit change in the corresponding predictor variable, holding other variables constant. The relationship between the linear predictor and the mean of the response variable is defined by the link function. For example, in a logistic regression (a type of GLM with a binomial family and logit link), the coefficients are on the log-odds scale. To interpret these coefficients meaningfully, especially for non-normal distributions, confidence intervals are essential.

The glm() function in R estimates these coefficients and their standard errors. The standard errors are crucial for constructing confidence intervals, as they quantify the precision of the coefficient estimates. A common method for calculating confidence intervals relies on the assumption that the coefficient estimates are approximately normally distributed, especially with large sample sizes. This allows us to use the Wald method.

flowchart TD
    A[GLM Model Fit in R] --> B{Extract Coefficients & Std. Errors}
    B --> C{Choose Confidence Interval Method}
    C --> D1[Wald Method (Normal Approximation)]
    C --> D2[Profile Likelihood Method]
    C --> D3[Bootstrapping]
    D1 --> E1[Calculate CI: Estimate +/- Z*Std.Error]
    D2 --> E2[Iteratively find bounds where log-likelihood drops by 1.92]
    D3 --> E3[Resample data, refit GLM, get empirical distribution of coefficients]
    E1 --> F[Interpret CI on Link Scale]
    E2 --> F
    E3 --> F
    F --> G{If needed: Transform CI to Response Scale}
    G --> H[Final Interpretation]

Workflow for obtaining confidence intervals from GLM in R.

Method 1: Wald Confidence Intervals (Normal Approximation)

The most straightforward way to get confidence intervals for GLM coefficients is using the Wald method, which assumes that the sampling distribution of the maximum likelihood estimates (MLEs) of the coefficients is approximately normal. This method is implemented by default in the confint() function for glm objects when method = "Wald" (though it's often the default and doesn't need to be explicitly stated).

For a 95% confidence interval, we typically use a Z-score of approximately 1.96 (for a two-tailed test). The interval is calculated as: Estimate ± Z * Standard Error.

Let's demonstrate with a logistic regression example.

# Load necessary library
library(stats)

# Simulate some data for a logistic regression
set.seed(123)
n <- 100
x <- rnorm(n)
y <- rbinom(n, 1, prob = plogis(0.5 + 1.5*x - 0.8*x^2))

# Fit a logistic regression model
model_logit <- glm(y ~ x + I(x^2), family = binomial(link = "logit"))

# View model summary
summary(model_logit)

# Get Wald confidence intervals (default for glm objects)
confint(model_logit)

# You can also manually calculate them for understanding
coefs <- coef(model_logit)
std_errs <- summary(model_logit)$coefficients[, "Std. Error"]
z_val <- qnorm(0.975) # For 95% CI

lower_ci <- coefs - z_val * std_errs
upper_ci <- coefs + z_val * std_errs

wald_ci_manual <- cbind(lower_ci, upper_ci)
colnames(wald_ci_manual) <- c("2.5 %", "97.5 %")
print(wald_ci_manual)

Fitting a logistic regression and obtaining Wald confidence intervals.

Method 2: Profile Likelihood Confidence Intervals

Profile likelihood confidence intervals are generally considered more reliable than Wald intervals, especially for GLMs and smaller sample sizes, because they do not rely on the assumption of normality of the coefficient estimates. Instead, they are derived by finding the range of parameter values for which the log-likelihood does not drop by more than a certain amount (e.g., 1.92 for a 95% CI, corresponding to a chi-squared distribution with 1 degree of freedom). This method directly explores the shape of the likelihood function.

The confint() function in R can also compute profile likelihood intervals by specifying method = "profile".

# Using the same model_logit from before

# Get profile likelihood confidence intervals
confint(model_logit, method = "profile")

# Compare with Wald intervals
wald_ci <- confint(model_logit, method = "Wald")
profile_ci <- confint(model_logit, method = "profile")

print("Wald CIs:")
print(wald_ci)
print("\nProfile Likelihood CIs:")
print(profile_ci)

Calculating profile likelihood confidence intervals for a GLM.

Transforming Confidence Intervals to the Response Scale

For many GLMs (e.g., logistic, Poisson), the coefficients and their confidence intervals are on the scale of the linear predictor (e.g., log-odds, log-counts). To make them more interpretable, especially for presentations, you might want to transform them back to the response scale (e.g., probabilities, counts). This involves applying the inverse of the link function to both the coefficient estimate and its confidence interval bounds.

For logistic regression, this means transforming log-odds to probabilities using the plogis() function (inverse of logit). For Poisson regression, it means transforming log-rates to rates using exp() (inverse of log).

Important: When transforming, you should transform the interval bounds directly, not the standard error. The interval on the transformed scale will generally not be symmetric around the transformed point estimate.

# Using the profile likelihood CIs for model_logit
profile_ci <- confint(model_logit, method = "profile")

# Transform the intercept and x coefficient CIs to probability scale
# (Note: This is for illustration. Interpreting individual coefficients on the response scale
# for non-linear models can be complex and often requires marginal effects or predictions.)

# For the intercept (representing probability when x=0)
intercept_ci_prob <- plogis(profile_ci["(Intercept)", ])
print(paste("Intercept CI on probability scale:", round(intercept_ci_prob[1], 3), "-", round(intercept_ci_prob[2], 3)))

# For the 'x' coefficient, transforming it directly to probability is not straightforward
# as it represents change in log-odds. Instead, we often look at odds ratios.

# Calculate Odds Ratios and their CIs for logistic regression
exp_coefs <- exp(coef(model_logit))
exp_profile_ci <- exp(profile_ci)

print("\nOdds Ratios:")
print(exp_coefs)
print("\nOdds Ratio CIs (from profile likelihood CIs):")
print(exp_profile_ci)

Transforming confidence intervals from log-odds to probability/odds ratio scale.