---
title: "Dual Change Score Model"
title-slide-attributes:
  data-background-image: lvmworkshop2012swirlbw.png
  data-background-size: contain
  data-background-opacity: "1"
author: "Rich Jones"
date: 2025-09-04
format: 
  revealjs:
    theme: white
    slide-number: true
    logo: logo-lvmworkshop-footer-2024-01-22.png
    self-contained: true
    code-copy: true
    code-line-numbers: true
    css: mplus-output.css
    fig-align: center
filters:
  - default-code-lang.lua
editor: source
---

## Where to learn more 

:::: {.columns}

::: {.column width="50%"}

![](figures/McArdle){height="2.5in"}

::: {style="font-size:0.5em;"}
McArdle, J. (2009). Latent variable modeling of differences and changes with longitudinal data. Annual Review of Psychology, 60, 577-605. 
:::
:::

::: {.column width="50%"}

![](figures/Serang){height="2.5in"}

::: {style="font-size:0.5em;"}
Serang, S., Grimm, K. J., & Zhang, Z. (2019). On the correspondence between the latent growth curve and latent change score models. Structural Equation Modeling: A Multidisciplinary Journal, 26(4), 623-635.
:::

:::

::::


## Dual Change Score Model

![](excalidraw/LVMWKSHP-DCSM-01.excalidraw.svg){height="30%"}

## Shape of the curve driven by $\beta$ 

::: {style="font-size:0.5em;"}
| $\beta$ | Equilibrium | Implication | Typical trajectory |
|-------------|-------------|------------------------|----------------| 
| $\beta>0$   | Exists but unstable| Proportional change, exponential growth or decay away from a equilibrium point ($-s/b$)| Moves away from equilibrium, monotone, no oscillation |
| $\beta=0$ | None (unless $s=0$)| Constant change, linear growth or decay | Straight line|
| $-1<\beta<0$ | Stable $y^* = -s/b$ | Proportional change, exponential monotone decay toward an equilibrium point| Monotone decay toward equilibrium |
| $-2<\beta<-1$ | Stable | Converges to equilibrium, oscillates around ($-s/b$)| Damped oscillation toward equilibrium |
| $\beta<-2$ | Unstable | Diverges with oscillation | Growing flips around equilibrium |

Equilibrium value for y is $y^* = -s/\beta$ and if $y$ ever reached this value the change would be 0 and $y$ would remain at this value. Related to an asymptote. Unstable equilibrium means that if $y$ is perturbed away from this value it will move further away. Stable means that if $y$ is perturbed away from this value it will move back toward it.
:::


## Baseline–Change Association

:::: {.columns}

::: {.column width="25%"}

![](excalidraw/LVMWKSHP-DCSM-03.excalidraw.svg){height="6in"}

:::

::: {.column width="75%"}
::: {style="font-size:0.5em;"}
In a Dual Chnge Score Model (DCSM), the baseline-change association has a dual nature: a between-person component ($\text{Cov}(y_0, s)$) and a within-person component ($\beta$). The random effects covariance $\text{Cov}(y_0, s)$ parameter captures the correlation between the baseline level ($y_0$) and constant change component ($s$). Although they are both random effects, this effect is at the "between-person" level in a multilevel framework. The covariance $\text{Cov}(y_0, s)$ tells us whether individuals who start higher also tend to have faster/slower long-run slopes. This is the between-person piece, describing population heterogeneity.
  
The dynamic feedback parameter ($\beta$) captures the structural regression of a person’s current level on their own subsequent change. It is interpreted as a within-person dynamic: does being higher now predict more or less change next time? Unlike the LGCM, where such dynamics are absorbed into the intercept–slope covariance, the DCSM models them explicitly. The recursion multiplier means that even though $\beta$ is fixed across individuals, the effect of $\beta$ on change varies as a function of the individual’s prior level of y.

In contrast, unconditional linear growth curve model (LGCM), the correlation between baseline and change is captured in one model parameter: the covariance between intercept and slope ($\text{Cov}(i,s)$). This is a between-persons effect.

**_However in practice, the constant change component and the proportional change parameter in the DCSM are often highly correlated, sometimes nearly collinear._** This means the dual nature of baseline–change correlation is not always empirically separable. More time points, smaller proportional change values, and constant change means near zero reduce the problem; few time points and extreme constant change means exacerbate it (Jacobucci R et al. Structural Equation Modeling: A Multidisciplinary Journal. 2019:26(6);924-930).

:::
:::

::::



## LDSM.EXE 

![](Zhang_ldsm_exe.png){height="35%"}

::: {style="font-size:0.5em;"}
LDSM.EXE is a Windows program that generates Mplus code for latent difference score models (LDSM), including dual change score models (DCSM).
::: 

