# Random Effects Model - forest, funnel, summary (interpret negative symbol as positive)
meta <- escalc(measure="SMCR",m1i = mydata$i.wb.m.t1, m2i=mydata$i.wb.m.t2, sd1i= mydata$i.wb.sd.t1, sd2i = mydata$i.wb.sd.t2, ni= mydata$i.sample.size.t1, ri=ri, data=mydata)

# "SMCR" for the standardized mean change using raw score standardization (Becker, 1988)
mc <- rma.uni(yi,vi, data=meta,  slab = meta$study.id)
summary(mc)
## 
## Random-Effects Model (k = 12; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc   
## -12.9434   25.8869   29.8869   30.6826   31.3869   
## 
## tau^2 (estimated amount of total heterogeneity): 0.5194 (SE = 0.2253)
## tau (square root of estimated tau^2 value):      0.7207
## I^2 (total heterogeneity / total variability):   99.80%
## H^2 (total variability / sampling variability):  502.71
## 
## Test for Heterogeneity:
## Q(df = 11) = 617.5099, p-val < .0001
## 
## Model Results:
## 
## estimate      se     zval    pval    ci.lb    ci.ub    
##  -0.5281  0.2099  -2.5160  0.0119  -0.9395  -0.1167  * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(mc)
## 
##        estimate    ci.lb     ci.ub 
## tau^2    0.5194   0.2627    1.9334 
## tau      0.7207   0.5126    1.3905 
## I^2(%)  99.8011  99.6075   99.9465 
## H^2    502.7144 254.7680 1868.4379
# PPIs improve wellbeing by about .53 on average (pvalue = .01), 95% CI = -94, -.12

#Creating vectors to keep figures consistent (this is for forest)
xl=c(-20, 20)
al=c(-10,10)
cx=1.5

forest(mc, alim = xl, xlim=al,header=TRUE)
text(-4.5, -1, pos=1, cex=.75, bquote(paste(
  "(estimate = -.53",
  ", se = ", .(mc$se), ", ",
  .(fmtp2(mc$pval)), "; ",
  I^2, " = ", .(fmtx(mc$I2, digits=1)), "%, ",
  tau^2, " = ", .(fmtx(mc$tau2, digits=2)), ")")))

#Creating vectors to keep figures consistent (this is for funnel)
xl=c(-4, 4)
al=c(-10,10)
cx=0.75

funnel(mc, xlim = xl, cex.lab=cx,cex.axis=cx,back='White',offset=1.5,shade="lightgrey",hlines="lightgrey", cex=cx, slab=mydata$study.id, label=FALSE, legend="topleft")

# categorical moderator
# frequency
m3 <- rma(yi, vi, data=meta,  slab = meta$study.id, mods = ~ frequency)
summary(m3)
## 
## Mixed-Effects Model (k = 12; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc   
##  -8.8069   17.6139   29.6139   29.2893  113.6139   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.6724 (SE = 0.3642)
## tau (square root of estimated tau^2 value):             0.8200
## I^2 (residual heterogeneity / unaccounted variability): 99.74%
## H^2 (unaccounted variability / sampling variability):   378.43
## R^2 (amount of heterogeneity accounted for):            0.00%
## 
## Test for Residual Heterogeneity:
## QE(df = 7) = 213.5241, p-val < .0001
## 
## Test of Moderators (coefficients 2:5):
## QM(df = 4) = 2.6930, p-val = 0.6104
## 
## Model Results:
## 
##                                   estimate      se     zval    pval    ci.lb 
## intrcpt                            -0.4061  0.4108  -0.9887  0.3228  -1.2112 
## frequencyMore than once per week   -1.0202  0.7276  -1.4023  0.1608  -2.4462 
## frequencyOnce                       0.1830  0.7115   0.2573  0.7970  -1.2115 
## frequencyVaried                    -0.0239  0.9176  -0.0261  0.9792  -1.8224 
## frequencyWeekly                     0.0108  0.6273   0.0172  0.9863  -1.2188 
##                                    ci.ub    
## intrcpt                           0.3990    
## frequencyMore than once per week  0.4057    
## frequencyOnce                     1.5776    
## frequencyVaried                   1.7745    
## frequencyWeekly                   1.2403    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(m3)
## 
##        estimate    ci.lb     ci.ub 
## tau^2    0.6724   0.2848    3.1303 
## tau      0.8200   0.5336    1.7693 
## I^2(%)  99.7357  99.3782   99.9431 
## H^2    378.4262 160.8316 1758.0524
# intercept reflects the estimated average change for daily activity use (this is the reference level)
# the other two coefficients estimate how much the average change for remaining frequencies differ from the reference level (daily)
# not significant, so frequency doesn't shape the effectiveness of the PPI

