Regression Models

Regression Models

Linear Regression

Stan

data {
  int<lower=0> N;
  int<lower=0> K;
  matrix[N, K] X;
  vector[N] y;
}
parameters {
  real alpha;
  vector[K] beta;
  real<lower=0> 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

data {
  int<lower=0> N;
  int<lower=0> K;
  matrix[N, K] X;
  array[N] int<lower=0,upper=1> 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

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

parameters {
  real alpha;
  vector[K] beta;
  real<lower=0> phi;  // Overdispersion
}
model {
  phi ~ exponential(1);
  y ~ neg_binomial_2_log(alpha + X * beta, phi);
}

Robust Regression (Student-t Errors)

Stan

parameters {
  real alpha;
  vector[K] beta;
  real<lower=0> sigma;
  real<lower=1> 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)

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<lower=0> sigma;
}
model {
  y ~ normal(Q * theta, sigma);
}
generated quantities {
  vector[K] beta = R_inv * theta;
}

Prior Recommendations