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)
