Determine individual consistency from mixed model?

Question

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.

Update based on Wolfgangs answer:
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.

Answer 1

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.