```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
```

## 1. A Behavioral Story: Counting Emergency Room Visits

Dr. Elena Rodriguez, a public health researcher, works with an urban hospital to understand patterns in emergency room visits. She's particularly interested in how social and behavioral factors influence the number of ER visits patients make in a year.

After analyzing patient records, she notices that the data doesn't look like typical linear regression data. Some patients have zero ER visits, others have one or two, and a few have many more. The distribution is skewed, with most values clustered near zero. The residuals from a linear regression model looked terrible—predictions included impossible negative visit counts!

This is where Generalized Linear Models (GLMs) become essential. GLMs extend ordinary linear regression to handle different types of outcome variables: counts, binary outcomes, proportions, and other non-normal distributions. For count data like ER visits, we use **Poisson regression** or **negative binomial regression**, two common GLM variants.

## 2. What are Generalized Linear Models?

Generalized Linear Models (GLMs) extend ordinary linear regression to accommodate response variables that don't follow a normal distribution. While standard linear regression assumes:

$$
Y_i \sim N(\mu_i, \sigma^2)
$$

GLMs allow the response variable to follow various distributions from the exponential family, including:
- **Poisson distribution** (count data)
- **Binomial distribution** (binary or proportion data)
- **Gamma distribution** (positive continuous data)
- **Negative binomial distribution** (overdispersed count data)

**The Three Components of GLMs:**

1. **Random component**: Specifies the probability distribution of the response variable
2. **Systematic component**: Linear combination of predictors ($X\beta$)
3. **Link function**: Connects the mean of the response to the linear predictors

## 3. The GLM Framework in Plain Form

The general GLM can be written as:

$$
g(\mu_i) = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \cdots + \beta_p X_{pi}
$$

where:
- $\mu_i = E(Y_i)$ is the expected value of the response
- $g(\cdot)$ is the link function
- $\beta_0, \beta_1, \ldots, \beta_p$ are regression coefficients

**For Poisson Regression (Count Data):**

We assume:

$$
Y_i \sim \text{Poisson}(\lambda_i)
$$

with the log link function:

$$
\log(\lambda_i) = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \cdots
$$

