library(metafor)
## Loading required package: Matrix
## Loading required package: metadat
## Loading required package: numDeriv
##
## Loading the 'metafor' package (version 4.4-0). For an
## introduction to the package please type: help(metafor)
library(ggplot2)
load data
dt = read.csv("dt/filtered_dt_2024.csv")
dt$author_name = paste(dt$Author, dt$Date, sep = " et al., ")
dt
smd_meta.psqi <- escalc(measure="SMD",
m1i=PSQI_concussion, m2i=PSQI_control,
sd1i=PSQI_concussion_sd, sd2i=PSQI_control_sd,
n1i=N_concussion,
n2i=N_controls,
slab = author_name, data=dt)
reml.psqi = rma(yi,vi,method = "REML", data = smd_meta.psqi)
summary(reml.psqi)
##
## Random-Effects Model (k = 6; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -3.5439 7.0879 11.0879 10.3067 17.0879
##
## tau^2 (estimated amount of total heterogeneity): 0.1749 (SE = 0.1712)
## tau (square root of estimated tau^2 value): 0.4182
## I^2 (total heterogeneity / total variability): 69.86%
## H^2 (total variability / sampling variability): 3.32
##
## Test for Heterogeneity:
## Q(df = 5) = 20.6706, p-val = 0.0009
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.8896 0.2138 4.1604 <.0001 0.4705 1.3087 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot
forest(reml.psqi)
text(-1.6, 18, "Author(s), Year", pos=4, cex=.8)
text( 1.6, 18, "Correlation [95% CI]", pos=2, cex=.8)
library(meta)
## Loading 'meta' package (version 7.0-0).
## Type 'help(meta)' for a brief overview.
## Readers of 'Meta-Analysis with R (Use R!)' should install
## older version of 'meta' package: https://tinyurl.com/dt4y5drs
dt.metacont = metacont(N_concussion, PSQI_concussion, PSQI_concussion_sd,
N_controls, PSQI_control, PSQI_control_sd,
comb.fixed = T, comb.random = T, studlab = author_name,
data = dt, sm = "SMD")
dt.metacont
## Number of studies: k = 6
## Number of observations: o = 520 (o.e = 191, o.c = 329)
##
## SMD 95%-CI z p-value
## Common effect model 0.7021 [0.4963; 0.9079] 6.69 < 0.0001
## Random effects model 0.8878 [0.4679; 1.3077] 4.14 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.1746 [0.0238; 1.1229]; tau = 0.4179 [0.1544; 1.0597]
## I^2 = 75.6% [45.1%; 89.2%]; H = 2.03 [1.35; 3.04]
##
## Test of heterogeneity:
## Q d.f. p-value
## 20.51 5 0.0010
##
## Details on meta-analytical method:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-Profile method for confidence interval of tau^2 and tau
## - Hedges' g (bias corrected standardised mean difference; using exact formulae)
plot
forest(dt.metacont, leftcols = c('studlab'))
funnel(dt.metacont)
dt_reformatted = dt[,c("author_name","PSQI_control","PSQI_control_sd","Mean.Age.in.Controls..years.","Sex.Female.n.....control","bias")]
colnames(dt_reformatted) <- c("author_name","PSQI","sd","age","female_ratio","bias")
dt_reformatted$concussed = FALSE
dt_reformatted_tmp = dt[,c("author_name","PSQI_concussion","PSQI_concussion_sd","Mean.Age.in.Concusssed..years.","Sex..Female.n.....concussed","bias")]
colnames(dt_reformatted_tmp) <- c("author_name","PSQI","sd","age","female_ratio","bias")
dt_reformatted_tmp$concussed = TRUE
dt_reformatted = cbind(dt_reformatted, dt_reformatted_tmp)
dt_reformatted$variance = dt_reformatted$sd^2
rma.da.adjusted <- rma(yi=PSQI, vi=variance, data = dt_reformatted, method = "DL",
mods = ~ age + female_ratio + bias)
rma.da.adjusted
##
## Mixed-Effects Model (k = 6; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0 (SE = 5.6262)
## tau (square root of estimated tau^2 value): 0
## I^2 (residual heterogeneity / unaccounted variability): 0.00%
## H^2 (unaccounted variability / sampling variability): 1.00
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 2) = 0.0309, p-val = 0.9847
##
## Test of Moderators (coefficients 2:4):
## QM(df = 3) = 0.8244, p-val = 0.8436
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 3.0951 13.4020 0.2309 0.8174 -23.1725 29.3626
## age 0.1593 0.6774 0.2351 0.8141 -1.1683 1.4868
## female_ratio 2.9145 6.9498 0.4194 0.6750 -10.7070 16.5359
## bias -1.0138 1.2156 -0.8340 0.4043 -3.3964 1.3688
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
regplot(rma.da.adjusted, mod = "age", pi=TRUE)
pdf("fig/rma_adjusted_age.pdf", width = 4, height = 4)
regplot(rma.da.adjusted, mod = "female_ratio")
pdf("fig/rma_adjusted_gender.pdf", width = 4, height = 4)
regplot(rma.da.adjusted, mod = "bias", pi=TRUE)