6
$\begingroup$

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:

Trace plots when sample size = 200 and number of sites = 10

The following are trace plots for scenario B:

Trace plots when sample size = 1000 and number of sites = 10

The following are trace plots for scenario C:

Trace plots when sample size = 1000 and number of sites = 50

The following are trace plots for when sample size = 200 and number of sites = 40. In response to Tommaso's comments below.

Trace plots when sample size = 200 and number of sites = 40.

$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – Tommaso Guerrini Mar 24 '17 at 10:50
  • $\begingroup$ 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 :) $\endgroup$ – Tommaso Guerrini Mar 24 '17 at 10:52
  • $\begingroup$ Does it make sense to you? $\endgroup$ – Tommaso Guerrini Mar 24 '17 at 10:59
  • 1
    $\begingroup$ 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. $\endgroup$ – AP30 Mar 24 '17 at 14:59
  • 1
    $\begingroup$ 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? $\endgroup$ – AP30 Mar 28 '17 at 11:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.