# Conditional intraclass correlation (ICC) from a linear mixed model as evidence for test-retest reliability?

In my experiment with two conditions (between-subjects design), participants completed a single-item scale three times: (1) before the experimental manipulation, (2) after the experimental manipulation, and (3) at the end of the study. I want to test the test-retest reliability of my single-item scale over three measurements.

The issue here is that the experimental manipulation is known to affect the response in the measure. I'm assuming that I need to control for any effect of the manipulation, so my guess is to use a linear mixed model (hierarchical linear model) with measurements nested within participants and include an experimental condition as a factor. This way, I can get an intraclass correlation (ICC) from this conditional model.

Would it be possible to interpret this conditional ICC as "a measure of test-retest reliability after controlling for the experimental effect"? Are there any existing studies using a similar approach?

Yes, you can do this and interpret it as you think. I have read about such an interpretation in the second chapter of Sophia Rabe-Hesketh and Anders Skrondal's Multilevel and Longitudinal Modeling using Stata book (Volume 1).

A more detailed explanation follows. Edit: I also added a simulation to demonstrate what is going on. Hat tip to Ariel Muldoon for a helpful blog post that aided me in creating this simulation.

In a random intercept model with no predictors, $$y_{ij} = \beta_0 + u_{0j} + \epsilon_{ij}$$ we get two variances, one for $$u_{0j}$$, which is $$\psi$$, and one for $$\epsilon_{ij}$$, which is $$\theta$$.

From these we can express between-subject dependence or reliability ($$\rho$$) as: $$\rho = \frac{\psi}{\psi+\theta}$$

In this equation, $$\psi$$ is the variance of subjects' true scores $$\beta_0 + u_{0j}$$ and $$\theta$$ is the measurement error variance, or squared standard error of measurement. $$\rho$$ becomes a test-retest reliability because of the repeated measurements.

In contrast to the Pearson correlation coefficient, $$\rho$$ is influenced by any linear transformations of measurements, which could include practice effects or experimentally-induced increases from time 1 to time 2. Thus, if you know of something in your data that induces linear changes, you must account for it in your mixed model.

In your case, you have a time-varying experimental manipulation (call it $$x_1$$). Including $$x_1$$ as a predictor in your random intercept model,

$$y_{ij} = \beta_0 + \beta_1x_1 +u_{0j} + \epsilon_{ij}$$

will (likely) have an effect on both $$\psi$$ and $$\theta$$. In so doing, the resulting estimates of $$\psi$$ and $$\theta$$ are no longer influenced by $$x_1$$, and you have an estimate of test-retest reliability robust to experimental effects.

Simulation

set.seed(807)

npart=1000 # number of particpants
ntime=3    # numer of observations (timepoints) per participant
mu=2.5     # mean value on the Likert item
sdp=1      # standard deviation of participant random effect (variance==1)
sd=.7071   # standard deviation of within participant (residual; variance = .5)

participant = rep(rep(1:npart, each = nobs),ntime)  # creating 1000 participants w/ 3 repeats
participant = participant[order(participant)]
time = rep(rep(1:ntime, each=1),1000)    # creating a time variable

parteff = rnorm(npart, 0, sdp)     # drawing from normal for participant deviation
parteff = rep(parteff, each=ntime) # ensuring participant effect is same for three observations

timeeff = rnorm(npart*ntime, 0, sd) # drawing from normal for within-participant residual

dat=data.frame(participant, time, parteff, timeeff) # create data frame

dat$resp = with(dat, mu + parteff + timeeff ) # creating response for each individual #Variance components model library(lme4) m1 <- lmer(resp ~ 1 + (1|participant), dat) summary(m1) # estimates close to simulated values Linear mixed model fit by REML ['lmerMod'] Formula: resp ~ 1 + (1 | participant) Data: dat REML criterion at convergence: 8523.8 Scaled residuals: Min 1Q Median 3Q Max -3.13381 -0.57238 0.01722 0.57846 2.84918 Random effects: Groups Name Variance Std.Dev. participant (Intercept) 1.0110 1.0055 Residual 0.5314 0.7289 Number of obs: 3000, groups: participant, 1000 Fixed effects: Estimate Std. Error t value (Intercept) 2.54142 0.03447 73.73 #Add treatment variable x1 which turns on at time 3 dat$trtmt = rep(c(0,0,1),1000)
b1 = .4 #average amount by which particpant's score increases b/c of treatment
x1 = runif(npart, .05, 1.5)