---

::: columns
::: {.column width="50%"}
::: {style="font-size:0.4em;"}
```
!*********************************************!
!This M+ latent difference score model scripts!
!are generated by LDSM generator.             !
!*********************************************!
MODEL:
!For the variable y
!Observed variables and latent level
ly1 by y1 @1;
ly2 by y2 @1;
ly3 by y3 @1;
!Autoregressive part
ly2 on ly1 @1;
ly3 on ly2 @1;
!Difference score on latent level
dy1 by ly2 @1;
dy2 by ly3 @1;
!Auto-proportion of difference score on level
dy1 on ly1 * (1);
dy2 on ly2 * (1);
!Model relationship between slope and ds
sy by dy1 @1;
sy by dy2 @1;
y0 by ly1 @1;
!Set the means and variance to be 0
[y1@0]; [ly1@0]; [dy1@0]; ly1@0; dy1@0;
[y2@0]; [ly2@0]; [dy2@0]; ly2@0; dy2@0;
[y3@0]; [ly3@0];ly3@0;
!Estimate means and variances for intercept and slope
[y0 sy];
y0 sy;
!Set all item-level residuals to be equal
y1* (2);
y2* (2);
y3* (2);
```
:::
:::
::: {.column width="50%"}
![](excalidraw/LVMWKSHP-DCSM-02.excalidraw.svg){height="50%"}
:::
:::

## DCSM of MSQTOT in EPESE

- MSQTOT assessed at baseline, +3 years, and +6 years
- Scores on 0-6 scale, higher score is more correct answers

--- 

{{< include mplus_engine.qmd >}}

::: {style="font-size:0.80em;"}
```{mplus} 
TITLE:    DCSM, MSQ sum scores from EPESE
          Input from LDSM.EXE
          y1, y2, and y3 are MSQTOT at wave 1, 4, and 7

DATA:     FILE = ex0201.dat;

VARIABLE: NAMES = y1 y2 y3 ; ! msqtot1 msqtot4 msqtot7 age black ;
          MISSING  = ALL (-9999) ;

OUTPUT:   TECH1 ;
          TECH4 ;
          
SAVEDATA: SAVE = FSCORES ;
          FILE = ex0399.dat ;

MODEL: 
          !For the variable y
          !Observed variables and latent level
             ly1 by y1 @1;
             ly2 by y2 @1;
             ly3 by y3 @1;
          !Autoregressive part
             ly2 on ly1 @1;
             ly3 on ly2 @1;
          !Difference score on latent level
             dy1 by ly2 @1;
             dy2 by ly3 @1;
          !Auto-proportion of difference score on level
             dy1 on ly1 * (1);
             dy2 on ly2 * (1);
          !Model relationship between slope and ds
             sy by dy1 @1;
             sy by dy2 @1;
             y0 by ly1 @1;
          ! Set the means and variance to be 0
             [y1@0]; [ly1@0]; [dy1@0]; ly1@0; dy1@0;
             [y2@0]; [ly2@0]; [dy2@0]; ly2@0; dy2@0;
             [y3@0]; [ly3@0];ly3@0;
          !Estimate means and variances for intercept and slope
             [y0 sy];
             y0 sy;
          !Set all residuals to be equal
             y1* (2);
             y2* (2);
             y3* (2);
``` 
::: 

::: {style="font-size:0.40em;"}
_Scroll down to see the entire Mplus output_
:::


<!-- copy output file to expected filename --> 
```{r}
#| echo: false
#| output: false
file.copy("formplus.out", "ex0399.out", overwrite = TRUE)
```

# Using R to explore the output and results 

- MplusAutomation package commands to extract information about the model and results
- Load the savedata file and plot the estimated latent trajectories

---

## Load the data in the SAVEDATA file

::: {style="font-size:0.5em;"}
I have a program to do this on my GitHub (rnj0nes/RNJmisc::runmplus_load_savedata.R)

This program will read the OUT file, find the name and contents of the SAVEDATA file, and read that file into R. 

Also `contents.r` is a program to provide a summary of a data frame, and `eme.r` will be discussed in the Appendix.

::: 

::: {style="font-size:0.95em;"}
```{r}
#| echo: true
library(dplyr)
library(tidyr)
library(ggplot2)
# RNJMisc programs
f1 <- "runmplus_load_savedata.R"
f2 <- "contents.r"
f3 <- "eme.r"
rnjmisc_url <- "https://raw.githubusercontent.com/rnj0nes/RNJmisc/master/"
devtools::source_url(paste0(rnjmisc_url, f1))
devtools::source_url(paste0(rnjmisc_url, f2))
devtools::source_url(paste0(rnjmisc_url, f3))
``` 
::: 


---