#continent
m4 <- rma(yi, vi, data=meta,  slab = meta$study.id, mods = ~ continent)
summary(m4)
## 
## Mixed-Effects Model (k = 12; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc   
##  -9.4584   18.9168   28.9168   29.3140   58.9168   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.5547 (SE = 0.2819)
## tau (square root of estimated tau^2 value):             0.7448
## I^2 (residual heterogeneity / unaccounted variability): 99.84%
## H^2 (unaccounted variability / sampling variability):   607.68
## R^2 (amount of heterogeneity accounted for):            0.00%
## 
## Test for Residual Heterogeneity:
## QE(df = 8) = 543.8498, p-val < .0001
## 
## Test of Moderators (coefficients 2:4):
## QM(df = 3) = 3.1864, p-val = 0.3638
## 
## Model Results:
## 
##                         estimate      se     zval    pval    ci.lb    ci.ub     
## intrcpt                  -1.2110  0.4422  -2.7383  0.0062  -2.0778  -0.3442  ** 
## continentAustralia        0.8496  0.8691   0.9776  0.3283  -0.8538   2.5529     
## continentEurope           0.9269  0.5245   1.7672  0.0772  -0.1011   1.9548   . 
## continentSouth America    0.7289  0.8672   0.8406  0.4006  -0.9707   2.4285     
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(m4)
## 
##        estimate    ci.lb     ci.ub 
## tau^2    0.5547   0.2458    2.4342 
## tau      0.7448   0.4958    1.5602 
## I^2(%)  99.8354  99.6295   99.9625 
## H^2    607.6848 269.8888 2663.4393
#no difference in wellbeing improvement by continent

#sample type
m5 <- rma(yi, vi, data=meta,  slab = meta$study.id, mods = ~ sample.type)
summary(m5)
## 
## Mixed-Effects Model (k = 12; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc   
## -11.3777   22.7554   30.7554   31.5443   40.7554   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.6459 (SE = 0.3093)
## tau (square root of estimated tau^2 value):             0.8037
## I^2 (residual heterogeneity / unaccounted variability): 99.81%
## H^2 (unaccounted variability / sampling variability):   527.83
## R^2 (amount of heterogeneity accounted for):            0.00%
## 
## Test for Residual Heterogeneity:
## QE(df = 9) = 374.8240, p-val < .0001
## 
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 0.4931, p-val = 0.7815
## 
## Model Results:
## 
##                          estimate      se     zval    pval    ci.lb    ci.ub    
## intrcpt                   -0.6839  0.3323  -2.0577  0.0396  -1.3353  -0.0325  * 
## sample.typeOlder adults    0.4370  0.6587   0.6634  0.5071  -0.8541   1.7281    
## sample.typeStudents        0.2250  0.5222   0.4309  0.6665  -0.7985   1.2486    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(m5)
## 
##        estimate    ci.lb     ci.ub 
## tau^2    0.6459   0.3002    2.6164 
## tau      0.8037   0.5479    1.6175 
## I^2(%)  99.8105  99.5932   99.9532 
## H^2    527.8295 245.8272 2135.0167
describeBy(meta$yi, meta$sample.type)
## 
##  Descriptive statistics by group 
## group: General/Community
##    vars n  mean  sd median trimmed  mad   min  max range  skew kurtosis   se
## X1    1 6 -0.74 1.2   -0.4   -0.74 0.29 -3.16 0.08  3.24 -1.26    -0.24 0.49
## ------------------------------------------------------------ 
## group: Older adults
##    vars n  mean   sd median trimmed  mad   min   max range skew kurtosis   se
## X1    1 2 -0.25 0.13  -0.25   -0.25 0.14 -0.34 -0.15  0.19    0    -2.75 0.09
## ------------------------------------------------------------ 
## group: Students
##    vars n  mean   sd median trimmed  mad   min   max range skew kurtosis   se
## X1    1 4 -0.46 0.08  -0.47   -0.46 0.08 -0.54 -0.36  0.18 0.15    -2.05 0.04
#wellbeing improvement strongest for general/community sample (compared to older adults and students)

