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

#Load other packages
library(tidyverse)
library(readxl)
library(skimr)

#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
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_post.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 present data: show.na=FALSE
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)

#Import dataset containing cluster variable
#NOTE only reads sheet 1, not subsequent sheets
mod <- read_excel("mod.xlsx")
mod

#Recode moderators/clusters
mod1 <- mutate(mod, 
                 Diagnosis1 = recode(Diagnosis, "1"="Dx Needed", "0"="Dx Not Needed"),
                 Child1 = recode(Child, "1"="Child", "0"="Adolescent"),
                 IPT1 = recode(IPT, "1"="IPT", "0"="CBT"))
#### Show the clusters
mod1

#Run fixed-effect TSSEM1 according to cluster labels
cluster1 <- tssem1(vector, n, method="FEM",
                   cluster=mod$IPT1)
#Run fixed-effect TSSEM1 according to dichotomous variable without labels
cluster1 <- tssem1(vector, n, method="FEM",
                   cluster=mod$IPT)

#Here is where I get error msg
sessionInfo()




summary(cluster1)


#Extract pooled corr matrix for each cluster
coef(cluster1)

#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

#Run fixed-effect TSSEM2 according to cluster labels
cluster2 <- tssem2(cluster1, Amatrix=A1, Smatrix=S1, diag.constraints=FALSE, intervals.type="LB", 
                  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(cluster2)

my.plots <- meta2semPlot(cluster2, manNames=c("Tx", "m1", "o1", "m3", "o3"))
## Setup two plots
layout(t(1:2))
semPaths(my.plots[[1]], whatLabels="est", nCharNodes=10, color="green")
title("IPT")
semPaths(my.plots[[2]], whatLabels="est", nCharNodes=10, color="yellow")
title("CBT")