Intergenerational Effects of Parental ACEs Through Mental Health and Parenting: A MASEM Study
Amy Y Liu
2025-09-13
R code modelled after Jak et al. (2021) https://doi.org/10.1002/jrsm.1498 and Bilici et al. (2025) https://doi.org/10.1017/rsm.2024.10
Data Management
Load packages
#Data management
library(readr)
library(haven)
library(psych)
library(dplyr)
library(tidyr)
# Analyses
library(lavaan)
library(metafor)
library(metaSEM)
# Plots
if (!require("qgraph")) install.packages('qgraph')
library(qgraph)
if (!require("semPlot")) install.packages('semPlot')
library(semPlot)
if (!require("semptools")) install.packages('semptools')
library(semptools)
## load extra functions (readstack) to read in data
source("http://www.suzannejak.nl/MASEM_functions.R")Read data
data <- read_csv("study coding cleaned.csv")Adjust dependent effect sizes
To deal with publications reporting on dependent samples, we adjust by using the simple average method, as guided by a simulation study (Bilici et al., 2025) https://doi.org/10.1017/rsm.2024.10
data_temp <- as.data.frame(data[,c(1:12)]) #include only STUDY_ID,N,and effect sizes columns
data_l <- gather(data_temp, Cell, effectsize, MH_ACE:SE_NPT, factor_key = T) #long data form
data_l <- data_l[complete.cases(data_l$effectsize),] #remove NAs
#sublist per STUDY_ID
for(i in 1:nrow(data_l)){
cor <- data_l[i, 4][!is.na(data_l[i, 4])]
data_l$variance[i] <- ((1-(cor^2))^2)/(data_l[i,]$N - 1)
#V_r = (1-r^2)^2 / n-1
}
final_d <- spread(data_l, key = "Cell", value = "effectsize")
sub <- split(final_d, final_d$STUDY_ID)
study <- matrix(NA, nrow = length(unique(data$STUDY_ID)), ncol = 12)
for(i in 1:length(sub)){
for(j in 4:13){
study[i,1] <- mean(sub[[i]][,1], na.rm = TRUE) #STUDY_ID
study[i,2] <- mean(sub[[i]][,2], na.rm = TRUE) #N
study[i,(j-1)] <- mean(sub[[i]][,j], na.rm = TRUE) # simple average
}
}
colnames(study)<- colnames(data)[1:12] #give new df column names
head(study)## STUDY_ID N MH_ACE PPT_ACE NPT_ACE SE_ACE PPT_MH NPT_MH
## [1,] 1 1188 0.2000000 NaN NaN NaN NaN NaN
## [2,] 2 126 0.3700000 NaN NaN NaN NaN NaN
## [3,] 3 417 0.3459234 NaN NaN NaN NaN NaN
## [4,] 4 1185 0.2990000 0.09400000 NaN NaN -0.135 NaN
## [5,] 5 130 0.1547689 NaN NaN -0.02100000 NaN NaN
## [6,] 6 295 0.2500000 0.01047176 0.1415346 0.01178511 0.030 0.4472494
## SE_MH NPT_PPT SE_PPT SE_NPT
## [1,] NaN NaN NaN NaN
## [2,] NaN NaN NaN NaN
## [3,] NaN NaN NaN NaN
## [4,] NaN NaN NaN NaN
## [5,] 0.2810138 NaN NaN NaN
## [6,] NaN -0.06227524 -0.1483161 NaNData summary
nvar <- 5
varnames <- c("ACE","MH","PPT","NPT","SE")
labels <- list(varnames,varnames)
## Create a list of correlation matrices from columns 3 to 12
cormatrices <- readstack(study[,3:12], no.var = nvar,
var.names = varnames, diag = FALSE)
## Create vector with sample sizes
n <- study[,"N"]
## Calculate the sampling covariance matrix of the correlations
my.df <- Cor2DataFrame(cormatrices, n, acov = "weighted")
## View number of studies per each pair of correlation
pattern.na(cormatrices, show.na=FALSE)## ACE MH PPT NPT SE
## ACE 47 35 18 21 21
## MH 35 35 11 12 16
## PPT 18 11 18 8 7
## NPT 21 12 8 21 8
## SE 21 16 7 8 21## View sample size per each pair of correlation
pattern.n(cormatrices, n=n)## ACE MH PPT NPT SE
## ACE 23206.85 17792.855 9394.0 14492.355 13432.855
## MH 17792.85 17792.855 8165.0 9766.855 10183.855
## PPT 9394.00 8165.000 9394.0 6346.000 5337.500
## NPT 14492.35 9766.855 6346.0 14492.355 9791.355
## SE 13432.85 10183.855 5337.5 9791.355 13432.855Testing for publication bias
study <- as.data.frame(study)
es_names <- colnames(study)[3:12]
for (var in es_names) {
# calculate standard error for each pair of correlation
study <- study %>%
mutate(
!!paste0("SE_", var) := ifelse(!is.na(.data[[var]]) & !is.na(N),
(1 - (.data[[var]]^2)) / sqrt(N - 2),
NA))
# calculate variance for each pair of correlation
study <- study %>%
mutate(
!!paste0("VAR_", var) := (.data[[paste0("SE_", var)]])^2)
tmp_data <- study%>%
select(r = all_of(var), se = all_of(paste0("SE_", var)), STUDY_ID) %>%
filter(!is.na(r), !is.na(se))
# generate rma objects for each pair of correlation
res <- rma(yi=r, sei=se, data = tmp_data, slab = tmp_data$STUDY_ID)
# trim and fill analysis, one-sided
tf <- trimfill(res, side = "left")
print(paste("trim-and-fill imputed estimate for", var))
print(tf$beta)
# create funnel plots with trim and fill
funnel(tf, at=c(-0.5,-0.3,-0.1,0,0.1,0.3,0.5), ylim =c(0,.15), xlab= "Effect Size (r)", main = paste("Funnel plot with trim and fill for path", var))
# conduct egger's test. Significant p values indicate asymmetry of funnel plot
egger_test <- regtest(res, model="rma")
print(paste("Egger's test p-value for", var, ":", egger_test$pval))
# create funnel plots for each pair of correlation
forest.rma(res, shade = T, main = paste("Forest Plot for", var), xlab = "Effect Size (r)")
# reference: https://wviechtb.github.io/metafor/reference/forest.rma.html#:~:text=predstyle=%22bar%22%20:%20the,as%20part%20of%20the%20annotations.
}## [1] "trim-and-fill imputed estimate for MH_ACE"
## [,1]
## intrcpt 0.2418828
## [1] "Egger's test p-value for MH_ACE : 0.750296360541277"
## [1] "trim-and-fill imputed estimate for PPT_ACE"
## [,1]
## intrcpt -0.07127955
## [1] "Egger's test p-value for PPT_ACE : 0.0588784543252634"
## [1] "trim-and-fill imputed estimate for NPT_ACE"
## [,1]
## intrcpt 0.1574161
## [1] "Egger's test p-value for NPT_ACE : 0.662783487071442"
## [1] "trim-and-fill imputed estimate for SE_ACE"
## [,1]
## intrcpt 0.117139
## [1] "Egger's test p-value for SE_ACE : 0.0790133371849875"
## [1] "trim-and-fill imputed estimate for PPT_MH"
## [,1]
## intrcpt -0.1145923
## [1] "Egger's test p-value for PPT_MH : 0.750943541254805"
## [1] "trim-and-fill imputed estimate for NPT_MH"
## [,1]
## intrcpt 0.3366817
## [1] "Egger's test p-value for NPT_MH : 0.836374492430556"
## [1] "trim-and-fill imputed estimate for SE_MH"
## [,1]
## intrcpt 0.2924895
## [1] "Egger's test p-value for SE_MH : 0.246051091655606"
## [1] "trim-and-fill imputed estimate for NPT_PPT"
## [,1]
## intrcpt -0.1737879
## [1] "Egger's test p-value for NPT_PPT : 0.58791854055823"
## [1] "trim-and-fill imputed estimate for SE_PPT"
## [,1]
## intrcpt -0.218553
## [1] "Egger's test p-value for SE_PPT : 0.337839510208832"
## [1] "trim-and-fill imputed estimate for SE_NPT"
## [,1]
## intrcpt 0.2942871
## [1] "Egger's test p-value for SE_NPT : 0.528844010620502"
Overall analysis with osmasem()
Specifying models
#lavaan syntax
model <-
'MH ~ b21*ACE
PPT ~ b31*ACE
NPT ~ b41*ACE
SE ~ b51*ACE + b52*MH + b53*PPT + b54*NPT
MH ~~ PPT
MH ~~ NPT
PPT ~~ NPT
ACE ~~ 1 * ACE
MH ~~ MH
PPT ~~ PPT
NPT ~~ NPT
SE ~~ SE'
# Create visualization of the specified model
plot(model, col="yellow")
#### Fit the model
# Convert the lavvan syntax into the RAM specification
RAM1 <- lavaan2RAM(model, obs.variables=varnames)
# Create the model implied correlation structure with implicit diagonal constraints
M0 <- create.vechsR(A0=RAM1$A, S0=RAM1$S)
# Create the heterogeneity variance-covariance matrix
T0 <- create.Tau2(RAM=RAM1, RE.type="Diag", Transform="expLog", RE.startvalues=0.05)
# Estimate CI's for three indirect effects
ind1 <- mxAlgebra(b21*b52, name="ACE_MH_SE")
ind2 <- mxAlgebra(b31*b53, name="ACE_PPT_SE")
ind3 <- mxAlgebra(b41*b54, name="ACE_NPT_SE")
# fit the model using One Stage MASEM
mx.fit0 <- osmasem(model.name="No moderator", Mmatrix=M0, Tmatrix=T0,
data=my.df, intervals.type = "LB",
mxModel.Args = list(ind1,ind2,ind3,mxCI(c("ACE_MH_SE","ACE_PPT_SE","ACE_NPT_SE")))
)
## View the results
#summary(mx.fit0)$parameters [, c(1, 5,6,12)] # View only name, estimates, se, and p value
summary(mx.fit0)## Summary of No moderator
##
## free parameters:
## name matrix row col Estimate Std.Error A z value Pr(>|z|)
## 1 b21 A0 MH ACE 0.27615779 0.02004107 13.779591 0.000000e+00
## 2 b31 A0 PPT ACE -0.06502425 0.02866833 -2.268156 2.331969e-02
## 3 b41 A0 NPT ACE 0.15001714 0.02833838 5.293780 1.198138e-07
## 4 b51 A0 SE ACE 0.07311381 0.03442705 2.123732 3.369260e-02
## 5 b52 A0 SE MH 0.18048413 0.03444196 5.240239 1.603685e-07
## 6 b53 A0 SE PPT -0.14137552 0.04886285 -2.893313 3.812015e-03
## 7 b54 A0 SE NPT 0.18933048 0.04312169 4.390609 1.130337e-05
## 8 MHWITHPPT S0 PPT MH -0.08832291 0.04675940 -1.888880 5.890785e-02
## 9 MHWITHNPT S0 NPT MH 0.27832137 0.04016769 6.928987 4.238609e-12
## 10 PPTWITHNPT S0 NPT PPT -0.15934592 0.11292387 -1.411092 1.582176e-01
## 11 Tau1_1 vecTau1 1 1 -4.60404402 0.32078587 -14.352391 0.000000e+00
## 12 Tau1_2 vecTau1 2 1 -4.66638134 0.49320787 -9.461287 0.000000e+00
## 13 Tau1_3 vecTau1 3 1 -4.36514407 0.41750283 -10.455364 0.000000e+00
## 14 Tau1_4 vecTau1 4 1 -4.20764369 0.37243847 -11.297554 0.000000e+00
## 15 Tau1_5 vecTau1 5 1 -3.96357434 0.66296658 -5.978543 2.251419e-09
## 16 Tau1_6 vecTau1 6 1 -4.18878016 0.52564096 -7.968900 1.554312e-15
## 17 Tau1_7 vecTau1 7 1 -5.07666454 0.56806939 -8.936698 0.000000e+00
## 18 Tau1_8 vecTau1 8 1 -2.32989357 0.52724410 -4.419004 9.915699e-06
## 19 Tau1_9 vecTau1 9 1 -4.96285852 0.81565983 -6.084471 1.168766e-09
## 20 Tau1_10 vecTau1 10 1 -5.04943637 0.61837333 -8.165676 2.220446e-16
##
## confidence intervals:
## lbound estimate ubound note
## No moderator.Amatrix[1,1] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Amatrix[2,1] 0.236209515 0.276157790 0.317243011
## No moderator.Amatrix[3,1] -0.127349358 -0.065024253 -0.008595234
## No moderator.Amatrix[4,1] 0.090830167 0.150017140 0.208365803
## No moderator.Amatrix[5,1] 0.003967792 0.073113808 0.144400930
## No moderator.Amatrix[1,2] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Amatrix[2,2] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Amatrix[3,2] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Amatrix[4,2] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Amatrix[5,2] 0.110417122 0.180484129 0.252702781
## No moderator.Amatrix[1,3] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Amatrix[2,3] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Amatrix[3,3] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Amatrix[4,3] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Amatrix[5,3] -0.250058598 -0.141375516 -0.037943578
## No moderator.Amatrix[1,4] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Amatrix[2,4] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Amatrix[3,4] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Amatrix[4,4] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Amatrix[5,4] 0.096492035 0.189330478 0.282223619
## No moderator.Amatrix[1,5] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Amatrix[2,5] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Amatrix[3,5] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Amatrix[4,5] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Amatrix[5,5] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Smatrix[1,1] 1.000000000 1.000000000 1.000000000 !!!
## No moderator.Smatrix[2,1] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Smatrix[3,1] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Smatrix[4,1] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Smatrix[5,1] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Smatrix[1,2] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Smatrix[2,2] 0.899426203 0.923736875 0.944070766
## No moderator.Smatrix[3,2] -0.191948155 -0.088322915 0.017889302
## No moderator.Smatrix[4,2] 0.194539236 0.278321365 0.366655906
## No moderator.Smatrix[5,2] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Smatrix[1,3] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Smatrix[2,3] -0.191948155 -0.088322915 0.017889302
## No moderator.Smatrix[3,3] 0.983759016 0.995771847 0.999908580
## No moderator.Smatrix[4,3] -0.411658499 -0.159345919 0.087499611
## No moderator.Smatrix[5,3] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Smatrix[1,4] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Smatrix[2,4] 0.194539236 0.278321365 0.366655906
## No moderator.Smatrix[3,4] -0.411658499 -0.159345919 0.087499611
## No moderator.Smatrix[4,4] 0.956608769 0.977494858 0.991256330
## No moderator.Smatrix[5,4] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Smatrix[1,5] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Smatrix[2,5] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Smatrix[3,5] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Smatrix[4,5] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Smatrix[5,5] 0.804755744 0.857132323 0.895313897
## No moderator.Tau2[1,1] 0.005356778 0.010011268 0.019066649
## No moderator.Tau2[2,1] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[3,1] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[4,1] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[5,1] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[6,1] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[7,1] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[8,1] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[9,1] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[10,1] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[1,2] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[2,2] 0.003447358 0.009406246 0.024998746
## No moderator.Tau2[3,2] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[4,2] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[5,2] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[6,2] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[7,2] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[8,2] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[9,2] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[10,2] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[1,3] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[2,3] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[3,3] 0.005678176 0.012712823 0.029715177
## No moderator.Tau2[4,3] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[5,3] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[6,3] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[7,3] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[8,3] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[9,3] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[10,3] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[1,4] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[2,4] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[3,4] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[4,4] 0.007398575 0.014881392 0.032293122
## No moderator.Tau2[5,4] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[6,4] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[7,4] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[8,4] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[9,4] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[10,4] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[1,5] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[2,5] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[3,5] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[4,5] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[5,5] 0.004959076 0.018995098 0.062259303
## No moderator.Tau2[6,5] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[7,5] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[8,5] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[9,5] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[10,5] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[1,6] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[2,6] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[3,6] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[4,6] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[5,6] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[6,6] 0.005448532 0.015164772 0.044905424
## No moderator.Tau2[7,6] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[8,6] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[9,6] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[10,6] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[1,7] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[2,7] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[3,7] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[4,7] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[5,7] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[6,7] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[7,7] 0.001930155 0.006240690 0.018620999
## No moderator.Tau2[8,7] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[9,7] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[10,7] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[1,8] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[2,8] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[3,8] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[4,8] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[5,8] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[6,8] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[7,8] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[8,8] 0.038777600 0.097306103 0.313370688
## No moderator.Tau2[9,8] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[10,8] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[1,9] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[2,9] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[3,9] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[4,9] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[5,9] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[6,9] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[7,9] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[8,9] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[9,9] 0.001360245 0.006992910 0.008435970
## No moderator.Tau2[10,9] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[1,10] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[2,10] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[3,10] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[4,10] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[5,10] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[6,10] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[7,10] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[8,10] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[9,10] 0.000000000 0.000000000 0.000000000 !!!
## No moderator.Tau2[10,10] 0.002122831 0.006412947 0.021204650
## No moderator.ACE_MH_SE[1,1] 0.030319162 0.049842098 0.071922634
## No moderator.ACE_PPT_SE[1,1] 0.001230016 0.009192837 0.022769664
## No moderator.ACE_NPT_SE[1,1] 0.014051527 0.028402817 0.047472139
##
## Model Statistics:
## | Parameters | Degrees of Freedom | Fit (-2lnL units)
## Model: 20 137 -180.8696
## Saturated: 65 92 NA
## Independence: 20 137 NA
## Number of observations/statistics: 23206.85/157
##
## Information Criteria:
## | df Penalty | Parameters Penalty | Sample-Size Adjusted
## AIC: -454.8696 -140.86959 -140.83336
## BIC: -1558.0214 20.17447 -43.38488
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2025-09-19 22:59:07
## Wall clock time: 128.2932 secs
## optimizer: SLSQP
## OpenMx version number: 2.21.13
## Need help? See help(mxSummary)Moderator analysis (Part 1)
Add moderator data to the dataset
## Get Study-level mean and standardize continuous moderators
data <- as.data.frame(data)
sub <- split(data, data$STUDY_ID)
for (i in seq_along(sub)) {
study[i,"QA"] <- mean(sub[[i]]$QA, na.rm = TRUE)
study[i,"ACE_TYPE"] <- mean(sub[[i]]$ACE_TYPE, na.rm = TRUE)
study[i,"PERCENT_MOM"] <- mean(sub[[i]]$PERCENT_MOM, na.rm = TRUE)
study[i,"DESIGN_LONG"] <- mean(sub[[i]]$DESIGN_LONG, na.rm = TRUE)
study[i,"HEALTHCARE"] <- mean(sub[[i]]$HEALTHCARE, na.rm = TRUE)
}
mean(study$QA); sd(study$QA)## [1] 4.003546## [1] 1.039375mean(study$ACE_TYPE); sd(study$ACE_TYPE)## [1] 11.44681## [1] 4.533971mean(study$PERCENT_MOM); sd(study$PERCENT_MOM)## [1] 0.9215186## [1] 0.2166399mean(study$DESIGN_LONG)## [1] 0.3829787mean(study$HEALTHCARE)## [1] 0.3829787# Scale continuous moderators
study$QA <- scale(study$QA)
study$ACE_TYPE <- scale(study$ACE_TYPE)
study$PERCENT_MOM <- scale(study$PERCENT_MOM)Continous moderator
Study quality
my.df$data <- data.frame(my.df$data, QA = study$QA,
check.names=FALSE)
Ax <- matrix (c(0,0,0,0,0,
"0*data.QA",0,0,0,0,
"0*data.QA",0,0,0,0,
"0*data.QA",0,0,0,0,
"0*data.QA","0*data.QA","0*data.QA","0*data.QA",0),
nrow=5,ncol=5,byrow = TRUE); Ax #Add moderation to all direct effects## [,1] [,2] [,3] [,4] [,5]
## [1,] "0" "0" "0" "0" "0"
## [2,] "0*data.QA" "0" "0" "0" "0"
## [3,] "0*data.QA" "0" "0" "0" "0"
## [4,] "0*data.QA" "0" "0" "0" "0"
## [5,] "0*data.QA" "0*data.QA" "0*data.QA" "0*data.QA" "0"M1 <- create.vechsR(A0=RAM1$A, S0=RAM1$S, Ax=Ax)
mx.fit1 <- osmasem(model.name = 'Study quality (z) as Moderator', Mmatrix = M1, Tmatrix = T0, data = my.df)
anova(mx.fit1, mx.fit0)## base comparison ep minus2LL df AIC
## 1 Study quality (z) as Moderator <NA> 27 -185.3351 130 -131.3351
## 2 Study quality (z) as Moderator No moderator 20 -180.8696 137 -140.8696
## diffLL diffdf p
## 1 NA NA NA
## 2 4.465558 7 0.7248606summary(mx.fit1)$parameters [, c(1, 5,6,12)] # View only name, estimates, se, and p value## name Estimate Std.Error Pr(>|z|)
## 1 b21 0.280631291 0.02039719 0.000000e+00
## 2 b31 -0.065822904 0.02763005 1.720516e-02
## 3 b41 0.149734675 0.02872826 1.867105e-07
## 4 b51 0.074072511 0.03339517 2.655071e-02
## 5 b52 0.174982250 0.03567298 9.334376e-07
## 6 b53 -0.135294548 0.05072082 7.643239e-03
## 7 b54 0.221347591 0.05511273 5.912654e-05
## 8 MHWITHPPT -0.090298642 0.04719617 5.571449e-02
## 9 MHWITHNPT 0.279274141 0.04006264 3.148370e-12
## 10 PPTWITHNPT -0.160329592 0.11311287 1.563572e-01
## 11 b21_1 -0.016971868 0.02016576 4.000018e-01
## 12 b31_1 0.018507934 0.02682068 4.901552e-01
## 13 b41_1 0.002354331 0.02713204 9.308519e-01
## 14 b51_1 -0.033960472 0.03910571 3.851604e-01
## 15 b52_1 -0.008088900 0.04003265 8.398717e-01
## 16 b53_1 -0.010512712 0.05748389 8.548914e-01
## 17 b54_1 -0.037135159 0.05445129 4.952462e-01
## 18 Tau1_1 -4.645562197 0.32691859 0.000000e+00
## 19 Tau1_2 -4.777706383 0.53590274 0.000000e+00
## 20 Tau1_3 -4.363272933 0.41876954 0.000000e+00
## 21 Tau1_4 -4.330879499 0.38434035 0.000000e+00
## 22 Tau1_5 -3.940673377 0.65633196 1.924129e-09
## 23 Tau1_6 -4.196779992 0.52712450 1.776357e-15
## 24 Tau1_7 -5.209255981 0.60280867 0.000000e+00
## 25 Tau1_8 -2.327546948 0.52719890 1.010420e-05
## 26 Tau1_9 -4.969363585 0.85243579 5.555886e-09
## 27 Tau1_10 -5.156574458 0.62281662 2.220446e-16ACEs type in the questionnaire
## Add ace to the dataset
my.df$data <- data.frame(my.df$data, ACE_TYPE = study$ACE_TYPE,
check.names=FALSE)
Ax <- matrix (c(0,0,0,0,0,
"0*data.ACE_TYPE",0,0,0,0,
"0*data.ACE_TYPE",0,0,0,0,
"0*data.ACE_TYPE",0,0,0,0,
"0*data.ACE_TYPE",0,0,0,0),
nrow=5,ncol=5,byrow = TRUE); Ax #Add moderation to all direct effect coming out from parental ACEs## [,1] [,2] [,3] [,4] [,5]
## [1,] "0" "0" "0" "0" "0"
## [2,] "0*data.ACE_TYPE" "0" "0" "0" "0"
## [3,] "0*data.ACE_TYPE" "0" "0" "0" "0"
## [4,] "0*data.ACE_TYPE" "0" "0" "0" "0"
## [5,] "0*data.ACE_TYPE" "0" "0" "0" "0"M1 <- create.vechsR(A0=RAM1$A, S0=RAM1$S, Ax=Ax)
mx.fit2 <- osmasem(model.name = 'ACEs type (z) as Moderator', Mmatrix = M1, Tmatrix = T0, data = my.df)
anova(mx.fit2,mx.fit0)## base comparison ep minus2LL df AIC diffLL
## 1 ACEs type (z) as Moderator <NA> 24 -182.8295 133 -134.8295 NA
## 2 ACEs type (z) as Moderator No moderator 20 -180.8696 137 -140.8696 1.959907
## diffdf p
## 1 NA NA
## 2 4 0.7431331summary(mx.fit2)$parameters [, c(1, 5,6,12)] # View only name, estimates, se, and p value## name Estimate Std.Error Pr(>|z|)
## 1 b21 0.275718297 0.02004524 0.000000e+00
## 2 b31 -0.065191795 0.02853882 2.235275e-02
## 3 b41 0.146031292 0.02855197 3.144427e-07
## 4 b51 0.074978558 0.03295969 2.291455e-02
## 5 b52 0.179264487 0.03435890 1.814418e-07
## 6 b53 -0.141289691 0.04897038 3.911564e-03
## 7 b54 0.191349296 0.04315320 9.242097e-06
## 8 MHWITHPPT -0.088647780 0.04686992 5.857664e-02
## 9 MHWITHNPT 0.278057407 0.04017993 4.506617e-12
## 10 PPTWITHNPT -0.158820715 0.11291201 1.595495e-01
## 11 b21_1 -0.001533683 0.02032769 9.398583e-01
## 12 b31_1 0.010571866 0.04367846 8.087504e-01
## 13 b41_1 -0.024189729 0.03647812 5.072474e-01
## 14 b51_1 0.040602276 0.02950184 1.687404e-01
## 15 Tau1_1 -4.601681798 0.32099448 0.000000e+00
## 16 Tau1_2 -4.678500633 0.49674749 0.000000e+00
## 17 Tau1_3 -4.400963041 0.42475861 0.000000e+00
## 18 Tau1_4 -4.330222678 0.38492457 0.000000e+00
## 19 Tau1_5 -3.959078330 0.66204520 2.230382e-09
## 20 Tau1_6 -4.187898467 0.52452962 1.332268e-15
## 21 Tau1_7 -5.072682568 0.57504953 0.000000e+00
## 22 Tau1_8 -2.330575578 0.52727435 9.868102e-06
## 23 Tau1_9 -4.965509461 0.81839950 1.300446e-09
## 24 Tau1_10 -5.054554152 0.61770133 2.220446e-16Percent Mother
my.df$data <- data.frame(my.df$data, PERCENT_MOM = study$PERCENT_MOM,
check.names=FALSE)
Ax <- matrix (c(0,0,0,0,0,
"0*data.PERCENT_MOM",0,0,0,0,
"0*data.PERCENT_MOM",0,0,0,0,
"0*data.PERCENT_MOM",0,0,0,0,
"0*data.PERCENT_MOM","0*data.PERCENT_MOM","0*data.PERCENT_MOM","0*data.PERCENT_MOM",0),
nrow=5,ncol=5,byrow = TRUE); Ax ## [,1] [,2] [,3]
## [1,] "0" "0" "0"
## [2,] "0*data.PERCENT_MOM" "0" "0"
## [3,] "0*data.PERCENT_MOM" "0" "0"
## [4,] "0*data.PERCENT_MOM" "0" "0"
## [5,] "0*data.PERCENT_MOM" "0*data.PERCENT_MOM" "0*data.PERCENT_MOM"
## [,4] [,5]
## [1,] "0" "0"
## [2,] "0" "0"
## [3,] "0" "0"
## [4,] "0" "0"
## [5,] "0*data.PERCENT_MOM" "0"M1 <- create.vechsR(A0=RAM1$A, S0=RAM1$S, Ax=Ax)
mx.fit3 <- osmasem(model.name = 'Percentage of Mother (z) as Moderator', Mmatrix = M1, Tmatrix = T0, data = my.df)
anova(mx.fit3,mx.fit0)## base comparison ep minus2LL df AIC
## 1 Percentage of Mother (z) as Moderator <NA> 27 -194.4541 130 -140.4541
## 2 Percentage of Mother (z) as Moderator No moderator 20 -180.8696 137 -140.8696
## diffLL diffdf p
## 1 NA NA NA
## 2 13.58455 7 0.05908351summary(mx.fit3)$parameters [, c(1, 5,6,12)] # View only name, estimates, se, and p value## name Estimate Std.Error Pr(>|z|)
## 1 b21 0.275549255 1.976726e-02 0.000000e+00
## 2 b31 -0.066084912 2.843008e-02 2.010022e-02
## 3 b41 0.145986938 2.912339e-02 5.367051e-07
## 4 b51 0.075477031 3.317536e-02 2.290035e-02
## 5 b52 0.174297642 3.642912e-02 1.713548e-06
## 6 b53 -0.264946169 5.028639e-02 1.373593e-07
## 7 b54 0.169442233 5.140485e-02 9.799158e-04
## 8 MHWITHPPT -0.090609064 4.680664e-02 5.289019e-02
## 9 MHWITHNPT 0.274923301 3.970169e-02 4.368284e-12
## 10 PPTWITHNPT -0.161969689 1.134193e-01 1.532742e-01
## 11 b21_1 -0.026586194 1.822989e-02 1.447344e-01
## 12 b31_1 -0.009714186 1.961387e-02 6.204086e-01
## 13 b41_1 -0.018753029 2.026742e-02 3.548204e-01
## 14 b51_1 -0.070224420 5.777281e-02 2.241651e-01
## 15 b52_1 0.012263617 5.693705e-02 8.294640e-01
## 16 b53_1 0.589824073 1.404526e-01 2.675624e-05
## 17 b54_1 0.089652547 9.126430e-02 3.259325e-01
## 18 Tau1_1 -4.642751768 3.228438e-01 0.000000e+00
## 19 Tau1_2 -4.713838524 4.912853e-01 0.000000e+00
## 20 Tau1_3 -4.349293151 4.135312e-01 0.000000e+00
## 21 Tau1_4 -4.514030147 3.838801e-01 0.000000e+00
## 22 Tau1_5 -3.964008939 6.573284e-01 1.634684e-09
## 23 Tau1_6 -4.238554661 5.382469e-01 3.330669e-15
## 24 Tau1_7 -5.241136393 5.669635e-01 0.000000e+00
## 25 Tau1_8 -2.331968928 5.271435e-01 9.698702e-06
## 26 Tau1_9 -33.019456917 1.199390e+04 9.978034e-01
## 27 Tau1_10 -5.085413499 6.106056e-01 0.000000e+00Categorical moderators
Whether recruited from primary healthcare settings e.g., prental care or primary pediatric
my.df$data <- data.frame(my.df$data, HEALTHCARE=study$HEALTHCARE,
check.names=FALSE)
Ax <- matrix(c(0,0,0,0,0,
"0*data.HEALTHCARE",0,0,0,0,
"0*data.HEALTHCARE",0,0,0,0,
"0*data.HEALTHCARE",0,0,0,0,
"0*data.HEALTHCARE","0*data.HEALTHCARE","0*data.HEALTHCARE","0*data.HEALTHCARE",0),
nrow=nvar, ncol=nvar, byrow=TRUE);Ax## [,1] [,2] [,3]
## [1,] "0" "0" "0"
## [2,] "0*data.HEALTHCARE" "0" "0"
## [3,] "0*data.HEALTHCARE" "0" "0"
## [4,] "0*data.HEALTHCARE" "0" "0"
## [5,] "0*data.HEALTHCARE" "0*data.HEALTHCARE" "0*data.HEALTHCARE"
## [,4] [,5]
## [1,] "0" "0"
## [2,] "0" "0"
## [3,] "0" "0"
## [4,] "0" "0"
## [5,] "0*data.HEALTHCARE" "0"M1 <- create.vechsR(A0=RAM1$A, S0=RAM1$S, Ax=Ax)
mx.fit4 <- osmasem(model.name="Recruitment site as moderator", Mmatrix=M1, Tmatrix=T0, data=my.df)
anova(mx.fit4,mx.fit0)## base comparison ep minus2LL df AIC diffLL
## 1 Recruitment site as moderator <NA> 27 -193.3918 130 -139.3918 NA
## 2 Recruitment site as moderator No moderator 20 -180.8696 137 -140.8696 12.5222
## diffdf p
## 1 NA NA
## 2 7 0.08464161summary(mx.fit4)$parameters [, c(1, 5,6,12)] # View only name, estimates, se, and p value## name Estimate Std.Error Pr(>|z|)
## 1 b21 0.28074870 0.02692133 0.000000e+00
## 2 b31 -0.08381021 0.03245417 9.811093e-03
## 3 b41 0.14058665 0.03311446 2.181293e-05
## 4 b51 0.09473954 0.04013904 1.826096e-02
## 5 b52 0.20075630 0.03471898 7.367667e-09
## 6 b53 -0.08824532 0.05916461 1.358245e-01
## 7 b54 0.21656299 0.03507198 6.624108e-10
## 8 MHWITHPPT -0.09019909 0.04557657 4.780845e-02
## 9 MHWITHNPT 0.27953079 0.03935331 1.219913e-12
## 10 PPTWITHNPT -0.15964381 0.11280107 1.569898e-01
## 11 b21_1 -0.01038499 0.04038505 7.970635e-01
## 12 b31_1 0.07183785 0.06019109 2.326750e-01
## 13 b41_1 0.03611623 0.06230219 5.621208e-01
## 14 b51_1 -0.05717589 0.07213193 4.279776e-01
## 15 b52_1 -0.04075425 0.05851850 4.861573e-01
## 16 b53_1 -0.13393478 0.08434248 1.122892e-01
## 17 b54_1 -0.15768787 0.07955914 4.747690e-02
## 18 Tau1_1 -4.59231153 0.32000387 0.000000e+00
## 19 Tau1_2 -4.76882178 0.50283049 0.000000e+00
## 20 Tau1_3 -4.37950946 0.41586748 0.000000e+00
## 21 Tau1_4 -4.27207020 0.36549828 0.000000e+00
## 22 Tau1_5 -4.02698525 0.67670148 2.666661e-09
## 23 Tau1_6 -4.24161731 0.53309498 1.776357e-15
## 24 Tau1_7 -5.62075670 0.75885976 1.292300e-13
## 25 Tau1_8 -2.33210676 0.52729940 9.745820e-06
## 26 Tau1_9 -5.62437611 1.32204320 2.096985e-05
## 27 Tau1_10 -5.82794752 0.72484323 8.881784e-16Study design. longitudinal = 1, cross-sectional/intervention study = 0
## Add moderator to the dataset
my.df$data <- data.frame(my.df$data, DESIGN_LONG=study$DESIGN_LONG,
check.names=FALSE)
Ax <- matrix(c(0,0,0,0,0,
"0*data.DESIGN_LONG",0,0,0,0,
"0*data.DESIGN_LONG",0,0,0,0,
"0*data.DESIGN_LONG",0,0,0,0,
"0*data.DESIGN_LONG","0*data.DESIGN_LONG","0*data.DESIGN_LONG","0*data.DESIGN_LONG",0),
nrow=nvar, ncol=nvar, byrow=TRUE)
Ax## [,1] [,2] [,3]
## [1,] "0" "0" "0"
## [2,] "0*data.DESIGN_LONG" "0" "0"
## [3,] "0*data.DESIGN_LONG" "0" "0"
## [4,] "0*data.DESIGN_LONG" "0" "0"
## [5,] "0*data.DESIGN_LONG" "0*data.DESIGN_LONG" "0*data.DESIGN_LONG"
## [,4] [,5]
## [1,] "0" "0"
## [2,] "0" "0"
## [3,] "0" "0"
## [4,] "0" "0"
## [5,] "0*data.DESIGN_LONG" "0"M1 <- create.vechsR(A0=RAM1$A, S0=RAM1$S, Ax=Ax)
mx.fit5 <- osmasem(model.name="Study design as moderator", Mmatrix=M1, Tmatrix=T0, data=my.df)
anova(mx.fit5,mx.fit0)## base comparison ep minus2LL df AIC diffLL
## 1 Study design as moderator <NA> 27 -194.4317 130 -140.4317 NA
## 2 Study design as moderator No moderator 20 -180.8696 137 -140.8696 13.56216
## diffdf p
## 1 NA NA
## 2 7 0.05953979summary(mx.fit5)$parameters [, c(1, 5,6,12)]## name Estimate Std.Error Pr(>|z|)
## 1 b21 0.3041404183 0.02696459 0.000000e+00
## 2 b31 -0.0860172314 0.03178938 6.813028e-03
## 3 b41 0.1505674995 0.03433784 1.160470e-05
## 4 b51 0.1576217624 0.04276348 2.278987e-04
## 5 b52 0.1650960712 0.04763602 5.286911e-04
## 6 b53 -0.1401594108 0.05961304 1.871530e-02
## 7 b54 0.2033885355 0.06080460 8.229644e-04
## 8 MHWITHPPT -0.0895524362 0.04539313 4.851618e-02
## 9 MHWITHNPT 0.2776269081 0.04026926 5.414558e-12
## 10 PPTWITHNPT -0.1586065888 0.11281476 1.597530e-01
## 11 b21_1 -0.0601823577 0.03926151 1.253106e-01
## 12 b31_1 0.0737411661 0.05711778 1.966908e-01
## 13 b41_1 0.0004637617 0.06121940 9.939558e-01
## 14 b51_1 -0.1531255359 0.05671785 6.938606e-03
## 15 b52_1 0.0214288891 0.06238440 7.312242e-01
## 16 b53_1 0.0117429833 0.09105156 8.973808e-01
## 17 b54_1 -0.0195278688 0.08035337 8.079861e-01
## 18 Tau1_1 -4.6566125145 0.32097885 0.000000e+00
## 19 Tau1_2 -4.8701257050 0.55992948 0.000000e+00
## 20 Tau1_3 -4.3507333032 0.42461197 0.000000e+00
## 21 Tau1_4 -4.7502771946 0.39834299 0.000000e+00
## 22 Tau1_5 -4.0320095600 0.67799391 2.731511e-09
## 23 Tau1_6 -4.1876876656 0.52559339 1.554312e-15
## 24 Tau1_7 -5.1527776115 0.59063555 0.000000e+00
## 25 Tau1_8 -2.3320479797 0.52722328 9.722064e-06
## 26 Tau1_9 -4.9826515867 0.80208490 5.227285e-10
## 27 Tau1_10 -5.0554568217 0.61169447 2.220446e-16Moderator analysis (Part 2)
To test the moderation of parent age, child age, we have to remove studies with missing values in these moderators
Child age
data <- read_csv("study coding cleaned_child.csv")## Rows: 51 Columns: 19
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (19): STUDY_ID, N, MH_ACE, PPT_ACE, NPT_ACE, SE_ACE, PPT_MH, NPT_MH, SE_...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.data_temp <- as.data.frame(data[,c(1:12)])
data_l <- gather(data_temp, Cell, effectsize, MH_ACE:SE_NPT, factor_key = T)
data_l <- data_l[complete.cases(data_l$effectsize),]
for(i in 1:nrow(data_l)){
cor <- data_l[i, 4][!is.na(data_l[i, 4])]
data_l$variance[i] <- ((1-(cor^2))^2)/(data_l[i,]$N - 1)
}
final_d <- spread(data_l, key = "Cell", value = "effectsize")
sub <- split(final_d, final_d$STUDY_ID)
study <- matrix(NA, nrow = length(unique(data$STUDY_ID)), ncol = 12)
for(i in 1:length(sub)){
for(j in 4:13){
study[i,1] <- mean(sub[[i]][,1], na.rm = TRUE)
study[i,2] <- mean(sub[[i]][,2], na.rm = TRUE)
study[i,(j-1)] <- mean(sub[[i]][,j], na.rm = TRUE)
}
}
colnames(study)<- colnames(data)[1:12]
head(study)## STUDY_ID N MH_ACE PPT_ACE NPT_ACE SE_ACE PPT_MH
## [1,] 1 1188.000 0.2000000 NaN NaN NaN NaN
## [2,] 3 417.000 0.3459234 NaN NaN NaN NaN
## [3,] 5 130.000 0.1547689 NaN NaN -0.02100000 NaN
## [4,] 6 295.000 0.2500000 0.01047176 0.1415346 0.04088650 0.03
## [5,] 7 1895.455 0.1765690 NaN 0.1000000 0.07246382 NaN
## [6,] 8 44.000 0.2350000 NaN 0.0940000 NaN NaN
## NPT_MH SE_MH NPT_PPT SE_PPT SE_NPT
## [1,] NaN NaN NaN NaN NaN
## [2,] NaN NaN NaN NaN NaN
## [3,] NaN 0.2810138 NaN NaN NaN
## [4,] 0.4472494 NaN -0.06227524 -0.06289462 NaN
## [5,] 0.3000000 0.1578109 NaN NaN 0.1518028
## [6,] 0.2470000 NaN NaN NaN NaNnvar <- 5
varnames <- c("ACE","MH","PPT","NPT","SE")
labels <- list(varnames,varnames)
cormatrices <- readstack(study[,3:12], no.var = nvar,
var.names = varnames, diag = FALSE)
n <- study[,"N"]
my.df <- Cor2DataFrame(cormatrices, n, acov = "weighted")
pattern.na(cormatrices, show.na=FALSE)## ACE MH PPT NPT SE
## ACE 35 25 14 17 21
## MH 25 25 9 10 16
## PPT 14 9 14 6 7
## NPT 17 10 6 17 8
## SE 21 16 7 8 21pattern.n(cormatrices, n=n)## ACE MH PPT NPT SE
## ACE 18886.021 13931.021 7637.167 13170.355 13592.021
## MH 13931.021 13931.021 6736.167 8903.855 10343.021
## PPT 7637.167 6736.167 7637.167 5834.000 5337.500
## NPT 13170.355 8903.855 5834.000 13170.355 9791.355
## SE 13592.021 10343.021 5337.500 9791.355 13592.021temp <- aggregate(CHILD_AGE ~ STUDY_ID, data = data, mean, na.rm = T)
mean(temp$CHILD_AGE); sd(temp$CHILD_AGE)## [1] 4.143183## [1] 4.364568temp$CHILD_AGE <- scale(temp$CHILD_AGE)
study <- merge(study, temp, by = "STUDY_ID")
my.df$data <- data.frame(my.df$data, CHILD_AGE = study$CHILD_AGE,
check.names=FALSE)
Ax <- matrix (c(0,0,0,0,0,
"0*data.CHILD_AGE",0,0,0,0,
"0*data.CHILD_AGE",0,0,0,0,
"0*data.CHILD_AGE",0,0,0,0,
"0*data.CHILD_AGE","0*data.CHILD_AGE","0*data.CHILD_AGE","0*data.CHILD_AGE",0),
nrow=5,ncol=5,byrow = TRUE); Ax## [,1] [,2] [,3]
## [1,] "0" "0" "0"
## [2,] "0*data.CHILD_AGE" "0" "0"
## [3,] "0*data.CHILD_AGE" "0" "0"
## [4,] "0*data.CHILD_AGE" "0" "0"
## [5,] "0*data.CHILD_AGE" "0*data.CHILD_AGE" "0*data.CHILD_AGE"
## [,4] [,5]
## [1,] "0" "0"
## [2,] "0" "0"
## [3,] "0" "0"
## [4,] "0" "0"
## [5,] "0*data.CHILD_AGE" "0"M1 <- create.vechsR(A0=RAM1$A, S0=RAM1$S, Ax=Ax)
mx.fit0 <- osmasem(model.name="No moderator", Mmatrix=M0, Tmatrix=T0,
data=my.df, intervals.type = "LB"
,
mxModel.Args = list(ind1, ind2,ind3, mxCI(c("ACE_MH_SE","ACE_PPT_SE","ACE_NPT_SE")))
)
mx.fit7 <- osmasem(model.name = 'Child Age (z) as Moderator', Mmatrix = M1, Tmatrix = T0, data = my.df)
anova(mx.fit7,mx.fit0)## base comparison ep minus2LL df AIC diffLL
## 1 Child Age (z) as Moderator <NA> 27 -166.9693 106 -112.9693 NA
## 2 Child Age (z) as Moderator No moderator 20 -158.8828 113 -118.8828 8.086516
## diffdf p
## 1 NA NA
## 2 7 0.3250269summary(mx.fit7)$parameters [, c(1,5,6,12)] # View only name, estimates, se, and p value## name Estimate Std.Error Pr(>|z|)
## 1 b21 0.283766535 0.02071866 0.000000e+00
## 2 b31 -0.091771400 0.03189179 4.007288e-03
## 3 b41 0.143281742 0.02818942 3.718722e-07
## 4 b51 0.058458861 0.03406493 8.614355e-02
## 5 b52 0.196686956 0.03379384 5.877661e-09
## 6 b53 -0.091615627 0.05461484 9.344747e-02
## 7 b54 0.207460734 0.04408464 2.526685e-06
## 8 MHWITHPPT -0.075323516 0.05866159 1.991298e-01
## 9 MHWITHNPT 0.296339431 0.04635229 1.624501e-10
## 10 PPTWITHNPT -0.162731374 0.14138572 2.497427e-01
## 11 b21_1 0.013466271 0.01926768 4.846116e-01
## 12 b31_1 -0.016646690 0.02789871 5.507189e-01
## 13 b41_1 0.001259066 0.02701032 9.628207e-01
## 14 b51_1 0.045610422 0.03421077 1.824602e-01
## 15 b52_1 -0.049283339 0.03235753 1.277369e-01
## 16 b53_1 -0.064881221 0.03672807 7.730701e-02
## 17 b54_1 -0.043241524 0.04694761 3.570195e-01
## 18 Tau1_1 -5.041348742 0.42851232 0.000000e+00
## 19 Tau1_2 -4.793851945 0.58159620 2.220446e-16
## 20 Tau1_3 -4.691830858 0.51774722 0.000000e+00
## 21 Tau1_4 -4.326907425 0.37517542 0.000000e+00
## 22 Tau1_5 -3.692748936 0.69607129 1.125925e-07
## 23 Tau1_6 -4.077599320 0.56743627 6.672440e-13
## 24 Tau1_7 -5.407765368 0.65333177 2.220446e-16
## 25 Tau1_8 -2.155592862 0.60404876 3.589291e-04
## 26 Tau1_9 -5.952118056 1.67751449 3.879077e-04
## 27 Tau1_10 -5.228974036 0.62966250 0.000000e+00Parent age
data <- read_csv("study coding cleaned_parent.csv")## Rows: 58 Columns: 19
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (19): STUDY_ID, N, MH_ACE, PPT_ACE, NPT_ACE, SE_ACE, PPT_MH, NPT_MH, SE_...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.data_temp <- as.data.frame(data[,c(1:12)])
data_l <- gather(data_temp, Cell, effectsize, MH_ACE:SE_NPT, factor_key = T)
data_l <- data_l[complete.cases(data_l$effectsize),]
for(i in 1:nrow(data_l)){
cor <- data_l[i, 4][!is.na(data_l[i, 4])]
data_l$variance[i] <- ((1-(cor^2))^2)/(data_l[i,]$N - 1)
}
final_d <- spread(data_l, key = "Cell", value = "effectsize")
sub <- split(final_d, final_d$STUDY_ID)
study <- matrix(NA, nrow = length(unique(data$STUDY_ID)), ncol = 12)
for(i in 1:length(sub)){
for(j in 4:13){
study[i,1] <- mean(sub[[i]][,1], na.rm = TRUE)
study[i,2] <- mean(sub[[i]][,2], na.rm = TRUE)
study[i,(j-1)] <- mean(sub[[i]][,j], na.rm = TRUE)
}
}
colnames(study)<- colnames(data)[1:12]
head(study)## STUDY_ID N MH_ACE PPT_ACE NPT_ACE SE_ACE PPT_MH NPT_MH
## [1,] 1 1188 0.2000000 NaN NaN NaN NaN NaN
## [2,] 2 126 0.3700000 NaN NaN NaN NaN NaN
## [3,] 3 417 0.3459234 NaN NaN NaN NaN NaN
## [4,] 4 1185 0.2990000 0.09400000 NaN NaN -0.135 NaN
## [5,] 5 130 0.1547689 NaN NaN -0.0210000 NaN NaN
## [6,] 6 295 0.2500000 0.01047176 0.1415346 0.0408865 0.030 0.4472494
## SE_MH NPT_PPT SE_PPT SE_NPT
## [1,] NaN NaN NaN NaN
## [2,] NaN NaN NaN NaN
## [3,] NaN NaN NaN NaN
## [4,] NaN NaN NaN NaN
## [5,] 0.2810138 NaN NaN NaN
## [6,] NaN -0.06227524 -0.06289462 NaNnvar <- 5
varnames <- c("ACE","MH","PPT","NPT","SE")
labels <- list(varnames,varnames)
cormatrices <- readstack(study[,3:12], no.var = nvar,
var.names = varnames, diag = FALSE)
n <- study[,"N"]
my.df <- Cor2DataFrame(cormatrices, n, acov = "weighted")
pattern.na(cormatrices, show.na=FALSE)## ACE MH PPT NPT SE
## ACE 42 32 14 18 20
## MH 32 32 9 11 15
## PPT 14 9 14 6 6
## NPT 18 11 6 18 8
## SE 20 15 6 8 20pattern.n(cormatrices, n=n)## ACE MH PPT NPT SE
## ACE 22174.021 17219.021 8711.167 13630.355 13481.021
## MH 17219.021 17219.021 7810.167 9363.855 10232.021
## PPT 8711.167 7810.167 8711.167 5834.000 5226.500
## NPT 13630.355 9363.855 5834.000 13630.355 9791.355
## SE 13481.021 10232.021 5226.500 9791.355 13481.021temp <- aggregate(PARENT_AGE ~ STUDY_ID, data = data, mean, na.rm = T)
mean(temp$PARENT_AGE); sd(temp$PARENT_AGE)## [1] 28.62071## [1] 2.344949temp$PARENT_AGE <- scale(temp$PARENT_AGE)
study <- merge(study, temp, by = "STUDY_ID")
my.df$data <- data.frame(my.df$data, PARENT_AGE = study$PARENT_AGE,
check.names=FALSE)
Ax <- matrix (c(0,0,0,0,0,
"0*data.PARENT_AGE",0,0,0,0,
"0*data.PARENT_AGE",0,0,0,0,
"0*data.PARENT_AGE",0,0,0,0,
"0*data.PARENT_AGE","0*data.PARENT_AGE","0*data.PARENT_AGE","0*data.PARENT_AGE",0),
nrow=5,ncol=5,byrow = TRUE); Ax## [,1] [,2] [,3]
## [1,] "0" "0" "0"
## [2,] "0*data.PARENT_AGE" "0" "0"
## [3,] "0*data.PARENT_AGE" "0" "0"
## [4,] "0*data.PARENT_AGE" "0" "0"
## [5,] "0*data.PARENT_AGE" "0*data.PARENT_AGE" "0*data.PARENT_AGE"
## [,4] [,5]
## [1,] "0" "0"
## [2,] "0" "0"
## [3,] "0" "0"
## [4,] "0" "0"
## [5,] "0*data.PARENT_AGE" "0"M1 <- create.vechsR(A0=RAM1$A, S0=RAM1$S, Ax=Ax)
mx.fit0 <- osmasem(model.name="No moderator", Mmatrix=M0, Tmatrix=T0,
data=my.df, intervals.type = "LB"
,
mxModel.Args = list(ind1, ind2,ind3, mxCI(c("ACE_MH_SE","ACE_PPT_SE","ACE_NPT_SE")))
)
mx.fit8 <- osmasem(model.name = 'Parent Age (z) as Moderator', Mmatrix = M1, Tmatrix = T0, data = my.df)
anova(mx.fit8,mx.fit0)## base comparison ep minus2LL df AIC diffLL
## 1 Parent Age (z) as Moderator <NA> 27 -173.3146 112 -119.3146 NA
## 2 Parent Age (z) as Moderator No moderator 20 -156.2352 119 -116.2352 17.07946
## diffdf p
## 1 NA NA
## 2 7 0.01689094summary(mx.fit8)$parameters [, c(1,5,6,12)] # View only name, estimates, se, and p value## name Estimate Std.Error Pr(>|z|)
## 1 b21 0.28166680 0.01835480 0.000000e+00
## 2 b31 -0.07248522 0.03321004 2.906311e-02
## 3 b41 0.13726816 0.02800639 9.519733e-07
## 4 b51 0.07329010 0.03460375 3.417664e-02
## 5 b52 0.17512570 0.03761145 3.221288e-06
## 6 b53 -0.16603319 0.05633847 3.208044e-03
## 7 b54 0.18207355 0.04631821 8.461857e-05
## 8 MHWITHPPT -0.08185906 0.06028737 1.745225e-01
## 9 MHWITHNPT 0.29567857 0.04101343 5.624390e-13
## 10 PPTWITHNPT -0.16418543 0.14300950 2.509382e-01
## 11 b21_1 -0.06767918 0.01888837 3.395234e-04
## 12 b31_1 0.02499744 0.03366403 4.577505e-01
## 13 b41_1 -0.02789970 0.02956852 3.453941e-01
## 14 b51_1 -0.03657532 0.03277315 2.644158e-01
## 15 b52_1 0.01824200 0.03435181 5.953947e-01
## 16 b53_1 -0.10371651 0.05701926 6.891550e-02
## 17 b54_1 -0.04212398 0.05140149 4.124961e-01
## 18 Tau1_1 -4.96676221 0.36701227 0.000000e+00
## 19 Tau1_2 -4.59446503 0.54305715 0.000000e+00
## 20 Tau1_3 -4.62713149 0.47716162 0.000000e+00
## 21 Tau1_4 -4.31100491 0.37772554 0.000000e+00
## 22 Tau1_5 -3.59967715 0.66258349 5.547960e-08
## 23 Tau1_6 -4.24812144 0.55182424 1.376677e-14
## 24 Tau1_7 -5.04851566 0.56865296 0.000000e+00
## 25 Tau1_8 -2.13861362 0.60365121 3.959076e-04
## 26 Tau1_9 -5.10328990 0.85304869 2.198554e-09
## 27 Tau1_10 -5.14330935 0.63728993 6.661338e-16