6
$\begingroup$

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?

$\endgroup$
6
$\begingroup$

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
$\endgroup$
2
$\begingroup$

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)
 ```
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$ – Erik Ruzek Jan 3 at 15:02

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.