##Install Packages to Run metaSEM Use 32 bit for R
install.packages('OpenMx', repos='http://openmx.psyc.virginia.edu/packages/')
install.packages(c('ellipse','MASS'))
install.packages(pkgs="e:/metaSEM_0.8-4.zip", repos=NULL)

##Run Example from Mike Cheung's Webpage - Code from http://courses.nus.edu.sg/course/psycwlm/internet/metaSEM/index.html
## Load the metaSEM library
library(metaSEM)

## LOAD DATA FILE INTO CONSOLE BY CLICKING AND DRAGGING FILE INTO CONSOLE

## READ DATA FILE AT THIS LOCATION
PUEUBIdata <- readFullMat(file = "C:\\Users\\Nicola\\Documents\\Dissertation_Process\\Proposal\\Article2\\Methods\\Analysis\\MetaAnalysis\\PUEUBI\\PUEUBIdata.dat")

## SEE DATA NAMED "PUEUBIdata"
PUEUBIdata

##WRITE MATRIX FOR SAMPLE SIZE OF EACH STUDY- Code modified from: http://www.r-tutor.com/r-introduction/matrix
PUEUBIn = matrix (c(425, 243, 102, 362, 109), nrow=1, ncol=5, byrow = TRUE)

## SEE MATRIX FOR SAMPLE SIZE OF EACH STUDY
PUEUBIn

##FIXED EFFECTS MASEM USING TSSEM - UNDER THE FIXED EFFECTS MODEL, IT IS ASSUMED THAT ALL POPULATION COVARIANCE MATRICES ARE THE SAME WHILE THERE ARE STUDY SPECIFIC COVARIANCE MATRICES UNDER THE RANDOM-EFFECTS MODEL. 
## TSSEM: STAGE 1 -  The tssem1() function is used to pool the correlation matrices with a fixed-effects model in the first stage analysis by specifying method="FEM" in the argument:
fixed1 <- tssem1(PUEUBIdata, PUEUBIn, method = "FEM")

## AFTER TSSEM1 RUNS ANALYSIS OF CORRELATION MATRIX (IN THIS CASE THE COVARIANCE MATRIX), OUTPUT RESULTS
summary(fixed1)

## OUTPUTS POOLED COVARIANCE MATRIX (THE PARAMETER ESTIMATES) EXTRACTED
coef(fixed1)

##RANDOM EFFECTS MASEM USING TSSEM
random1 <- tssem1(PUEUBIdata, PUEUBIn, method = "REM", RE.type = "Diag")
summary(random1)

## Drop variable x3 (modified syntax from: http://www.statmethods.net/management/subset.html)
newdata <- PUEUBIdata[c(-3)]

## CREATES "A MATRIX" WHICH SPECIFIES THE ASYMMETRIC PATHS (ASYMTOTIC MATRIX)
## OPTION 1 - ONE FACTOR
Lambda <- matrix(c(".3*Alpha_x1", ".3*Alpha_x2", ".3*Alpha_x3",rep(0,3)), ncol = 1, nrow = 3)
A1 <- rbind( cbind(matrix(0, ncol = 3, nrow = 3), Lambda), matrix(0, ncol = 4, nrow = 1) )
dimnames(A1) <- list(c("x1", "x2", "x3", "Alpha"), c("x1", "x2", "x3", "Alpha"))
A1
A1 <- as.mxMatrix(A1)

## CREATES "S MATRIX" WHICH SPECIFIES THE SYMMETRIC VARIANCE/COVARIANCE MATRIX
Phi <- matrix(c(1), ncol=1, nrow=1)
S1 <- bdiagMat(list(Diag(c(".3*e1", ".3*e2", ".3*e3")), Phi))
dimnames(S1) <- list(c("x1", "x2", "x3", "Alpha"), c("x1", "x2", "x3", "Alpha"))
S1
S1 <- as.mxMatrix(S1)

## CREATES "F MATRIX" WHICH SPECIFIES A SELECTION MATRIX USED TO OBSERVED VARIABLES
F1 <- create.Fmatrix(c(1, 1, 1, 0), name = "F1")

## OPTION 2 - THREE FACTORS (WITH NO PATHS BETWEEN THREE FACTORS)
Lambda <- create.mxMatrix(c ("0", "1*Alpha_x2", "0", ".5*Beta_x1", ".5*Beta_x2", "0"), ncol=2, nrow=3, as.mxMatrix=FALSE)
A1 <- rbind( cbind(matrix(0, ncol = 3, nrow = 3), Lambda), matrix(0,  ncol = 5, nrow = 2) )
dimnames(A1) <- list(c("x1", "x2", "x3", "Alpha", "Beta"), c("x1", "x2", "x3", "Alpha", "Beta"))
A1
A1 <- as.mxMatrix(A1)

## CREATES "S MATRIX" WHICH SPECIFIES THE SYMMETRIC VARIANCE/COVARIANCE MATRIX
Phi <- matrix(c("1", "0", "0", "1"), ncol=2, nrow=2)
S1 <- bdiagMat(list(Diag(c(".3*e1", ".3*e2", ".3*e3")), Phi))
dimnames(S1) <- list(c("x1", "x2", "x3", "Alpha", "Beta"), c("x1", "x2", "x3", "Alpha", "Beta"))
S1
S1 <- as.mxMatrix(S1)

## CREATES "F MATRIX" WHICH SPECIFIES A SELECTION MATRIX USED TO OBSERVED VARIABLES
F1 <- create.Fmatrix(c(1, 1, 1, 0, 0), name = "F1")

## TSSEM: STAGE 2 
random2 <- tssem2(random1, Amatrix = A1, Smatrix = S1, Fmatrix = F1, model.name = "TSSEM2 PUEUBI")

## AFTER TSSEM2 RUNS ANALYSIS, OUTPUT RESULTS 
summary(random2)

#### SUB-GROUP ANALYSIS: Fixed-effects model: with LMS as cluster - Modified from http://openmx.psyc.virginia.edu/thread/2360
## SET UP CLUSTER
LMS <- c("Proprietary", "Proprietary", "Self", "Self", "Proprietary")

## First stage analysis
fixed1.cluster <- tssem1(PUEUBIdata, PUEUBIn, method = "FEM", cluster=LMS)
summary(fixed1.cluster)

## Second stage analysis
fixed2.cluster <- tssem2(fixed1.cluster, Amatrix = A1, Smatrix = S1, Fmatrix = F1)
summary(fixed2.cluster)