## Demonstrate discrepancy between lavaan and OpenMx 
## in the computation of the standardized covariance 
## among residuals of endogenous variables in a path model
## (using a data set that comes with lavaan)

library(lavaan)
library(OpenMx)

# lavaan 
#
pathl <- '
 t1 ~ x1 + x2
 t2 ~ x1 + x2
 t3 ~ x1 + x2 + t1 + t2
 x1 ~~ x2
 t1 ~~ t2
'
fitl <- sem(pathl,data=Demo.growth,missing="ML")
summary(fitl,standardized=T,fit.measures=T)
standardizedsolution(fitl) -> stdl
parametertable(fitl) -> unsl

# OpenMx
#
mvars <- c("x1","x2","t1","t2","t3")
dataRaw <- mxData(observed=Demo.growth,type="raw")
varPaths <- mxPath(from=mvars,arrows=2,free=T,values=rep(1,5),labels=paste("r",mvars,sep="_"))
regPaths1 <- mxPath(from=mvars[1:2],to="t1",arrows=1,free=T,values=c(1,1)    ,labels=c("b_x1t1","b_x2t1"))
regPaths2 <- mxPath(from=mvars[1:2],to="t2",arrows=1,free=T,values=c(1,1)    ,labels=c("b_x1t2","b_x2t2"))
regPaths3 <- mxPath(from=mvars[1:4],to="t3",arrows=1,free=T,values=c(1,1,1,1),labels=c("b_x1t3","b_x2t3","b_c1t3","b_t2t3"))
covPaths1 <- mxPath(from=mvars[1],to=mvars[2],arrows=2,free=T,values=0.5,labels="c_x1x2")
covPaths2 <- mxPath(from=mvars[3],to=mvars[4],arrows=2,free=T,values=0.5,labels="c_t1t2")
means    <- mxPath(from="one",to=mvars,arrows=1,free=T,values=rep(1,5),labels=paste("m",mvars,sep="_"))
uniRegMod <- mxModel(model="Demo",type="RAM",dataRaw,manifestVars=mvars,varPaths,regPaths1,regPaths2,regPaths3,covPaths1,covPaths2,means)
uniRegFit <- mxRun(uniRegMod)
summary(uniRegFit)
mxStandardizeRAMpaths(uniRegFit,SE=T) -> stdP

# Models are identical; AIC, BIC, logLilkelihood match

# all parameters match, unstandardized AND standardized (x1~x2 and t1~t2 not picked up here due to different order of row/col)
#
merge(stdl,stdP,by.x=c("lhs","rhs"),by.y=c("row","col")) -> stdboth
stdboth$discrep <- stdboth$est.std - stdboth$Std.Value 
stdboth 
#
merge(unsl,stdP,by.x=c("lhs","rhs"),by.y=c("row","col")) -> unsboth
unsboth$discrep <- unsboth$est - unsboth$Raw.Value 
unsboth 

# MATCH (standardized covariance of exogenous residuals)
stdl[stdl$lhs=="x1" & stdl$rhs=="x2",]
stdP[stdP$col=="x1" & stdP$row=="x2",]

# MISMATCH (standardized covariance of endogenous residuals)
stdl[stdl$lhs=="t1" & stdl$rhs=="t2",]
stdP[stdP$col=="t1" & stdP$row=="t2",]