This ensures $\lambda_i > 0$ (counts can't be negative). Exponentiating both sides:

$$
\lambda_i = \exp(\beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \cdots)
$$

**Interpreting Coefficients:**

For Poisson regression, $\exp(\beta_j)$ represents the multiplicative change in the expected count for a one-unit increase in $X_j$.

## 4. Simulating a Dataset

Let's create data similar to Dr. Rodriguez's study. We'll simulate ER visits based on age and chronic disease status.

```{r simulate-data}
# Set seed for reproducibility
set.seed(456)

# Sample size
n <- 200

# Simulate predictors
age <- round(rnorm(n, mean = 45, sd = 15))
age <- pmax(age, 18)  # Minimum age 18

chronic_disease <- rbinom(n, size = 1, prob = 0.3)  # 30% have chronic disease

# Generate Poisson counts
# Log-linear model: log(lambda) = beta0 + beta1*age + beta2*chronic_disease
beta0 <- 0.2
beta1 <- 0.01  # Age effect
beta2 <- 0.8   # Chronic disease effect

log_lambda <- beta0 + beta1 * age + beta2 * chronic_disease
lambda <- exp(log_lambda)

# Simulate ER visits from Poisson distribution
er_visits <- rpois(n, lambda = lambda)

# Create dataframe
er_data <- data.frame(
  age = age,
  chronic_disease = factor(chronic_disease, levels = c(0, 1), 
                           labels = c("No", "Yes")),
  er_visits = er_visits
)

# Display first few rows
head(er_data)

# Summary statistics
summary(er_data)

# Distribution of ER visits
table(er_data$er_visits)
```

## 5. Visualizing the Relationship

Let's visualize the distribution of ER visits and its relationship with our predictors.

```{r visualize-distribution}
library(ggplot2)

# Histogram of ER visits
ggplot(er_data, aes(x = er_visits)) +
  geom_histogram(binwidth = 1, fill = "steelblue", color = "black", alpha = 0.7) +
  labs(
    title = "Distribution of ER Visits",
    x = "Number of ER Visits",
    y = "Frequency"
  ) +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(color = "gray90", linetype = "dashed"),
    panel.grid.minor = element_blank(),
    axis.line.x = element_line(color = "black"),
    axis.line.y = element_line(color = "black"),
    panel.border = element_blank(),
    axis.line.x.top = element_blank(),
    axis.line.y.right = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold")
  )
```

The distribution is right-skewed with many low values and few high values—typical of count data.

```{r visualize-relationships}
# ER visits by chronic disease
ggplot(er_data, aes(x = chronic_disease, y = er_visits, fill = chronic_disease)) +
  geom_boxplot(alpha = 0.6) +
  geom_jitter(width = 0.2, alpha = 0.3) +
  labs(
    title = "ER Visits by Chronic Disease Status",
    x = "Chronic Disease",
    y = "Number of ER Visits"
  ) +
  scale_fill_brewer(palette = "Set2") +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(color = "gray90", linetype = "dashed"),
    panel.grid.minor = element_blank(),
    axis.line.x = element_line(color = "black"),
    axis.line.y = element_line(color = "black"),
    panel.border = element_blank(),
    axis.line.x.top = element_blank(),
    axis.line.y.right = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "none"
  )

# ER visits vs age
ggplot(er_data, aes(x = age, y = er_visits, color = chronic_disease)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "glm", method.args = list(family = "poisson"), se = FALSE) +
  labs(
    title = "ER Visits by Age and Chronic Disease Status",
    x = "Age (years)",
    y = "Number of ER Visits",
    color = "Chronic Disease"
  ) +
  scale_color_brewer(palette = "Set1") +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(color = "gray90", linetype = "dashed"),
    panel.grid.minor = element_blank(),
    axis.line.x = element_line(color = "black"),
    axis.line.y = element_line(color = "black"),
    panel.border = element_blank(),
    axis.line.x.top = element_blank(),
    axis.line.y.right = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold")
  )
```

## 6. Fitting the Poisson GLM in R

We use the `glm()` function with `family = poisson` for Poisson regression:

```{r fit-poisson}
# Fit Poisson GLM
poisson_model <- glm(er_visits ~ age + chronic_disease, 
                     family = poisson(link = "log"), 
                     data = er_data)

# Display summary
summary(poisson_model)
```

## 7. Interpreting the Results

Let's extract and interpret the key values:

```{r extract-poisson-results}
# Extract coefficients
coef_summary <- summary(poisson_model)$coefficients
print(coef_summary)

# Extract specific coefficients
beta0_hat <- coef(poisson_model)[1]
beta_age <- coef(poisson_model)[2]
beta_chronic <- coef(poisson_model)[3]

cat("\nIntercept:", round(beta0_hat, 4), "\n")
cat("Age coefficient:", round(beta_age, 4), "\n")
cat("Chronic disease coefficient:", round(beta_chronic, 4), "\n")

# Calculate exponentiated coefficients (rate ratios)
exp_coefs <- exp(coef(poisson_model))
cat("\nExponentiated coefficients (Rate Ratios):\n")
print(round(exp_coefs, 4))
```

**Interpretation on the Log Scale:**

- **Intercept** ($\beta_0$) = `r round(beta0_hat, 3)`: The log expected count of ER visits for a 0-year-old without chronic disease (not meaningful; we need to center age for practical interpretation)

- **Age coefficient** ($\beta_1$) = `r round(beta_age, 4)`: For each additional year of age, the log expected count of ER visits increases by `r round(beta_age, 4)`

- **Chronic disease coefficient** ($\beta_2$) = `r round(beta_chronic, 3)`: Having a chronic disease increases the log expected count of ER visits by `r round(beta_chronic, 3)`

**Interpretation on the Count Scale (Rate Ratios):**

Exponentiating the coefficients gives us rate ratios (also called incidence rate ratios):

- **Age**: $\exp(\beta_1)$ = `r round(exp_coefs[2], 4)` means for each additional year of age, the expected number of ER visits is multiplied by `r round(exp_coefs[2], 4)`, or increases by approximately `r round((exp_coefs[2] - 1) * 100, 2)`%

- **Chronic disease**: $\exp(\beta_2)$ = `r round(exp_coefs[3], 3)` means having a chronic disease multiplies the expected ER visits by `r round(exp_coefs[3], 3)`, or increases them by approximately `r round((exp_coefs[3] - 1) * 100, 1)`%

## 8. Making Predictions

Let's predict ER visits for specific patient profiles:

```{r predictions}
# Create new data for prediction
new_patients <- data.frame(
  age = c(30, 30, 60, 60),
  chronic_disease = factor(c("No", "Yes", "No", "Yes"), levels = c("No", "Yes"))
)

# Predict expected counts
predictions <- predict(poisson_model, newdata = new_patients, type = "response")
new_patients$predicted_visits <- predictions

print(new_patients)
```

For example:
- A 30-year-old without chronic disease: `r round(predictions[1], 2)` expected ER visits
- A 30-year-old with chronic disease: `r round(predictions[2], 2)` expected ER visits
- A 60-year-old without chronic disease: `r round(predictions[3], 2)` expected ER visits
- A 60-year-old with chronic disease: `r round(predictions[4], 2)` expected ER visits

## 9. Checking Model Fit: Deviance and Overdispersion

**Deviance:**

The deviance measures how well the model fits the data. For Poisson regression, we check the residual deviance against the residual degrees of freedom.

```{r check-deviance}
# Extract deviance information
residual_deviance <- poisson_model$deviance
residual_df <- poisson_model$df.residual
deviance_ratio <- residual_deviance / residual_df

cat("Residual Deviance:", round(residual_deviance, 2), "\n")
cat("Residual DF:", residual_df, "\n")
cat("Deviance / DF:", round(deviance_ratio, 3), "\n")
```

If the deviance/DF ratio is close to 1, the Poisson model fits well. A ratio substantially greater than 1 indicates **overdispersion**—more variance than the Poisson distribution allows.

**Testing for Overdispersion:**

```{r test-overdispersion}
# Dispersion parameter estimate
dispersion <- sum(residuals(poisson_model, type = "pearson")^2) / residual_df
cat("Dispersion parameter:", round(dispersion, 3), "\n")

if (dispersion > 1.5) {
  cat("Warning: Overdispersion detected. Consider negative binomial regression.\n")
} else {
  cat("No severe overdispersion detected.\n")
}
```

## 10. Negative Binomial Regression (for Overdispersion)

If overdispersion is detected, we can fit a negative binomial model, which includes an extra parameter to account for excess variance:

```{r fit-negbin}
library(MASS)

# Fit negative binomial GLM
nb_model <- glm.nb(er_visits ~ age + chronic_disease, data = er_data)

# Display summary
summary(nb_model)

# Compare with Poisson
cat("\nPoisson AIC:", AIC(poisson_model), "\n")
cat("Negative Binomial AIC:", AIC(nb_model), "\n")
```

The negative binomial model includes a dispersion parameter (theta). A lower AIC indicates better model fit. If the negative binomial AIC is substantially lower, use that model instead.

## 11. Visualizing the Fitted Model

```{r visualize-fitted}
# Create prediction grid
age_seq <- seq(min(er_data$age), max(er_data$age), length.out = 100)
pred_data <- expand.grid(
  age = age_seq,
  chronic_disease = factor(c("No", "Yes"), levels = c("No", "Yes"))
)

# Get predictions with confidence intervals
pred_data$predicted <- predict(poisson_model, newdata = pred_data, type = "response")

# Plot
ggplot(er_data, aes(x = age, y = er_visits, color = chronic_disease)) +
  geom_point(alpha = 0.4) +
  geom_line(data = pred_data, aes(y = predicted), size = 1.2) +
  labs(
    title = "Fitted Poisson Regression Model",
    x = "Age (years)",
    y = "Expected ER Visits",
    color = "Chronic Disease"
  ) +
  scale_color_brewer(palette = "Set1") +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(color = "gray90", linetype = "dashed"),
    panel.grid.minor = element_blank(),
    axis.line.x = element_line(color = "black"),
    axis.line.y = element_line(color = "black"),
    panel.border = element_blank(),
    axis.line.x.top = element_blank(),
    axis.line.y.right = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold")
  )
```

## 12. Model Diagnostics

```{r diagnostics}
# Residuals vs fitted
plot(poisson_model, which = 1, main = "Residuals vs Fitted")

# Q-Q plot
plot(poisson_model, which = 2, main = "Q-Q Plot")

# Influential observations
plot(poisson_model, which = 5, main = "Residuals vs Leverage")
```

For GLMs, we examine:
- **Residuals vs Fitted**: Look for patterns suggesting poor fit
- **Q-Q Plot**: Assess normality of deviance residuals
- **Influential Points**: Check for observations with high leverage or Cook's distance

## 13. Reporting Results in APA Style

Here's how to report Poisson regression results:

> A Poisson regression was conducted to examine the relationship between age, chronic disease status, and the number of emergency room visits per year. The model provided a good fit to the data, with a residual deviance of `r round(residual_deviance, 2)` on `r residual_df` degrees of freedom (dispersion parameter = `r round(dispersion, 2)`).
>
> Age was a significant predictor of ER visits, *b* = `r round(beta_age, 4)`, *z* = `r round(coef_summary[2, 3], 2)`, *p* `r ifelse(coef_summary[2, 4] < .001, "< .001", paste("=", round(coef_summary[2, 4], 3)))`. Specifically, each additional year of age was associated with a `r round((exp_coefs[2] - 1) * 100, 2)`% increase in the expected number of ER visits (rate ratio = `r round(exp_coefs[2], 3)`).
>
> Chronic disease status was also a significant predictor, *b* = `r round(beta_chronic, 3)`, *z* = `r round(coef_summary[3, 3], 2)`, *p* `r ifelse(coef_summary[3, 4] < .001, "< .001", paste("=", round(coef_summary[3, 4], 3)))`. Patients with chronic diseases had `r round((exp_coefs[3] - 1) * 100, 1)`% more ER visits compared to those without chronic diseases (rate ratio = `r round(exp_coefs[3], 2)`), holding age constant.

## 14. Practical Considerations

**When to Use GLMs:**
- When the outcome variable is not normally distributed
- For count data (Poisson or negative binomial)
- For binary outcomes (logistic regression)
- For rates or proportions

**Model Selection:**
- **Poisson**: When variance equals the mean
- **Negative Binomial**: When variance exceeds the mean (overdispersion)
- Compare models using AIC or likelihood ratio tests

**Common Mistakes:**
- Using linear regression for count data (can predict negative counts)
- Ignoring overdispersion in Poisson models
- Misinterpreting coefficients (they're on the log scale)
- Not checking model assumptions and diagnostics

**Extensions:**
- **Zero-inflated models**: When there are excess zeros
- **Hurdle models**: Separate processes for zero vs. non-zero counts
- **Offset terms**: When exposure varies (e.g., person-years)

## 15. Conclusion

Generalized Linear Models extend regression to handle non-normal response variables. For count data like ER visits, Poisson regression provides a natural framework that respects the discrete, non-negative nature of counts.

In our analysis, we found that both age and chronic disease significantly predict ER visit frequency. Older patients and those with chronic diseases make more ER visits, with chronic disease having a particularly strong effect (more than doubling expected visits).

GLMs are powerful tools in public health and psychology research, enabling analysis of diverse outcome types while maintaining the interpretability of regression coefficients. Mastering GLMs opens doors to analyzing binary outcomes (logistic regression), survival times (Cox regression), and many other data types common in health sciences.
