Bayesian Regression Models for Non-Normal Data


Bayesian Regression Models for Non-Normal Data

Our last post covered how to do bayesian regression for normally distributed data in R using STAN. In this post, we’ll take a look at how to fit a regression model adapted to non-normal model in STAN using two common distributions seen in empirical data, namely, binomial and negative-binomial. As mentioned before in Bayesian modelling we use a set of sampling methods known as Markov Chain Monte Carlo (MCMC), we define a statistical model and find probabilistic estimates for the parameters.

Logistic Regression

The likelihood of a binary outcome, such as pass or fail, is estimated using logistic models (but this model can be extended to include more than two outcomes). This is accomplished by using the logit function to convert a standard regression. The main parameter that we focus on here describes odds of the outcome derived from probabilities and transformed using a logit function. The following is the code for a logistic regression model with one predictor and an intercept.

library(tidyverse)
library(kableExtra)
library(arm) 
library(emdbook) 
library(rstan)
library(rstanarm) 

set.seed(1234)
 x1 = rnorm(10000)           

 z = 1 + 2*x1       
 pr = 1/(1+exp(-z))         
 y = rbinom(10000,1,pr)     

 df = data.frame(y=y,x1=x1)
 

 
head(df)%>%kableExtra::kable()
y x1
0 -1.2070657
0 0.2774292
1 1.0844412
0 -2.3456977
1 0.4291247
1 0.5060559
 glm( y~x1,data=df,family="binomial")%>%summary()%>%pluck(coefficients)%>%kableExtra::kable()
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.036922 0.0294863 35.16621 0
x1 1.979352 0.0417219 47.44162 0
   model_stan = "
   data {
  int<lower=0> N;
  vector[N] x;
  int<lower=0,upper=1> y[N];
}
parameters {
  real alpha;
  real beta;
}
model {
  y ~ bernoulli_logit(alpha + beta * x);
}
   "
writeLines(model_stan, con = "model_stan.stan")
   cat(model_stan)
## 
##    data {
##   int<lower=0> N;
##   vector[N] x;
##   int<lower=0,upper=1> y[N];
## }
## parameters {
##   real alpha;
##   real beta;
## }
## model {
##   y ~ bernoulli_logit(alpha + beta * x);
## }
## 
stan_data <- list(
  N = 10000,
  x = df$x1,
  y = df$y
)

fit_rstan <- rstan::stan(
  file = "model_stan.stan",
  data = stan_data
)
fit_rstan
## Inference for Stan model: anon_model.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##           mean se_mean   sd     2.5%      25%      50%      75%    97.5% n_eff
## alpha     1.04    0.00 0.03     0.98     1.02     1.04     1.06     1.10  1936
## beta      1.98    0.00 0.04     1.90     1.95     1.98     2.01     2.06  1804
## lp__  -4339.45    0.03 1.05 -4342.29 -4339.80 -4339.12 -4338.74 -4338.49  1479
##       Rhat
## alpha    1
## beta     1
## lp__     1
## 
## Samples were drawn using NUTS(diag_e) at Tue Jun 21 14:03:59 2022.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

This also can be expanded to several predictors. Below you can find the code as well as some recommendations for making sense of priors.

set.seed(1234)
 x1 = rnorm(10000)           
 x2 = rnorm(10000)
 x3 = rnorm(10000)
 z = 1 + 2*x1 + 3*x2 - 1*x3       
 pr = 1/(1+exp(-z))         
 y = rbinom(10000,1,pr)     

 df2 = data.frame(y=y,x1=x1,x2=x2,x3=x3)
 

 
