# Load the data load("data for OpenMX.rda") load("n for OpenMX.rda") load("grade for OpenMX.rda") # Library library(metaSEM) library(beepr) #Reminder to check if the process is complete. mxOption(NULL, 'Number of Threads', parallel::detectCores()-2) # Calculate the sampling covariance matrix of the correlations my.df <- Cor2DataFrame(data, n, acov = "weighted") # Adding a moderator my.df$data <- data.frame(my.df$data, grade=grade, check.names=FALSE) # check.names has to be False # Proposed model model0 <- "## Factor loadings DECOD =~ d1*FLUEN + d2*ACCU LANG =~ l1*VOC + l2*LC META =~ m1*PA + m2*RAN + m3*MOR + m4*OTH ## Factor correlations META ~~ metaWITHlang*LANG LANG ~~ langWITHdecod*DECOD ## RC is modeled by LANG and DECOD RC ~ lang2rc*LANG + decod2rc*DECOD ## RC is modeled by META DECOD ~ meta2decod*META ## Errors FLUEN ~~ FLUEN_err*FLUEN ACCU ~~ ACCU_err*ACCU VOC ~~ VOC_err*VOC LC ~~ LL_err*LC RC ~~ RC_err*RC PA ~~ PA_err*PA RAN ~~ RAN_err*RAN MOR ~~ MOR_err*MOR OTH ~~ OTH_err*OTH DECOD ~~ 1*DECOD LANG ~~ 1*LANG META ~~ 1*META" RAM0 <- lavaan2RAM(model0, obs.variables=c("FLUEN","ACCU","VOC","LC","RC","PA","RAN","MOR","OTH")) ## Create heterogeneity variances T0 <- create.Tau2(RAM=RAM0, RE.type="Diag", Transform="expLog", RE.startvalues=0.05) ## With no moderator M0 <- create.vechsR(A0=RAM0$A, S0=RAM0$S, F0=RAM0$F) fit0 <- osmasem(model.name="With no moderator", Mmatrix=M0, Tmatrix=T0, data=my.df) summary(fit0) plot(fit0, color="green") beep(1) ## 1st step: Replace the A matrix with the moderator "grade" Ax <- matrix(c(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,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,0,0,0,0,0,0,0,"0*data.grade","0*data.grade",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, 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,0,0,0,0,0,0,0,0,0,0,"0*data.grade", 0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0), byrow=TRUE, nrow=12, ncol=12) M1 <- create.vechsR(A0=RAM0$A, S0=RAM0$S, F0=RAM0$F, Ax=Ax) ## With moderator fit1 <- osmasem(model.name="With moderator", Mmatrix=M1, Tmatrix=T0, data=my.df) summary(fit1) plot(fit1, color="blue") coef(fit1) beep(2) ## R2 for RC Sres <- fit0$mx.fit$algebras$Smatrix$result Sres <- as.matrix(diag(Sres)) rownames(RAM0$A) dimnames(Sres) <- list(rownames(RAM0$A),'(residual) variance') Sres R2forRC <- 1-Sres[5] round(R2forRC,3) ## Model explained 52.7% variance in RC. ## Get the estimated coefficients coef(fit0) coef(fit1) ## R2 for moderation osmasemR2(fit1, fit0) ## Anova anova(fit1, fit0) ## Get the estimated A0 and A1 A0 <- mxEval(A0, fit1$mx.fit) A0 A1 <- mxEval(A1, fit1$mx.fit) A1 ## Compute the estimated A matrix at -1SD (-1) of the standardized grade A0 - A1 ## Compute the estimated A matrix at 0 (mean) of the standardized grade A0 ## Compute the estimated A matrix at +1SD (+1) of the standardized grade A0 + A1 ## When completed beep(3)