# Computing Cronbach's Alpha and Intraclass Correlation coefficient from an lmer model

I am trying to calculate the Intraclass Correlation for a rater study using R and the library lme4 and the function lmer. The data has the following design: The same 6 raters (at least 4) are rating 25 horses live and all raters are rating a subset of 10 horses on video.

The two way random model applied:

m1 <- lmer(Score ~ -1 + (1|HorseID) + (1|RaterID) + Time, data=mydata)


The absolute agreement ICC is calculated using the estimated coefficients:

xVars <- function(model) {
exvars = lme4::VarCorr(model)
vars = c(exvars$HorseID[1,1], exvars$RaterID[1,1], attr(exvars,"sc")^2)
names(vars) <- c('item var', 'judge var', 'residual var')
vars }

# helper function for ICC(k) variations

icck <- function(variances, k=1) {
icc = variances[1] / (variances[1] + (variances[2] + variances[3]) / k)
names(icc) = c(paste('ICC', k, sep=''))
icc }


And the ICC as:

ICC.m1 <- icck(xVars(m1))


• Confidence Intervals for the ICC
• Cronbach's alpha

But I can't figure out at smart way of doing so? Help would be greatly appreciated!!

• You don't need a mixed-effect model to compute Cronbach alpha, only the variance of the total score and that of individual items. However, the average measures ICC for the 2-way mixed model (ICC(2,k)) should be close enough. Re. ICC, you can use the bootstrap or asymptotic formula, e.g., Comparison of confidence interval methods for an intra-class correlation coefficient (ICC). (but there are many other approaches that were suggested in the past, especially for rater reliability).
– chl
Commented Oct 25, 2020 at 9:32

### Setup

Actually you can now get the ICC fairly easy from the model using the performance package. You can also get the alpha coefficient quickly with the psych package. Here below I load the requisite libraries for my example:

#### Load Libraries ####
library(lmerTest)
library(performance)
library(psych)
library(tidyverse)


The dataset I will use for the ICC example comes from cake dataset of the nlme package, which measures breakage angles from different cakes. Since this data doesn't have multiple-item scales from what I recall, I've used the BFI dataset from the psych package separately for that example. You can inspect the data below:

#### Load Datasets ####
cake <- lme4::cake
head(bfi) # lots of items for alpha use


### ICC Example

From there I can fit a mixed model with the cake dataset and get the ICC from the saved model:

#### Run Model ####
model <- lmer(angle ~ recipe * temperature
+ (1|recipe:replicate),
cake,
REML= FALSE)

#### Run ICC ####
icc(model)


Which gives me this quick output:

# Intraclass Correlation Coefficient



### Reliability Example

You don't need to fit a mixed model to get the reliability coefficients and the previous example didn't have multiple item scales from what I recall, so from here I use the BFI data. Since I'm not interested in clumping all the items together, I will only select the A items from the BFI dataset, then run the reliability function on it directly:

#### Check Only A Items ####
a.items <- bfi %>%
select(contains("A"),
-age,
-education)

#### Run Reliability ####
reliability(a.items) # reliability measures


From there I get a host of reliability metrics, but you can see alpha listed among them:

Measures of reliability
reliability(keys = a.items)
omega_h alpha omega.tot  Uni r.fit fa.fit max.split
All_items    0.64  0.71      0.75 0.89   0.9   0.99      0.75
min.split mean.r med.r n.items
All_items      0.62   0.33  0.34       5


For the alpha CI, you can also directly use the alpha function, which only brings back alpha metrics. I specify it explicitly here because sometimes it gives an error otherwise with specific packages loaded:

a.items %>%
psych::alpha()


Which gives me what I want:

Some items ( A1 ) were negatively correlated with the total scale and
probably should be reversed.
To do this, run the function again with the 'check.keys=TRUE' option
Reliability analysis
Call: psych::alpha(x = .)

raw_alpha std.alpha G6(smc) average_r  S/N   ase mean   sd
0.43      0.46    0.53      0.15 0.85 0.016  4.2 0.74
median_r
0.32

95% confidence boundaries
lower alpha upper
Feldt      0.4  0.43  0.46
Duhachek   0.4  0.43  0.46

Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r  S/N alpha se  var.r
A1      0.72      0.73    0.67     0.398 2.64   0.0087 0.0065
A2      0.28      0.30    0.39     0.097 0.43   0.0219 0.1098
A3      0.18      0.21    0.31     0.061 0.26   0.0249 0.1015
A4      0.25      0.31    0.44     0.099 0.44   0.0229 0.1607
A5      0.21      0.24    0.36     0.072 0.31   0.0238 0.1311
med.r
A1 0.376
A2 0.081
A3 0.081
A4 0.105
A5 0.095

Item statistics
n raw.r std.r r.cor r.drop mean  sd
A1 2784 0.066 0.024 -0.39  -0.31  2.4 1.4
A2 2773 0.630 0.666  0.58   0.37  4.8 1.2
A3 2774 0.724 0.742  0.72   0.48  4.6 1.3
A4 2781 0.686 0.661  0.50   0.37  4.7 1.5
A5 2784 0.700 0.719  0.64   0.45  4.6 1.3

Non missing response frequency for each item
1    2    3    4    5    6 miss
A1 0.33 0.29 0.14 0.12 0.08 0.03 0.01
A2 0.02 0.05 0.05 0.20 0.37 0.31 0.01
A3 0.03 0.06 0.07 0.20 0.36 0.27 0.01
A4 0.05 0.08 0.07 0.16 0.24 0.41 0.01
A5 0.02 0.07 0.09 0.22 0.35 0.25 0.01