Psychometrics of Intensive Longitudinal Measures of Emotional States

Overview

This tutorial will walk through the examples provided in Chapter 7 of Bolger & Laurenceau’s Intensive Longitudinal Methods (2013). The book chapter contains MPlus, SPSS, and SAS code, so we will translate those methods into R. The data set used in this tutorial consists of 50 individuals measured on 4 items over 10 time points.

Generalizability Theory

Generalizability Theory (or G Theory) can be used to answer the question, “Can within-subject changes be measured reliably?” To answer this question, we must specify the dimensions of generalizability. In this example, the dimensions are time points, persons, items.

The analysis uses an ANOVA with random effects for each of the specified dimensions.

Preliminaries

Loading libraries used in this script.

library(lme4)      #multilevel models
library(psych)     #psychometrics basic functions
library(tidyverse) #ggplot2, dplyr, tidyr, readr, purr, tibble, stringr, forcats

We use the example data from Bolger & Laurenceau (2013) Chapter 7.

Loading the data

filepath <- "https://raw.githubusercontent.com/The-Change-Lab/collaborations/refs/heads/main/Bolger_Laurenceau_2013/BL2013_data_psychometrics.csv"
psychometrics <- read.csv(filepath)
head(psychometrics, 10)
##    person time item y
## 1     301    1    1 2
## 2     301    1    2 2
## 3     301    1    3 3
## 4     301    1    4 4
## 5     301    2    1 2
## 6     301    2    2 3
## 7     301    2    3 3
## 8     301    2    4 2
## 9     301    3    1 3
## 10    301    3    2 3
psych::describe(psychometrics)
##        vars    n   mean    sd median trimmed   mad min max range  skew kurtosis
## person    1 2000 438.68 58.10  449.0  443.65 65.23 301 521   220 -0.61    -0.51
## time      2 2000   5.50  2.87    5.5    5.50  3.71   1  10     9  0.00    -1.23
## item      3 2000   2.50  1.12    2.5    2.50  1.48   1   4     3  0.00    -1.36
## y         4 1802   2.45  1.07    2.0    2.40  1.48   1   5     4  0.26    -0.69
##          se
## person 1.30
## time   0.06
## item   0.03
## y      0.03

Plotting the data

In this example, a \(m = 4\) item scale was completed on \(K = 10\) occasions by \(N = 50\) persons.

#pivoting to wider
datawide <- psychometrics %>%
  pivot_wider(names_from = item,
              values_from = y)
#column renaming
names(datawide) <- c("person", "time", "item1", "item2", "item3", "item4")

#describe the 4 items
describe(datawide)
##        vars   n   mean    sd median trimmed   mad min max range  skew kurtosis
## person    1 500 438.68 58.14  449.0  443.65 65.23 301 521   220 -0.61    -0.52
## time      2 500   5.50  2.88    5.5    5.50  3.71   1  10     9  0.00    -1.23
## item1     3 450   2.44  0.97    2.0    2.41  1.48   1   5     4  0.20    -0.54
## item2     4 453   2.77  1.14    3.0    2.76  1.48   1   5     4  0.03    -0.75
## item3     5 447   2.27  1.06    2.0    2.19  1.48   1   5     4  0.37    -0.75
## item4     6 452   2.30  1.04    2.0    2.23  1.48   1   5     4  0.33    -0.69
##          se
## person 2.60
## time   0.13
## item1  0.05
## item2  0.05
## item3  0.05
## item4  0.05
#plotting by person
datawide %>%
  ggplot(aes(x=time,y=item1, group=person)) +
  geom_line(aes(x=time,y=item1), color="blue") +
  geom_line(aes(x=time,y=item2), color="red") +
  geom_line(aes(x=time,y=item3), color="green") +
  geom_line(aes(x=time,y=item4), color="orange") +
  xlab("Time") +
  ylab("Rating") +
  scale_x_continuous(breaks=seq(1,10,by=1)) +
  theme_bw() +
  facet_wrap( ~ person)
