Meta-Analysis (Hands-on Training)
GAYI KOMI SELASSI
2025-10-28
Load necessary libraries
library(metafor)## Warning: package 'metafor' was built under R version 4.4.3## Loading required package: Matrix## Loading required package: metadat## Warning: package 'metadat' was built under R version 4.4.3## Loading required package: numDeriv##
## Loading the 'metafor' package (version 4.8-0). For an
## introduction to the package please type: help(metafor)library(meta)## Warning: package 'meta' was built under R version 4.4.3## Loading 'meta' package (version 8.2-0).
## Type 'help(meta)' for a brief overview.##Data for meta-analysis
hypertension_data <- data.frame(
Author = c("Mahapatra et al., (2021)", "Santosh et al., (2019)",
"Basavanagowdappa et al., (2016)", "Boppana et al., (2023)",
"Annie Caroline et al., (2021)", "Bharatia et al., (2016)",
"Gupta et al., (2019)", "Mandal et al., (2014)",
"Umamaheswari et al., (2019)","Sagarad et al.,(2013)","Garg et al.,(2018)"),
Year = c(2021, 2019, 2016, 2023, 2021, 2016, 2019, 2014, 2019, 2013, 2018), # Pub. years
Events = c(31, 67, 248, 100, 130, 922, 525, 70, 18, 30, 4), # Hypertension cases
N = c(275, 169, 1537, 200, 200, 4725, 3073, 300, 100, 990, 42) # Total sample size
)Sorting the data by Year in ascending order
# Sorting
hypertension_data <- hypertension_data[order(hypertension_data$Year), ]
# Display the sorted data
print(hypertension_data)## Author Year Events N
## 10 Sagarad et al.,(2013) 2013 30 990
## 8 Mandal et al., (2014) 2014 70 300
## 3 Basavanagowdappa et al., (2016) 2016 248 1537
## 6 Bharatia et al., (2016) 2016 922 4725
## 11 Garg et al.,(2018) 2018 4 42
## 2 Santosh et al., (2019) 2019 67 169
## 7 Gupta et al., (2019) 2019 525 3073
## 9 Umamaheswari et al., (2019) 2019 18 100
## 1 Mahapatra et al., (2021) 2021 31 275
## 5 Annie Caroline et al., (2021) 2021 130 200
## 4 Boppana et al., (2023) 2023 100 200Step 1: Calculate observed proportions
Calculate prevalence (proportions) and their variances
ies <- escalc(measure = "PR", xi = Events, ni = N, data = hypertension_data)
View(ies)Step 2: Perform random-effects meta-analysis
Using DerSimonian-Laird (DL) method for random-effects model
pes <- rma(yi, vi, data = ies, method = "DL")
print(pes, digits = 2)##
## Random-Effects Model (k = 11; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.01 (SE = 0.01)
## tau (square root of estimated tau^2 value): 0.11
## I^2 (total heterogeneity / total variability): 98.91%
## H^2 (total variability / sampling variability): 92.00
##
## Test for Heterogeneity:
## Q(df = 10) = 920.04, p-val < .01
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.24 0.03 7.33 <.01 0.18 0.31 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1confint(pes)##
## estimate ci.lb ci.ub
## tau^2 0.0114 0.0165 0.1090
## tau 0.1070 0.1284 0.3302
## I^2(%) 98.9131 99.2422 99.8847
## H^2 92.0042 131.9623 867.6348##Step 3: Identifying outliers with residuals
stud_res <- rstudent(pes)
abs_z <- abs(stud_res$z)
outliers <- stud_res[order(-abs_z), ]
print(outliers)##
## resid se z
## 10 0.4462 0.1036 4.3058
## 1 -0.2343 0.0855 -2.7392
## 11 0.2809 0.1123 2.5018
## 6 0.1671 0.1165 1.4343
## 5 -0.1620 0.1217 -1.3313
## 9 -0.1449 0.1165 -1.2435
## 3 -0.0919 0.1231 -0.7467
## 7 -0.0818 0.1312 -0.6237
## 8 -0.0699 0.1194 -0.5854
## 4 -0.0549 0.1295 -0.4240
## 2 -0.0117 0.1157 -0.1008##Step 4: Leave-One-Out Analysis
l1o <- leave1out(pes)
# Extract estimates and confidence intervals
yi <- l1o$estimate # Leave-one-out effect sizes
sei <- l1o$se # Standard errors
ci.lb <- yi - 1.96 * sei # Lower bounds of 95% CI
ci.ub <- yi + 1.96 * sei # Upper bounds of 95% CI# Create Forest Plot for Leave-One-Out Analysis (based on prevalence)
forest(yi, sei = sei,
slab = hypertension_data$Author,
xlab = "Summary Proportions Leaving Out Each Study",
refline = pes$b, # Overall summary proportion
digits = 4, # Number of digits
alim = c(min(ci.lb), max(ci.ub)), # Adjust x-axis limits
cex = 0.8, # Text size
col = "blue", # Point color
main = "Leave-One-Out Sensitivity Analysis (Prevalence)")
# Calculate I-squared values for Leave-One-Out Analysis
isq <- sapply(1:nrow(hypertension_data), function(i) {
temp_pes <- rma(yi = ies$yi[-i], vi = ies$vi[-i], method = "DL") # Subset yi and vi correctly
return(temp_pes$I2)
})# Create Forest Plot for Leave-One-Out Analysis (based on I-squared)
forest(isq, sei = rep(0, length(isq)),
slab = hypertension_data$Author,
xlab = "I-squared Values Leaving Out Each Study",
refline = pes$I2, # Overall I-squared value
digits = 2, # Number of digits
alim = c(min(isq), max(isq)), # Adjust x-axis limits
cex = 0.8, # Text size
col = "red", # Point color
main = "Leave-One-Out Sensitivity Analysis (I-squared)")
####################Step 5: Baujat Plot (Heterogeneity Diagnostics)####################
baujat(pes, main = "Baujat Plot for Heterogeneity Diagnostics")
####################Step 6: Influence Diagnostics######################################
inf <- influence(pes)
print(inf)##
## rstudent dffits cook.d cov.r tau2.del QE.del hat weight dfbs inf
## 1 -2.7392 -0.8213 0.3868 0.6625 0.0066 335.0862 0.0965 9.6485 -0.8056
## 2 -0.1008 -0.0333 0.0011 1.1122 0.0116 904.8254 0.0919 9.1946 -0.0333
## 3 -0.7467 -0.2581 0.0791 1.3024 0.0136 913.5870 0.0960 9.5997 -0.2591
## 4 -0.4240 -0.1610 0.0343 1.4413 0.0151 786.3395 0.0965 9.6455 -0.1621
## 5 -1.3313 -0.3991 0.1604 1.0974 0.0115 919.1107 0.0820 8.2034 -0.3990
## 6 1.4343 0.4441 0.1901 1.0525 0.0110 872.8535 0.0861 8.6088 0.4446
## 7 -0.6237 -0.2266 0.0696 1.4762 0.0155 892.0304 0.0963 9.6347 -0.2282
## 8 -0.5854 -0.1805 0.0329 1.1046 0.0116 918.8868 0.0857 8.5687 -0.1804
## 9 -1.2435 -0.4041 0.1701 1.1481 0.0119 918.1161 0.0938 9.3756 -0.4044
## 10 4.3058 1.3641 1.4492 0.8441 0.0087 688.4011 0.0880 8.7992 1.3713 *
## 11 2.5018 0.7833 0.5565 0.9884 0.0103 814.8919 0.0872 8.7212 0.7852 *# Plot the Influence Diagnostics
par(mar = c(5, 5, 4, 2) + 0.1)
plot(inf, cex = 0.8) # Reduce text size to 80% of default
# Add a custom title to the plot
mtext("Influence Diagnostics for Hypertension Meta-Analysis", side = 3, line = 4, cex = 1.2, font = 2)
####################Step 7: Remove Outliers (if needed, for example 10,11)############
outlier_removed <- ies[-c(10,11), ]
pes_updated <- rma(yi, vi, data = outlier_removed, method = "DL")
print(pes_updated, digits = 2)##
## Random-Effects Model (k = 9; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.01 (SE = 0.01)
## tau (square root of estimated tau^2 value): 0.09
## I^2 (total heterogeneity / total variability): 98.62%
## H^2 (total variability / sampling variability): 72.59
##
## Test for Heterogeneity:
## Q(df = 8) = 580.71, p-val < .01
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.17 0.03 5.79 <.01 0.11 0.23 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1####################Step 8: Forest Plot for Overall Meta-Analysis##################### # Forest plot for all studies
overall_meta <- metaprop(
event = hypertension_data$Events,
n = hypertension_data$N,
studlab = hypertension_data$Author,
sm = "PR", method.ci = "NAsm", method.tau = "DL"
)Generate forest plot for overall meta-analysis
forest(
overall_meta,
leftcols = c("studlab", "event", "n", "w.random"), # Study details and weights
leftlabs = c("Study", "Cases", "Total", "Weight"),
rightcols = c("effect", "ci"), # Effect size and confidence intervals
rightlabs = c("Prevalence", "95% CI"),
comb.fixed = FALSE, # Disable common-effect (fixed-effect) model
comb.random = TRUE, # Display random-effects model only
col.square = "blue", # Color for individual study markers
col.diamond = "red", # Color for random-effect model (diamond)
col.diamond.lines = "black", # Diamond outline
type.random = "diamond", # Random-effect model as a diamond
prediction = FALSE, # Show prediction interval
print.Q = TRUE, # Display Q statistic for heterogeneity
print.I2 = TRUE, # Display I-squared heterogeneity
print.tau2 = TRUE, # Display tau-squared heterogeneity
overall = TRUE, # Include overall random-effect summary
fixed = FALSE, # Completely suppress fixed-effects model
digits = 3, # Precision for numbers
xlim = c(0, 1), # Adjust x-axis range
main = "Forest Plot: Random Effect Model Only"
)## Warning: Use argument 'common' instead of 'fixed' (deprecated).
Generate forest plot for meta-analysis after removing outliers
forest(
metaprop(
event = hypertension_data$Events[-c(10,11)], # Remove outlier events
n = hypertension_data$N[-c(10,11)], # Remove outlier sample sizes
studlab = hypertension_data$Author[-c(10,11)], # Remove outlier study labels
sm = "PR", method.ci = "NAsm", method.tau = "DL"
),
leftcols = c("studlab", "event", "n", "w.random"), # Study details and weights
leftlabs = c("Study", "Cases", "Total", "Weight"),
rightcols = c("effect", "ci"), # Effect size and confidence intervals
rightlabs = c("Prevalence", "95% CI"),
comb.fixed = FALSE, # Disable common-effect (fixed-effect) model
comb.random = TRUE, # Display random-effects model only
col.square = "blue", # Color for individual study markers
col.diamond = "red", # Color for random-effect model (diamond)
col.diamond.lines = "black", # Diamond outline
type.random = "diamond", # Random-effect model as a diamond
prediction = FALSE, # Show prediction interval
print.Q = TRUE, # Display Q statistic for heterogeneity
print.I2 = TRUE, # Display I-squared heterogeneity
print.tau2 = TRUE, # Display tau-squared heterogeneity
overall = TRUE, # Include overall random-effect summary
common = FALSE, # Completely suppress fixed-effects model
digits = 3, # Precision for numbers
xlim = c(0, 1), # Adjust x-axis range
main = "Forest Plot: Random Effect Model After Removing Outlier (Study 5)"
)
#####################Step 9: Funnel Plot for Publication Bias###########################
funnel(pes, main = "Funnel Plot for Hypertension Data")
#####################Step 10: Trim-and-Fill Analysis for Publication Bias###################
pes_trimfill <- trimfill(pes)
print(pes_trimfill)##
## Estimated number of missing studies on the left side: 0 (SE = 1.9507)
##
## Random-Effects Model (k = 11; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.0114 (SE = 0.0087)
## tau (square root of estimated tau^2 value): 0.1070
## I^2 (total heterogeneity / total variability): 98.91%
## H^2 (total variability / sampling variability): 92.00
##
## Test for Heterogeneity:
## Q(df = 10) = 920.0419, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.2439 0.0333 7.3287 <.0001 0.1787 0.3091 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1funnel(pes_trimfill, main = "Trim-and-Fill Funnel Plot")
#####################Step 11: Egger’s Regression Test for Small-Study Effects##############
egger_test <- regtest(pes, model = "lm", predictor = "sei")
print(egger_test)##
## Regression Test for Funnel Plot Asymmetry
##
## Model: weighted regression with multiplicative dispersion
## Predictor: standard error
##
## Test for Funnel Plot Asymmetry: t = 1.7516, df = 9, p = 0.1138
## Limit Estimate (as sei -> 0): b = 0.0801 (CI: -0.0173, 0.1775)#####################Step 12: Rank Correlation Test######################################
rank_corr <- ranktest(pes)
print(rank_corr)##
## Rank Correlation Test for Funnel Plot Asymmetry
##
## Kendall's tau = -0.0182, p = 1.0000####################Step 13: Perform meta-regression with sample size as a moderator#################
meta_reg_sample <- rma(yi, vi, mods = ~ N, data = ies, method = "DL")
print(meta_reg_sample, digits = 2)##
## Mixed-Effects Model (k = 11; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.01 (SE = 0.01)
## tau (square root of estimated tau^2 value): 0.12
## I^2 (residual heterogeneity / unaccounted variability): 98.79%
## H^2 (unaccounted variability / sampling variability): 82.47
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 9) = 742.21, p-val < .01
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 1.59, p-val = 0.21
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.28 0.05 6.16 <.01 0.19 0.37 ***
## N -0.00 0.00 -1.26 0.21 -0.00 0.00
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Summary of meta-regression
summary(meta_reg_sample)##
## Mixed-Effects Model (k = 11; tau^2 estimator: DL)
##
## logLik deviance AIC BIC AICc
## 1.8330 64.2517 2.3341 3.5278 5.7627
##
## tau^2 (estimated amount of residual heterogeneity): 0.0137 (SE = 0.0110)
## tau (square root of estimated tau^2 value): 0.1171
## I^2 (residual heterogeneity / unaccounted variability): 98.79%
## H^2 (unaccounted variability / sampling variability): 82.47
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 9) = 742.2057, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 1.5878, p-val = 0.2076
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.2787 0.0453 6.1558 <.0001 0.1900 0.3674 ***
## N -0.0000 0.0000 -1.2601 0.2076 -0.0001 0.0000
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1####################Step 13a: Perform meta-regression with COVID period as a moderator#################
Add a classification column for Pre- and Post-COVID years
hypertension_data$COVID_Period <- ifelse(hypertension_data$Year >= 2020, "Post-COVID", "Pre-COVID")Convert COVID_Period to a factor and set Pre-COVID as the reference level
hypertension_data$COVID_Period <- factor(hypertension_data$COVID_Period, levels = c("Pre-COVID", "Post-COVID"))Add COVID_Period to the ies dataframe
ies$COVID_Period <- hypertension_data$COVID_PeriodPerform meta-regression with COVID_Period as a moderator
meta_reg_covid <- rma(yi, vi, mods = ~ COVID_Period, data = ies, method = "DL")
print(meta_reg_covid, digits = 2)##
## Mixed-Effects Model (k = 11; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.01 (SE = 0.01)
## tau (square root of estimated tau^2 value): 0.10
## I^2 (residual heterogeneity / unaccounted variability): 98.90%
## H^2 (unaccounted variability / sampling variability): 90.59
## R^2 (amount of heterogeneity accounted for): 7.86%
##
## Test for Residual Heterogeneity:
## QE(df = 9) = 815.31, p-val < .01
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 10.35, p-val < .01
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.18 0.04 4.84 <.01 0.11 0.25 ***
## COVID_PeriodPost-COVID 0.23 0.07 3.22 <.01 0.09 0.37 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Summary of meta-regression
summary(meta_reg_covid)##
## Mixed-Effects Model (k = 11; tau^2 estimator: DL)
##
## logLik deviance AIC BIC AICc
## 4.2323 59.4531 -2.4645 -1.2709 0.9640
##
## tau^2 (estimated amount of residual heterogeneity): 0.0105 (SE = 0.0083)
## tau (square root of estimated tau^2 value): 0.1027
## I^2 (residual heterogeneity / unaccounted variability): 98.90%
## H^2 (unaccounted variability / sampling variability): 90.59
## R^2 (amount of heterogeneity accounted for): 7.86%
##
## Test for Residual Heterogeneity:
## QE(df = 9) = 815.3117, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 10.3473, p-val = 0.0013
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.1811 0.0374 4.8377 <.0001 0.1078 0.2545 ***
## COVID_PeriodPost-COVID 0.2324 0.0722 3.2167 0.0013 0.0908 0.3740 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1####################Step 13b: Perform bivariate meta-regression with COVID period and sample size#################
meta_reg_bivariate <- rma(yi, vi, mods = ~ COVID_Period + N, data = ies, method = "DL")
print(meta_reg_bivariate, digits = 2)##
## Mixed-Effects Model (k = 11; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.01 (SE = 0.01)
## tau (square root of estimated tau^2 value): 0.10
## I^2 (residual heterogeneity / unaccounted variability): 98.48%
## H^2 (unaccounted variability / sampling variability): 65.88
## R^2 (amount of heterogeneity accounted for): 6.50%
##
## Test for Residual Heterogeneity:
## QE(df = 8) = 527.03, p-val < .01
##
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 10.29, p-val < .01
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.19 0.05 3.76 <.01 0.09 0.29 ***
## COVID_PeriodPost-COVID 0.22 0.08 2.88 <.01 0.07 0.38 **
## N -0.00 0.00 -0.27 0.79 -0.00 0.00
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Summary of bivariate meta-regression
summary(meta_reg_bivariate)##
## Mixed-Effects Model (k = 11; tau^2 estimator: DL)
##
## logLik deviance AIC BIC AICc
## 4.3333 59.2510 -0.6666 0.9250 6.0001
##
## tau^2 (estimated amount of residual heterogeneity): 0.0107 (SE = 0.0092)
## tau (square root of estimated tau^2 value): 0.1035
## I^2 (residual heterogeneity / unaccounted variability): 98.48%
## H^2 (unaccounted variability / sampling variability): 65.88
## R^2 (amount of heterogeneity accounted for): 6.50%
##
## Test for Residual Heterogeneity:
## QE(df = 8) = 527.0311, p-val < .0001
##
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 10.2862, p-val = 0.0058
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.1902 0.0506 3.7608 0.0002 0.0911 0.2894 ***
## COVID_PeriodPost-COVID 0.2248 0.0781 2.8789 0.0040 0.0718 0.3779 **
## N -0.0000 0.0000 -0.2690 0.7880 -0.0001 0.0000
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1####################Step 14: Bubble Plot for Meta-Regression with Sample Size########################
Bubble plot for Sample Size with modified CI and PI lines
bubble_sample <- predict(meta_reg_sample, newmods = ies$N)
y_range_sample <- range(c(bubble_sample$ci.lb, bubble_sample$ci.ub, ies$yi))
plot(ies$N, ies$yi,
xlab = "Sample Size (N)",
ylab = "Observed Effect Size (Log Proportions)",
main = "Bubble Plot for Meta-Regression (Sample Size)",
pch = 21, bg = "lightblue", cex = 1.5, ylim = y_range_sample)
# Add meta-regression line
lines(sort(ies$N), bubble_sample$pred[order(ies$N)], col = "blue", lwd = 2)
# Add 95% confidence intervals (solid black lines)
lines(sort(ies$N), bubble_sample$ci.lb[order(ies$N)], col = "black", lty = 1, lwd = 1.5)
lines(sort(ies$N), bubble_sample$ci.ub[order(ies$N)], col = "black", lty = 1, lwd = 1.5)
# Add legend
legend("topright",
legend = c("Meta-Regression Line", "95% CI"),
col = c("blue", "black"),
lty = c(1, 1),
lwd = c(2, 1.5))
####################Step 14a: Bubble Plot for Meta-Regression with COVID Period###########
Predict values for bubble plot with COVID Period
bubble_covid <- predict(meta_reg_covid, newmods = as.numeric(ies$COVID_Period == "Post-COVID"))Create the bubble plot
plot(as.numeric(ies$COVID_Period), ies$yi,
xlab = "COVID Period (1 = Pre-COVID, 2 = Post-COVID)",
ylab = "Observed Effect Size (Log Proportions)",
main = "Bubble Plot for Meta-Regression (COVID Period)",
pch = 21, bg = "lightgreen", cex = 1.5,
ylim = range(c(bubble_covid$ci.lb, bubble_covid$ci.ub)))
# Add meta-regression line
lines(as.numeric(ies$COVID_Period), bubble_covid$pred, col = "blue", lwd = 2)
# Add 95% confidence intervals (solid black lines)
lines(as.numeric(ies$COVID_Period), bubble_covid$ci.lb, col = "black", lty = 1, lwd = 1.5)
lines(as.numeric(ies$COVID_Period), bubble_covid$ci.ub, col = "black", lty = 1, lwd = 1.5)
# Add legend
legend("topleft",
legend = c("Meta-Regression Line", "95% CI"),
col = c("blue", "black"),
lty = c(1, 1),
lwd = c(2, 1.5))
Step 15: Distribution of the True Effect
Ensure that pes\(b and pes\)tau2 are numeric and properly defined
if (!is.null(pes$b) && !is.null(pes$tau2) && pes$tau2 > 0) {
true_effects <- seq(
from = as.numeric(pes$b - 4 * sqrt(pes$tau2)),
to = as.numeric(pes$b + 4 * sqrt(pes$tau2)),
length.out = 1000
)
} else {
stop("Error: Invalid 'pes' object or missing 'tau2' value.")
}Density of the distribution
density_values <- dnorm(true_effects, mean = as.numeric(pes$b), sd = sqrt(as.numeric(pes$tau2)))Plot the distribution
plot(true_effects, density_values, type = "l", lwd = 2, col = "blue",
xlab = "True Effect Size (Log Proportions)",
ylab = "Density",
main = "Distribution of the True Effect")
# Add vertical line for the estimated overall effect size
abline(v = as.numeric(pes$b), col = "red", lwd = 2, lty = 2)
# Add legend
legend("topright", legend = c("True Effect Distribution", "Overall Effect Size"),
col = c("blue", "red"), lty = c(1, 2), lwd = 2)