M1 STAGE 1 (x=v1, m=v2, y=v3)

data <- read.table(file = "Data.dat", header = TRUE)

head(data)

nvar <- 3
varnames <- c("x","m","y") 
labels <- list(varnames,varnames)

cordat <- list() 
for (i in 1:nrow(data)){   
     cordat[[i]] <- vec2symMat(as.matrix(data[i,3:5]), 
     diag = FALSE) 
     dimnames(cordat[[i]]) <- labels}
for (i in 1 :length (cordat)){
for (j in 1 :nrow (cordat[[i]])){
if ( sum ( is.na (cordat[[i]][j,])) == nvar-1)
{cordat[[i]][j,j] <- NA}}}
for (i in 1 :length (cordat)){
for (j in 1 :nrow (cordat[[i]])){
for (k in 1 : nvar){
if ( is.na (cordat[[i]][j,k]) == TRUE
&is.na (cordat[[i]][j,j]) != TRUE
&is.na (cordat[[i]][k,k]) != TRUE){
if ( sum ( is.na (cordat[[i]])[j,]) >sum ( is.na (cordat[[i]])[k,]))
{cordat[[i]][k,k] <- NA}
if ( sum ( is.na (cordat[[i]])[j,]) <=sum ( is.na (cordat[[i]])[k,]))
{cordat[[i]][j,j] <- NA}
}}}}

stage1random <- tssem1 (Cov=cordat, n=data$N, method="REM")
stage1random <- rerun (stage1random, autofixtau2 = TRUE)
summary (stage1random)

M1 STAGE 2

A<- create.mxMatrix ( c (0,0,0,
"0.1*b_mx",0,0,
"0.1*b_yx","0.1*b_ym",0),
type = "Full",
nrow = 3,
ncol = 3,
byrow = TRUE,
name = "A",
dimnames = list (varnames,varnames))

dimnames(A)[[1]] <- dimnames(A)[[2]] <- c("x","m","y")
A

S <- create.mxMatrix (
c (1,
0, ".1*p_m",
0,0,".1*p_y"),
type="Symm", byrow = TRUE, name="S",
dimnames = list (varnames,varnames))

dimnames(S)[[1]] <- dimnames(S)[[2]] <- c("x","m","y")
S

stage2 <- tssem2 (stage1random, Amatrix=A, Smatrix=S,
diag.constraints=FALSE, intervals="LB",
mx.algebras=list(Ind=mxAlgebra(b_mx* b_ym,name="Ind")))

summary(stage2)
mxEval(Smatrix, stage2$mx.fit)