```{r}
#| echo: true
#| results: markup
# Read SAVEDATA from LDSM
res <- runmplus_load_savedata("ex0399.out")
df <- res$data 
```

--- 

`contents` of `df` (the SAVEDATA file)

::: columns
::: {.column width="75%"}
::: {style="font-size:0.5em;"}
```{r}
#| echo: false
#| results: asis
contents(df)
``` 
:::
:::
::: {.column width="25%"}
::: {style="font-size:0.4em;"}
y1-y3 are the observed MSQTOT scores at waves 1, 4, and 7. Everything else is a latent variable or factor score estimate, or a derivation of a factor score (the standard errors, XXX`_se`). Notice the sample sizes show missing data by death and attrition for the observed variables, but not for the latent variables (which are estimated for all persons).

The `ly` variables are the estimated true scores at each time point. Their `_se` variables are estimates of the standard errors of measurement for these factor scores, and are analogous to the standard deviations of the posterior distribution of the factor scores given the model and the data, where the factor score estimates are the means of these posterior distributions.

The `dy` variables are estimated latent difference scores, `y0` is the estimated intercept, and `sy` is the estimated slope factor score.
:::
:::
:::


--- 

## Plots

::: {style="font-size:0.5em;"}
We will make two plots. The first is a spaghetti plot of the estimated true scores (`ly1`, `ly2`, `ly3`) over time, with the mean curve overlaid. The second plot is a histogram of the person × occasion residuals (the differences between the observed scores and the estimated true scores).

Both of these will require data in "long" format, so we will reshape the data frame.
:::

::: {style="font-size:0.7em;"}
```{r}
#| echo: true
#| results: markup
# Reshape wide to long
# df$ly1, ly2, and ly3 are the expected values 
df_long <- df %>%
   select(y1:y3, ly1, ly2, ly3) %>% 
   mutate(id = row_number()) %>%         # step 1: row identifier
   pivot_longer(
      cols = c(y1:y3, ly1:ly3),          # step 2: reshape
      names_to = c(".value", "obs"),     # split into 'y'/'ly' and 'obs'
      names_pattern = "([a-z]+)([0-9]+)"
   ) %>%
   arrange(id, obs)
head(df_long)
```
::: 

---

::: {style="font-size:0.5em;"}
With the reshaped data we can easily compute correlations and R-squared values.
:::

::: {style="font-size:0.8em;"}
```{r}
#| echo: true
#| results: markup
vars <- c("y", "ly")
cor_matrix <- cor(df_long[vars], use = "pairwise.complete.obs")
print(format(round(cor_matrix, 2), nsmall = 2, trim = TRUE))
```


```{r}
#| echo: true
#| results: markup
r_sq <- cor_matrix["y", "ly"]^2
cat("R-squared:", round(r_sq, 2), "\n")
``` 
::: 

---

::: {style="font-size:0.5em;"}
In the spaghetti plot we will overlay the mean curve from the data (in red) and the expected mean curve implied by the estimated parameters of the model (in blue). We can compute the expected means from the estimated parameters in the Mplus output. This provides a check on my understanding of the model, as these two curves should overlap exactly. In this code I use model parameters copied from the Mplus output, but it would have been better to extract them using MplusAutomation functions. I will demonstrate that in an appendix.
:::


::: {style="font-size:0.8em;"}
```{r}
#| echo: true
#| results: asis
# ensure obs is numeric (pivot_longer often makes it character)
df_long <- df_long %>% mutate(obs = as.integer(obs))

# Extract from Mplus output 
mean_y0    <- 4.782  # mean of "y0" from Mplus output
beta       <- 0.617  # beta from Mplus output (Dy on Ly)
mean_sy    <- -3.103 # mean of "sy" from Mplus output
# Three different ways to compute the implied means all are equal
# implied_y1 <- mean_y0 # mean of "Y0" from Mplus output
# implied_y2 <- implied_y1 + beta*implied_y1 + mean_sy
# another way  
# implied_y3 <- implied_y2 * (1 + beta) + mean_sy
# closed form showing proportional (~exponential) and linear change
implied_y1 <- (1 + beta)^(1-1) * (mean_y0 + mean_sy/beta) - (mean_sy/beta)
implied_y2 <- (1 + beta)^(2-1) * (mean_y0 + mean_sy/beta) - (mean_sy/beta)
implied_y3 <- (1 + beta)^(3-1) * (mean_y0 + mean_sy/beta) - (mean_sy/beta)

# expected means data
expected_df <- tibble::tibble(
   obs = c(1L, 2L, 3L),
   mu  = c(implied_y1, implied_y2 , implied_y3)
)

# mean curve by obs
mean_curve <- df_long %>%
   group_by(obs) %>%
   summarize(ly_bar = mean(ly, na.rm = TRUE), .groups = "drop")
```
::: {style="font-size:0.30em;"}
_Scroll down to see the entire code chunk_
:::

