--- name: regression-models description: Bayesian regression models including linear, logistic, Poisson, negative binomial, and robust regression with Stan and JAGS implementations. --- # Regression Models ## Linear Regression ### Stan ```stan data { int N; int K; matrix[N, K] X; vector[N] y; } parameters { real alpha; vector[K] beta; real sigma; } model { alpha ~ normal(0, 10); beta ~ normal(0, 5); sigma ~ exponential(1); y ~ normal(alpha + X * beta, sigma); } generated quantities { array[N] real y_rep; for (n in 1:N) y_rep[n] = normal_rng(alpha + X[n] * beta, sigma); } ``` ### JAGS ``` model { for (i in 1:N) { y[i] ~ dnorm(mu[i], tau) mu[i] <- alpha + inprod(X[i,], beta[]) } alpha ~ dnorm(0, 0.001) for (k in 1:K) { beta[k] ~ dnorm(0, 0.001) } tau ~ dgamma(0.001, 0.001) sigma <- 1\/sqrt(tau) } ``` ## Logistic Regression ### Stan ```stan data { int N; int K; matrix[N, K] X; array[N] int y; } parameters { real alpha; vector[K] beta; } model { alpha ~ normal(0, 2.5); beta ~ normal(0, 2.5); y ~ bernoulli_logit(alpha + X * beta); } ``` ### JAGS ``` model { for (i in 1:N) { y[i] ~ dbern(p[i]) logit(p[i]) <- alpha + inprod(X[i,], beta[]) } alpha ~ dnorm(0, 0.4) # SD ≈ 1.58 for (k in 1:K) { beta[k] ~ dnorm(0, 0.4) } } ``` ## Poisson Regression ### Stan ```stan model { alpha ~ normal(0, 5); beta ~ normal(0, 2.5); y ~ poisson_log(alpha + X * beta); } ``` ### JAGS ``` model { for (i in 1:N) { y[i] ~ dpois(lambda[i]) log(lambda[i]) <- alpha + inprod(X[i,], beta[]) } } ``` ## Negative Binomial (Overdispersed Counts) ### Stan ```stan parameters { real alpha; vector[K] beta; real phi; \/\/ Overdispersion } model { phi ~ exponential(1); y ~ neg_binomial_2_log(alpha + X * beta, phi); } ``` ## Robust Regression (Student-t Errors) ### Stan ```stan parameters { real alpha; vector[K] beta; real sigma; real nu; \/\/ Degrees of freedom } model { nu ~ gamma(2, 0.1); \/\/ Prior on df y ~ student_t(nu, alpha + X * beta, sigma); } ``` ## QR Decomposition (For Correlated Predictors) ```stan transformed data { matrix[N, K] Q = qr_thin_Q(X) * sqrt(N - 1.0); matrix[K, K] R = qr_thin_R(X) \/ sqrt(N - 1.0); matrix[K, K] R_inv = inverse(R); } parameters { vector[K] theta; real sigma; } model { y ~ normal(Q * theta, sigma); } generated quantities { vector[K] beta = R_inv * theta; } ``` ## Prior Recommendations | Parameter | Weakly Informative | Reference | |-----------|-------------------|-----------| | Intercept | `normal(0, 10)` | Scale of outcome | | Coefficients | `normal(0, 2.5)` | Gelman et al. | | SD (sigma) | `exponential(1)` | Half-normal alternative | | Logistic coef | `normal(0, 2.5)` | ~4 logit units = extreme |