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
