Please also provide the following information in addition to your question:

- Operating System: Mac OS X El Capitan Version 10.11.6
- brms Version: 2.3.1

Background:

I am conducting a study on methylmercury (MeHg) contamination of nestling Red-winged blackbirds (*Agelaius phoeniceus*) and Great-tailed grackles (*Quiscalus mexicanus*). I collected blood samples from nestlings to measure MeHg concentration in parts per billion, and I attempted to collect blood samples when the nestlings were 5-6 days old, 8-9 days old, and 11-12 days old. However, it was common to only get one sample per individual due to predation, nest failure, or finding the nest right before fledge.

For blackbirds, I collected 424 blood samples from 243 nestlings in 84 nests from 14 different ponds. For grackles, I collected 97 blood samples from 40 nestlings in 17 nests from 2 different ponds.

The structure of my data is:

- 1-3 blood samples per individual nestling
- Individual nestlings in a nest
- Nests in a pond

Question:

I’m going through the model building process, and I am starting with simple models and working my way toward more complex models. I am currently working on a model to see how species predicts MeHg contamination in nestling birds. My predictor is gtgr, with 0 = blackbirds and 1 = grackles. I’ve looked at a varying intercepts model, and I am now working on a varying intercepts, varying slope model.

I first tested the model with the lme4 package:

model3_lmer <- lmer(MeHg ~ gtgr + (gtgr|pond_index/nest_index/ind_index),

data = rwbl_filter)

I got the following error message:

Error: number of observations (=521) <= number of random effects (=560) for term (gtgr | ind_index:(nest_index:pond_index)); the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable

So I modified the code for lme4 so it would assume the intercept and slope were uncorrelated:

model3_lmer <- lmer(MeHg ~ gtgr + (gtgr||pond_index/nest_index/ind_index),

data = rwbl_filter)

Now the model ran and here’s the summary output:

Linear mixed model fit by REML [‘lmerMod’]

Formula: MeHg ~ gtgr + ((1 | pond_index/nest_index/ind_index) + (0 + gtgr |

pond_index/nest_index/ind_index))

Data: rwbl_filter

REML criterion at convergence: 3420.6

Scaled residuals:

Min 1Q Median 3Q Max

-5.8062 -0.4411 -0.0656 0.3616 6.0925

Random effects:

Groups Name Variance Std.Dev.

ind_index…nest_index.pond_index. gtgr 0.00 0.000

ind_index…nest_index.pond_index…1 (Intercept) 0.00 0.000

nest_index.pond_index gtgr 13.29 3.645

nest_index.pond_index.1 (Intercept) 83.08 9.115

pond_index gtgr 15.10 3.885

pond_index.1 (Intercept) 116.29 10.784

Residual 22.26 4.718

Number of obs: 521, groups:

ind_index:(nest_index:pond_index), 280; nest_index:pond_index, 105; pond_index, 14

Fixed effects:

Estimate Std. Error t value

(Intercept) 21.317 3.195 6.671

gtgr 10.415 4.890 2.130

Correlation of Fixed Effects:

(Intr)

gtgr -0.044

I then moved on to run the model with brms, and here’s the code I used for that:

model3_fm <- bf(MeHg ~ gtgr + (gtgr|pond_index/nest_index/ind_index))

get_prior(model3_fm, data = rwbl_filter)

priors_3 <- c(prior(“student_t(3,0,10)”, class = “b”),

prior(“lkj(1)”, class = cor),

prior(“student_t(3,20,10)”, class = “Intercept”),

prior(“student_t(3,0,10)”, class = “sd”),

prior(“student_t(3,0,10)”, class = “sigma”))

model3 <- brm(model3_fm, data = rwbl_filter, prior = priors_3,

iter = 2000, chains = 4, cores = 4,

control = list(adapt_delta = 0.999))

I will note that I am getting 1-4 divergent iterations even with the adapt_delta set to 0.999. Here’s the summary output of the model:

Family: gaussian

Links: mu = identity; sigma = identity

Formula: MeHg ~ gtgr + (gtgr | pond_index/nest_index/ind_index)

Data: rwbl_filter (Number of observations: 521)

Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;

total post-warmup samples = 4000

Group-Level Effects:

~pond_index (Number of levels: 14)

Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat

sd(Intercept) 11.03 2.52 7.15 16.96 1344 1.00

sd(gtgr) 7.70 7.03 0.25 25.77 4000 1.00

cor(Intercept,gtgr) -0.09 0.58 -0.97 0.93 4000 1.00

~pond_index:nest_index (Number of levels: 105)

Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat

sd(Intercept) 9.24 0.79 7.83 10.90 1164 1.00

sd(gtgr) 4.72 4.27 0.16 16.43 389 1.01

cor(Intercept,gtgr) -0.02 0.51 -0.91 0.91 1082 1.01

~pond_index:nest_index:ind_index (Number of levels: 280)

Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat

sd(Intercept) 0.36 0.27 0.01 1.02 2311 1.00

sd(gtgr) 0.76 0.60 0.03 2.21 4000 1.00

cor(Intercept,gtgr) -0.12 0.58 -0.96 0.93 4000 1.00

Population-Level Effects:

Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat

Intercept 21.04 3.12 15.00 27.40 1318 1.00

gtgr 7.38 5.94 -4.92 19.01 2325 1.00

Family Specific Parameters:

Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat

sigma 4.73 0.16 4.43 5.06 4000 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample

is a crude measure of effective sample size, and Rhat is the potential

scale reduction factor on split chains (at convergence, Rhat = 1).

Warning message:

There were 4 divergent transitions after warmup. Increasing adapt_delta above 0.999 may help.

See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

One of the things I noticed was that the estimated rho for each level was close to zero but there was a lot of uncertainty associated with this estimate. Lme4 also estimated low correlation between the varying intercepts and varying slopes. Would it make sense and be justifiable to run the model again assuming that there is no correlation between the intercepts and slopes? Would that make the model computationally easier and would I get fewer divergent iterations?

I’m not as familiar with the pros and cons of using models with uncorrelated varying intercepts and slopes, and I’m having a hard time finding literature on the subject. I would very much appreciate if anyone has advice on the matter.

Thank you, and sorry for the long winded question.

Tom