Lavaan - Example 5.1: Cross-Lagged Model with Observed Variables
Loading Libraries
library(MplusAutomation)
## Version: 0.7-3
## We work hard to write this free software. Please help us get credit by citing:
##
## Hallquist, M. N. & Wiley, J. F. (2018). MplusAutomation: An R Package for Facilitating Large-Scale Latent Variable Analyses in Mplus. Structural Equation Modeling, 25, 621-638. doi: 10.1080/10705511.2017.1402334.
##
## -- see citation("MplusAutomation").
library(texreg)
## Warning: package 'texreg' was built under R version 3.6.2
## Version: 1.37.5
## Date: 2020-06-17
## Author: Philip Leifeld (University of Essex)
##
## Consider submitting praise using the praise or praise_interactive functions.
## Please cite the JSS article in your publications -- see citation("texreg").
library(datasets)
library(lavaan)
## Warning: package 'lavaan' was built under R version 3.6.2
## This is lavaan 0.6-7
## lavaan is BETA software! Please report any bugs.
library(semPlot)
## Registered S3 methods overwritten by 'huge':
## method from
## plot.sim BDgraph
## print.sim BDgraph
library(readtext)
Preparing data for R
socex1.1 <- read.table("./Lista M6/Newsom_Ch5/Data/socex1.dat",header=FALSE)
Reading Lavaan Sintax
lav<-readtext::readtext("./Lista M6/Newsom_Ch5/R/ex5-1.R")
cat(lav$text)
## #title: Newsom Longitudinal SEM Chapter 5, Example 5.1
##
## library(lavaan)
## setwd("<your directory path")
## socex1.1 <- read.table ("socex1.dat", header=FALSE)
##
##
## names(socex1.1) = c("w1vst1", "w1vst2", "w1vst3", "w2vst1", "w2vst2",
## "w2vst3", "w3vst1", "w3vst2", "w3vst3", "w1unw1", "w1unw2", "w1unw3",
## "w2unw1", "w2unw2", "w2unw3", "w3unw1", "w3unw2", "w3unw3", "w1dboth",
## "w1dsad", "w1dblues", "w1ddep", "w2dboth", "w2dsad","w2dblues", "w2ddep",
## "w3dboth", "w3dsad", "w3dblues", "w3ddep", "w1marr2", "w1happy", "w1enjoy",
## "w1satis", "w1joyful", "w1please", "w2happy", "w2enjoy", "w2satis", "w2joyful",
## "w2please", "w3happy", "w3enjoy", "w3satis", "w3joyful", "w3please", "w1lea",
## "w2lea", "w3lea")
##
## #+++++++++++++++++++++++++++++++
## # Chapter 5, Example 5.1
## #+++++++++++++++++++++++++++++++
##
## attach(socex1.1)
##
## #define variables
##
## socex1.1$w1unw = (w1unw1 + w1unw2 + w1unw3)/3
## socex1.1$w2unw = (w2unw1 + w2unw2 + w2unw3)/3
## socex1.1$w1posaff = (w1happy + w1enjoy + w1satis + w1joyful + w1please)/5
## socex1.1$w2posaff = (w2happy + w2enjoy + w2satis + w2joyful + w2please)/5
##
## model5.2a <- ' #cross-lagged panel model with observed variables;
## w2unw ~ w1posaff + w1unw
## w2posaff ~ w1unw + w1posaff
## w1unw ~~ w1posaff
## w2unw ~~ w2posaff
## '
##
## fitmodel5.2a <- sem(model5.2a, data=socex1.1, information="expected")
## summary(fitmodel5.2a, fit.measures=TRUE, standardized=TRUE, rsquare=TRUE)
Naming Dataframe
names(socex1.1) = c("w1vst1", "w1vst2", "w1vst3", "w2vst1", "w2vst2",
"w2vst3", "w3vst1", "w3vst2", "w3vst3", "w1unw1", "w1unw2", "w1unw3",
"w2unw1", "w2unw2", "w2unw3", "w3unw1", "w3unw2", "w3unw3", "w1dboth",
"w1dsad", "w1dblues", "w1ddep", "w2dboth", "w2dsad","w2dblues", "w2ddep",
"w3dboth", "w3dsad", "w3dblues", "w3ddep", "w1marr2", "w1happy", "w1enjoy",
"w1satis", "w1joyful", "w1please", "w2happy", "w2enjoy", "w2satis", "w2joyful",
"w2please", "w3happy", "w3enjoy", "w3satis", "w3joyful", "w3please", "w1lea",
"w2lea", "w3lea")
Preparing Data
#+++++++++++++++++++++++++++++++
# Chapter 5, Example 5.1
#+++++++++++++++++++++++++++++++
# Makes socex1.1 default
attach(socex1.1)
#define variables
socex1.1$w1unw = (w1unw1 + w1unw2 + w1unw3)/3
socex1.1$w2unw = (w2unw1 + w2unw2 + w2unw3)/3
socex1.1$w1posaff = (w1happy + w1enjoy + w1satis + w1joyful + w1please)/5
socex1.1$w2posaff = (w2happy + w2enjoy + w2satis + w2joyful + w2please)/5
Lavaan Model
lavaan.model <-
'
#Cross-lagged panel model with observed variables;
# Regressions ON (~)
w2unw ~ w1posaff + w1unw
w2posaff ~ w1unw + w1posaff
# Co-variation WITH (~~)
w1unw ~~ w1posaff
w2unw ~~ w2posaff
'
Model Fit
lavaan.fit <- sem(lavaan.model, data=socex1.1, information="expected", estimator="ML")
Model Summary
summary(lavaan.fit, fit.measures=TRUE, standardized=TRUE, rsquare=TRUE)
## lavaan 0.6-7 ended normally after 21 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of free parameters 10
##
## Number of observations 574
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Model Test Baseline Model:
##
## Test statistic 273.360
## Degrees of freedom 6
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.000
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -2387.273
## Loglikelihood unrestricted model (H1) -2387.273
##
## Akaike (AIC) 4794.547
## Bayesian (BIC) 4838.073
## Sample-size adjusted Bayesian (BIC) 4806.327
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.000
## P-value RMSEA <= 0.05 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## w2unw ~
## w1posaff -0.199 0.050 -3.997 0.000 -0.199 -0.155
## w1unw 0.344 0.038 9.011 0.000 0.344 0.349
## w2posaff ~
## w1unw -0.029 0.030 -0.968 0.333 -0.029 -0.036
## w1posaff 0.490 0.039 12.436 0.000 0.490 0.465
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## w1posaff ~~
## w1unw -0.093 0.022 -4.222 0.000 -0.093 -0.179
## .w2unw ~~
## .w2posaff -0.043 0.018 -2.377 0.017 -0.043 -0.100
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .w2unw 0.546 0.032 16.941 0.000 0.546 0.835
## .w2posaff 0.343 0.020 16.941 0.000 0.343 0.777
## w1posaff 0.398 0.023 16.941 0.000 0.398 1.000
## w1unw 0.675 0.040 16.941 0.000 0.675 1.000
##
## R-Square:
## Estimate
## w2unw 0.165
## w2posaff 0.223
Fit Measures - Clean Output
lavaan.fit.measures<-as.data.frame(lavaan::fitmeasures(lavaan.fit, fit.measures = c("rmsea", "rmsea.ci.lower","rmsea.ci.upper","cfi","tli" )))
names(lavaan.fit.measures)[1]<-"Fit Measures"
print(lavaan.fit.measures)
## Fit Measures
## rmsea 0
## rmsea.ci.lower 0
## rmsea.ci.upper 0
## cfi 1
## tli 1
Parameters Estimate - Clean Output
fit.parameter<-lavaan::parameterestimates(lavaan.fit,output ="text", standardized =T,rsquare=T)
fit.parameter
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## w2unw ~
## w1posaff -0.199 0.050 -3.997 0.000 -0.296 -0.101
## w1unw 0.344 0.038 9.011 0.000 0.269 0.419
## w2posaff ~
## w1unw -0.029 0.030 -0.968 0.333 -0.089 0.030
## w1posaff 0.490 0.039 12.436 0.000 0.413 0.567
## Std.lv Std.all
##
## -0.199 -0.155
## 0.344 0.349
##
## -0.029 -0.036
## 0.490 0.465
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## w1posaff ~~
## w1unw -0.093 0.022 -4.222 0.000 -0.136 -0.050
## .w2unw ~~
## .w2posaff -0.043 0.018 -2.377 0.017 -0.079 -0.008
## Std.lv Std.all
##
## -0.093 -0.179
##
## -0.043 -0.100
##
## Variances:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## .w2unw 0.546 0.032 16.941 0.000 0.482 0.609
## .w2posaff 0.343 0.020 16.941 0.000 0.303 0.383
## w1posaff 0.398 0.023 16.941 0.000 0.352 0.444
## w1unw 0.675 0.040 16.941 0.000 0.596 0.753
## Std.lv Std.all
## 0.546 0.835
## 0.343 0.777
## 0.398 1.000
## 0.675 1.000
##
## R-Square:
## Estimate
## w2unw 0.165
## w2posaff 0.223
Basic Gray Plot
semPaths(lavaan.fit,whatLabel = "est", style="ram", what="path")
Colored and Weighted Plot
semPaths(lavaan.fit,whatLabel = "est", style="ram", what="est")
Updating Previous Lavaan Model
lavaan.model <-
'
#Cross-lagged panel model with observed variables;
# Regressions ON (~)
w2unw ~ w1posaff
w2posaff ~ w1unw + w1posaff
# Co-variation WITH (~~)
w1unw ~~ w1posaff
w2unw ~~ w2posaff
'
Updated Model Fit
lavaan.fit <- sem(lavaan.model, data=socex1.1, information="expected", estimator="MLR")
Updated Fit Measures - Clean Output
lavaan.fit.measures<-as.data.frame(lavaan::fitmeasures(lavaan.fit, fit.measures = c("rmsea", "rmsea.ci.lower","rmsea.ci.upper","cfi","tli" )))
names(lavaan.fit.measures)[1]<-"Fit Measures"
print(lavaan.fit.measures)
## Fit Measures
## rmsea 0.3613232
## rmsea.ci.lower 0.2950712
## rmsea.ci.upper 0.4323809
## cfi 0.7197105
## tli -0.6817368
Updated Parameters Estimate - Clean Output
fit.parameter<-lavaan::parameterestimates(lavaan.fit,output ="text", standardized =T,rsquare=T)
fit.parameter
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## w2unw ~
## w1posaff -0.279 0.054 -5.186 0.000 -0.384 -0.173
## w2posaff ~
## w1unw -0.002 0.028 -0.074 0.941 -0.058 0.053
## w1posaff 0.496 0.038 12.900 0.000 0.421 0.572
## Std.lv Std.all
##
## -0.279 -0.217
##
## -0.002 -0.003
## 0.496 0.471
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## w1posaff ~~
## w1unw -0.093 0.021 -4.458 0.000 -0.134 -0.052
## .w2unw ~~
## .w2posaff -0.049 0.018 -2.673 0.008 -0.085 -0.013
## Std.lv Std.all
##
## -0.093 -0.179
##
## -0.049 -0.106
##
## Variances:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## .w2unw 0.623 0.035 17.931 0.000 0.555 0.691
## .w2posaff 0.344 0.020 17.326 0.000 0.305 0.382
## w1posaff 0.398 0.025 16.117 0.000 0.350 0.446
## w1unw 0.675 0.034 19.559 0.000 0.607 0.742
## Std.lv Std.all
## 0.623 0.953
## 0.344 0.778
## 0.398 1.000
## 0.675 1.000
##
## R-Square:
## Estimate
## w2unw 0.047
## w2posaff 0.222
Updated Basic Gray Plot
semPaths(lavaan.fit,whatLabel = "est", style="ram", what="path")
Updated Colored and Weighted Plot
semPaths(lavaan.fit,whatLabel = "est", style="ram", what="est")
Mplus - Example 5.1: Cross-Lagged Model with Observed Variables
Reading original .inp file
inp<-readtext::readtext("./Lista M6/Newsom_Ch5/MPLUS/ex5-1a.inp")
cat(inp$text)
## title: Newsom Longitudinal SEM Chapter 5, Example 5.1,
## Cross-lagged Panel Models;
##
## data: file=socex1.dat; format=free;
##
## define:
## ! create equally weighted composites for unwanted advice
## ! and positive affect;
## w1unw = (w1unw1 + w1unw2 + w1unw3)/3;
## w2unw = (w2unw1 + w2unw2 + w2unw3)/3;
## w1posaff = (w1happy + w1enjoy + w1satis + w1joyful + w1please)/5;
## w2posaff = (w2happy + w2enjoy + w2satis + w2joyful + w2please)/5;
##
##
## variable:
## names=
## w1vst1 w1vst2 w1vst3
## w2vst1 w2vst2 w2vst3
## w3vst1 w3vst2 w3vst3
## w1unw1 w1unw2 w1unw3
## w2unw1 w2unw2 w2unw3
## w3unw1 w3unw2 w3unw3
## w1dboth w1dsad w1dblues w1ddep
## w2dboth w2dsad w2dblues w2ddep
## w3dboth w3dsad w3dblues w3ddep
## w1marr2
## w1happy w1enjoy w1satis w1joyful w1please
## w2happy w2enjoy w2satis w2joyful w2please
## w3happy w3enjoy w3satis w3joyful w3please
## w1lea w2lea w3lea;
##
## usevariables=
## w1posaff w2posaff w1unw w2unw;
##
## analysis: type=general;
## information=expected; model=nomeanstructure; !remove if means desired;
##
## model:
## !cross-lagged and autoregressive effects;
## w2posaff on w1posaff w1unw;
## w2unw on w1unw w1posaff;
##
## output: sampstat stdyx ;
Writing a Mplus Object using R
mplus.object<-mplusObject(
# Basic Mplus .inp
TITLE = NULL,
DATA = NULL,
VARIABLE = NULL,
DEFINE = NULL,
MONTECARLO = NULL,
MODELPOPULATION = NULL,
MODELMISSING = NULL,
ANALYSIS = NULL,
# Model
MODEL =
"w2posaff on w1posaff w1unw;
w2unw on w1unw w1posaff;",
MODELINDIRECT = NULL,
MODELCONSTRAINT = NULL,
MODELTEST = NULL,
MODELPRIORS = NULL,
# Output
OUTPUT = NULL,
SAVEDATA = NULL,
PLOT = NULL,
# R notation
usevariables = NULL, rdata =socex1.1, autov = TRUE, imputed = FALSE)
## R variables selected automatically as any variable name that
## occurs in the MODEL, VARIABLE, or DEFINE section.
## If any issues, suggest explicitly specifying USEVARIABLES.
## A starting point may be:
## USEVARIABLES = w1unw w2unw w1posaff w2posaff;
Creating a complete Mplus Object as a list in R
mplus.model<-mplusModeler(
# Sintaxe R
mplus.object,
# Dados formato R para formato .dat
dataout ="dados.dat",
# Sintaxe Mplus
modelout ="sintaxe.inp",
# Run one time
run=1)
## Wrote model to: sintaxe.inp
## File(s) with md5 hash matching data found, using
## dados_aab16c8f476c462166e677cdbfb8033f.dat
##
## Running model: sintaxe.inp
## System command: cd "." && "/Applications/Mplus/mplus" "sintaxe.inp"
## Reading model: sintaxe.out
Summary Basic Model
summary(mplus.model)
## Estimated using ML
## Number of obs: 574, number of (free) parameters: 9
##
## Model: Chi2(df = 0) = 0, p = 0
## Baseline model: Chi2(df = 5) = 254.661, p = 0
##
## Fit Indices:
##
## CFI = 1, TLI = 1, SRMR = 0
## RMSEA = 0, 90% CI [0, 0], p < .05 = 0
## AIC = 2308.202, BIC = 2347.376
Results Fit Indices
mplus.model.fit<-data.frame(
c(
mplus.model$results$summaries$RMSEA_Estimate,
mplus.model$results$summaries$RMSEA_90CI_LB,
mplus.model$results$summaries$RMSEA_90CI_UB,
mplus.model$results$summaries$CFI,
mplus.model$results$summaries$TLI)
)
names(mplus.model.fit)[1]<-"Fit Indices"
rownames(mplus.model.fit)<-c("RMSE Estimate","RMSEA 90 CI Lower","RMSEA 90 CI Upper","CFI","TLI")
mplus.model.fit
## Fit Indices
## RMSE Estimate 0
## RMSEA 90 CI Lower 0
## RMSEA 90 CI Upper 0
## CFI 1
## TLI 1
Results Parameters
mplus.model$results$parameters$unstandardized
## paramHeader param est se est_se pval
## 1 W2POSAFF.ON W1POSAFF 0.490 0.039 12.436 0.000
## 2 W2POSAFF.ON W1UNW -0.029 0.030 -0.968 0.333
## 3 W2UNW.ON W1UNW 0.344 0.038 9.011 0.000
## 4 W2UNW.ON W1POSAFF -0.199 0.050 -3.997 0.000
## 5 W2POSAFF.WITH W2UNW -0.043 0.018 -2.377 0.017
## 6 Intercepts W2UNW 1.905 0.183 10.417 0.000
## 7 Intercepts W2POSAFF 1.497 0.145 10.325 0.000
## 8 Residual.Variances W2UNW 0.546 0.032 16.941 0.000
## 9 Residual.Variances W2POSAFF 0.343 0.020 16.941 0.000
Updated Version
mplus.model.update<-update(
# Previous model
mplus.model,
# Changed Parameters
ANALYSIS = ~ + "ESTIMATOR = MLR;",
MODEL = ~
"w2posaff on w1posaff w1unw;
w2unw on w1posaff;"
)
## R variables selected automatically as any variable name that
## occurs in the MODEL, VARIABLE, or DEFINE section.
## If any issues, suggest explicitly specifying USEVARIABLES.
## A starting point may be:
## USEVARIABLES = w1unw w2unw w1posaff w2posaff;
Updated Running model
mplus.model.update<-mplusModeler(
# Sintaxe R
mplus.model.update,
# Dados formato R para formato .dat
dataout ="dados.dat",
# Sintaxe Mplus
modelout ="sintaxe.inp",
# Run one time
run=1)
## Wrote model to: sintaxe.inp
## File(s) with md5 hash matching data found, using
## dados_aab16c8f476c462166e677cdbfb8033f.dat
##
## Running model: sintaxe.inp
## System command: cd "." && "/Applications/Mplus/mplus" "sintaxe.inp"
## Reading model: sintaxe.out
Update Summary Models
summary(mplus.model.update)
## Estimated using MLR
## Number of obs: 574, number of (free) parameters: 8
##
## Model: Chi2(df = 1) = 77.57, p = 0
## Baseline model: Chi2(df = 5) = 267.438, p = 0
##
## Fit Indices:
##
## CFI = 0.708, TLI = -0.459, SRMR = 0.113
## RMSEA = 0.365, 90% CI [0.299, 0.436], p < .05 = 0
## AIC = 2382.14, BIC = 2416.961
Update Results Fit Indices
mplus.model.update.fit<-data.frame(
c(
mplus.model.update$results$summaries$RMSEA_Estimate,
mplus.model.update$results$summaries$RMSEA_90CI_LB,
mplus.model.update$results$summaries$RMSEA_90CI_UB,
mplus.model.update$results$summaries$CFI,
mplus.model.update$results$summaries$TLI)
)
names(mplus.model.update.fit)[1]<-"Fit Indices"
names(mplus.model.update.fit)[1]<-"Fit Indices"
rownames(mplus.model.update.fit)<-c("RMSE Estimate","RMSEA 90 CI Lower","RMSEA 90 CI Upper","CFI","TLI")
mplus.model.update.fit
## Fit Indices
## RMSE Estimate 0.365
## RMSEA 90 CI Lower 0.299
## RMSEA 90 CI Upper 0.436
## CFI 0.708
## TLI -0.459
Update Results Parameters
mplus.model.update$results$parameters$unstandardized
## paramHeader param est se est_se pval
## 1 W2POSAFF.ON W1POSAFF 0.496 0.039 12.885 0.000
## 2 W2POSAFF.ON W1UNW -0.002 0.030 -0.070 0.944
## 3 W2UNW.ON W1POSAFF -0.279 0.054 -5.186 0.000
## 4 W2POSAFF.WITH W2UNW -0.049 0.020 -2.524 0.012
## 5 Intercepts W2UNW 2.828 0.168 16.825 0.000
## 6 Intercepts W2POSAFF 1.424 0.147 9.697 0.000
## 7 Residual.Variances W2UNW 0.623 0.035 17.931 0.000
## 8 Residual.Variances W2POSAFF 0.344 0.020 17.322 0.000