## Warning: Removed 16 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 17 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_line()`).

The Generalizability Theory Model (Variance Decomposition)

To invoke the G-theory perspective we specify a three-way, crossed random effects ANOVA model (time * person * item). For person j, responding at time point i to item k, the model for the observations \(Y_{ijk}\) is

\[Y_{ijk} + \mu +T_{i} + P_{j} + I_{k} + (TP)_{ij} + (TI)_{jk} + (PI)_{jk} + (TPI)_{ijk} + e_{ijk}\]

In the case of completely balanced data, where each person provides scores for each item for each time point, the variances associated with each random effect sum to the total variance. Thus the model provides for a complete variance decomposition.

\[\sigma^2_{Y_{ijk}} = \sigma^2_{T_{i}} + \sigma^2_{P_{j}} + \sigma^2_{I_{k}} + \sigma^2_{(TP)_{ij}} + \sigma^2_{(TI)_{ik}} + \sigma^2_{(PI)_{jk}} + [\sigma^2_{(TPI)_{ijk}} + \sigma^2_{e_{ijk}}]\]

Note that the final two variances are bracketed to indicate that they cannot be separated.

Implementation of Variance Decomposition Model

The lme4 package contains the function lmer() which fits linear mixed-effects models. It is used here to engage the variance decomposition. Basically, this is an anova.

The summary() of our model will show us the random effects (i.e.,variances) associated with each dimension of generalizability.

The model is specified as an intercept-only model with random effects for all the dimensions of generalizability.

Restructuring the variables so that they are explicitly categorical variables using the factor() function.

psychometrics <- psychometrics %>%
  mutate(person = factor(person),
         time = factor(time),
         item = factor(item))

Specification of the ANOVA with random effects model in lmer() is …

#specifying model
model1 <- lmer(y ~  1 + 
                 (1|person) + 
                 (1|time) + 
                 (1|item) + 
                 (1|person:time) + 
                 (1|person:item) + 
                 (1|time:item), 
               data=psychometrics)
