Determine individual consistency from mixed model?

To investigate the difference between three tasks I tested 30 individuals on each of the three tasks and ran the linear model as follows:

model <- lme(behaviour ~ task, random = ~ 1|subject, data=dat, method="ML")


With contrasts using glht to investigate if task 1 was different from task 2 and task 3.

To furthermore investigate if individuals are consistent in their behaviour across the three different tasks I now simply test this by running three correlation tests such as cor.test(behaviour.task1, behaviour.task2). However, I was wondering if it would not possible to get this information directly from the model, i.e. by looking at the intercept.

I have found a couple papers that have done something similar like this by reporting if the intercept is significant, however they don't explain how they performed their analyses nor can I see how showing that the intercept is significantly different from 0 in a mixed model means that individuals were consistent.

First, to calculate the overall ICC of the model I ran

as.numeric(VarCorr(model)[1,1]) / (as.numeric(VarCorr(model)[1,1]) + model$sigma^2)  This shows the ICC = 0.336. To determine if the consistency of individuals is significant, I calculated if the ICC is different from 0 as follows: model0 <- gls(behaviour ~ task, data=dat, method="ML") anova(model0, model)  The output shows: Model df AIC BIC logLik Test L.Ratio p-value model0 1 4 1211.739 1221.648 -601.8695 model 2 5 1205.301 1217.688 -597.6505 1 vs 2 8.438081 0.0037  Thus the ICC is significantly larger than 0 (P < 0.005). Now, as I expect the tasks will differently affect the consistency of individuals, I run a new model as follows: model2 <- gls(behaviour ~ task, correlation = corSymm(form = ~ 1 | subject), data=dat, method="ML") summary(model2)  And get the following output: Generalized least squares fit by maximum likelihood Model: t.notcovered ~ session Data: cot.sum Subset: maleSub AIC BIC logLik 1203.922 1221.263 -594.961 Correlation Structure: General Formula: ~1 | rat Parameter estimate(s): Correlation: 1 2 2 0.233 3 0.453 0.584  thus showing that the correlation between task 1 and 2 is 0.233, between task 1 and 3 is 0.453 and between task 2 and 3 is 0.584. Now to see if this model provides a better fit than the model that assumes a single common correlation I ran the following: anova(model, model2) Model df AIC BIC logLik Test L.Ratio p-value model 1 5 1205.301 1217.688 -597.6505 model2 2 7 1203.922 1221.264 -594.9610 1 vs 2 5.378838 0.0679  showing that the more complex model does not yield a significantly better fit, suggesting that the consistency does not differ significantly between the tasks, but as we know from the simpler model, individuals were overall consistent in their behaviour. I added the above information with the idea that other people with a similar question may learn from this as well. 1 Answer You can estimate the intraclass correlation coefficient (ICC) (wikipedia link). It tells you how correlated the behavioral responses are for the same individual. It is defined as: $$ICC = \tau^2 / (\tau^2 + \sigma^2),$$ where$\tau^2$is the "intercept" variance and$\sigma^2$is the residual variance. Substituting the estimated values into the equation then gives you the estimated ICC. Here is an example (with some simulated data): library(nlme) set.seed(1231233) dat <- data.frame(subject=rep(1:30, each=3), task=factor(1:3)) dat$behaviour <- rep(rnorm(30, 0, 1), each=3) + c(2,4,5) + rnorm(30*3, 0, 1)

res <- lme(behaviour ~ task, random = ~ 1|subject, data=dat, method="ML")
as.numeric(VarCorr(res)[1,1]) / (as.numeric(VarCorr(res)[1,1]) + res$sigma^2)  Running this yields 0.4714821, which is quite close to the true ICC for these data (since I simulated the data with an intercept variance of 1 and a residual variance of 1, the true ICC is:$1 / (1 + 1) = 0.5$). Note that a zero intercept variance implies that the responses for the same individual are uncorrelated (since the ICC is then equal to 0). Therefore, we can test whether the ICC is larger than 0 by testing whether the intercept variance is larger than zero. This can be done as follows: res0 <- gls(behaviour ~ task, data=dat, method="ML") anova(res0, res)  with output:  Model df AIC BIC logLik Test L.Ratio p-value res0 1 4 345.1997 355.1989 -168.5998 res 2 5 328.8656 341.3647 -159.4328 1 vs 2 18.33407 <.0001  So clearly the model with the intercept variance component fits better - which implies that the ICC is significantly larger than 0. Finally, if you want to let the correlation between different task pairs to differ, you can fit the model: res2 <- gls(behaviour ~ task, correlation = corSymm(form = ~ 1 | subject), data=dat, method="ML") summary(res2)  with output: Generalized least squares fit by maximum likelihood Model: behaviour ~ task Data: dat AIC BIC logLik 329.6571 347.1557 -157.8285 Correlation Structure: General Formula: ~1 | subject Parameter estimate(s): Correlation: 1 2 2 0.387 3 0.394 0.633 Coefficients: Value Std.Error t-value p-value (Intercept) 2.185489 0.2925499 7.470483 0 task2 1.785381 0.3239530 5.511233 0 task3 2.993919 0.3220204 9.297298 0 [...]  So, .387 is the estimated correlation between tasks 1 and 2, .394 is the estimated correlation between tasks 1 and 3, and .633 is the estimated correlation between tasks 2 and 3. Especially the last value seems quite a bit different. But does that model really provide a better fit than the model that assumes a single common correlation? You can test that with: anova(res, res2) Model df AIC BIC logLik Test L.Ratio p-value res 1 5 328.8656 341.3647 -159.4328 res2 2 7 329.6571 347.1557 -157.8285 1 vs 2 3.208543 0.201  And the answer is no. The more complex model does not yield a significantly better fit. So, for these data, the ICC of .47 seems to provide a reasonable estimate of the correlation in the responses across the 3 tasks. • Thanks, I did not know about ICC. Calculating ICC according to your explanation I get a value of 0.35. How do I know the significance/effect size of this effect? Secondly, as I test individuals across three tasks and expect consistency to differ between tasks 1 and 2 compared to 2 and 3 or 1 and 3, is there a way to do some kind of contrast on the ICC or would it in the end be best to do a correlation test instead? Thanks again. Dec 1, 2013 at 19:26 • The ICC is a correlation, so you can interpret it that way. Whether a correlation of .35 is small, medium, or large for your particular application is something you will have to think about. As for testing its significance - you can just test whether the intercept variance is larger than zero. That is the same as testing whether the ICC is larger than zero. You can also allow the correlation to differ for the different task pairs. I'll update my answer to show you how that can be done. Dec 1, 2013 at 19:42 • thanks @Wolfgang, I updated my question to show the steps in relation to my dataset in the hope it may help others. Please let me know if my interpretation was wrong. Thanks again for your great help. Dec 2, 2013 at 10:08 • Sorry, an extra question: if I run a mixed model with an interaction like lme(behaviour ~ task*sex, random = ~ 1|subject, data=dat, method="ML") would it be possible to get separate ICC per sex and within tasks per sex and then to compare these, i.e. to investigate if males and females differed in their consistency across the tasks? Dec 2, 2013 at 10:29 • Your interpretation seems right. The model allowing for different correlations between the three task pairs does not fit significantly better strictly speaking ($p = .07\$), but its AIC is lower, which suggests that it may be a better approximation to the true model (then again, its BIC is higher). I would actually tentatively conclude that there may be evidence that the three task pairs correlate differently. I'll get to your follow-up question later (this part is more tricky). Dec 11, 2013 at 11:15