My last post demonstrated a dual change model for one variable, now I want to demonstrate a bivariate dual change model. A SEM path diagram for a bivariate dual change model is below, taken from Wang, Zhou, and Zhang (2016)
(If you do not have access to that link you can view a similar path diagram in Jones, King, Gilrane, McCausland, Cortina, & Grimm, 2016)
Essentially, we have two dual change processes and a coupling parameter from the latent true score on one variable to the latent change score on the other.
The DGP
\[\begin{equation} y_t = constant_y + (1 + proportion_y)*y_{t-1} + coupling_{xy}*x_{t-1} + e \end{equation}\]
where \(constant_y\) is the change factor (or latent slope) on \(y\), \(proportion_y\) is the proportion change factor, and \(coupling_xy\) is the coupling parameter relating \(x\) to \(y\). The DGP for \(x\) is
\[\begin{equation} x_t = constant_x + (1 + proportion_x)*x_{t-1} + coupling_{yx}*y_{t-1} + e \end{equation}\]
where the terms are similar but now applied to values of \(x\). The true values used in the DGP are:
\[\begin{align} y_t &= 0.5 + (1 + -0.32)y_{t-1} + 0.4x_{t-1} + e \\ x_t &= 0.5 + (1 + 0.22)x_{t-1} - 0.4y_{t-1} + e \end{align}\]
with initial values for both \(x\) and \(y\) sampled from \(N\) ~ (10, 1).
people <- 700
time <- 6
x_cause_y <- 0.4
y_cause_x <- -0.4
const_x <- 0.5
const_y <- 0.5
prop_x <- 0.22
prop_y <- -0.32
df_mat <- matrix(, ncol = 4, nrow = people*time)
count <- 0
for(i in 1:people){
unob_het_y <- rnorm(1, 0, 3)
unob_het_x <- rnorm(1, 0, 3)
for(j in 1:time){
count <- count + 1
if(j == 1){
df_mat[count, 1] <- i
df_mat[count, 2] <- j
df_mat[count, 3] <- rnorm(1, 10, 1)
df_mat[count, 4] <- rnorm(1, 10, 1)
}else{
df_mat[count, 1] <- i
df_mat[count, 2] <- j
df_mat[count, 3] <- const_x + (1+prop_x)*df_mat[count - 1, 3] + y_cause_x*df_mat[count - 1, 4] + unob_het_x + rnorm(1,0,1)
df_mat[count, 4] <- const_y + (1+prop_y)*df_mat[count - 1, 4] + x_cause_y*df_mat[count - 1, 3] + unob_het_y + rnorm(1,0,1)
}
}
}
library(tidyverse)
library(ggplot2)
library(reshape2)
df <- data.frame(df_mat)
names(df) <- c('id', 'time', 'x', 'y')
Values of \(y\) over time.
random_nums <- sample(c(1:700), 6)
df_sample <- df %>%
filter(id %in% random_nums)
ggplot(df, aes(x = time, y = y, group = id)) +
geom_point(color = 'grey85') +
geom_line(color = 'grey85') +
geom_point(data = df_sample, aes(x = time, y = y, group = id)) +
geom_line(data = df_sample, aes(x = time, y = y, group = id))
Values of \(x\) over time.
plot_single_response <- function(y_axis){
plot_it <- ggplot(df, aes(x = time, y = !!y_axis, group = id)) +
geom_point(color = 'grey85') +
geom_line(color = 'grey85') +
geom_point(data = df_sample, aes(x = time, y = !!y_axis, group = id)) +
geom_line(data = df_sample, aes(x = time, y = !!y_axis, group = id))
return(plot_it)
}
plot_single_response(quo(x))
Three randomly selected individuals with \(x\) and \(y\) plotted simultaneously.
three_cases <- df %>%
filter(id == 4 | id == 500 | id == 322) %>%
gather(x, y, key = 'variable', value = 'response')
ggplot(three_cases, aes(x = time, y = response, color = variable)) +
geom_point() +
geom_line() +
facet_wrap(~id)
Dual Change Model on Y
df_wide_y <- df %>%
select(id, time, y) %>%
reshape(idvar = 'id', timevar = 'time', direction = 'wide')
library(lavaan)
dual_change_y_string <- [1672 chars quoted with ''']
dc_y_model <- sem(dual_change_y_string, data = df_wide_y)
summary(dc_y_model, fit.measures = T)
lavaan 0.6-9 ended normally after 67 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 16
Number of equality constraints 9
Number of observations 700
Model Test User Model:
Test statistic 6082.258
Degrees of freedom 20
P-value (Chi-square) 0.000
Model Test Baseline Model:
Test statistic 7974.233
Degrees of freedom 15
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.238
Tucker-Lewis Index (TLI) 0.429
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -11526.640
Loglikelihood unrestricted model (H1) -8485.511
Akaike (AIC) 23067.280
Bayesian (BIC) 23099.138
Sample-size adjusted Bayesian (BIC) 23076.911
Root Mean Square Error of Approximation:
RMSEA 0.658
90 Percent confidence interval - lower 0.644
90 Percent confidence interval - upper 0.672
P-value RMSEA <= 0.05 0.000
Standardized Root Mean Square Residual:
SRMR 2.366
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Latent Variables:
Estimate Std.Err z-value P(>|z|)
ly1 =~
y.1 1.000
ly2 =~
y.2 1.000
ly3 =~
y.3 1.000
ly4 =~
y.4 1.000
ly5 =~
y.5 1.000
ly6 =~
y.6 1.000
cy2 =~
ly2 1.000
cy3 =~
ly3 1.000
cy4 =~
ly4 1.000
cy5 =~
ly5 1.000
cy6 =~
ly6 1.000
l_intercept =~
ly1 1.000
l_slope =~
cy2 1.000
cy3 1.000
cy4 1.000
cy5 1.000
cy6 1.000
Regressions:
Estimate Std.Err z-value P(>|z|)
ly2 ~
ly1 1.000
ly3 ~
ly2 1.000
ly4 ~
ly3 1.000
ly5 ~
ly4 1.000
ly6 ~
ly5 1.000
cy2 ~
ly1 (prop) 0.367 0.020 18.281 0.000
cy3 ~
ly2 (prop) 0.367 0.020 18.281 0.000
cy4 ~
ly3 (prop) 0.367 0.020 18.281 0.000
cy5 ~
ly4 (prop) 0.367 0.020 18.281 0.000
cy6 ~
ly5 (prop) 0.367 0.020 18.281 0.000
Covariances:
Estimate Std.Err z-value P(>|z|)
l_intercept ~~
l_slope -2.346 0.250 -9.399 0.000
.cy2 ~~
.cy3 0.000
.cy4 0.000
.cy5 0.000
.cy6 0.000
.cy3 ~~
.cy4 0.000
.cy5 0.000
.cy6 0.000
.cy4 ~~
.cy5 0.000
.cy6 0.000
.cy5 ~~
.cy6 0.000
Intercepts:
Estimate Std.Err z-value P(>|z|)
l_slope -4.638 0.214 -21.714 0.000
l_intercept 11.437 0.116 98.207 0.000
.ly1 0.000
.ly2 0.000
.ly3 0.000
.ly4 0.000
.ly5 0.000
.ly6 0.000
.cy2 0.000
.cy3 0.000
.cy4 0.000
.cy5 0.000
.cy6 0.000
.y.1 0.000
.y.2 0.000
.y.3 0.000
.y.4 0.000
.y.5 0.000
.y.6 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
l_ntrcp 6.704 0.496 13.517 0.000
l_slope 1.953 0.144 13.533 0.000
.ly1 0.000
.ly2 0.000
.ly3 0.000
.ly4 0.000
.ly5 0.000
.ly6 0.000
.cy2 0.000
.cy3 0.000
.cy4 0.000
.cy5 0.000
.cy6 0.000
.y.1 (rs_v) 6.341 0.169 37.417 0.000
.y.2 (rs_v) 6.341 0.169 37.417 0.000
.y.3 (rs_v) 6.341 0.169 37.417 0.000
.y.4 (rs_v) 6.341 0.169 37.417 0.000
.y.5 (rs_v) 6.341 0.169 37.417 0.000
.y.6 (rs_v) 6.341 0.169 37.417 0.000Code to change the \(y\)’s in the string to \(x\)’s without manually deleting and inserting \(x\) into the string above. All you have to do is paste the string into a .txt document and save the file as “y_file.txt”
Dual Change Model on X
df_wide_x <- df %>%
select(id, time, x) %>%
reshape(idvar = 'id', timevar = 'time', direction = 'wide')
library(lavaan)
dual_change_x_string <- [1672 chars quoted with ''']
dc_x_model <- sem(dual_change_x_string, data = df_wide_x)
summary(dc_x_model, fit.measures = T)
lavaan 0.6-9 ended normally after 72 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 16
Number of equality constraints 9
Number of observations 700
Model Test User Model:
Test statistic 4093.422
Degrees of freedom 20
P-value (Chi-square) 0.000
Model Test Baseline Model:
Test statistic 11625.807
Degrees of freedom 15
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.649
Tucker-Lewis Index (TLI) 0.737
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -10086.040
Loglikelihood unrestricted model (H1) -8039.329
Akaike (AIC) 20186.080
Bayesian (BIC) 20217.937
Sample-size adjusted Bayesian (BIC) 20195.711
Root Mean Square Error of Approximation:
RMSEA 0.539
90 Percent confidence interval - lower 0.526
90 Percent confidence interval - upper 0.553
P-value RMSEA <= 0.05 0.000
Standardized Root Mean Square Residual:
SRMR 0.806
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Latent Variables:
Estimate Std.Err z-value P(>|z|)
lx1 =~
x.1 1.000
lx2 =~
x.2 1.000
lx3 =~
x.3 1.000
lx4 =~
x.4 1.000
lx5 =~
x.5 1.000
lx6 =~
x.6 1.000
cx2 =~
lx2 1.000
cx3 =~
lx3 1.000
cx4 =~
lx4 1.000
cx5 =~
lx5 1.000
cx6 =~
lx6 1.000
l_intercept =~
lx1 1.000
l_slope =~
cx2 1.000
cx3 1.000
cx4 1.000
cx5 1.000
cx6 1.000
Regressions:
Estimate Std.Err z-value P(>|z|)
lx2 ~
lx1 1.000
lx3 ~
lx2 1.000
lx4 ~
lx3 1.000
lx5 ~
lx4 1.000
lx6 ~
lx5 1.000
cx2 ~
lx1 (prop) 0.145 0.004 34.642 0.000
cx3 ~
lx2 (prop) 0.145 0.004 34.642 0.000
cx4 ~
lx3 (prop) 0.145 0.004 34.642 0.000
cx5 ~
lx4 (prop) 0.145 0.004 34.642 0.000
cx6 ~
lx5 (prop) 0.145 0.004 34.642 0.000
Covariances:
Estimate Std.Err z-value P(>|z|)
l_intercept ~~
l_slope -1.043 0.250 -4.177 0.000
.cx2 ~~
.cx3 0.000
.cx4 0.000
.cx5 0.000
.cx6 0.000
.cx3 ~~
.cx4 0.000
.cx5 0.000
.cx6 0.000
.cx4 ~~
.cx5 0.000
.cx6 0.000
.cx5 ~~
.cx6 0.000
Intercepts:
Estimate Std.Err z-value P(>|z|)
l_slope -3.559 0.123 -29.026 0.000
l_intercept 10.347 0.075 138.246 0.000
.lx1 0.000
.lx2 0.000
.lx3 0.000
.lx4 0.000
.lx5 0.000
.lx6 0.000
.cx2 0.000
.cx3 0.000
.cx4 0.000
.cx5 0.000
.cx6 0.000
.x.1 0.000
.x.2 0.000
.x.3 0.000
.x.4 0.000
.x.5 0.000
.x.6 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
l_ntrcp 2.752 0.201 13.690 0.000
l_slope 9.980 0.570 17.502 0.000
.lx1 0.000
.lx2 0.000
.lx3 0.000
.lx4 0.000
.lx5 0.000
.lx6 0.000
.cx2 0.000
.cx3 0.000
.cx4 0.000
.cx5 0.000
.cx6 0.000
.x.1 (rs_v) 2.105 0.056 37.417 0.000
.x.2 (rs_v) 2.105 0.056 37.417 0.000
.x.3 (rs_v) 2.105 0.056 37.417 0.000
.x.4 (rs_v) 2.105 0.056 37.417 0.000
.x.5 (rs_v) 2.105 0.056 37.417 0.000
.x.6 (rs_v) 2.105 0.056 37.417 0.000Bivariate Dual Change Model
bi_dc_string <- [3578 chars quoted with ''']
df_both <- df %>%
reshape(idvar = 'id', timevar = 'time', direction = 'wide')
bi_dc_model <- sem(bi_dc_string, data = df_both)
summary(bi_dc_model, fit.measures = T)
lavaan 0.6-9 ended normally after 115 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 46
Number of equality constraints 26
Number of observations 700
Model Test User Model:
Test statistic 1593.425
Degrees of freedom 70
P-value (Chi-square) 0.000
Model Test Baseline Model:
Test statistic 23686.764
Degrees of freedom 66
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.936
Tucker-Lewis Index (TLI) 0.939
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -15278.191
Loglikelihood unrestricted model (H1) -14481.478
Akaike (AIC) 30596.382
Bayesian (BIC) 30687.404
Sample-size adjusted Bayesian (BIC) 30623.900
Root Mean Square Error of Approximation:
RMSEA 0.176
90 Percent confidence interval - lower 0.169
90 Percent confidence interval - upper 0.184
P-value RMSEA <= 0.05 0.000
Standardized Root Mean Square Residual:
SRMR 0.097
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Latent Variables:
Estimate Std.Err z-value P(>|z|)
ly1 =~
y.1 1.000
ly2 =~
y.2 1.000
ly3 =~
y.3 1.000
ly4 =~
y.4 1.000
ly5 =~
y.5 1.000
ly6 =~
y.6 1.000
cy2 =~
ly2 1.000
cy3 =~
ly3 1.000
cy4 =~
ly4 1.000
cy5 =~
ly5 1.000
cy6 =~
ly6 1.000
l_intercept =~
ly1 1.000
l_slope =~
cy2 1.000
cy3 1.000
cy4 1.000
cy5 1.000
cy6 1.000
lx1 =~
x.1 1.000
lx2 =~
x.2 1.000
lx3 =~
x.3 1.000
lx4 =~
x.4 1.000
lx5 =~
x.5 1.000
lx6 =~
x.6 1.000
cx2 =~
lx2 1.000
cx3 =~
lx3 1.000
cx4 =~
lx4 1.000
cx5 =~
lx5 1.000
cx6 =~
lx6 1.000
lx_intercept =~
lx1 1.000
lx_slope =~
cx2 1.000
cx3 1.000
cx4 1.000
cx5 1.000
cx6 1.000
Regressions:
Estimate Std.Err z-value P(>|z|)
ly2 ~
ly1 1.000
ly3 ~
ly2 1.000
ly4 ~
ly3 1.000
ly5 ~
ly4 1.000
ly6 ~
ly5 1.000
cy2 ~
ly1 (prop) -0.315 0.004 -71.219 0.000
cy3 ~
ly2 (prop) -0.315 0.004 -71.219 0.000
cy4 ~
ly3 (prop) -0.315 0.004 -71.219 0.000
cy5 ~
ly4 (prop) -0.315 0.004 -71.219 0.000
cy6 ~
ly5 (prop) -0.315 0.004 -71.219 0.000
lx2 ~
lx1 1.000
lx3 ~
lx2 1.000
lx4 ~
lx3 1.000
lx5 ~
lx4 1.000
lx6 ~
lx5 1.000
cx2 ~
lx1 (prpx) 0.222 0.002 97.772 0.000
cx3 ~
lx2 (prpx) 0.222 0.002 97.772 0.000
cx4 ~
lx3 (prpx) 0.222 0.002 97.772 0.000
cx5 ~
lx4 (prpx) 0.222 0.002 97.772 0.000
cx6 ~
lx5 (prpx) 0.222 0.002 97.772 0.000
cy2 ~
lx1 (xy) 0.405 0.002 184.920 0.000
cy3 ~
lx2 (xy) 0.405 0.002 184.920 0.000
cy4 ~
lx3 (xy) 0.405 0.002 184.920 0.000
cy5 ~
lx4 (xy) 0.405 0.002 184.920 0.000
cy6 ~
lx5 (xy) 0.405 0.002 184.920 0.000
cx2 ~
ly1 (yx) -0.433 0.005 -93.914 0.000
cx3 ~
ly2 (yx) -0.433 0.005 -93.914 0.000
cx4 ~
ly3 (yx) -0.433 0.005 -93.914 0.000
cx5 ~
ly4 (yx) -0.433 0.005 -93.914 0.000
cx6 ~
ly5 (yx) -0.433 0.005 -93.914 0.000
Covariances:
Estimate Std.Err z-value P(>|z|)
l_intercept ~~
l_slope 0.392 0.141 2.786 0.005
.cy2 ~~
.cy3 0.000
.cy4 0.000
.cy5 0.000
.cy6 0.000
.cy3 ~~
.cy4 0.000
.cy5 0.000
.cy6 0.000
.cy4 ~~
.cy5 0.000
.cy6 0.000
.cy5 ~~
.cy6 0.000
lx_intercept ~~
lx_slope 0.015 0.125 0.123 0.902
.cx2 ~~
.cx3 0.000
.cx4 0.000
.cx5 0.000
.cx6 0.000
.cx3 ~~
.cx4 0.000
.cx5 0.000
.cx6 0.000
.cx4 ~~
.cx5 0.000
.cx6 0.000
.cx5 ~~
.cx6 0.000
l_intercept ~~
lx_intercept -0.154 0.050 -3.095 0.002
lx_slope 0.028 0.136 0.207 0.836
l_slope ~~
lx_intercept -0.206 0.131 -1.578 0.115
lx_slope 0.098 0.333 0.293 0.770
Intercepts:
Estimate Std.Err z-value P(>|z|)
l_slope 0.214 0.119 1.790 0.073
l_intercept 10.004 0.046 218.589 0.000
.ly1 0.000
.ly2 0.000
.ly3 0.000
.ly4 0.000
.ly5 0.000
.ly6 0.000
.cy2 0.000
.cy3 0.000
.cy4 0.000
.cy5 0.000
.cy6 0.000
.y.1 0.000
.y.2 0.000
.y.3 0.000
.y.4 0.000
.y.5 0.000
.y.6 0.000
lx_slope 0.550 0.118 4.644 0.000
lx_intercept 10.074 0.043 236.997 0.000
.lx1 0.000
.lx2 0.000
.lx3 0.000
.lx4 0.000
.lx5 0.000
.lx6 0.000
.cx2 0.000
.cx3 0.000
.cx4 0.000
.cx5 0.000
.cx6 0.000
.x.1 0.000
.x.2 0.000
.x.3 0.000
.x.4 0.000
.x.5 0.000
.x.6 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
l_ntr 0.990 0.076 13.059 0.000
l_slp 8.651 0.481 17.991 0.000
.ly1 0.000
.ly2 0.000
.ly3 0.000
.ly4 0.000
.ly5 0.000
.ly6 0.000
.cy2 0.000
.cy3 0.000
.cy4 0.000
.cy5 0.000
.cy6 0.000
.y.1 (res_vr) 0.738 0.020 37.356 0.000
.y.2 (res_vr) 0.738 0.020 37.356 0.000
.y.3 (res_vr) 0.738 0.020 37.356 0.000
.y.4 (res_vr) 0.738 0.020 37.356 0.000
.y.5 (res_vr) 0.738 0.020 37.356 0.000
.y.6 (res_vr) 0.738 0.020 37.356 0.000
lx_nt 1.028 0.065 15.705 0.000
lx_sl 8.406 0.459 18.296 0.000
.lx1 0.000
.lx2 0.000
.lx3 0.000
.lx4 0.000
.lx5 0.000
.lx6 0.000
.cx2 0.000
.cx3 0.000
.cx4 0.000
.cx5 0.000
.cx6 0.000
.x.1 (rs_vrx) 0.438 0.012 36.455 0.000
.x.2 (rs_vrx) 0.438 0.012 36.455 0.000
.x.3 (rs_vrx) 0.438 0.012 36.455 0.000
.x.4 (rs_vrx) 0.438 0.012 36.455 0.000
.x.5 (rs_vrx) 0.438 0.012 36.455 0.000
.x.6 (rs_vrx) 0.438 0.012 36.455 0.000Bo\(^2\)m =)