#load metaSEM and semPlot
library(metaSEM)
library(semPlot)

#set working directory(library)
setwd ("D:/Research/MY MS/4 In Prep/Study 1/Analyses/Data/Aug 2018/CME")

#read a text file, enter working file name and location of external file to load
#csv doesn't seem to be read correctly as matrix
vector <- read.table ("v_pre_posti.txt")
vector
#read vector form, ie., corr in a straight line 1 row per matrix
vector <- readStackVec ("v_pre_posti.txt")
vector


#give column and row headings, enter columns first, check vector
  vector <- lapply(vector,function(x)
  {dimnames(x) <- list(c("Tx", "m1", "o1", "m3", "o3"),
                       c("Tx", "m1", "o1", "m3", "o3"))
   x})
  vector 

#input ns
n <- read.table("vectormaxn.txt")
n

## Show the missing data: show.na=TRUE
pattern.na(vector, show.na=FALSE)
## Display the accumulative sample sizes for each correlation
pattern.n(vector, n)

#test if matrix is positive definite
is.pd(vector)

random1 <- tssem1(vector, n, method="REM", RE.type="Diag")
#summary of REM
options("scipen"=100, "digits"=20)
summary(random1)
#OpenMx status = 0 at first run

#Rerun analysis
#random1 <-rerun(random1)
#summary(random1)

#To extract estimated ave corr matrix in matrix form
## Select the fixed effects and convert it into a correlation matrix
vec2symMat(coef(random1, select="fixed"),diag=FALSE)

#Create A matrix for fit to model
A2 <- create.mxMatrix(c(0, 0,  0, 0, 0,
                        0, 0, 0, 0, 0, 
                        0, 0, 0, 0, 0, 
                        "0.1*Txm3", "0.1*m1m3", 0, 0, 0, 
                        "0.1*Txo3", 0, "0.1*o1o3", "0.1*m3o3", 0),
                        type="Full", byrow=TRUE, ncol=5, nrow=5, as.mxMatrix=FALSE)
## Label matrix for inspecting the model.
dimnames(A2)[[1]] <- dimnames(A2)[[2]] <- c("Tx", "m1", "o1", "m3", "o3") 
A2

#Create S matrix-symmetric variance covariance matrix
S1 <- create.mxMatrix(c(1, 
                        0, 1,
                        0, "0.1*r_m1o1", 1,
                        0, 0, 0, "0.1*var_m3",
                        0, 0, 0, 0, "0.1*var_o3"), 
                      byrow=TRUE, type="Symm", as.mxMatrix=FALSE)
## label matrix variables, useful for inspecting the model.
dimnames(S1)[[1]] <- dimnames(S1)[[2]] <- c("Tx", "m1", "o1", "m3", "o3")
S1

## Second stage analysis
random2 <- tssem2(random1, Amatrix=A2, Smatrix=S1, diag.constraints=FALSE,
                  intervals.type="LB", model.name="TSSEM2 test",
                  mx.algebras=list( medE=mxAlgebra(Txm3*m3o3, name="medE"),
                                    dirE=mxAlgebra(Txo3, name="dirE"), 
                                    totE=mxAlgebra(Txo3 + Txm3*m3o3, name="totE") ))

summary(random2)
#openmx 0 at first run, no NA bounds

#rerun
#random2 <- rerun(random2)
#summary(random2)

#Plot the model and label the parameters for checking
## Convert the model to semPlotModel object
my.plot <- meta2semPlot(random2, manNames=c("Tx", "m1", "o1", "m3", "o3"))
                        
## Plot the model with labels
## The labels are overlapped. We may modify it by using layout="spring"
semPaths(my.plot, whatLabels="path", nCharEdges=10, nCharNodes=10,
         layout="spring", color="yellow", edge.label.cex=0.8)
#Plot the parameter estimates in the figure
## Plot the parameter estimates
semPaths(my.plot, whatLabels="est", nCharNodes=10, layout="spring",
         color="green", edge.label.cex=1.2)