## boundary (singular) fit: see help('isSingular')
#model summary
summary(model1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ 1 + (1 | person) + (1 | time) + (1 | item) + (1 | person:time) +  
##     (1 | person:item) + (1 | time:item)
##    Data: psychometrics
## 
## REML criterion at convergence: 4046.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2721 -0.5311 -0.0105  0.5044  3.8614 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  person:time (Intercept) 0.255272 0.50524 
##  person:item (Intercept) 0.190111 0.43602 
##  person      (Intercept) 0.361774 0.60148 
##  time:item   (Intercept) 0.004983 0.07059 
##  time        (Intercept) 0.000000 0.00000 
##  item        (Intercept) 0.048596 0.22044 
##  Residual                0.299542 0.54730 
## Number of obs: 1802, groups:  
## person:time, 455; person:item, 200; person, 50; time:item, 40; time, 10; item, 4
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   2.4340     0.1456   16.71
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

Using the VarCorr() function, we can extract and save each variance value from the summary table.

#looking at VarCorr Object
VarCorr(model1)
##  Groups      Name        Std.Dev.
##  person:time (Intercept) 0.505244
##  person:item (Intercept) 0.436017
##  person      (Intercept) 0.601477
##  time:item   (Intercept) 0.070587
##  time        (Intercept) 0.000000
##  item        (Intercept) 0.220444
##  Residual                0.547304

We then extract each of the standard deviations, square them into variances and place into separate objects for easy use in later calculations.

#extracting and naming each component (be careful with naming)
(personTime <- VarCorr(model1)[[1]][1,1]) #person:time 
## [1] 0.2552719
(personItem <- VarCorr(model1)[[2]][1,1]) #person:item 
## [1] 0.1901106
(person <- VarCorr(model1)[[3]][1,1]) #person 
## [1] 0.3617743
(timeItem <- VarCorr(model1)[[4]][1,1]) #time:item 
## [1] 0.004982546
(time <- VarCorr(model1)[[5]][1,1]) #time 
## [1] 0
(item <- VarCorr(model1)[[6]][1,1]) #item 
## [1] 0.04859551
(residual <- sigma(model1)^2) #residual
## [1] 0.2995419

Great! - We can now use these components to calcualte a variety of different generalizability coefficients (i.e., reliabilities).

Generalizability Coefficients

Using the variance decomposition, we can examine the reliability of multi-item repeated measures scales. In this example, a \(m = 4\) item scale was completed on \(K = 10\) occasions by \(N = 50\) persons.

Reliability of Within-person Change

Back to our initial question: Are there reliable within-person differences in change over time?

We use the following formula to calculate the reliability coefficient: \[Rc = \frac{\sigma^2_{person:time}}{\left(\sigma^2_{person:time} + \frac{(\sigma^2_{person:time:item} + \sigma^2_{e})}{m}\right)}\]

where \(m\) refers to the number of items. In our case, \(m = 4\).

We cannot separate the \(\sigma^2_{person:time:item}\) from the \(\sigma^2_{e}\). Both are included in the residual error term, \(\sigma^2_{error}\).

#setting k
m <- 4

#Calculating Rc
(Rc <- personTime/(personTime + residual/m))
## [1] 0.7731825

This coefficient represents the degree to which these repeated measures are adequate and systematic. Using the same interpretation rules of thumb as Cronbach’s alpha, we can determine that four items can capture within person change reliably, \(R_{c} = 0.77\).

Reliability of Between-person Differences

Another reliability coefficient of interest is \(R_{1F}\), which is calculated as

\[R_{1F} = \frac{\sigma^2_{person} + [\frac{\sigma^2_{person:item}}{m}]}{\sigma^2_{person} + [\frac{\sigma^2_{person:item}}{m}] + [\frac{\sigma^2_{error}}{m}]}\]

where \(m\) refers to the number of items. In our case, \(m = 4\).

Our calculation is

m <- 4
(R1f <- (person + (personItem/m))/(person + (personItem/m) + (residual/m)))
## [1] 0.8453379

The reliability coefficient \(R_{1F}\) is the expected between-person reliability estimate for one fixed day, a kind of average of day-specific Cronbach’s alphas across occasions. Using the same interpretation rules of thumb as Cronbach’s alpha, we can determine that four items can capture between-person differences on any given day reliably, \(R_{1F} = 0.85\).

Reliability of Between-person Differences (Across Entire Set of Repeated Measures)

Another reliability coefficient of interest is \(R_{KF}\), which is calculated as

\[R_{KR} = \frac{\sigma^2_{person} + [\frac{\sigma^2_{person:item}}{m}]}{\sigma^2_{person} + [\frac{\sigma^2_{person:item}}{m}] + [\frac{\sigma^2_{error}}{Km}]}\] where \(K\) refers to the number of occasions. and \(m\) refers to the number of items. In our case, \(K = 10\) and \(m = 4\).

Our calculation is

K <- 10
m <- 4
(RKR <- (person + (personItem/m))/(person + (personItem/m) + (residual/(K*m))))
## [1] 0.9820328

The reliability coefficient \(R_{KRF}\) is the expected between-person reliability estimate averaging across all days and items. Using the same interpretation rules of thumb as Cronbach’s alpha, we can determine that four items can capture between-person differences across the 10 repeated measures reliably, \(R_{KR} = 0.98\).

Conclusion Notes

Expanding Classical Test Theory notions of reliability, the Generalizability Theory approach provides opportunity to consider the extent to which measurement generalizes across many different facets (including all 3-dimensions of the Data Box!).

Note: The book chapter also contains MPlus code to calculate the same metrics. To see this example in SAS, MPlus, SPSS, and another R example, please visit [http://www.intensivelongitudinal.com/ch7/ch7index.html]. The lavaan package in R can be used to do the structural equation modeling, and some basic forms of multilevel structural equation modeling.

Thank you for playing!!