Estimate reliability using ICC

In our example where we got 30 dogs and weighed them, and then weighed them again, we might have data like the below. We are wanting to separate out the variance within each dog’s observed weights and attribute it to measurement error. The variance between dogs is taken as the true variance. So Dougal first got 20.6, and then 25.9, so if we take his average of 23.25 we can think of ‘within-dougal-variance’ as how his measurements vary around his average (due to measurement error). And Dougal’s average 23.25 differs from Rufus’s average of 29.7, and so on - these values have a variance which we take to reflect Dougal truly varying from Rufus and Daisy etc.

dogdata

  dog    measurement_1 measurement_2
  <chr>          <dbl>         <dbl>
1 dougal          20.6         25.9 
2 rufus           28.1         31.3 
3 daisy           10.7         9.79
4 bella           19.8         17.2 
5 shell           13.4         8.88
6 coco            20.8         19.6 
. ...             ...          ...
. ...             ...          ...

The ICC() function from the psych package will calculate the ICC for us.

library(psych)
ICC(dogdata[,2:3])

Call: ICC(x = dogdata[, 2:3])

Intraclass correlation coefficients 
                         type  ICC  F df1 df2       p lower bound upper bound
Single_raters_absolute   ICC1 0.90 19  29  30 1.0e-12        0.81        0.95
Single_random_raters     ICC2 0.90 19  29  29 2.1e-12        0.81        0.95
Single_fixed_raters      ICC3 0.90 19  29  29 2.1e-12        0.80        0.95
Average_raters_absolute ICC1k 0.95 19  29  30 1.0e-12        0.89        0.98
Average_random_raters   ICC2k 0.95 19  29  29 2.1e-12        0.89        0.98
Average_fixed_raters    ICC3k 0.95 19  29  29 2.1e-12        0.89        0.98

 Number of subjects = 30     Number of Judges =  2
See the help file for a discussion of the other 4 McGraw and Wong estimates,

The ICC function doesn’t know what kind of reliability we want to measure, because it doesn’t know how exactly we got our measurement multiple times. So to cover all its bases, the ICC function spits out a whole bunch of different types of ICC. Here’s what each type of ICC means:

ICC1: Each dog is measured at a different time point and the time points are selected at random.
ICC2: Dogs are all measured at a random sample of k time points.
ICC3: Dogs are all measured at a specific fixed set of k time points.

Then for each of these, we also get the equivalent for if there’s more than one rating for each dog at each time point.

In the official terminology that you’ll see on the help pages of the ICC function, the same ideas will be phrased more generally as:

ICC1: Each target is rated by a different judge and the judges are selected at random.
ICC2: A random sample of k judges rate each target. The measure is one of absolute agreement in the ratings.
ICC3: A fixed set of k judges rate each target. There is no generalization to a larger population of judges.

This all uses the language of “judges” and “targets”, which in our case corresponds to “timepoints” and “dogs”. There’s a nice paper (Koo & Li 2016) here explaining this in more depth and with some flowcharts to help us think about which one we want. In our case, we have a random sample of timepoints (it’s not like we measured all dogs before and after some event - we don’t expect any systematic differences between timepoints), so ICC2 is probably what we want.

Rabbit hole: doesn’t ICC come up in mixed models!?

It did! We already saw ICC! and as it happens, it’s the same thing.
This box shows how we can use the multilevel model machinery to reconstruct the same ICC value computed above.

If we pivot our data to be long:

doglong <- dogdata |> 
  pivot_longer(measurement_1:measurement_2,
               names_to = "time", values_to = "score")

doglong

   dog    time          score
   <chr>  <chr>         <dbl>
 1 dougal measurement_1 20.6 
 2 dougal measurement_2 25.9 
 3 rufus  measurement_1 28.1 
 4 rufus  measurement_2 31.3 
 5 daisy  measurement_1 10.7 
 . ...    ...           ...
 . ...    ...           ...

We can then use a multilevel/mixed effects model to regress the ratings on to a model with just an intercept and by-dog and by-timepoint random intercepts:

doglong$score=ifelse(doglong$time=="measurement_2",doglong$score+4,doglong$score)
mod <- lmer(score ~ 1 + (1 | dog) + (1 | time), 
            data = doglong)
summary(mod)

...
Random effects:
 Groups   Name        Variance  Std.Dev.  
 dog      (Intercept) 48.78     6.985
 time     (Intercept) 2.413e-09 4.912e-5
 Residual             5.292     2.300
Number of obs: 60, groups:  dog, 30; time, 2
...

From the random effect variances, we can work out the ICC! This is because the estimated dog-level variability in intercepts can be taken as an estimate of “true variance”!

48.78 / (48.78 + 2.413e-09 + 5.292)

[1] 0.902

This is the same number that we got for the ICC2 from ICC() above.