#### Approach 3 - MODEL 3
#### Estimate a Random-Effects pooled correlation matrix with with those variables present in the research model
#### Fit the Random-Effects stage 2 model using this pooled correlation matrix 

### Step 1: Import Libraries
library("OpenMx")
library("lavaan")
library("metaSEM")
library("semPlot")

## Step 2: Load Data.  
## DATA FILE FOR MODEL 3
data3 <- read.csv("/Users/srikanthparameswaran/Desktop/TSSEM Project 2/Approach 3/Model 3 Data.csv")

### Step 3: Setup the correlation matrices
nvar <- 9
varnames <- c("X","Y","A","B","C","M","D","E","F")
labels <- list(varnames,varnames)
cordat <- list()

for (i in 1:nrow(data3)){
cordat[[i]] <- vec2symMat(as.matrix(data3[i,2:37]),diag = FALSE)
dimnames(cordat[[i]]) <- labels}

data3$data<-cordat

### Step 4: Put NA on diagonal if variable is missing

for (i in 1:length(data3$data)){
    for (j in 1:nrow(data3$data[[i]])){ 
        if (sum(is.na(data3$data[[i]][j,]))==nvar-1) 
         {data3$data[[i]][j,j] <- NA}
    }}

### Step 5: Put NA on diagonal for variable with least correlations

for (i in 1:length(data3$data)){
for (j in 1:nrow(data3$data[[i]])){
    for (k in 1:nvar){              
       if (is.na(data3$data[[i]][j,k])==TRUE
          &is.na(data3$data[[i]][j,j])!=TRUE
          &is.na(data3$data[[i]][k,k])!=TRUE){

if(sum(is.na(data3$data[[i]])[j,])>sum(is.na(data3$data[[i]])[k,]))
{data3$data[[i]][k,k] = NA}
if(sum(is.na(data3$data[[i]])[j,])<=sum(is.na(data3$data[[i]])[k,]))
{data3$data[[i]][j,j] = NA} 
}}}}

# Display processed data
data3$data

pattern.na(data3$data, show.na = FALSE)
pattern.n(data3$data, data3$Corr_Sample)

### Step 6: Check Positive-Definiteness. 
### Everything should be TRUE because I am using the already processed file.
is.pd(data3$data)
length(is.pd(data3$data))

### Step 7: Random Effects Stage 1 Estimation
random_stage1<- tssem1(data3$data, data3$Corr_Sample, method="REM", RE.type="Diag")
summary(random_stage1)
vec2symMat(coef(random_stage1, select="fixed"), diag = FALSE)

pooled <- vec2symMat(coef(random_stage1, select="fixed"), diag = FALSE)
dimnames(pooled)<- list(varnames,varnames)
pooled

### Step 8: Random Effects Stage 2

### Step 8a: Model setup in Lavaan. Converting Lavaan to RAM
model1 <- 'X ~  A2X*A + M2X*M + B2X*B + C2X*C + D2X*D + E2X*E + F2X*F
	   M ~  A2M*A + B2M*B
	   Y ~  A2Y*A + M2Y*M + B2Y*B + C2Y*C + D2Y*D + E2Y*E + F2Y*F + X2Y*X'

plot(model3, col="yellow")

RAM1 <- lavaan2RAM(model1, obs.variables=c("X","Y","A","B","C","M","D","E","F"))

RAM1

RAM1$S 

## Step 8b: All variances except the DVs were set to 1. 
RAM1$S[3,3] <- 1 
RAM1$S[4,4] <- 1  
RAM1$S[5,5] <- 1   
RAM1$S[7,7] <- 1 
RAM1$S[8,8] <- 1 
RAM1$S[9,9] <- 1 

RAM1$S

RAM1

### Step 8c: Random Effects Stage 2 Model estimation

random_stage2 <- tssem2(random_stage1, Amatrix=RAM1$A, Smatrix=RAM1$S )
summary(random_stage2)

my.plot1 <- meta2semPlot(random_stage2)
semPaths(my.plot1, whatLabels="est", layout = "tree2", sizeMan=12, edge.label.cex=2,
color="yellow", edge.color = "black", weighted=FALSE)

(Smatrix <- diag(mxEval(Smatrix, random_stage2$mx.fit)))
1 - Smatrix