rm(list=ls())  # clearing memory

## Load the library
library("metaSEM")


##	create names for matrix rows/columns
mtrxnames <- c("AT" , "SN" , "PBC" , "IN" , "PA")


godin1986 <- as.matrix(read.csv("godin1986.csv", header=FALSE))
rownames(godin1986) <-mtrxnames
colnames(godin1986) <-mtrxnames

latimer2005 <- as.matrix(read.csv("latimer2005.csv", header=FALSE))
rownames(latimer2005) <-mtrxnames
colnames(latimer2005) <-mtrxnames

Kosma2007 <- as.matrix(read.csv("Kosma2007.csv", header=FALSE))
rownames(Kosma2007) <-mtrxnames
colnames(Kosma2007) <-mtrxnames

Ellis2009 <- as.matrix(read.csv("Ellis2009.csv", header=FALSE))
rownames(Ellis2009) <-mtrxnames
colnames(Ellis2009) <-mtrxnames

latimer2004paramild <- as.matrix(read.csv("latimer2004paramild.csv", header=FALSE))
rownames(latimer2004paramild) <-mtrxnames
colnames(latimer2004paramild) <-mtrxnames

latimer2004paramoderate <- as.matrix(read.csv("latimer2004paramoderate.csv", header=FALSE))
rownames(latimer2004paramoderate) <-mtrxnames
colnames(latimer2004paramoderate) <-mtrxnames

latimer2004paraintense <- as.matrix(read.csv("latimer2004paraintense.csv", header=FALSE))
rownames(latimer2004paraintense) <-mtrxnames
colnames(latimer2004paraintense) <-mtrxnames

latimer2004tetramild <- as.matrix(read.csv("latimer2004tetramild.csv", header=FALSE))
rownames(latimer2004tetramild) <-mtrxnames
colnames(latimer2004tetramild) <-mtrxnames

latimer2004tetramoderate <- as.matrix(read.csv("latimer2004tetramoderate.csv", header=FALSE))
rownames(latimer2004tetramoderate) <-mtrxnames
colnames(latimer2004tetramoderate) <-mtrxnames

latimer2004tetraintense <- as.matrix(read.csv("latimer2004tetraintense.csv", header=FALSE))
rownames(latimer2004tetraintense) <-mtrxnames
colnames(latimer2004tetraintense) <-mtrxnames

Haegele2017 <- as.matrix(read.csv("Haegele2017.csv", header=FALSE))
rownames(Haegele2017) <-mtrxnames
colnames(Haegele2017) <-mtrxnames

Kosma2009hr <- as.matrix(read.csv("Kosma2009hr.csv", header=FALSE))
rownames(Kosma2009hr) <-mtrxnames
colnames(Kosma2009hr) <-mtrxnames

Kosma2009MET <- as.matrix(read.csv("Kosma2009MET.csv", header=FALSE))
rownames(Kosma2009MET) <-mtrxnames
colnames(Kosma2009MET) <-mtrxnames

Arbour2010 <- as.matrix(read.csv("Arbour2010.csv", header=FALSE))
rownames(Arbour2010) <-mtrxnames
colnames(Arbour2010) <-mtrxnames


#	14 studies
dataL <- list(godin1986 = godin1986  ,  latimer2005 = latimer2005 , Kosma2007 = Kosma2007 , Ellis2009 = Ellis2009 , latimer2004paramild = latimer2004paramild , latimer2004paramoderate = latimer2004paramoderate , latimer2004paraintense = latimer2004paraintense, latimer2004tetramild = latimer2004tetramild, latimer2004tetramoderate = latimer2004tetramoderate, latimer2004tetraintense = latimer2004tetraintense , Haegele2017 = Haegele2017 , Kosma2009hr = Kosma2009hr , Kosma2009MET = Kosma2009MET , Arbour2010 = Arbour2010)

nL <- c(36, 104, 223, 223, 54, 54, 54, 70, 70, 70, 209, 141, 141, 574 )



head(dataL)	
head(godin1986)
head(latimer2005)
head(Kosma2007)
head(Ellis2009)
head(latimer2004paramild)
head(latimer2004paramoderate)
head(latimer2004paraintense)
head(latimer2004tetramild)
head(latimer2004tetramoderate)
head(latimer2004tetraintense)
head(Haegele2017)
head(Kosma2009hr)
head(Kosma2009MET)
head(Arbour2010)

#  test for positive definite matrix
is.pd(dataL)

my.df<-Cor2DataFrame(dataL, nL, acov = "weighted")
pattern.na(my.df, show.na=FALSE, type="osmasem")


model4 <- 'IN ~ b41*AT + b42*SN + b43*PBC
           PA ~ b54*IN + b53*PBC 
           AT ~~ p21*SN
           AT ~~ p31*PBC
           SN ~~ p32*PBC
           AT ~~ 1*AT
           SN ~~ 1*SN
           PBC ~~ 1*PBC
           IN ~~ p44*IN
           PA ~~ p44*PA'

plot(model4, color="yellow") 

#Creating A1 matrix for 1-factor
A1  mtrxnames1 <- c("AT" , "SN" , "PBC" , "IN" , "PA")
A1 <- create.mxMatrix(c(0, 0, 0, 0, 0,
                        		0, 0, 0, 0, 0,
				0, 0, 0, 0, 0, 
                        		"0*b41", "0*b42", "0*b43", 0, 0,
                        		0, 0, "0*b53", "0*b54", 0),
                      type="Full", byrow=TRUE, ncol=5, nrow=5,
                      as.mxMatrix=FALSE)
dimnames(A1)[[1]] <- dimnames(A1)[[2]] <- c("AT","SN","PBC", "IN","PA")
A1

#Creating S1 matrix
S1 <- create.mxMatrix(c(1, 
                       	        "0*p21", 1, 
                        	         "0*p31", "0*p32", 1,
                       	         0, 0, 0, "0*p44",
                       	         0, 0, 0, 0, "0*p44"),
                      byrow=TRUE, type="Symm", as.mxMatrix=FALSE)
dimnames(S1)[[1]] <- dimnames(S1)[[2]] <- c("AT","SN","PBC", "IN","PA")
S1


#Creating F1 matrix for 1-factor
F1 <- create.Fmatrix(c(1, 1, 1, 1, 1, 0), as.mxMatrix=FALSE)
dimnames(F1) <- list(c("AT","SN","PBC","IN","PA"))
F1


RAM4 <- lavaan2RAM(model4, obs.variables=c("AT", "SN", "PBC", "IN", "PA"))
RAM4

## Create mean structure
M0 <- create.vechsR(A0=RAM4$A, S0=RAM4$S)

## Create heterogeneity variances
T0 <- create.Tau2(RAM=RAM4, RE.type="Diag")

mx.fit0 <- osmasem(model.name="No moderator", Mmatrix=M0, Tmatrix=T0, data=my.df)

## Get the chi-square statistic on the proposed model
summary(mx.fit0, Saturated=TRUE)