#code obtained from: #tutorial #https://cran.r-project.org/web/packages/metaSEM/vignettes/Examples.html#two-stage-structural-equation-modeling-tssem # download datasets from github: https://github.com/mikewlcheung/metasem/tree/master/data #Helpful walkthrough - Mathias Harrer example: resilience and depression #https://bookdown.org/MathiasHarrer/Doing_Meta_Analysis_in_R/mediation.html #------------- #PACKAGE INFO install.packages('metaSEM') install.packages("diagram") #TSSEM functions # stage 1 = tssem1() # method = "FEM" (fixed effects) # method = "REM" (random effects) # categorical variables: use argument 'cluster' to perform fixed-effects TSSEM on each group # stage 2 = tssem2() # automatically handles whether fixed or random effects model is used #!! when a pooled correlation matrix is used IT IS IMPORTANT TO ENSURE THAT THE DIAGONALS of the model-implied correlation matrix are fixed at 1 # diag.constraint = TRUE (imposes this nonlinear constraint) # in this case, SEs are not reported in OpenMx since they are not accurate, get likelihood based confidence intervals instead # intervals = "LB" (likelihood-based CIs) # Reticular Action Model (RAM) formulation, which specifies 3 matrices, A,S,F: # 1. A: asymmetric path (regression coefficients) from independent variables to the dependent variables # 2. S: symmmetric paths (variances and covariances) of variables # 3. F: used to select the two observed variables library(metaSEM) library(diagram) library(readr) #----------------------------------------------------------------------------------------- #Organise data #----------------------------------------------------------------------------------------- rm(list=ls()) #source("./TSSEM_file_info.R") datapath <- "~/path_to_data" setwd(datapath) #med_var <- read_csv("./DBAS_extracted_data_130320_singlemediators.csv",skip_empty_rows = TRUE) med_var <- read_csv("./DBAS_extracted_data_130320_singlemediators-noEidelman.csv",skip_empty_rows = TRUE) #med_var <- read_csv("./Arousal_extracted_data_130320_singlemediator.csv",skip_empty_rows = TRUE) #med_var <- read_csv("./Arousal_extracted_data_290720_multimediator.csv",skip_empty_rows = TRUE) mediator <- "Change in DBAS" 'mediator <- "Change in arousal"' outcome <- "Change in sleep outcome" med_var <- med_var[rowSums(is.na(med_var)) != ncol(med_var),] #delete empty rows med_var$data <- as.matrix(med_var[c(2:4)]) #read just A, B, C path columns #assign data structure med_var$data <- lapply(split(med_var$data,1:nrow(med_var)), function(x) {mat <- matrix(1, ncol=3, nrow=3); mat[upper.tri(mat, diag=FALSE)] <- x; mat[lower.tri(mat)] <- t(mat)[lower.tri(mat)]; mat}) #add variable names med_var$data <- lapply(med_var$data, function(x, var.names) {dimnames(x) <- list(var.names, var.names); x}, var.names=c("Treatment", "Outcome", "Mediator")) #add study names names(med_var$data) <- as.matrix(med_var[c(1)]) #add sample sizes med_var$n <- as.matrix(med_var[c(5)]) #----------------------------------------------------------------------------------------- #STAGE 1 TSSEM #----------------------------------------------------------------------------------------- med1 <- tssem1(med_var$data, med_var$n, method = "REM") med1 <- rerun(med1) #option to rerun summary(med1) #The pooled correlation coefficients (fixed effects) and the variance components #(the random effects) can be extracted by the use of the coef() command via the #select argument. ffx <- coef(med1,"fixed") rfx <- coef(med1,"random") #----------------------------------------------------------------------------------------- #STAGE 2 TSSEM #----------------------------------------------------------------------------------------- # Build "A" matrix A <- matrix(c(0 , 0 , 0 , "0.2*Tx_Out", 0 , 0 , "0.2*Tx_Med" , "0.2*Med_Out", 0), ncol = 3, nrow=3, byrow=TRUE) # Set column and row labels dimnames(A)[[1]] <- dimnames(A)[[2]] <- c("Treatment", "Outcome", "Mediator") A <- as.mxMatrix(A) #We use starting values of 0.1 # Build "S" matrix S <- Diag(c(1, "0.1*ErrVarE", "0.1*ErrVarD")) # Set column and row labels dimnames(S)[[1]] <- dimnames(S)[[2]] <- c("Treatment", outcome, mediator) ##ADDED CODE to convert path coefficients using implied R>>>>>>>> A1_impliedR <- impliedR(Amatrix = A$values, Smatrix = S) ##SUBSTITUTING 'b' and 'c' path estimates A$values[3] <- A1_impliedR$SigmaObs[3] A$values[6] <- A1_impliedR$SigmaObs[6] ##>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> #Model fitting med2 <- tssem2(med1, Amatrix = A, Smatrix = S, intervals.type = "LB", diag.constraints = TRUE, mx.algebras=list(indirect=mxAlgebra(Amatrix[3,2]*Amatrix[3,1],name="indirect"), direct=mxAlgebra(Amatrix[2,1],name="direct"), total=mxAlgebra(direct+indirect,name="total"))) # Rerun med2 <- rerun(med2) summary(med2) med_var_coef <- coef(med2) #save coefficients to structure coef(med2,random) #----------------------------------------------------------------------------------------- #PLOT #----------------------------------------------------------------------------------------- coef1 <- toString(round(med_var_coef[c(2)],digits=2)) coef2 <- toString(round(med_var_coef[c(3)],digits=2)) coef3 <- toString(round(med_var_coef[c(1)],digits=2)) data <- c(0, coef1, 0, 0, 0, 0, coef2, coef3, 0) M<- matrix (nrow=3, ncol=3, byrow = TRUE, data=data) plot<- plotmat (M, pos=c(1,2), name= c(mediator, "Treatment", outcome), box.type = "rect", box.size = 0.12, box.prop=0.5, curve=0, shadow.size = 0, dtext = .5, prefix = "") c_prime <- med_var_coef[c(2)] - med_var_coef[c(3)]*med_var_coef[c(1)]