head(df2)%>%kableExtra::kable()
y x1 x2 x3
0 -1.2070657 -1.8168975 -1.6878627
1 0.2774292 0.6271668 -0.9552011
1 1.0844412 0.5180921 -0.6480572
0 -2.3456977 0.1409218 0.2610342
1 0.4291247 1.4572719 -1.2196940
1 0.5060559 -0.4935965 -1.5501888
 glm( y~x1+x2+x3,data=df2,family="binomial")%>%summary()%>%pluck(coefficients)%>%kableExtra::kable()
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.0108545 0.0367817 27.48253 0
x1 1.9471836 0.0494298 39.39289 0
x2 2.9598129 0.0641890 46.11088 0
x3 -0.9448316 0.0372918 -25.33619 0
model_stan2 ="


data {
  int<lower=0> N;   // number of observations
  int<lower=0> K;   // number of predictors
  matrix[N, K] X;   // predictor matrix
  int<lower=0,upper=1> y[N];      // outcome vector
}
parameters {
  real alpha;           // intercept
  vector[K] beta;       // coefficients for predictors
}
model {
  y ~ bernoulli_logit(alpha + X * beta); 
}

"

writeLines(model_stan2, con = "model_stan2.stan")
cat(model_stan2)
## 
## 
## 
## data {
##   int<lower=0> N;   // number of observations
##   int<lower=0> K;   // number of predictors
##   matrix[N, K] X;   // predictor matrix
##   int<lower=0,upper=1> y[N];      // outcome vector
## }
## parameters {
##   real alpha;           // intercept
##   vector[K] beta;       // coefficients for predictors
## }
## model {
##   y ~ bernoulli_logit(alpha + X * beta); 
## }
predictors <- df2[,2:4]

stan_data2 <- list(
  N = 10000,
  K = 3,
  X = predictors,
  y = df2$y
)


fit_rstan2 <- rstan::stan(
  file = "model_stan2.stan",
  data = stan_data2
)
## Trying to compile a simple C file
fit_rstan2
## Inference for Stan model: anon_model.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##             mean se_mean   sd     2.5%      25%      50%      75%    97.5%
## alpha       1.01    0.00 0.04     0.94     0.99     1.01     1.04     1.08
## beta[1]     1.95    0.00 0.05     1.85     1.92     1.95     1.98     2.05
## beta[2]     2.96    0.00 0.06     2.84     2.92     2.96     3.00     3.09
## beta[3]    -0.95    0.00 0.04    -1.02    -0.97    -0.95    -0.92    -0.87
## lp__    -3006.39    0.03 1.38 -3009.82 -3007.07 -3006.08 -3005.37 -3004.68
##         n_eff Rhat
## alpha    2845    1
## beta[1]  2338    1
## beta[2]  2221    1
## beta[3]  2160    1
## lp__     1788    1
## 
## Samples were drawn using NUTS(diag_e) at Tue Jun 21 14:05:38 2022.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

Negative Binomial

The Poisson distribution, which assumes that the variance is equal to the mean, is a popular choice for modelling count data. The data are said to be overdispersed when the variance exceeds the mean, the Negative Binomial distribution can be used in this case. Note that the negative binomial distribution is a probability distribution of discrete random variables.

Let’s say we have a response variable y that has a negative binomial distribution and is influenced by a collection of k explanatory variables X. The negative binomial distribution to model this data has two parameters: The μ or the expected value that need to be positive so a log link function can be used to map the linear predictor (the explanatory variables times the regression parameters) and ϕ which is the overdispersion parameter, where a small value means a large deviation from a Poisson distribution, as ϕ gets larger the negative binomial looks more and more like a Poisson distribution.

Let’s simulate some data and fit a STAN model to them:

N<-100000
df3 <-data.frame(x1=runif(N,-2,2),x2=runif(N,-2,2))
#the model
X<-model.matrix(~x1*x2,df3)
K<-dim(X)[2] #number of regression params
#the regression slopes
betas<-runif(K,-1,1)
#the overdispersion for the simulated data
phi<-5
#simulate the response
y_nb<-rnbinom(N,size=phi,mu=exp(X%*%betas))
hist(y_nb)

