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