I am trying to estimate the intraclass correlation coefficient (ICC; or repeatability, R, as termed in ecology) of a nonlinear parameter at the individual level after fitting a model with brms. My data are spring molt progression rates of individual bighorn sheep. I have multiple molt estimates per individual per year (≈7 per id-year) which allow an estimation of a random intercept of year nested within id for the inflection point of a nonlinear logistic growth curve (see phi1 below). A simplified version of my code is the following:
# Define model priors prior1 = prior(normal(1, 0.5), nlpar = "phi2") + prior(normal(175, 8), nlpar = "phi1") + prior(cauchy(2,3), class="sd", group="id", coef="Intercept", nlpar="phi1") + prior(cauchy(8,4), class="sd", group="id:year", coef="Intercept", nlpar="phi1") + prior(student_t(3, 0, 2.5), class="sigma") # Fit model final <- brm(bf(moltProp ~ 1 / (1 + exp((phi1-jj)/exp(phi2))), phi2 ~ 1, phi1 ~ 1 + (1|id/year), nl=TRUE), data = molt, prior = prior1, iter = 8000, warmup = 4000, cores = 4, chains = 4) # jj is Julian date
To estimate an ICC at the individual level for a linear multilevel model with (1|ID) I would simply use the the standard deviation estimated from the ID random intercept (sd_id__Intercept) and the model residuals (sigma):
hyp <- "sd_id__Intercept^2 / (sd_id__Intercept^2 + sigma^2) = 0" (hyp <- hypothesis(final, hyp, class = NULL))
However, in the case of the nonlinear model, (1|ID/Year) is estimated for phi1 only, whereas the sigma that is returned in the summary is at the moltprop (response variable) level. Am I correct in thinking that sd_id__phi1_Intercept and sd_id:year__phi1_Intercept, respectively represent among- and within-individual variation for phi1? If so, could I estimate an ICC for phi1 using the following code?
hyp <- "sd_id__phi1_Intercept^2 / (sd_id__phi1_Intercept^2 + sd_id:year__phi1_Intercept^2) = 0" (hyp <- hypothesis(final, hyp, class = NULL))
If not, how could an ICC be estimated for phi1 at the individual level? In other words, is there a way to correctly extract or estimate within-individual variance for phi1 only?
Data to run code are shared as a CSV file.
molt.csv (46.8 KB)
- Operating System: Windows
- brms Version: 2.15.0
Thank you very much!