# Why does increasing number of observations in linear mixed model cause Bayesian modelling approach to fail?

I have a fairly good understanding of the theory behind Bayesian modeling and I have started to attempt some practical modeling using jags in R. I have been following examples online to get used to the syntax. However in the below example I simply increase the number of observations and it has an unexpected effect, inducing bias in the results. I am trying to figure out why.

I have been using the resource here. In the third example, they attempt to create a linear mixed model: b_mass ~ b_length + r.e(site). Body mass (of a snake) has a linear relationship with the continuous variable length, and a random intercept for the site at which the snake was found.

Questions:
1. Why does the Bayesian model perform so poorly when the number of observations increases?

2. Is the ratio of observations to number of random effects important in a mixed model? It doesn't seem to be when fitting a standard linear mixed effects model under a frequentist framework.

3. What do trace plots with lots of outliers/long tails indicate? And what adjustments should be made for them?

Details:

SCENARIO (A):
The code to generate the data and fit the model, taken directly from the website above is:

## First create the data
set.seed(42)

samplesize <- 200
nsites <- 10 # Number of sites
b_length <- sort(rnorm(samplesize)) # Body length (explanatory variable)
sites <- sample(1:10, samplesize, replace = T) # Site (grouping variable)
table(sites)

int_true_mean <- 45 # True mean intercept
int_true_sigma <- 10 # True standard deviation of random intercepts
int_true_sites <- rnorm(n = nsites, mean = int_true_mean, sd = int_true_sigma) # True intercept of each site

# Intercept of each snake individual (depending on the site where it was captured):
sitemat <- matrix(0, nrow = samplesize, ncol = nsites)
for (i in 1:nrow(sitemat)) sitemat[i, sites[i]] <- 1
int_true <- sitemat %*% int_true_sites

slope_true <- 10 # True slope (coefficient of length)
mu <- int_true + slope_true * b_length # True means of normal distributions
sigma <- 5 # True standard deviation of normal distribution (random variation)

b_mass <- rnorm(samplesize, mean = mu, sd = sigma) # Body mass (response variable)

snakes3 <- data.frame(b_length = b_length, b_mass = b_mass, site = sites)
head(snakes3)