MASS::glm.nb(y_nb ~ X[,2:K])%>%summary()%>%pluck(coefficients)%>%kableExtra::kable()
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.0788902 0.0039192 20.12938 0
X[, 2:K]x1 -0.0815968 0.0032703 -24.95109 0
X[, 2:K]x2 0.7941334 0.0031740 250.20341 0
X[, 2:K]x1:x2 0.4301282 0.0026264 163.77308 0
data {
  int N; //the number of observations
  int K; //the number of columns in the model matrix
  int y[N]; //the response
  matrix[N,K] X; //the model matrix
}
parameters {
  vector[K] beta; //the regression parameters
  real phi; //the overdispersion parameters
}
transformed parameters {
  vector[N] mu;//the linear predictor
  mu <- exp(X*beta); //using the log link 
}
model {  
  beta[1] ~ cauchy(0,10); //prior for the intercept following Gelman 2008

  for(i in 2:K)
   beta[i] ~ cauchy(0,2.5);//prior for the slopes following Gelman 2008
  
  y ~ neg_binomial_2(mu,phi);
}
generated quantities {
 vector[N] y_rep;
 for(n in 1:N){
  y_rep[n] <- neg_binomial_2_rng(mu[n],phi); //posterior draws to get posterior predictive checks
 }
}

stan_mod = "data {
  int N; //the number of observations
  int K; //the number of columns in the model matrix
  int y[N]; //the response
  matrix[N,K] X; //the model matrix
}
parameters {
  vector[K] beta; //the regression parameters
  real phi; //the overdispersion parameters
}
transformed parameters {
  vector[N] mu;//the linear predictor
  mu <- exp(X*beta); //using the log link 
}
model {  
  beta[1] ~ cauchy(0,10); //prior for the intercept following Gelman 2008

  for(i in 2:K)
   beta[i] ~ cauchy(0,2.5);//prior for the slopes following Gelman 2008
  
  y ~ neg_binomial_2(mu,phi);
}
"
writeLines(stan_mod, con = "stan_mod.stan")

cat(stan_mod)
## data {
##   int N; //the number of observations
##   int K; //the number of columns in the model matrix
##   int y[N]; //the response
##   matrix[N,K] X; //the model matrix
## }
## parameters {
##   vector[K] beta; //the regression parameters
##   real phi; //the overdispersion parameters
## }
## transformed parameters {
##   vector[N] mu;//the linear predictor
##   mu <- exp(X*beta); //using the log link 
## }
## model {  
##   beta[1] ~ cauchy(0,10); //prior for the intercept following Gelman 2008
## 
##   for(i in 2:K)
##    beta[i] ~ cauchy(0,2.5);//prior for the slopes following Gelman 2008
##   
##   y ~ neg_binomial_2(mu,phi);
## }
stan_data3 <- list(
  N = N,
  K = K,
  X = X,
  y = y_nb
)
fit_rstan3 <- rstan::stan(
  file = "stan_mod.stan",
  data = stan_data3
)
fit_rstan3%>%summary()
## $summary
##                     mean      se_mean           sd          2.5%           25%
## beta[1]     7.899069e-02 6.825217e-05 0.0039010484  7.162531e-02  7.635923e-02
## beta[2]    -8.157657e-02 5.689044e-05 0.0032236243 -8.786688e-02 -8.376046e-02
## beta[3]     7.941039e-01 5.618629e-05 0.0031476117  7.880475e-01  7.919520e-01
## beta[4]     4.301036e-01 4.552321e-05 0.0026265744  4.248578e-01  4.283603e-01
## phi         4.908958e+00 2.039012e-03 0.0831883094  4.745765e+00  4.854542e+00

Conclusion

This tutorial provided only a quick overview of how to fit logistic and negative binomial regression models with the Bayesian software STAN using the rstan library/API and to extract a collection of useful summaries from the models. Future postings will address the question of outliers and the use of robust linear models.

References