# Attempting to compare Bayesian and Frequentist mixed effects models

This question may be better suited for stack overflow (happy to move it if deemed too off topic).

I am currently in the process of learning Bayesian analysis using stan in R as my software. Currently as part of this process I am trying to build regression models with uninformative priors as my understanding was this will lead to very similar point estimates of parameters as Frequentist modelling (for example the standard linear model I implemented had nearly identical parameter estimates in both cases).

The main goal of this is so that I can tell if I am programming the model correctly.

I am currently trying to extend this approach to mixed effects models however I can not get even remotely comparable parameter estimates between the two implementations.

I guess my question boils down to two parts:

a) Will uninformative priors in Bayesian mixed models (heirarchical models ?) lead to similar point estimates as Frequentist appraochs ? and if yes

b) Where am I going wrong in my code / approach ?

APPENDIX:

data:

machines <- structure(list(Worker = structure(c(4L, 4L, 4L, 2L, 2L, 2L, 5L,
5L, 5L, 3L, 3L, 3L, 6L, 6L, 6L, 1L, 1L, 1L, 4L, 4L, 4L, 2L, 2L,
2L, 5L, 5L, 5L, 3L, 3L, 3L, 6L, 6L, 6L, 1L, 1L, 1L, 4L, 4L, 4L,
2L, 2L, 2L, 5L, 5L, 5L, 3L, 3L, 3L, 6L, 6L, 6L, 1L, 1L, 1L), .Label = c("6",
"2", "4", "1", "3", "5"), class = c("ordered", "factor")), Machine = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L), .Label = c("A", "B", "C"), class = "factor"),
score = c(52, 52.8, 53.1, 51.8, 52.8, 53.1, 60, 60.2, 58.4,
51.1, 52.3, 50.3, 50.9, 51.8, 51.4, 46.4, 44.8, 49.2, 62.1,
62.6, 64, 59.7, 60, 59, 68.6, 65.8, 69.7, 63.2, 62.8, 62.2,
64.8, 65, 65.4, 43.7, 44.2, 43, 67.5, 67.2, 66.9, 61.5, 61.7,
62.3, 70.8, 70.6, 71, 64.1, 66.2, 64, 72.1, 72, 71.1, 62,
61.4, 60.5)), .Names = c("Worker", "Machine", "score"), row.names = c(NA,
-54L), class = "data.frame")


model:

$Y_{ijk} = \beta_j + b_i + \epsilon_{ijk}$
$b_i \sim N(0 , \sigma_1^2)$
$\epsilon_{ijk} \sim N(0 ,\sigma^2)$

where
Y = score
$\beta_j$ = Machine (A / B / C) - The fixed effects
$b_i$ = Worker ( 1 / 2 / 3 / 4 / 5 / 6) - The random effects
k = replicate indicator

Frequentist Approach

library(lme4)
fm1 <- lmer(
formula = score ~ Machine -1 + ( 1 | Worker ),
data = machines
)
summary(fm1)
ranef(fm1)


Bayesian Approach

note in addition I apply priors

$Ga( 1 / \sigma_1 | 5 , 5 )$
$Ga( 1 / \sigma | 5 , 5 )$
$N(\beta_j | 0, 100)$

stan model:

data {
int n;          // number of observations
int k_fixed;    // number of fixed effects
int k_random;   // number of random effects

matrix[n,k_fixed] fixed_data;
matrix[n,k_random] random_data;

matrix[n, 1] raw_y; //the response variable

real a1;  // prior a for lam1
real b1;  // prior b for lam1

real a2;  // prior a for lam2
real b2;  // prior b for lam2

real m;   // prior mu for fixed effects
real v;   // prior variance for fixed effects

}

transformed data{
vector[n] y;
y = col(raw_y, 1);  // convert y matrix to vector
}

parameters {
real<lower=0> lambda_one;
real<lower=0> lambda_two;

vector[k_fixed] beta;
vector[k_random] zeta;
}

transformed parameters {
vector[n] mu;
mu = fixed_data * beta + random_data * zeta;
}

model {

target += gamma_lpdf( lambda_one | a1 , b1);

target += gamma_lpdf( lambda_two | a2 , b2);

target += normal_lpdf( zeta | 0 , sqrt(1/lambda_one) ) ;

target += normal_lpdf( beta | m , sqrt(v) );

target += normal_lpdf( y | mu ,  sqrt(1 / lambda_two ));
}