## Create list for jags to work on, include number of sites as an additional element in the list:
Nsites <- length(levels(as.factor(snakes3$site))) jagsdata_s3 <- with(snakes3, list(b_mass = b_mass, b_length = b_length, site = site, N = length(b_mass), Nsites = Nsites)) ## Define model lm3_jags <- function(){ # Likelihood: for (i in 1:N){ b_mass[i] ~ dnorm(mu[i], tau) # tau is precision (1 / variance) mu[i] <- alpha + a[site[i]] + beta * b_length[i] # Random intercept for site } # Priors: alpha ~ dnorm(0, 0.01) # intercept sigma_a ~ dunif(0, 100) # standard deviation of random effect (variance between sites) tau_a <- 1 / (sigma_a * sigma_a) # convert to precision for (j in 1:Nsites){ a[j] ~ dnorm(0, tau_a) # random intercept for each site } beta ~ dnorm(0, 0.01) # slope sigma ~ dunif(0, 100) # standard deviation of fixed effect (variance within sites) tau <- 1 / (sigma * sigma) # convert to precision } ## Define initial values init_values <- function(){ list(alpha = rnorm(1), sigma_a = runif(1), beta = rnorm(1), sigma = runif(1)) } # Define parameters to kept in results params <- c("alpha", "beta", "sigma", "sigma_a") # Fit the model fit_lm3 <- jags(data = jagsdata_s3, inits = init_values, parameters.to.save = params, model.file = lm3_jags, n.chains = 3, n.iter = 20000, n.burnin = 5000, n.thin = 2, DIC = F) # I have reduced n.thin to 2, so output slightly different from website  The output is: ## Output fit_lm3 mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff alpha 42.583 4.481 30.875 41.144 43.476 45.208 48.110 1.004 3600 beta 9.893 0.349 9.208 9.662 9.895 10.127 10.582 1.001 19000 sigma 4.755 0.248 4.300 4.581 4.745 4.915 5.272 1.001 22000 sigma_a 8.803 4.386 4.610 6.292 7.689 9.845 19.898 1.002 2400  The correct answers should be: alpha (average intercept) = 45. beta (coefficient for length)= 10. sigma (residual sd)= 5. sigma_a (sd of intercepts) = 10. The results are quite good, 42.583 is close to 45, 9.893 is close to 10, 4.755 is close to 5, and 8.803 is close to 10. SCENARIO (B): However I thought it would be nice to get ever more accurate results, so I increased the number of observations to 1000. samplesize <- 1000  The rest of the code remains the same, and the output is now:  mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff alpha 33.587 8.982 11.181 28.906 35.602 40.078 45.933 1.001 6100 beta 9.839 0.164 9.519 9.727 9.838 9.950 10.165 1.001 22000 sigma 5.265 0.118 5.038 5.184 5.263 5.343 5.503 1.001 22000 sigma_a 19.635 9.749 8.960 12.965 16.852 23.377 45.792 1.001 5300  The coefficient for length and residual variation estimates are still good (beta and sigma), but the estimates for variation of intercepts, and the mean of the intercepts are now very poor. SCENARIO (C): I thought this might have to do with the ratio of #observations/#random intercepts, so I increased the number of sites to 50, and the model performs well again. However I'm not sure why this is the case. samplesize <- 1000 nsites <- 50 # Number of sites b_length <- sort(rnorm(samplesize)) # Body length (explanatory variable) sites <- sample(1:50, samplesize, replace = T) # Site (grouping variable)  Output:  mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff alpha 43.398 1.533 40.296 42.395 43.420 44.436 46.353 1.001 22000 beta 9.836 0.162 9.522 9.726 9.836 9.945 10.156 1.001 22000 sigma 5.013 0.114 4.794 4.936 5.011 5.089 5.241 1.001 22000 sigma_a 10.816 1.151 8.832 10.009 10.728 11.517 13.362 1.001 22000  I have attached trace plots below for scenario (A). Clearly there are some issues as the densities for alpha and sigma_a have long tails, and many 'outlier' values. I would like to post more but my reputation is not high enough. The trace plots for scenario (B) look similar (poor) whereas for scenario (C) they look good. This does not fit in with the results. Finally, I fit a standard linear mixed effects model in all these scenarios, and it performed consistently well across all three. The following are trace plots for scenario A: The following are trace plots for scenario B: The following are trace plots for scenario C: The following are trace plots for when sample size = 200 and number of sites = 40. In response to Tommaso's comments below. • Very very very nice question! Let me guess: little data, small number of sites chosen ($A$) works well, large data, small number of sites chosen ($B$) works bad, large data, high number of sites ($C$) works well. I guess in model$B$(and$A$) you are trying to model with a few sites data coming from many different sites, so that your random intercept coefficient is trying to average between different sites. This to me can be seen in the long positive tail for sigm_a, which may mean that you need a higher standard deviation over the random effects – Tommaso Guerrini Mar 24 '17 at 10:50 • Instead, when you increase the number of sites the variation induced by different sites is already caught by a different intercept, so that you don't need to allow for a greater variance.. Try model$A\$ with more sites and show the posterior distribution for alpha and sigma please :) – Tommaso Guerrini Mar 24 '17 at 10:52
• Does it make sense to you? – Tommaso Guerrini Mar 24 '17 at 10:59
• In model B (and A) I am trying to model with a few sites data that only comes from a few sites. The number of sites the data is generated by is exactly the number of sites I try to model it with. You are correct in understanding of the issue, but i'm not entirely sure I understand your solution. Please expand... I have done as you asked (scenario A but with 40 sites instead of 10, this means on average about 5 snakes per site). The estimates are good again: alpha 42.167 beta 9.846 sigma 5.283 sigma_a 9.630 And I will upload trace plots above. – AP30 Mar 24 '17 at 14:59
• Tommaso, thank you again, but I am not sure I follow your logic. My most important question is 1) When the number of observations is increased, but the number of sites stays the same (A -> B), do we predict the variation in the random effects so poorly. In this scenario, do we not have more information about each site as there are five times as many snakes in each? Why should the predictive performance drop so quickly? – AP30 Mar 28 '17 at 11:59