::: 


--- 

::: {style="font-size:0.5em;"}
Code for the spaghetti plot using ggplot2.
:::

```{r}
#| echo: true
#| results: false
spagplot <- ggplot(df_long, aes(x = obs, y = ly, group = id)) +
   geom_line(alpha = 0.01, linewidth = 0.3) + # spaghetti
   geom_point(alpha = 0.01, size = 0.6) +     # optional points
   
   # data mean curve (red)
   geom_line(
      data = mean_curve,
      aes(x = obs, y = ly_bar, group = 1),
      linewidth = 3.1, color = "red", inherit.aes = FALSE
   ) +
   
   # expected mean curve (blue)
   geom_line(
      data = expected_df,
      aes(x = obs, y = mu, group = 1),
      linewidth = 1.1, color = "blue", inherit.aes = FALSE
   ) +
   geom_point(
      data = expected_df,
      aes(x = obs, y = mu),
      size = 0.6 , color = "blue", inherit.aes = FALSE
   ) +
   
   scale_x_continuous(breaks = sort(unique(df_long$obs))) +
   labs(x = "obs", y = "ly", title = "Spaghetti plot of ly by obs") +
   theme_minimal()
``` 
::: {style="font-size:0.40em;"}
_Scroll down to see the entire code chunk_
:::

--- 

```{r}
#| echo: false
#| results: asis
#| fig-width: 6
#| fig-height: 4
#| fig-align: center
print(spagplot)
```

::: {style="font-size:0.5em;"}
Each line is one persons estimated trajectory of the "true score" (ly) over time. The thick red line is the mean of the estimated true scores at each time point. The blue line is the expected mean trajectory implied by the estimated parameters of the model. The fact that these overlap exactly provides a check on my understanding of the model. The observed and model-implied means are identical because there are no covariates and there is no missing data.
:::



--- 

::: {style="font-size:0.5em;"}
Code for the plot of person × occasion residuals.
::: 

```{r}
#| echo: true
#| results: false
df_long <- df_long %>%
   mutate(diff = y - ly)   # new variable: difference

residhist <- ggplot(df_long, aes(x = diff)) +
   geom_histogram(bins = 30, color = "white") +
   labs(x = "Difference (y - ly)", y = "Count",
        title = "Histogram of person × occasion residuals") +
   theme_minimal()
``` 

--- 

```{r}
#| echo: false
#| results: asis
#| fig-width: 6
#| fig-height: 4
print(residhist)
``` 

--- 

## Appendix: Extracting parameters from Mplus output

::: {style="font-size:0.5em;"}
Here we show how to extract the estimated parameters from the Mplus output using MplusAutomation functions. This is better than copying and pasting numbers from the output, as it is less error-prone and can be automated.

Even better would be to run Mplus from R using MplusAutomation functions, and use the H5RESULTS to get parameter estimates, but that is for another time.

Assume we have already run the model and have the output file `ex0399.out`.
::: 


::: {style="font-size:0.8em;"}
```{r}
#| echo: true
library(MplusAutomation)
# read the Mplus output file
out <- MplusAutomation::readModels("ex0399.out")
# extract and view the parameter estimates
params <- out$parameters$unstandardized
head(params)
```
::: 

---

::: {style="font-size:0.5em;"}
Let's just look at the parameters that are really estimates, defined by having a standard error.
::: 

::: {style="font-size:0.8em;"}
```{r}
#| echo: true
# let's at those entries in params with SE > 0
params |> dplyr::filter(se>0) |> head(n=sum(params$se > 0))
```
::: 

::: {style="font-size:0.5em;"}
Pull the mean of y0 into an object in R:
::: 

::: {style="font-size:0.8em;"}
```{r}
#| echo: true
mean_y0 <- params |> dplyr::filter(paramHeader == "Means" & param == "Y0") |> pull(est)
mean_y0
```
::: 


--- 

::: {style="font-size:0.5em;"}
That code is not hard but I do have a helper function to make it easier. This function is in my RNJmisc repo on GitHub, and it's called `eme.r`.
::: 

::: {style="font-size:0.7em;"}
```{r}
#| echo: true
# this creates an object in the R environment called means_y0 (if 
# there is a matching entry in `params`. Note the case insensitivity
# of the function. Syntax is eme(params_object, "paramHeader", "param")
eme(params, "means", "y0")
means_y0$est
means_y0$se
```

```{r}
#| echo: true
#| message: true
# and if your guess at the name is not close you'll get a message
# and some help
eme(params, "y0_with","sy")
```
:::

---

(fin)

