install.packages("metaSEM")
require("metaSEM")
install.packages("readxl")
library(readxl)
dataset <- read_excel("c:/users/use/Documents/R/paper6_2/paper6_2 data/masem_data_add.xlsx")
head(dataset)
## make list of cormatrices (cordat), NA on diagonal
nvar <- 8
varnames <- c("RND","CEO","ABS","COM","TEC","ENV","REL","PER")
labels <- list(varnames, varnames)
cordat <- list()
for (i in 1:nrow(dataset)){	
 cordat[[i]] <- vec2symMat(as.matrix(dataset[i,3:30]),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}
   }}
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}	
       }}}}
head(cordat)
dataset$N

#### stage 1 
stage1random <- tssem1(Cov=cordat, n=dataset$N, method="REM", RE.type="Diag")

# pattern of corelation matrices
pattern.na(cordat,show.na=FALSE)

summary(stage1random)

# returns models via maTryHard()
random1_rerun <- rerun(stage1random)



### stage 2
require("OpenMx")
library(OpenMx)
A <- create.mxMatrix(
 c( 0,"0.1*CEO2RND",0,0,0, "0.1*ENV2RND",0,0,
    0,0,0,0,0,0,0,0,
    0,"0.1*CEO2ABS",0,0,0,0,"0.1*REL2ABS",0,
    0,"0.1*CEO2COM",0,0,"0.1*TEC2COM","0.1*ENV2COM","0.1*REL2COM",0,
    "0.1*RND2TEC",0,"0.1*ABS2TEC",0,0,0,0,0,
    0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,
    0,0,0,"0.1*COM2PER",0,0,0,0),
   type = "Full",
   nrow = 8,
   ncol = 8,
   byrow = TRUE,
   name = "A",
   dimnames = list(varnames,varnames))
dimnames(A)[[1]] <-dimnames(A)[[2]] <- c("RND","CEO","ABS","COM","TEC","ENV","REL","PER")
A

S <- create.mxMatrix(
 c("1*var_RND",
   0,1,
   0,0,"1*var_ABS",
   0,0,0,"1*var_COM",
   0,0,0,0,"1*var_TEC",
   0,"0.1*cor_CEOENV",0,0,0,1,
   0,"0.1*cor_CEOREL",0,0,0,"0.1*cor_ENVREL",1,
   0,0,0,0,0,0,0,"1*var_PER"),
 type = "Symm",
 byrow = TRUE, 
 name="S",
 dimnames = list(varnames,varnames))
S

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

summary(stage2)