1
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – crazjo Dec 1 '13 at 19:26
  • $\begingroup$ 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. $\endgroup$ – Wolfgang Dec 1 '13 at 19:42
  • $\begingroup$ 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. $\endgroup$ – crazjo Dec 2 '13 at 10:08
  • $\begingroup$ 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? $\endgroup$ – crazjo Dec 2 '13 at 10:29
  • 1
    $\begingroup$ 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). $\endgroup$ – Wolfgang Dec 11 '13 at 11:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.