#age
m6 <- rma(yi, vi, data=meta,  slab = meta$study.id, mods = ~ age.mean)
summary(m6)
## 
## Mixed-Effects Model (k = 12; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc   
## -12.2825   24.5649   30.5649   31.4727   34.5649   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.5873 (SE = 0.2670)
## tau (square root of estimated tau^2 value):             0.7664
## I^2 (residual heterogeneity / unaccounted variability): 99.82%
## H^2 (unaccounted variability / sampling variability):   554.20
## R^2 (amount of heterogeneity accounted for):            0.00%
## 
## Test for Residual Heterogeneity:
## QE(df = 10) = 617.1709, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## QM(df = 1) = 0.0710, p-val = 0.7899
## 
## Model Results:
## 
##           estimate      se     zval    pval    ci.lb   ci.ub    
## intrcpt    -0.6693  0.5623  -1.1902  0.2340  -1.7714  0.4329    
## age.mean    0.0037  0.0138   0.2664  0.7899  -0.0233  0.0307    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(m6)
## 
##        estimate    ci.lb     ci.ub 
## tau^2    0.5873   0.2851    2.2725 
## tau      0.7664   0.5339    1.5075 
## I^2(%)  99.8196  99.6290   99.9533 
## H^2    554.2021 269.5299 2141.6029
#no age differences
## RESOURCES ##
# https://training.cochrane.org/handbook/current/chapter-15

# estimate of the average true standardized mean difference is .53
# Q-test suggests that the underlying true standardized mean differences are heterogeneous
# can put weighted=false in the rma.uni code, which would give us an estimate reflecting unweighted average of 
# the true standardized mean differences of these studies

# for continuous outcome measures, we can present summary results for studies using natural units of measurement 
# or as minimal important differences when all studies use the same scale (e.g. SWL)
# our studies measure the same construct but with different scales, so need to find a way to interpret the 
# standardized mean difference

# resources for interpreting figures
# https://s4be.cochrane.org/blog/2023/06/27/how-to-read-a-funnel-plot/
# https://www.google.com/url?sa=i&url=https%3A%2F%2Fs4be.cochrane.org%2Fblog%2F2016%2F07%2F11%2Ftutorial-read-forest-plot%2F&psig=AOvVaw3QnJDUY9cHQ1Gc4fa1jKN0&ust=1739462356615000&source=images&cd=vfe&opi=89978449&ved=0CBQQjRxqFwoTCMCbzrHAvosDFQAAAAAdAAAAABAE
# https://www.google.com/url?sa=i&url=https%3A%2F%2Fbookdown.org%2FMathiasHarrer%2FDoing_Meta_Analysis_in_R%2Fforest.html&psig=AOvVaw3QnJDUY9cHQ1Gc4fa1jKN0&ust=1739462356615000&source=images&cd=vfe&opi=89978449&ved=0CBQQjRxqFwoTCMCbzrHAvosDFQAAAAAdAAAAABAJ

### forest plot
# line is of a null effect (i.e., 0 - no change in wellbeing)
# black box is the point estimate of the effect size (smc)
# any black box that crosses the null effect line is n.s. (none in this case)
# diamond (at R.E.M. line) is the point estimate and CIs when you combine and 
# average all individual studies together 
# I2 tells us heterogeneity of the studies (how consistent they are) 
# reflects variability of true effect sizes
# over 50% means studies may be inconsistent due to a reason other than chance
# tau2 is the estimated standard deviation of the distribution of effect sizes
# quantifies proportion of total variation attributable to heterogenity

### funnel plot
# 95% confidence interval is the triangle lines
# y axis is study precision (higher up the dot, higher precision) 
# x axis is the estimated effect size (standardized mean change)
# 95% of studies would be expected to be within the triangle but because the 
# rule of thumb is to only interpret a funnel plot's asymmetry with 10 or more studies