#create list with files from Kalokerinos Study
kalok_files = paste0("Dataset ", 1:11, " public.csv.xlsx")
# add Dataset 12
kalok_files <- append(kalok_files, "results_wd.xlsx" )
#add Dataset 13
kalok_files <- append(kalok_files, "results_LS.xlsx")
#read data
K_data <- lapply(kalok_files, read_excel)
# adding column with name of dataset
for( i in 1:13){
K_data[[i]]$Dataset <- NA
}
for( i in 1:13){
K_data[[i]]$Dataset <- rep(i, 7)
}
K_data <- bind_rows(K_data)
#add author names
K_data$Dataset[K_data$Dataset == 1] <- "Van Ryckeghem et al."
K_data$Dataset[K_data$Dataset == 2] <- "Dejonckheere et al."
K_data$Dataset[K_data$Dataset == 3] <- "Kalokerinos et al."
K_data$Dataset[K_data$Dataset == 4] <- "Grommisch et al."
K_data$Dataset[K_data$Dataset == 5] <- "Erbas et al."
K_data$Dataset[K_data$Dataset == 6] <- "Brans et al."
K_data$Dataset[K_data$Dataset == 7] <- "Holland et al."
K_data$Dataset[K_data$Dataset == 8] <- "Koval et al."
K_data$Dataset[K_data$Dataset == 9] <- "Koval et al."
K_data$Dataset[K_data$Dataset == 10] <- "Dejonckheere et al."
K_data$Dataset[K_data$Dataset == 11] <- "Kalokerinos et al."
K_data$Dataset[K_data$Dataset == 12] <- "Denissen et al."
K_data$Dataset[K_data$Dataset == 13] <- "Own diary data"
K_data <- bind_rows(K_data)
### seperate columns in value, est.error and CI
K_data <- K_data %>%
tidyr::separate(negemo_b_neuro,
c("neg_b_neuro", "est.error_n_b_neuro", "KIn_b_neuro"), sep = " ")
K_data <-K_data %>%
tidyr::separate(negemo_b_neuro_sigma,
c("neg_b_neuro_sigma", "est.error_n_b_neuro_sigma", "KIn_b_neuro_sigma"), sep = " ")
K_data <-K_data %>%
tidyr::separate(posemo_b_neuro,
c("pos_b_neuro", "est.error_p_b_neuro", "KIp_b_neuro"), sep = " ")
K_data <-K_data %>%
tidyr::separate(posemo_b_neuro_sigma,
c("pos_b_neuro_sigma", "est.error_p_b_neuro_sigma", "KIp_b_neuro_sigma"), sep = " ")
K_data <-K_data %>%
tidyr::separate(b_neg_sigma_sex,
c("b_neg_sigma_sex", "est.error_b_neg_sigma_sex", "KIb_neg_sigma_sex"), sep = " ")
K_data <-K_data %>%
tidyr::separate(b_pos_sigma_sex,
c("b_pos_sigma_sex", "est.error_b_pos_sigma_sex", "KIb_pos_sigma_sex"), sep = " ")
# save Table
library("writexl")
write_xlsx(K_data,"~/new_results2.xlsx")
## Filter Dataset for plots
K_b_sigma <- K_data %>% filter(model== "model3")
K_RSD <- K_data %>% filter(model == "RSD")
K_RSD_w <- K_data %>% filter(model == "RSD_weight")
K_gender <- K_data %>% filter(model == "gender")
K_SD <- K_data %>% filter(model == "SD")
apatheme = theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
panel.border = element_blank(),
text=element_text(family='Arial'),
legend.title=element_blank(),
legend.position=c(.7,.8),
axis.line.x = element_line(color='black'),
axis.line.y = element_line(color='black'))
# 1. b sigma
K_b_sigma$est.error_n_b_neuro_sigma <- as.numeric( K_b_sigma$est.error_n_b_neuro_sigma)
K_b_sigma$neg_b_neuro_sigma <- as.numeric( K_b_sigma$neg_b_neuro_sigma)
mod_b_sigma <- rma(yi= neg_b_neuro_sigma, sei= est.error_n_b_neuro_sigma, data= K_b_sigma )
summary(mod_b_sigma)
##
## Random-Effects Model (k = 13; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 17.5255 -35.0511 -31.0511 -30.0813 -29.7177
##
## tau^2 (estimated amount of total heterogeneity): 0.0013 (SE = 0.0012)
## tau (square root of estimated tau^2 value): 0.0354
## I^2 (total heterogeneity / total variability): 48.02%
## H^2 (total variability / sampling variability): 1.92
##
## Test for Heterogeneity:
## Q(df = 12) = 20.1979, p-val = 0.0634
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0974 0.0155 6.2937 <.0001 0.0670 0.1277 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest(mod_b_sigma, slab =c(paste0(1:13,".", " ", K_b_sigma$Dataset)),
mlab = (paste("Random Effects Model Study 1-13 ")))
par(bg="white")
K_b_sigma$est.error_p_b_neuro_sigma <- as.numeric( K_b_sigma$est.error_p_b_neuro_sigma)
K_b_sigma$pos_b_neuro_sigma <- as.numeric( K_b_sigma$pos_b_neuro_sigma)
mod_b_sigma2 <- rma(yi= pos_b_neuro_sigma, sei= est.error_p_b_neuro_sigma, data= K_b_sigma )
summary(mod_b_sigma2)
##
## Random-Effects Model (k = 13; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 16.6627 -33.3254 -29.3254 -28.3556 -27.9920
##
## tau^2 (estimated amount of total heterogeneity): 0.0023 (SE = 0.0015)
## tau (square root of estimated tau^2 value): 0.0480
## I^2 (total heterogeneity / total variability): 74.31%
## H^2 (total variability / sampling variability): 3.89
##
## Test for Heterogeneity:
## Q(df = 12) = 53.6705, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0427 0.0173 2.4727 0.0134 0.0089 0.0766 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pos_sigma <- forest(mod_b_sigma2, slab =c(paste0(1:13,".", " ", K_b_sigma$Dataset)),
xlim=c(-0.6,0.6), at=seq(-0.2,0.3,by=.1),
mlab = (paste("Random Effects Model Study 1-13 ")))
pos_sigma
## $xlim
## [1] -0.6 0.6
##
## $alim
## [1] -0.2 0.3
##
## $at
## [1] -0.2 -0.1 0.0 0.1 0.2 0.3
##
## $ylim
## [1] -1.5 16.0
##
## $rows
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13
##
## $cex
## [1] 1
##
## $cex.lab
## [1] 1
##
## $cex.axis
## [1] 1
##
## $ilab.xpos
## NULL
##
## $ilab.pos
## NULL
##
## $textpos
## [1] -0.6 0.6
# b mean
K_b_sigma$est.error_n_b_neuro <- as.numeric( K_b_sigma$est.error_n_b_neuro)
K_b_sigma$neg_b_neuro <- as.numeric( K_b_sigma$neg_b_neuro)
mod_b <- rma(yi= neg_b_neuro, sei= est.error_n_b_neuro, data= K_b_sigma )
forest(mod_b, slab =c(paste0(1:13,".", " ", K_b_sigma$Dataset)),
mlab = (paste("Random Effects Model Study 1-13 ")))
summary(mod_b)
##
## Random-Effects Model (k = 13; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 11.4488 -22.8975 -18.8975 -17.9277 -17.5642
##
## tau^2 (estimated amount of total heterogeneity): 0.0063 (SE = 0.0038)
## tau (square root of estimated tau^2 value): 0.0791
## I^2 (total heterogeneity / total variability): 78.79%
## H^2 (total variability / sampling variability): 4.72
##
## Test for Heterogeneity:
## Q(df = 12) = 78.9647, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.2458 0.0272 9.0461 <.0001 0.1925 0.2990 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
K_b_sigma$est.error_p_b_neuro <- as.numeric( K_b_sigma$est.error_p_b_neuro)
K_b_sigma$pos_b_neuro <- as.numeric( K_b_sigma$pos_b_neuro)
mod_b2 <- rma(yi= pos_b_neuro, sei= est.error_p_b_neuro, data= K_b_sigma )
forest(mod_b2, slab =c(paste0(1:13,".", " ", K_b_sigma$Dataset)),
mlab = (paste("Random Effects Model Study 1-13 ")))
summary(mod_b2)
##
## Random-Effects Model (k = 13; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 10.4089 -20.8179 -16.8179 -15.8480 -15.4845
##
## tau^2 (estimated amount of total heterogeneity): 0.0052 (SE = 0.0034)
## tau (square root of estimated tau^2 value): 0.0723
## I^2 (total heterogeneity / total variability): 72.39%
## H^2 (total variability / sampling variability): 3.62
##
## Test for Heterogeneity:
## Q(df = 12) = 45.8999, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.2389 0.0262 -9.1255 <.0001 -0.2902 -0.1876 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
K_RSD$neg_b_neuro_sigma<- as.numeric( K_RSD$neg_b_neuro_sigma)
K_RSD$est.error_n_b_neuro_sigma<- as.numeric( K_RSD$est.error_n_b_neuro_sigma)
K_RSD$pos_b_neuro<- as.numeric( K_RSD$pos_b_neuro)
K_RSD$est.error_p_b_neuro<- as.numeric( K_RSD$est.error_p_b_neuro)
#negative
mod_RSD <- rma(yi= neg_b_neuro_sigma, vi= est.error_n_b_neuro_sigma, data= K_RSD )
summary(mod_RSD)
##
## Random-Effects Model (k = 13; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 8.5811 -17.1622 -13.1622 -12.1924 -11.8289
##
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0112)
## tau (square root of estimated tau^2 value): 0
## I^2 (total heterogeneity / total variability): 0.00%
## H^2 (total variability / sampling variability): 1.00
##
## Test for Heterogeneity:
## Q(df = 12) = 0.6464, p-val = 1.0000
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0069 0.0486 0.1420 0.8871 -0.0883 0.1021
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest(mod_RSD, slab =c(paste0(1:13,".", " ", K_b_sigma$Dataset)),xlim = c(-1.68, 1.58), at=seq(-0.6,0.6,by=.2),
mlab = (paste("Random Effects Model Study 1-13 ")))
#positive
mod_RSDp <- rma(yi= pos_b_neuro, vi= est.error_p_b_neuro, data= K_RSD )
summary(mod_RSDp)
##
## Random-Effects Model (k = 13; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 8.5439 -17.0879 -13.0879 -12.1181 -11.7546
##
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0109)
## tau (square root of estimated tau^2 value): 0
## I^2 (total heterogeneity / total variability): 0.00%
## H^2 (total variability / sampling variability): 1.00
##
## Test for Heterogeneity:
## Q(df = 12) = 0.8759, p-val = 1.0000
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0383 0.0479 0.7987 0.4244 -0.0556 0.1321
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest(mod_RSDp, slab =c(paste0(1:13,".", " ", K_b_sigma$Dataset)),
mlab = (paste("Random Effects Model Study 1-13 ")))
K_RSD_w$neg_b_neuro_sigma<- as.numeric( K_RSD_w$neg_b_neuro_sigma)
K_RSD_w$est.error_n_b_neuro_sigma<- as.numeric( K_RSD_w$est.error_n_b_neuro_sigma)
K_RSD_w$pos_b_neuro<- as.numeric( K_RSD_w$pos_b_neuro)
K_RSD_w$est.error_p_b_neuro<- as.numeric( K_RSD_w$est.error_p_b_neuro)
mod_RSD_w <- rma(yi= neg_b_neuro_sigma, sei= est.error_n_b_neuro_sigma, data= K_RSD_w )
summary(mod_RSD_w)
##
## Random-Effects Model (k = 13; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 20.4499 -40.8997 -36.8997 -35.9299 -35.5664
##
## tau^2 (estimated amount of total heterogeneity): 0.0008 (SE = 0.0005)
## tau (square root of estimated tau^2 value): 0.0278
## I^2 (total heterogeneity / total variability): 66.37%
## H^2 (total variability / sampling variability): 2.97
##
## Test for Heterogeneity:
## Q(df = 12) = 37.0590, p-val = 0.0002
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0327 0.0104 3.1383 0.0017 0.0123 0.0531 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest(mod_RSD_w, slab =c(paste0(1:13,".", " ", K_b_sigma$Dataset)),xlim = c(-1.68, 1.58), at=seq(-0.6,0.6,by=.2),
mlab = (paste("Random Effects Model Study 1-13 ")))
mod_RSD_w2 <- rma(yi= pos_b_neuro, sei= est.error_p_b_neuro, data= K_RSD_w )
summary(mod_RSD_w2)
##
## Random-Effects Model (k = 13; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 16.7812 -33.5623 -29.5623 -28.5925 -28.2290
##
## tau^2 (estimated amount of total heterogeneity): 0.0025 (SE = 0.0012)
## tau (square root of estimated tau^2 value): 0.0505
## I^2 (total heterogeneity / total variability): 92.07%
## H^2 (total variability / sampling variability): 12.60
##
## Test for Heterogeneity:
## Q(df = 12) = 143.9779, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0305 0.0154 1.9799 0.0477 0.0003 0.0608 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest(mod_RSD_w2, slab =c(paste0(1:13,".", " ", K_b_sigma$Dataset)), xlim = c(-0.4, 0.31), at=seq(-0.15,0.15,by=.05),
mlab = (paste("Random Effects Model Study 1-13 ")))
K_gender$neg_b_neuro <- as.numeric(K_gender$neg_b_neuro)
K_gender$est.error_n_b_neuro <- as.numeric(K_gender$est.error_n_b_neuro)
K_gender$neg_b_neuro_sigma <- as.numeric(K_gender$neg_b_neuro_sigma)
K_gender$est.error_n_b_neuro_sigma <- as.numeric(K_gender$est.error_n_b_neuro_sigma)
K_gender$b_neg_sigma_sex <- as.numeric(K_gender$b_neg_sigma_sex)
K_gender$est.error_b_neg_sigma_sex <- as.numeric(K_gender$est.error_b_neg_sigma_sex)
mod_sex2 <- rma(yi= neg_b_neuro_sigma, sei= est.error_n_b_neuro_sigma, data= K_gender )
## Warning: Studies with NAs omitted from model fitting.
summary(mod_sex2)
##
## Random-Effects Model (k = 10; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 15.1072 -30.2144 -26.2144 -25.8200 -24.2144
##
## tau^2 (estimated amount of total heterogeneity): 0.0019 (SE = 0.0010)
## tau (square root of estimated tau^2 value): 0.0434
## I^2 (total heterogeneity / total variability): 92.35%
## H^2 (total variability / sampling variability): 13.06
##
## Test for Heterogeneity:
## Q(df = 9) = 158.3208, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0927 0.0148 6.2460 <.0001 0.0636 0.1218 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest(mod_sex2, slab =c(paste0(1:13,".", " ", K_b_sigma$Dataset)),
mlab = (paste("Random Effects Model Study 1-13* ")))
mod_sex3 <- rma(yi= neg_b_neuro, sei= est.error_n_b_neuro, data= K_gender )
## Warning: Studies with NAs omitted from model fitting.
forest(mod_sex3, slab =c(paste0(1:13,".", " ", K_b_sigma$Dataset)),
mlab = (paste("Random Effects Model Study 1-13 ")))
summary(mod_sex3)
##
## Random-Effects Model (k = 10; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 9.8948 -19.7896 -15.7896 -15.3952 -13.7896
##
## tau^2 (estimated amount of total heterogeneity): 0.0046 (SE = 0.0034)
## tau (square root of estimated tau^2 value): 0.0678
## I^2 (total heterogeneity / total variability): 75.43%
## H^2 (total variability / sampling variability): 4.07
##
## Test for Heterogeneity:
## Q(df = 9) = 48.8704, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.2636 0.0274 9.6038 <.0001 0.2098 0.3174 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
noC <- Censo %>% filter(model== "modelnoCensoring")
apatheme = theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
panel.border = element_blank(),
text=element_text(family='Arial'),
legend.title=element_blank(),
legend.position=c(.7,.8),
axis.line.x = element_line(color='black'),
axis.line.y = element_line(color='black'))
names(noC)
## [1] "model" "negemo_b_neuro" "est.error_b_neg"
## [4] "KI_b_neg" "negemo_b_neuro_sigma" "est.error_b_neg_sigma"
## [7] "KI_b_neg_sigma" "posemo_b_neuro" "posemo_b_neuro_sigma"
## [10] "Dataset"
noC$negemo_b_neuro_sigma<- as.numeric(noC$negemo_b_neuro_sigma)
noC$Dataset <- as.numeric(noC$Dataset)
noC$est.error_b_neg_sigma <- as.numeric(noC$est.error_b_neg_sigma)
noC$negemo_b_neuro<- as.numeric(noC$negemo_b_neuro)
noC$Dataset <- as.numeric(noC$Dataset)
noC$est.error_b_neg <- as.numeric(noC$est.error_b_neg)
# b sigma without censoring
mod_b_sigma <- rma(yi= negemo_b_neuro_sigma, sei= est.error_b_neg_sigma, data= noC )
summary(mod_b_sigma)
##
## Random-Effects Model (k = 13; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 13.2993 -26.5985 -22.5985 -21.6287 -21.2652
##
## tau^2 (estimated amount of total heterogeneity): 0.0043 (SE = 0.0027)
## tau (square root of estimated tau^2 value): 0.0653
## I^2 (total heterogeneity / total variability): 73.59%
## H^2 (total variability / sampling variability): 3.79
##
## Test for Heterogeneity:
## Q(df = 12) = 50.8640, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.1675 0.0232 7.2244 <.0001 0.1220 0.2129 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest(mod_b_sigma, slab =c(paste0("Study ", noC$Dataset)),
mlab = (paste("RE Model Study 1-13 ")))
#b mean without censoring
mod_b_mean <- rma(yi= negemo_b_neuro, sei= est.error_b_neg, data= noC )
summary(mod_b_mean)
##
## Random-Effects Model (k = 13; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 12.5662 -25.1324 -21.1324 -20.1626 -19.7991
##
## tau^2 (estimated amount of total heterogeneity): 0.0054 (SE = 0.0031)
## tau (square root of estimated tau^2 value): 0.0735
## I^2 (total heterogeneity / total variability): 80.30%
## H^2 (total variability / sampling variability): 5.08
##
## Test for Heterogeneity:
## Q(df = 12) = 88.9269, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.2153 0.0246 8.7547 <.0001 0.1671 0.2635 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest(mod_b_mean, slab =c(paste0("Study ", noC$Dataset)),
mlab = (paste("RE Model Study 1-13 ")))
CensoP <- Censo %>%
tidyr::separate(posemo_b_neuro_sigma,
c("pos_b_neuro_sigma", "est.error_b_pos_sigma", "KI_b_pos_sigma"), sep = " ")
CensoP <- CensoP %>%
tidyr::separate(posemo_b_neuro,
c("pos_b_neuro", "est.error_b_pos", "KI_b_pos"), sep = " ")
write_xlsx(CensoP,"~/CensoP.xlsx")
noCp <- CensoP %>% filter(model== "modelnoCensoring")
names(noCp)
## [1] "model" "negemo_b_neuro" "est.error_b_neg"
## [4] "KI_b_neg" "negemo_b_neuro_sigma" "est.error_b_neg_sigma"
## [7] "KI_b_neg_sigma" "pos_b_neuro" "est.error_b_pos"
## [10] "KI_b_pos" "pos_b_neuro_sigma" "est.error_b_pos_sigma"
## [13] "KI_b_pos_sigma" "Dataset"
noCp$pos_b_neuro_sigma<- as.numeric(noCp$pos_b_neuro_sigma)
noCp$Dataset <- as.numeric(noCp$Dataset)
noCp$est.error_b_pos_sigma <- as.numeric(noCp$est.error_b_pos_sigma)
#b sigma without censoring
mod_b_sigmaP <- rma(yi= pos_b_neuro_sigma, sei= est.error_b_pos_sigma, data= noCp )
summary(mod_b_sigmaP)
##
## Random-Effects Model (k = 13; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 14.5194 -29.0388 -25.0388 -24.0690 -23.7055
##
## tau^2 (estimated amount of total heterogeneity): 0.0025 (SE = 0.0017)
## tau (square root of estimated tau^2 value): 0.0504
## I^2 (total heterogeneity / total variability): 75.19%
## H^2 (total variability / sampling variability): 4.03
##
## Test for Heterogeneity:
## Q(df = 12) = 58.4155, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0450 0.0184 2.4506 0.0143 0.0090 0.0810 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest(mod_b_sigmaP, slab =c(paste0("Study ", noCp$Dataset)),
mlab = (paste("RE Model Study 1-13 ")))
#Exclude study 11
quant <- quant[-c(21:22),]
apatheme = theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
panel.border = element_blank(),
text=element_text(family='Arial'),
legend.title=element_blank(),
legend.position=c(.7,.8),
axis.line.x = element_line(color='black'),
axis.line.y = element_line(color='black'))
names(quant)
## [1] "model" "b_neuro"
## [3] "est.error_b_neuro" "KI_b_neuro"
## [5] "b_neuro_sigma" "est.error_b_neuro_sigma"
## [7] "KI_b_neuro_sigma" "negemo_sigma"
## [9] "Study"
quant$b_neuro_sigma <- as.numeric(quant$b_neuro_sigma)
quant$Study <- as.numeric(quant$Study)
quant$est.error_b_neuro_sigma <- as.numeric(quant$est.error_b_neuro_sigma)
quant$b_neuro <- as.numeric(quant$b_neuro)
quant$est.error_b_neuro <- as.numeric(quant$est.error_b_neuro)
quant2 <- quant %>% filter(model== "BCSLM_q")
#b sigma two part model C:
mod_b_sigma <- rma(yi= b_neuro_sigma, sei= est.error_b_neuro_sigma, data= quant2 )
summary(mod_b_sigma)
##
## Random-Effects Model (k = 12; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 16.1752 -32.3504 -28.3504 -27.5546 -26.8504
##
## tau^2 (estimated amount of total heterogeneity): 0.0012 (SE = 0.0012)
## tau (square root of estimated tau^2 value): 0.0349
## I^2 (total heterogeneity / total variability): 47.32%
## H^2 (total variability / sampling variability): 1.90
##
## Test for Heterogeneity:
## Q(df = 11) = 18.2142, p-val = 0.0767
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0990 0.0159 6.2336 <.0001 0.0679 0.1301 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest(mod_b_sigma, slab =c(paste0("Study ", quant2$Study)),
mlab = (paste("RE Model Study 1-13 ")))
#b mean two part model C:
mod_b_mean <- rma(yi= b_neuro, sei= est.error_b_neuro, data= quant2 )
summary(mod_b_mean)
##
## Random-Effects Model (k = 12; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 9.8937 -19.7875 -15.7875 -14.9917 -14.2875
##
## tau^2 (estimated amount of total heterogeneity): 0.0071 (SE = 0.0043)
## tau (square root of estimated tau^2 value): 0.0840
## I^2 (total heterogeneity / total variability): 78.80%
## H^2 (total variability / sampling variability): 4.72
##
## Test for Heterogeneity:
## Q(df = 11) = 59.3091, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.2438 0.0294 8.3031 <.0001 0.1862 0.3013 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest(mod_b_mean, slab =c(paste0("Study ", quant2$Study)),
mlab = (paste("RE Model Study 1-13 ")))
#exclude study 2 and 3
quant <- quant[-c(3:6),]
quant$b_neuro_sigma <- as.numeric(quant$b_neuro_sigma)
quant$Study <- as.numeric(quant$Study)
quant$est.error_b_neuro_sigma <- as.numeric(quant$est.error_b_neuro_sigma)
quant$b_neuro <- as.numeric(quant$b_neuro)
quant$est.error_b_neuro <- as.numeric(quant$est.error_b_neuro)
quant2 <- quant %>% filter(model== "BCSLM_q")
#b sigma two part model C:
mod_b_sigma <- rma(yi= b_neuro_sigma, sei= est.error_b_neuro_sigma, data= quant2 )
summary(mod_b_sigma)
##
## Random-Effects Model (k = 11; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 14.7635 -29.5270 -25.5270 -24.9218 -23.8127
##
## tau^2 (estimated amount of total heterogeneity): 0.0019 (SE = 0.0015)
## tau (square root of estimated tau^2 value): 0.0440
## I^2 (total heterogeneity / total variability): 66.34%
## H^2 (total variability / sampling variability): 2.97
##
## Test for Heterogeneity:
## Q(df = 10) = 35.3863, p-val = 0.0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0382 0.0176 2.1704 0.0300 0.0037 0.0727 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest(mod_b_sigma, slab =c(paste0("Study ", quant2$Study)),
mlab = (paste("RE Model Study 1-13 ")))
#b mean two part model C:
mod_b_mean <- rma(yi= b_neuro, sei= est.error_b_neuro, data= quant2 )
summary(mod_b_mean)
##
## Random-Effects Model (k = 11; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 8.4083 -16.8167 -12.8167 -12.2115 -11.1024
##
## tau^2 (estimated amount of total heterogeneity): 0.0059 (SE = 0.0040)
## tau (square root of estimated tau^2 value): 0.0767
## I^2 (total heterogeneity / total variability): 74.87%
## H^2 (total variability / sampling variability): 3.98
##
## Test for Heterogeneity:
## Q(df = 10) = 42.1501, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.1678 0.0290 -5.7784 <.0001 -0.2247 -0.1109 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest(mod_b_mean, slab =c(paste0("Study ", quant2$Study)),
mlab = (paste("RE Model Study 1-13 ")))