Uncorrelated Varying Intercept and Slope Models


#1

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. 1-3 blood samples per individual nestling
  2. Individual nestlings in a nest
  3. 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


#2

I looked over my post, and I realize copying and pasting summary outputs doesn’t transition well into the post. Anyone have tips on how do to this so the summary outputs are more legible? Sorry, I’m new to posting things on this forum, as well as github.


#3

Use three backticks both before and after a chunk of code or output. Or you can click on the </> icon or press Control+Shift+C.


#4

To be honest, I don’t think 1-4 divergent transitions will be a relevant problem for this kind of model (please don’t tell Michael, I said that), but I can’t give you any guarantee of course. You may want to specify stronger priors on the varying effects SDs (see ?set_prior).

Personally, I only ignore correlations when it is not computationally feasible and with only 521 it seems this isn’t a problem. So I would probably just include the correlations regardless.


#5

Thanks for the input, Paul!

In case other people want to see my model code and outputs, here is the lme4 code I used for an uncorrelated intercept and slope model:

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

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

And here is the brms code I ran for model with correlated intercepts and slopes, and its respective outputs:

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))
 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).

To reiterate the question, what are people’s thoughts when comes the correlation between intercepts and slopes? I know it’s preferred to keep the correlation, but is it justified to assume the two are uncorrelated when the model is estimating low correlation and it might be computationally easier to run?