library(dplyr)
dat <- dat %>% mutate(resp2=case_when
(time==3 ~ (mu+b1*x1+parteff+timeeff),
TRUE ~ resp))
glimpse(dat)

#run m1 without covariate for trtmt
m2 <- lmer(resp2 ~ 1 + (1|participant), dat)
summary(m2)

Linear mixed model fit by REML ['lmerMod']
Formula: resp2 ~ 1 + (1 | participant)
Data: dat

REML criterion at convergence: 8659.9

Scaled residuals:
Min       1Q   Median       3Q      Max
-2.72238 -0.56861  0.01894  0.57177  3.10610

Random effects:
Groups      Name        Variance Std.Dev.
participant (Intercept) 1.0070   1.0035
Residual                0.5669   0.7529
Number of obs: 3000, groups:  participant, 1000

Fixed effects:
Estimate Std. Error t value
(Intercept)  2.64169    0.03458   76.39

#add trtmt as a fixed effect predictor
m3 <- lmer(resp2 ~ 1 + trtmt + (1|participant), dat)
summary(m3)

Linear mixed model fit by REML ['lmerMod']
Formula: resp2 ~ 1 + trtmt + (1 | participant)
Data: dat

REML criterion at convergence: 8546.7

Scaled residuals:
Min       1Q   Median       3Q      Max
-3.06878 -0.57650  0.02712  0.57887  2.89709

Random effects:
Groups      Name        Variance Std.Dev.
participant (Intercept) 1.0178   1.0088
Residual                0.5346   0.7311
Number of obs: 3000, groups:  participant, 1000

Fixed effects:
Estimate Std. Error t value
(Intercept)  2.53746    0.03585   70.78
trtmt        0.31270    0.02832   11.04

Correlation of Fixed Effects:
(Intr)
trtmt -0.263

> texreg::screenreg(c(m1, m2, m3))

======================================================================
Model 1       Model 2       Model 3
----------------------------------------------------------------------
(Intercept)                       2.54 ***      2.64 ***      2.54 ***
(0.03)        (0.03)        (0.04)
trtmt                                                         0.31 ***
(0.03)
----------------------------------------------------------------------
AIC                            8529.83       8665.86       8554.72
BIC                            8547.85       8683.88       8578.75
Log Likelihood                -4261.92      -4329.93      -4273.36
Num. obs.                      3000          3000          3000
Num. groups: participant       1000          1000          1000
Var: participant (Intercept)      1.01          1.01          1.02
Var: Residual                     0.53          0.57          0.53
======================================================================
*** p < 0.001; ** p < 0.01; * p < 0.05


This post really helped me and I wanted to thank you. In case other users ran to the same issue I had - I am adding a slight change to the simulation above. The only thing here is that this shows that Pearson corr for two times measurements is exactly the same as $$\rho$$. Nothing special - only nice to see the numbers match :) Also, ever so slight correction in the participant vector to make this work.

Cheers

Nitzan

set.seed(807)

npart=1000 # number of particpants
ntime=2   # numer of observations (timepoints) per participant
mu=2.5     # mean value on the Likert item
sdp=1      # standard deviation of participant random effect (variance==1)
sd=.7071   # standard deviation of within participant (residual; variance = .5)

participant = rep(rep(1:npart, each = nobs),ntime)  # creating 1000 participants w/ 3 repeats
participant = participant[order(participant)]
time        = rep(rep(1:ntime, each=1),1000)        # creating a time variable

parteff = rnorm(npart, 0, sdp)     # drawing from normal for participant deviation
parteff = rep(parteff, each=ntime) # ensuring participant effect is same for three observations

timeeff = rnorm(npart*ntime, 0, sd) # drawing from normal for within-participant residual

dat=data.frame(participant, time, parteff, timeeff) # create data frame

dat\$resp = with(dat, mu + parteff + timeeff ) # creating response for each individual

#Variance components model
library(lme4)

m1 <- lmer(resp ~ 1 + (1|participant), dat)
summary(m1) # estimates close to simulated values

#calculate pearson corr
library(reshape2)
df.wide   <-dcast(dat,participant~time,mean,value.var='resp')[,-1]
cor(df.wide)

#get the same from the HLM fit
print(VarCorr(m1))
.95478^2/(.95478^2+0.74685^2)
$$$$
`
• Thank you, @nitzan shahar! (+1) A good catch on the code for creating participants. I updated in the code in my post accordingly. And yes, nice to see the confirmation from the Pearson correlation once you reshape the data to wide. – Erik Ruzek Jan 3 at 15:02