Beta Regression

# Beta Regression Example 1: Exam Pass Rates
# Proportion outcome: class pass rates

library(betareg)     # Beta regression
library(ggplot2)     # Visualization
library(dplyr)       # Data manipulation
library(lmtest)      # Likelihood ratio tests

set.seed(2025)
n <- 90

# Simulate classroom data
data <- data.frame(
  class_size = sample(10:40, n, replace=TRUE),
  teacher_exp = sample(0:30, n, replace=TRUE)
)

# True model: logit(pass_rate) = 0.5 - 0.04*class_size + 0.03*teacher_exp
# Smaller classes → higher pass rates
# More experienced teachers → higher pass rates
logit_mu <- 0.5 - 0.04*data$class_size + 0.03*data$teacher_exp
mu <- plogis(logit_mu)  # Inverse logit: (0,1)

# Generate beta-distributed pass rates with precision φ=30
phi <- 30  # Higher φ = less variability
alpha <- mu * phi
beta_param <- (1 - mu) * phi
data$pass_rate <- rbeta(n, alpha, beta_param)

# Ensure no exact 0s or 1s (requirement for beta regression)
data$pass_rate <- pmax(0.001, pmin(0.999, data$pass_rate))

cat("=== Data Summary ===")
summary(data)
cat("\nPass rate range:", range(data$pass_rate), "\n")
cat("Any exact 0s or 1s?", any(data$pass_rate == 0 | data$pass_rate == 1), "\n")

# === STEP 1: Visualize Data ===
p1 <- ggplot(data, aes(x=class_size, y=pass_rate)) +
  geom_point(alpha=0.6, size=2) +
  geom_smooth(method="loess", se=TRUE, color="blue") +
  labs(title="Pass Rate vs. Class Size",
       x="Class Size", y="Pass Rate(Proportion)") +
  ylim(0, 1) +
  theme_classic()

p2 <- ggplot(data, aes(x=teacher_exp, y=pass_rate)) +
  geom_point(alpha=0.6, size=2) +
  geom_smooth(method="loess", se=TRUE, color="darkgreen") +
  labs(title="Pass Rate vs. Teacher Experience",
       x="Teacher Experience(years)", y="Pass Rate(Proportion)") +
  ylim(0, 1) +
  theme_classic()

library(gridExtra)
grid.arrange(p1, p2, ncol=2)

# Histogram of pass rates
ggplot(data, aes(x=pass_rate)) +
  geom_histogram(bins=20, fill="steelblue", color="black", alpha=0.7) +
  labs(title="Distribution of Pass Rates",
       x="Pass Rate", y="Frequency") +
  theme_classic()

# === STEP 2: WRONG APPROACH - Linear Regression (OLS) ===
# Common mistake: treating proportions as normal continuous outcome
model_ols <- lm(pass_rate ~ class_size + teacher_exp, data=data)
summary(model_ols)

cat("\n=== OLS Problems ===")
# Check for predictions outside [0,1]
pred_ols <- predict(model_ols)
cat("OLS predictions < 0:", sum(pred_ols < 0), "\n")
cat("OLS predictions > 1:", sum(pred_ols > 1), "\n")
cat("Range of OLS predictions:", round(range(pred_ols), 3), "\n")

# Check heteroscedasticity (common with proportions)
library(lmtest)
bptest(model_ols)
cat("Breusch-Pagan p-value suggests heteroscedasticity\n")

# Residual plot shows heteroscedasticity
par(mfrow=c(2,2))
plot(model_ols)
par(mfrow=c(1,1))

# === STEP 3: CORRECT APPROACH - Beta Regression ===
model_beta <- betareg(pass_rate ~ class_size + teacher_exp, data=data,
                      link="logit")  # Default: logit link
summary(model_beta)

cat("\n=== Beta Regression Results ===")
print(summary(model_beta))

# Coefficients interpretation (on logit scale)
coef_beta <- coef(model_beta)
cat("\nCoefficients(logit scale):")
print(coef_beta)

# === STEP 4: Coefficient Interpretation ===
# On logit scale: β represents change in log-odds
# For substantive interpretation: marginal effects at means

# Average marginal effects (AME)
library(margins)
margins_beta <- margins(model_beta)
summary(margins_beta)

cat("\n=== Average Marginal Effects ===")
print(summary(margins_beta))

# Manual calculation for interpretation
# At mean values:
mean_class <- mean(data$class_size)
mean_exp <- mean(data$teacher_exp)

# Predict at mean
pred_mean <- predict(model_beta, 
                     newdata=data.frame(class_size=mean_class, 
                                       teacher_exp=mean_exp),
                     type="response")

# Predict with class_size + 1
pred_class_plus1 <- predict(model_beta,
                           newdata=data.frame(class_size=mean_class+1,
                                             teacher_exp=mean_exp),
                           type="response")

cat("\n=== Substantive Interpretation ===")
cat("At mean class size(", round(mean_class, 1), ") and mean experience(",
    round(mean_exp, 1), " years):\n")
cat("Predicted pass rate:", round(pred_mean, 3), "\n")
cat("\nEffect of 1-student increase in class size:\n")
cat("New predicted pass rate:", round(pred_class_plus1, 3), "\n")
cat("Absolute change:", round(pred_class_plus1 - pred_mean, 4), 
    "(~", round((pred_class_plus1 - pred_mean)*100, 2), "percentage points)\n")

# === STEP 5: Model Diagnostics ===

# Residual plots
par(mfrow=c(2,2))

