1 Aggregate Results

#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)

2 Data Preparation

### 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") 

3 Plot Data

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'))

3.1 BCLSM

3.1.1 Meta-Analysis b sigma negative

# 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")

3.1.2 Meta-Analysis b sigma positive

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

3.1.3 Meta-Analysis b mean negative

# 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

3.1.4 Meta-Analysis b mean positive

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

3.2 RVI

3.2.1 Meta-Analysis unweighted RVI negative

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 ")))

3.2.2 Meta-Analysis unweighted RVI positive

#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 ")))

3.3 RVI with weights

3.3.1 Meta-Analysis weighted RVI negative

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 ")))

3.3.2 Meta-Analysis weighted RVI positive

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 ")))

3.4 Gender Effects

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

3.5 Model comparison

3.5.1 censoring vs. no censoring

3.5.1.1 no censoring (negative Emotion)

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 "))) 

3.5.1.2 no censoring (positive Emotion)

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 "))) 

3.5.2 model C (two part)

3.5.2.1 negative Emotion

#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 "))) 

3.5.2.2 positive Emotion

#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 ")))