rm(list=ls())
library(lavaan)
library(MASS)
library(lme4)
library(lmerTest)
library(ggplot2)set.seed(1)
# stabilisco il sample size
N = 5000
# genero i dati progressivamente per ogni wave
x0 = rnorm(N,0,1)
x1 = rnorm(N,0,1) + x0*0.8
x2 = rnorm(N,0,1) + x1*0.8
x3 = rnorm(N,0,1) + x2*0.8
# metto tutto in un dataset
df = data.frame(x0,x1,x2,x3)model = "
x1 ~ x0
x2 ~ x1
x3 ~ x2
"
fit = sem(model=model, data=df)
summary(fit, standardized=T)## lavaan 0.6.16 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 6
##
## Number of observations 5000
##
## Model Test User Model:
##
## Test statistic 1.989
## Degrees of freedom 3
## P-value (Chi-square) 0.575
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## x1 ~
## x0 0.798 0.014 58.085 0.000 0.798 0.635
## x2 ~
## x1 0.814 0.011 75.818 0.000 0.814 0.731
## x3 ~
## x2 0.794 0.010 80.527 0.000 0.794 0.751
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .x1 0.996 0.020 50.000 0.000 0.996 0.597
## .x2 0.960 0.019 50.000 0.000 0.960 0.465
## .x3 1.002 0.020 50.000 0.000 1.002 0.435set.seed(1)
# stabilisco il sample size
N = 5000
# genero le intercette individuali degli N soggetti
X_int = rnorm(N,0,2)
# genero i "residui-within" autoregressivi progressivamente per ogni wave
wx0 = rnorm(N,0,1)
wx1 = rnorm(N,0,1) + wx0*0.5
wx2 = rnorm(N,0,1) + wx1*0.5
wx3 = rnorm(N,0,1) + wx2*0.5
# sommo tutto (intercette, e residui-within) generando le variabili osservate
x0 = X_int*1 + wx0*1
x1 = X_int*1 + wx1*1
x2 = X_int*1 + wx2*1
x3 = X_int*1 + wx3*1
# metto tutte le variabili osservate in un dataset
df = data.frame(x0,x1,x2,x3)model = "
# DEFINISCO RANDOM INTERCEPT COME LATENTE
X_intercept =~ 1*x0 + 1*x1 + 1*x2 + 1*x3
# RESIDUI WITHIN-SUBJECT
wx0 =~ 1*x0
wx1 =~ 1*x1
wx2 =~ 1*x2
wx3 =~ 1*x3
# PARTE AUTOREGRESSIVA (SU RESIDUI WITHIN)
wx1 ~ wx0
wx2 ~ wx1
wx3 ~ wx2
# SPECIFICO VALORE MEDIO DELL'INTERCETTA
X_intercept ~ 1
# SPECIFICO VARIANZA DELL'INTERCETTA
X_intercept ~~ X_intercept
# SPECIFICO VARIANZE DEI RESIDUI
wx0 ~~ wx0
wx1 ~~ wx1
wx2 ~~ wx2
wx3 ~~ wx3
"
fit = lavaan(model=model, data=df)
summary(fit, standardized=T)## lavaan 0.6.16 ended normally after 49 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 9
##
## Number of observations 5000
##
## Model Test User Model:
##
## Test statistic 11.263
## Degrees of freedom 5
## P-value (Chi-square) 0.046
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## X_intercept =~
## x0 1.000 2.049 0.897
## x1 1.000 2.049 0.877
## x2 1.000 2.049 0.870
## x3 1.000 2.049 0.871
## wx0 =~
## x0 1.000 1.012 0.443
## wx1 =~
## x1 1.000 1.123 0.481
## wx2 =~
## x2 1.000 1.161 0.493
## wx3 =~
## x3 1.000 1.154 0.491
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## wx1 ~
## wx0 0.527 0.034 15.661 0.000 0.475 0.475
## wx2 ~
## wx1 0.515 0.026 19.872 0.000 0.498 0.498
## wx3 ~
## wx2 0.511 0.021 23.936 0.000 0.514 0.514
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## X_intercept -0.007 0.031 -0.214 0.831 -0.003 -0.003
## .x0 0.000 0.000 0.000
## .x1 0.000 0.000 0.000
## .x2 0.000 0.000 0.000
## .x3 0.000 0.000 0.000
## wx0 0.000 0.000 0.000
## .wx1 0.000 0.000 0.000
## .wx2 0.000 0.000 0.000
## .wx3 0.000 0.000 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## X_intercept 4.200 0.102 41.243 0.000 1.000 1.000
## wx0 1.024 0.057 18.092 0.000 1.000 1.000
## .wx1 0.977 0.031 31.726 0.000 0.774 0.774
## .wx2 1.013 0.029 35.224 0.000 0.752 0.752
## .wx3 0.979 0.027 36.921 0.000 0.736 0.736
## .x0 0.000 0.000 0.000
## .x1 0.000 0.000 0.000
## .x2 0.000 0.000 0.000
## .x3 0.000 0.000 0.000set.seed(1)
# stabilisco il sample size
N = 5000
# genero le intercette e le slope individuali degli N soggetti
X_int = rnorm(N,0,2)
X_slo = rnorm(N,0,.5) + X_int*0.1 # slope correlata all'intercetta
cor(X_int, X_slo)## [1] 0.3792832# genero i "residui-within" autoregressivi progressivamente per ogni wave
wx0 = rnorm(N,0,1)
wx1 = rnorm(N,0,1) + wx0*0.5
wx2 = rnorm(N,0,1) + wx1*0.5
wx3 = rnorm(N,0,1) + wx2*0.5
# sommo tutto (intercette, slope, e residui-within) generando le variabili osservate
x0 = X_int*1 + X_slo*0 + wx0*1
x1 = X_int*1 + X_slo*1 + wx1*1
x2 = X_int*1 + X_slo*2 + wx2*1
x3 = X_int*1 + X_slo*3 + wx3*1
# metto tutte le variabili osservate in un dataset
df = data.frame(x0,x1,x2,x3)model = "
# DEFINISCO RANDOM INTERCEPT E SLOPE COME LATENTI
X_intercept =~ 1*x0 + 1*x1 + 1*x2 + 1*x3
X_slope =~ 0*x0 + 1*x1 + 2*x2 + 3*x3
# RESIDUI WITHIN-SUBJECT
wx0 =~ 1*x0
wx1 =~ 1*x1
wx2 =~ 1*x2
wx3 =~ 1*x3
# PARTE AUTOREGRESSIVA (SU RESIDUI WITHIN)
wx1 ~ wx0
wx2 ~ wx1
wx3 ~ wx2
# SPECIFICO CORRELAZIONE INTERCETTA-SLOPE
X_intercept ~~ X_slope
# SPECIFICO VALORI MEDI DELL'INTERCETTA e DELLA SLOPE
X_intercept ~ 1
X_slope ~ 1
# SPECIFICO VARIANZE DELL'INTERCETTA e DELLA SLOPE
X_intercept ~~ X_intercept
X_slope ~~ X_slope
# SPECIFICO VARIANZE DEI RESIDUI
wx0 ~~ wx0
wx1 ~~ wx1
wx2 ~~ wx2
wx3 ~~ wx3
"
fit = lavaan(model=model, data=df)
summary(fit, standardized=T)## lavaan 0.6.16 ended normally after 146 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 12
##
## Number of observations 5000
##
## Model Test User Model:
##
## Test statistic 5.282
## Degrees of freedom 2
## P-value (Chi-square) 0.071
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## X_intercept =~
## x0 1.000 2.126 0.926
## x1 1.000 2.126 0.828
## x2 1.000 2.126 0.735
## x3 1.000 2.126 0.650
## X_slope =~
## x0 0.000 0.000 0.000
## x1 1.000 0.493 0.192
## x2 2.000 0.986 0.341
## x3 3.000 1.479 0.452
## wx0 =~
## x0 1.000 0.865 0.377
## wx1 =~
## x1 1.000 1.081 0.421
## wx2 =~
## x2 1.000 1.246 0.431
## wx3 =~
## x3 1.000 1.417 0.433
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## wx1 ~
## wx0 0.436 0.145 3.007 0.003 0.349 0.349
## wx2 ~
## wx1 0.589 0.078 7.575 0.000 0.511 0.511
## wx3 ~
## wx2 0.687 0.115 5.997 0.000 0.604 0.604
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## X_intercept ~~
## X_slope 0.330 0.081 4.062 0.000 0.315 0.315
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## X_intercept -0.024 0.032 -0.728 0.467 -0.011 -0.011
## X_slope 0.004 0.010 0.442 0.659 0.009 0.009
## .x0 0.000 0.000 0.000
## .x1 0.000 0.000 0.000
## .x2 0.000 0.000 0.000
## .x3 0.000 0.000 0.000
## wx0 0.000 0.000 0.000
## .wx1 0.000 0.000 0.000
## .wx2 0.000 0.000 0.000
## .wx3 0.000 0.000 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## X_intercept 4.521 0.332 13.636 0.000 1.000 1.000
## X_slope 0.243 0.049 4.967 0.000 1.000 1.000
## wx0 0.747 0.318 2.347 0.019 1.000 1.000
## .wx1 1.026 0.061 16.810 0.000 0.878 0.878
## .wx2 1.147 0.125 9.151 0.000 0.739 0.739
## .wx3 1.276 0.137 9.283 0.000 0.635 0.635
## .x0 0.000 0.000 0.000
## .x1 0.000 0.000 0.000
## .x2 0.000 0.000 0.000
## .x3 0.000 0.000 0.000set.seed(1)
# stabilisco il sample size
N = 50000
# genero le intercette e le slope individuali degli N soggetti
X_int = rnorm(N,0,2)
X_slo = rnorm(N,0,.5) + X_int*0.1 # slope di X correlata all'intercetta di X
Y_int = rnorm(N,0,1.5) + X_int*0.2 # intercetta di Y correlata all'intercetta di X
Y_slo = rnorm(N,0,.8)
cor(X_int, X_slo)## [1] 0.3690963cor(X_int, Y_int)## [1] 0.2521154# genero i "residui-within" autoregressivi e cross-lagged (y -> x) progressivamente per ogni wave
wx0 = rnorm(N,0,1)
wy0 = rnorm(N,0,1)
wx1 = rnorm(N,0,1) + wx0*0.5 + wy0*0.2
wy1 = rnorm(N,0,1) + wy0*0.4
wx2 = rnorm(N,0,1) + wx1*0.5 + wy1*0.2
wy2 = rnorm(N,0,1) + wy1*0.4
wx3 = rnorm(N,0,1) + wx2*0.5 + wy2*0.2
wy3 = rnorm(N,0,1) + wy2*0.4
# sommo tutto (intercette, slope, e residui-within) generando le variabili osservate
x0 = X_int*1 + X_slo*0 + wx0*1
x1 = X_int*1 + X_slo*1 + wx1*1
x2 = X_int*1 + X_slo*2 + wx2*1
x3 = X_int*1 + X_slo*3 + wx3*1
y0 = Y_int*1 + Y_slo*0 + wy0*1
y1 = Y_int*1 + Y_slo*1 + wy1*1
y2 = Y_int*1 + Y_slo*2 + wy2*1
y3 = Y_int*1 + Y_slo*3 + wy3*1
# metto tutte le variabili osservate in un dataset
df = data.frame(x0,x1,x2,x3, y0,y1,y2,y3)model = "
# DEFINISCO INTERCETTE E SLOPE LATENTI
X_intercept =~ 1*x0 + 1*x1 + 1*x2 + 1*x3
X_slope =~ 0*x0 + 1*x1 + 2*x2 + 3*x3
Y_intercept =~ 1*y0 + 1*y1 + 1*y2 + 1*y3
Y_slope =~ 0*y0 + 1*y1 + 2*y2 + 3*y3
# RESIDUI WITHIN-SUBJECT
wx0 =~ 1*x0
wx1 =~ 1*x1
wx2 =~ 1*x2
wx3 =~ 1*x3
wy0 =~ 1*y0
wy1 =~ 1*y1
wy2 =~ 1*y2
wy3 =~ 1*y3
# PARTE AUTOREGRESSIVA (SU RESIDUI WITHIN) CON CONSTRAINTS TEMPORALI
wx1 ~ arx*wx0 + clyx*wy0
wx2 ~ arx*wx1 + clyx*wy1
wx3 ~ arx*wx2 + clyx*wy2
wy1 ~ ary*wy0 + clxy*wx0
wy2 ~ ary*wy1 + clxy*wx1
wy3 ~ ary*wy2 + clxy*wx2
# SPECIFICO CORRELAZIONI INTERCETTA-SLOPE E ANCHE TRA VARIABILI
X_intercept ~~ X_slope
Y_intercept ~~ Y_slope
X_intercept ~~ Y_intercept
X_slope ~~ Y_slope
# SPECIFICO VALORI MEDI DELL'INTERCETTA e DELLA SLOPE
X_intercept ~ 1
X_slope ~ 1
Y_intercept ~ 1
Y_slope ~ 1
# SPECIFICO VARIANZE DELL'INTERCETTA e DELLA SLOPE
X_intercept ~~ X_intercept
X_slope ~~ X_slope
Y_intercept ~~ Y_intercept
Y_slope ~~ Y_slope
# SPECIFICO VARIANZE DEI RESIDUI
wx0 ~~ wx0
wx1 ~~ wx1
wx2 ~~ wx2
wx3 ~~ wx3
wy0 ~~ wy0
wy1 ~~ wy1
wy2 ~~ wy2
wy3 ~~ wy3
"
fit = lavaan(model=model, data=df)
summary(fit, standardized=T)## lavaan 0.6.16 ended normally after 111 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 32
## Number of equality constraints 8
##
## Number of observations 50000
##
## Model Test User Model:
##
## Test statistic 143.436
## Degrees of freedom 20
## P-value (Chi-square) 0.000
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## X_intercept =~
## x0 1.000 2.011 0.894
## x1 1.000 2.011 0.792
## x2 1.000 2.011 0.703
## x3 1.000 2.011 0.624
## X_slope =~
## x0 0.000 0.000 0.000
## x1 1.000 0.528 0.208
## x2 2.000 1.056 0.369
## x3 3.000 1.585 0.492
## Y_intercept =~
## y0 1.000 1.505 0.814
## y1 1.000 1.505 0.736
## y2 1.000 1.505 0.607
## y3 1.000 1.505 0.492
## Y_slope =~
## y0 0.000 0.000 0.000
## y1 1.000 0.788 0.385
## y2 2.000 1.575 0.635
## y3 3.000 2.363 0.772
## wx0 =~
## x0 1.000 1.009 0.449
## wx1 =~
## x1 1.000 1.157 0.456
## wx2 =~
## x2 1.000 1.206 0.422
## wx3 =~
## x3 1.000 1.218 0.378
## wy0 =~
## y0 1.000 1.073 0.581
## wy1 =~
## y1 1.000 1.101 0.539
## wy2 =~
## y2 1.000 1.118 0.451
## wy3 =~
## y3 1.000 1.128 0.369
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## wx1 ~
## wx0 (arx) 0.517 0.012 44.205 0.000 0.451 0.451
## wy0 (clyx) 0.204 0.004 48.821 0.000 0.190 0.190
## wx2 ~
## wx1 (arx) 0.517 0.012 44.205 0.000 0.496 0.496
## wy1 (clyx) 0.204 0.004 48.821 0.000 0.187 0.187
## wx3 ~
## wx2 (arx) 0.517 0.012 44.205 0.000 0.512 0.512
## wy2 (clyx) 0.204 0.004 48.821 0.000 0.188 0.188
## wy1 ~
## wy0 (ary) 0.433 0.011 38.738 0.000 0.422 0.422
## wx0 (clxy) -0.005 0.004 -1.112 0.266 -0.004 -0.004
## wy2 ~
## wy1 (ary) 0.433 0.011 38.738 0.000 0.426 0.426
## wx1 (clxy) -0.005 0.004 -1.112 0.266 -0.005 -0.005
## wy3 ~
## wy2 (ary) 0.433 0.011 38.738 0.000 0.429 0.429
## wx2 (clxy) -0.005 0.004 -1.112 0.266 -0.005 -0.005
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## X_intercept ~~
## X_slope 0.390 0.014 27.178 0.000 0.368 0.368
## Y_intercept ~~
## Y_slope 0.040 0.012 3.422 0.001 0.033 0.033
## X_intercept ~~
## Y_intercept 0.802 0.018 43.733 0.000 0.265 0.265
## X_slope ~~
## Y_slope 0.006 0.003 2.081 0.037 0.015 0.015
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## X_intercept -0.008 0.010 -0.841 0.400 -0.004 -0.004
## X_slope 0.000 0.003 0.020 0.984 0.000 0.000
## Y_intercept 0.009 0.008 1.112 0.266 0.006 0.006
## Y_slope -0.004 0.004 -0.993 0.321 -0.005 -0.005
## .x0 0.000 0.000 0.000
## .x1 0.000 0.000 0.000
## .x2 0.000 0.000 0.000
## .x3 0.000 0.000 0.000
## .y0 0.000 0.000 0.000
## .y1 0.000 0.000 0.000
## .y2 0.000 0.000 0.000
## .y3 0.000 0.000 0.000
## wx0 0.000 0.000 0.000
## .wx1 0.000 0.000 0.000
## .wx2 0.000 0.000 0.000
## .wx3 0.000 0.000 0.000
## wy0 0.000 0.000 0.000
## .wy1 0.000 0.000 0.000
## .wy2 0.000 0.000 0.000
## .wy3 0.000 0.000 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## X_intercept 4.043 0.058 69.728 0.000 1.000 1.000
## X_slope 0.279 0.007 39.777 0.000 1.000 1.000
## Y_intercept 2.266 0.040 56.928 0.000 1.000 1.000
## Y_slope 0.620 0.007 86.647 0.000 1.000 1.000
## wx0 1.018 0.050 20.224 0.000 1.000 1.000
## .wx1 1.018 0.012 86.789 0.000 0.761 0.761
## .wx2 1.023 0.015 68.126 0.000 0.704 0.704
## .wx3 1.017 0.016 63.078 0.000 0.685 0.685
## wy0 1.152 0.037 31.374 0.000 1.000 1.000
## .wy1 0.997 0.011 86.943 0.000 0.822 0.822
## .wy2 1.023 0.016 63.453 0.000 0.819 0.819
## .wy3 1.038 0.015 67.248 0.000 0.816 0.816
## .x0 0.000 0.000 0.000
## .x1 0.000 0.000 0.000
## .x2 0.000 0.000 0.000
## .x3 0.000 0.000 0.000
## .y0 0.000 0.000 0.000
## .y1 0.000 0.000 0.000
## .y2 0.000 0.000 0.000
## .y3 0.000 0.000 0.000