# Pearson residuals vs. fitted
resid_pearson <- residuals(model_beta, type="pearson")
fitted_beta <- fitted(model_beta)

plot(fitted_beta, resid_pearson, 
     main="Pearson Residuals vs. Fitted",
     xlab="Fitted Values", ylab="Pearson Residuals")
abline(h=0, col="red", lty=2)

# Q-Q plot
qqnorm(resid_pearson, main="Q-Q Plot(Pearson Residuals)")
qqline(resid_pearson, col="red")

# Residuals vs. class_size
plot(data$class_size, resid_pearson,
     main="Residuals vs. Class Size",
     xlab="Class Size", ylab="Pearson Residuals")
abline(h=0, col="red", lty=2)

# Residuals vs. teacher_exp
plot(data$teacher_exp, resid_pearson,
     main="Residuals vs. Teacher Experience",
     xlab="Teacher Experience", ylab="Pearson Residuals")
abline(h=0, col="red", lty=2)

par(mfrow=c(1,1))

# === STEP 6: Pseudo R-squared ===
cat("\n=== Model Fit ===")
cat("Pseudo R² (McFadden):", round(model_beta$pseudo.r.squared, 3), "\n")
cat("Log-likelihood:", round(logLik(model_beta), 2), "\n")
cat("AIC:", round(AIC(model_beta), 2), "\n")
cat("BIC:", round(BIC(model_beta), 2), "\n")

# Compare to null model
model_null <- betareg(pass_rate ~ 1, data=data)
lrtest(model_null, model_beta)

# === STEP 7: Predictions and Visualization ===

# Create prediction grid
class_range <- seq(10, 40, length.out=50)
exp_range <- seq(0, 30, length.out=50)
pred_grid <- expand.grid(class_size=class_range, teacher_exp=mean(data$teacher_exp))
pred_grid$predicted <- predict(model_beta, newdata=pred_grid, type="response")

# Plot predictions
ggplot(pred_grid, aes(x=class_size, y=predicted)) +
  geom_line(color="blue", size=1.2) +
  geom_point(data=data, aes(x=class_size, y=pass_rate), alpha=0.4) +
  labs(title="Beta Regression: Predicted Pass Rate vs. Class Size",
       subtitle=paste("Holding teacher experience at mean(", round(mean(data$teacher_exp), 1), " years)"),
       x="Class Size", y="Pass Rate") +
  ylim(0, 1) +
  theme_classic()

# Predicted vs. Observed
data$predicted_beta <- fitted(model_beta)

ggplot(data, aes(x=predicted_beta, y=pass_rate)) +
  geom_point(alpha=0.6, size=2) +
  geom_abline(slope=1, intercept=0, color="red", linetype="dashed") +
  labs(title="Predicted vs. Observed Pass Rates",
       x="Predicted Pass Rate(Beta Regression)",
       y="Observed Pass Rate") +
  xlim(0, 1) + ylim(0, 1) +
  theme_classic()

# === STEP 8: Compare OLS vs. Beta Predictions ===
data$predicted_ols <- fitted(model_ols)

ggplot(data, aes(x=class_size)) +
  geom_point(aes(y=pass_rate), alpha=0.5, size=2) +
  geom_line(aes(y=predicted_beta, color="Beta Regression"), size=1.2) +
  geom_line(aes(y=predicted_ols, color="OLS"), size=1.2, linetype="dashed") +
  scale_color_manual(values=c("Beta Regression"="blue", "OLS"="red")) +
  labs(title="Beta Regression vs. OLS Predictions",
       x="Class Size", y="Pass Rate", color="Model") +
  ylim(0, 1) +
  theme_classic()

# === STEP 9: Specific Predictions ===
new_classrooms <- data.frame(
  class_size = c(15, 25, 35),
  teacher_exp = c(5, 15, 25)
)

preds <- predict(model_beta, newdata=new_classrooms, type="response")
se_preds <- predict(model_beta, newdata=new_classrooms, type="response", se.fit=TRUE)

cat("\n=== Predictions for Specific Classrooms ===")
result_df <- cbind(new_classrooms, 
                   Predicted_Pass_Rate = round(preds, 3))
print(result_df)

# === APA-Style Reporting ===
cat("\n=== APA-Style Report ===")
cat("Beta regression was used to model exam pass rates(proportions in (0,1))\n")
cat("as a function of class size and teacher experience(N=", n, " classrooms).\n")
cat("Beta regression is appropriate for proportions, avoiding issues of linear\n")
cat("regression(predictions outside [0,1], heteroscedasticity). A logit link\n")
cat("function was used to map proportions to the real line.\n")
cat("\n")
cat("The model was significant(LR test vs. null: p < .001), with pseudo-R²=",
    round(model_beta$pseudo.r.squared, 2), ".\n")
cat("Class size had a significant negative effect(β=", round(coef_beta[2], 3),
    ", z=", round(summary(model_beta)$coefficients$mean[2,3], 2), ", p<.001),\n")
cat("indicating larger classes associated with lower pass rates. At mean teacher\n")
cat("experience, each additional student decreased pass rate by approximately\n")
cat(round(abs((pred_class_plus1 - pred_mean)*100), 2), "percentage points.\n")
cat("\n")
cat("Teacher experience had a significant positive effect(β=", round(coef_beta[3], 3),
    ", z=", round(summary(model_beta)$coefficients$mean[3,3], 2), ", p<.001).\n")
cat("Diagnostic plots showed good model fit with no systematic residual patterns.\n")
cat("Beta regression provided superior fit compared to OLS(which produced\n")
cat("out-of-range predictions and heteroscedastic residuals).\n")