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

nvar <- 6 
varnames <- c("v1 ","v2","v3","v4","v5","v6") 
labels <- list(varnames,varnames)
cordat <- list()
for (i in 1:nrow(data)){   
     cordat[[i]] <- vec2symMat(as.matrix(data[i,3:17]), 
     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} 
}}

cordat[[13]]

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}   
}}}}

stage1fixed  <- tssem1(Cov=cordat, n=data$N, method="REM", RE.type="Zero")
summary(stage1fixed)

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

A<-create.mxMatrix (c(0,0,0,0,0,0
,0,0,0,0,0,0,
0,0,0,0,0,0,
"0.1*b41","0.1*b42",0,0,0,0,
0,0,"0.1*b53","0.1*b54",0,0,
0,0,"0.1*b63","0.1*b64",0,0),
type = "Full", 
      nrow = 6, 
      ncol = 6, 
      byrow = TRUE, 
      name = "A", 
      dimnames = list(varnames,varnames))

S <- create.mxMatrix( 
c(1, 
".1*p21",1,
".1*p31",".1*p32",1, 
0,0,0,"1*p44",
0,0,0,0,"1*p55",
0,0,0,0,0,"1*p66"),  
type="Symm", byrow = TRUE, name="S", 
      dimnames = list(varnames,varnames))

stage2 <- tssem2(stage1random, Amatrix=A, Smatrix=S, 
                 diag.constraints=TRUE, intervals="LB")

summary(stage2)