R code:

library(rstan)
library(tidyverse)

### Use model.matrix to build a design matrix for stan to use
Y <- machines %>% select( score)
x_fixed  <- model.matrix( score ~ Machine -1 , machines)
x_random <- model.matrix( score ~ Worker -1 , machines)

dat <- list(
n = nrow(machines) ,
k_fixed = ncol(x_fixed) ,
k_random = ncol(x_random) ,

fixed_data = x_fixed,
random_data = x_random,

raw_y = Y ,

a1 = 5,
b1 = 5,

a2 = 5,
b2 = 5,

m = 0 ,
v = 100
)

fit <- stan(
file="mixed_effects.stan" ,
data = dat ,
iter = 20000,
warmup = 5000,
chains = 1
)

samples <- rstan::extract(fit)

(1/samples$lambda_one) %>% mean (1/samples$lambda_two ) %>%  mean
samples$beta %>% apply(MARGIN = 2 , mean) samples$zeta %>% apply(MARGIN = 2 , mean)

• No, the estimates will not always be the same. The differences are especially pronounced in multilevel models for the variance estimates. Paradoxically, with small data, the noninformative prior is some times the most informative prior. I would use a half cauchy(0, SD), where SD us determined in reference to the SD of y. I usually go a bit smaller than that. – D_Williams Aug 16 '17 at 1:05
• Thank you for your comment @D_Williams - Indeed my priors were actually over informative - Using the half cauchy instead got me point estimates very similar to that of the frequentist approach (also with help from the answer below) – gowerc Aug 16 '17 at 21:36

I can't comment your code because I don't understand it. I proceed as follows and I get results similar to the frequentist ones.

stancode <- 'data {
int<lower=1> N; // number of observations
vector[N] y; // observations
int<lower=1> J;
int<lower=1> I;
int<lower=1,upper=J> jj[N];
int<lower=1,upper=I> ii[N];
}
parameters {
real<lower=0> sigmab;
real<lower=0> sigmaw;
vector[I] b;
vector[J] beta;
}
transformed parameters {
vector[N] mu;
for(k in 1:N){
mu[k] = beta[jj[k]] + b[ii[k]];
}
}
model {
y ~ normal(mu, sigmaw); // sampling distribution of the observations
sigmaw ~ cauchy(0, 5); // prior on the within standard deviation
b ~ normal(0, sigmab);
sigmab ~ cauchy(0,5);
beta ~ normal(0, 100);
}'

stanmodel <- stan_model(model_code = stancode, model_name="stanmodel")

standata <- list(N=nrow(machines),
y = machines$score, J = nlevels(machines$Machine),
jj = as.integer(machines$Machine), I = nlevels(machines$Worker),
ii = as.integer(machines\$Worker))

samples <- sampling(stanmodel, data = standata, #init=inits,
iter = 10000, warmup = 1000, chains = 4)

library(coda)
codasamples <- do.call(mcmc.list,
plyr::alply(
rstan::extract(
samples,
pars=c("sigmab", "sigmaw", "beta[1]",
"beta[2]", "beta[3]"),
permuted=FALSE), 2, mcmc))

> summary(codasamples)

Iterations = 1:9000
Thinning interval = 1
Number of chains = 4
Sample size per chain = 9000

1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:

Mean     SD Naive SE Time-series SE
sigmab   5.901 2.3905 0.012599       0.031054
sigmaw   3.233 0.3445 0.001816       0.002652
beta[1] 52.205 2.7963 0.014738       0.043209
beta[2] 60.167 2.7898 0.014704       0.042541
beta[3] 66.114 2.7925 0.014718       0.043094

2. Quantiles for each variable:

2.5%    25%    50%    75%  97.5%
sigmab   3.090  4.399  5.390  6.809 11.554
sigmaw   2.646  2.988  3.201  3.445  3.988
beta[1] 46.793 50.629 52.238 53.834 57.592
beta[2] 54.787 58.588 60.204 61.783 65.495
beta[3] 60.701 64.526 66.143 67.752 71.424

• Hi @Stéphane Laurent - Thank you for this; it was exactly what I needed and seems to generate the expected numbers. Additionally I was able to use your working code to help fix my own. For reference the issue with my code was that my priors were way too informative ! Using a cauchy prior as your have (and as suggest above in a comment) and increasing the variance of beta I was able to get matching values - Thank you !!!. – gowerc Aug 16 '17 at 20:07