# Load the readxl library to read Excel files
library(readxl)
# Specify the file path
file_path <- "C:/Users/avira/OneDrive/Documents/Research/SRMA Project/SRMA 5 (vaccine)/SRMA 5.1 analysis/DES SRMA 5.xlsx"
# Read the data from the specified Excel file
data <- read_excel(file_path)
# Print the names of the columns in the dataset
print(colnames(data))
## [1] "S.no."
## [2] "DOI"
## [3] "STUDY TITLE"
## [4] "AUTHOR"
## [5] "PUBLICATION YEAR"
## [6] "STUDY DESIGN"
## [7] "Place"
## [8] "Public/private intervention"
## [9] "SAMPLE SIZE"
## [10] "Study Population"
## [11] "Target population"
## [12] "Intervention details"
## [13] "Intervention group"
## [14] "Control group"
## [15] "Total vaccinated population (single dose)"
## [16] "Total vaccinated population (two doses)"
## [17] "Total vaccinated population (Three doses)"
## [18] "Total vaccinated population (%)"
## [19] "Delivery mode of Intervention"
## [20] "I1"
## [21] "I2"
## [22] "Overall_effect_size"
## [23] "STATISTICAL SIGNIFICANCE"
## [24] "DIRECTION OF EFFECT"
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
library(readxl)
# Load the data
file_path <- "C:/Users/avira/OneDrive/Documents/Research/SRMA Project/SRMA 5 (vaccine)/SRMA 5.1 analysis/DES SRMA 5.xlsx"
data <- read_excel(file_path)
# Visualization for Overall Effect Size
ggplot(data, aes(x = "", y = Overall_effect_size)) +
geom_boxplot(fill = "lightblue", colour = "darkblue") +
geom_violin(fill = "lightblue", alpha = 0.4) +
geom_jitter(width = 0.1, colour = "darkblue", size = 1.5, alpha = 0.9) +
labs(title = "Overall Effect Size", y = "Overall Effect Size", x = "") +
theme_minimal()
# Visualization for Infection in Intervention Group (I1)
ggplot(data, aes(x = "", y = I1)) +
geom_boxplot(fill = "lightgreen", colour = "darkgreen") +
geom_violin(fill = "lightgreen", alpha = 0.4) +
geom_jitter(width = 0.1, colour = "darkgreen", size = 1.5, alpha = 0.9) +
labs(title = "Infection in Intervention Group (I1)", y = "I1", x = "") +
theme_minimal()
# Visualization for Infection in Control Group (I2)
ggplot(data, aes(x = "", y = I2)) +
geom_boxplot(fill = "lightcoral", colour = "darkred") +
geom_violin(fill = "lightcoral", alpha = 0.4) +
geom_jitter(width = 0.1, colour = "darkred", size = 1.5, alpha = 0.9) +
labs(title = "Infection in Control Group (I2)", y = "I2", x = "") +
theme_minimal()
# Visualization for Sample Size
ggplot(data, aes(x = "", y = `SAMPLE SIZE`)) +
geom_boxplot(fill = "lightyellow", colour = "goldenrod") +
geom_violin(fill = "lightyellow", alpha = 0.4) +
geom_jitter(width = 0.1, colour = "goldenrod", size = 1.5, alpha = 0.9) +
labs(title = "Sample Size Distribution", y = "Sample Size", x = "") +
theme_minimal()
# Print all plots sequentially without saving
library(dplyr)
library(knitr)
library(kableExtra)
##
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
##
## group_rows
library(readxl)
library(tidyr)
# Load the data
file_path <- "C:/Users/avira/OneDrive/Documents/Research/SRMA Project/SRMA 5 (vaccine)/SRMA 5.1 analysis/DES SRMA 5.xlsx"
data <- read_excel(file_path)
# Define numerical columns
numerical_cols <- c("SAMPLE SIZE", "Total vaccinated population (single dose)",
"Total vaccinated population (two doses)", "Total vaccinated population (Three doses)",
"Total vaccinated population (%)", "I1", "I2", "Overall_effect_size")
# Select only numerical data
numerical_data <- data %>%
select(all_of(numerical_cols))
# Compute descriptive statistics for the selected numerical data
descr_stats <- numerical_data %>%
summarise(across(everything(), list(
mean = ~mean(., na.rm = TRUE),
sd = ~sd(., na.rm = TRUE),
min = ~min(., na.rm = TRUE),
median = ~median(., na.rm = TRUE),
max = ~max(., na.rm = TRUE),
iqr = ~IQR(., na.rm = TRUE)
), .names = "{col}_{fn}"))
# Convert the wide format summary to a long format for better display
descr_stats_long <- pivot_longer(descr_stats, cols = everything(), names_to = "Measure", values_to = "Value")
# Format the values to two decimal places and convert to character to prevent scientific notation
descr_stats_long$Value <- sprintf("%.2f", as.numeric(descr_stats_long$Value))
# Create the table using kable
table <- descr_stats_long %>%
kable("html", escape = FALSE, align = 'lcc', digits = 2, col.names = c("Measure", "Value")) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) %>%
add_header_above(c("Descriptive Summary Statistics" = 2))
# Print the table in the R console
print(table)
## <table class="table table-striped table-hover table-condensed" style="color: black; width: auto !important; margin-left: auto; margin-right: auto;">
## <thead>
## <tr><th style="border-bottom:hidden;padding-bottom:0; padding-left:3px;padding-right:3px;text-align: center; " colspan="2"><div style="border-bottom: 1px solid #ddd; padding-bottom: 5px; ">Descriptive Summary Statistics</div></th></tr>
## <tr>
## <th style="text-align:left;"> Measure </th>
## <th style="text-align:center;"> Value </th>
## </tr>
## </thead>
## <tbody>
## <tr>
## <td style="text-align:left;"> SAMPLE SIZE_mean </td>
## <td style="text-align:center;"> 16804.78 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> SAMPLE SIZE_sd </td>
## <td style="text-align:center;"> 38738.27 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> SAMPLE SIZE_min </td>
## <td style="text-align:center;"> 27.00 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> SAMPLE SIZE_median </td>
## <td style="text-align:center;"> 5792.50 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> SAMPLE SIZE_max </td>
## <td style="text-align:center;"> 371108.00 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> SAMPLE SIZE_iqr </td>
## <td style="text-align:center;"> 18201.75 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (single dose)_mean </td>
## <td style="text-align:center;"> 16395.58 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (single dose)_sd </td>
## <td style="text-align:center;"> 39161.86 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (single dose)_min </td>
## <td style="text-align:center;"> 46.00 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (single dose)_median </td>
## <td style="text-align:center;"> 4597.50 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (single dose)_max </td>
## <td style="text-align:center;"> 241288.00 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (single dose)_iqr </td>
## <td style="text-align:center;"> 15749.25 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (two doses)_mean </td>
## <td style="text-align:center;"> 12181.71 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (two doses)_sd </td>
## <td style="text-align:center;"> 19813.99 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (two doses)_min </td>
## <td style="text-align:center;"> 0.00 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (two doses)_median </td>
## <td style="text-align:center;"> 6214.00 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (two doses)_max </td>
## <td style="text-align:center;"> 129720.00 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (two doses)_iqr </td>
## <td style="text-align:center;"> 13816.75 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (Three doses)_mean </td>
## <td style="text-align:center;"> 7202.00 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (Three doses)_sd </td>
## <td style="text-align:center;"> 11357.86 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (Three doses)_min </td>
## <td style="text-align:center;"> 39.00 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (Three doses)_median </td>
## <td style="text-align:center;"> 3508.00 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (Three doses)_max </td>
## <td style="text-align:center;"> 56302.00 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (Three doses)_iqr </td>
## <td style="text-align:center;"> 5546.50 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (%)_mean </td>
## <td style="text-align:center;"> 37797.46 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (%)_sd </td>
## <td style="text-align:center;"> 310933.71 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (%)_min </td>
## <td style="text-align:center;"> 22.00 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (%)_median </td>
## <td style="text-align:center;"> 4370.00 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (%)_max </td>
## <td style="text-align:center;"> 3711008.00 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Total vaccinated population (%)_iqr </td>
## <td style="text-align:center;"> 16474.75 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> I1_mean </td>
## <td style="text-align:center;"> 0.28 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> I1_sd </td>
## <td style="text-align:center;"> 0.64 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> I1_min </td>
## <td style="text-align:center;"> 0.00 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> I1_median </td>
## <td style="text-align:center;"> 0.04 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> I1_max </td>
## <td style="text-align:center;"> 3.96 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> I1_iqr </td>
## <td style="text-align:center;"> 0.28 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> I2_mean </td>
## <td style="text-align:center;"> 1.28 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> I2_sd </td>
## <td style="text-align:center;"> 4.86 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> I2_min </td>
## <td style="text-align:center;"> 0.00 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> I2_median </td>
## <td style="text-align:center;"> 0.11 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> I2_max </td>
## <td style="text-align:center;"> 51.10 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> I2_iqr </td>
## <td style="text-align:center;"> 0.61 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Overall_effect_size_mean </td>
## <td style="text-align:center;"> 0.66 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Overall_effect_size_sd </td>
## <td style="text-align:center;"> 0.23 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Overall_effect_size_min </td>
## <td style="text-align:center;"> 0.14 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Overall_effect_size_median </td>
## <td style="text-align:center;"> 0.67 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Overall_effect_size_max </td>
## <td style="text-align:center;"> 1.00 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Overall_effect_size_iqr </td>
## <td style="text-align:center;"> 0.40 </td>
## </tr>
## </tbody>
## </table>
# Load required libraries
library(metafor)
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
## Loading required package: metadat
## Loading required package: numDeriv
##
## Loading the 'metafor' package (version 4.6-0). For an
## introduction to the package please type: help(metafor)
library(readxl)
library(dplyr)
# Set the file path
file_path <- "C:/Users/avira/OneDrive/Documents/Research/SRMA Project/SRMA 5 (vaccine)/SRMA 5.1 analysis/DES SRMA 5.xlsx"
# Read the data from Excel file
data <- read_excel(file_path)
# Prepare the data for meta-analysis
# Assuming 'I1' and 'I2' are the incidence in the intervention and control groups, respectively
# and 'Overall_effect_size' is the log of the effect size
meta_data <- data %>%
select(I1, I2, Overall_effect_size) %>%
mutate(
log_effect_size = as.numeric(Overall_effect_size), # Ensure it's numeric
variance = 1 / (I1 + 1) + 1 / (I2 + 1) # Approximate variance of log effect size
)
# Run a random-effects meta-analysis
meta_results <- rma(yi = log_effect_size, vi = variance, data = meta_data, method = "REML")
# Print the results of the meta-analysis
print(summary(meta_results))
##
## Random-Effects Model (k = 142; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -162.2725 324.5450 328.5450 334.4425 328.6319
##
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.1342)
## 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 = 141) = 5.2313, p-val = 1.0000
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.6654 0.0982 6.7752 <.0001 0.4729 0.8579 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Generate a forest plot
forest(meta_results, main="Forest Plot of Vaccine Efficacy")
# Load required libraries
library(metafor)
library(readxl)
# Set the file path
file_path <- "C:/Users/avira/OneDrive/Documents/Research/SRMA Project/SRMA 5 (vaccine)/SRMA 5.1 analysis/DES SRMA 5.xlsx"
# Read the data from Excel file
data <- read_excel(file_path)
# Prepare the data for meta-analysis
# Assuming 'Overall_effect_size' is already the log of the effect size and the dataset includes 'variance'
meta_data <- data %>%
select(I1, I2, Overall_effect_size) %>%
mutate(
log_effect_size = as.numeric(Overall_effect_size),
variance = 1 / (I1 + 1) + 1 / (I2 + 1) # Approximate variance of log effect size
)
# Run a random-effects meta-analysis
meta_results <- rma(yi = log_effect_size, vi = variance, data = meta_data, method = "REML")
# Set up plotting area with smaller margins to utilize space better
par(mar = c(0.1, 4, 1, 2)) # You may need to adjust these values
# Create the forest plot with smaller text and point sizes
forest(meta_results,
cex = 0.5, # Adjust text size
pch = 19, # Type of point to use for the mean estimates
cex.lab = 0.5, # Adjust label text size
cex.axis = 0.5, # Adjust axis text size
mgp = c(1, 0.5, 0), # Adjust axis title, label and line position
mar = c(0.1, 0.1, 0.1, 0.1)) # Adjust margins of the plot
# Restore default plotting area settings
par(mar = c(5, 4, 4, 2) + 0.1)
# Assuming the metafor package and your data have been loaded and prepared as per previous examples.
# Run the meta-analysis if not already done
meta_results <- rma(yi = log_effect_size, vi = variance, data = meta_data, method = "REML")
# Create a Galbraith plot
galbraith(meta_results, main = "Galbraith Plot to Investigate Heterogeneity")
# Load required libraries
library(metafor)
library(readxl)
library(dplyr)
# Set the file path
file_path <- "C:/Users/avira/OneDrive/Documents/Research/SRMA Project/SRMA 5 (vaccine)/SRMA 5.1 analysis/DES SRMA 5.xlsx"
# Read the data from Excel file
data <- read_excel(file_path)
# Prepare the data for meta-analysis
# Assuming 'I1' and 'I2' are the incidences in the intervention and control groups, respectively
# and 'Overall_effect_size' is the log of the effect size
meta_data <- data %>%
select(`I1`, `I2`, `Overall_effect_size`) %>%
mutate(
log_effect_size = as.numeric(`Overall_effect_size`),
variance = 1 / (`I1` + 1) + 1 / (`I2` + 1) # Approximate variance of log effect size
)
# Run a random-effects meta-analysis
meta_results <- rma(yi = log_effect_size, vi = variance, data = meta_data, method = "REML")
# Create the funnel plot
funnel(meta_results, xlab = "Log Effect Size", ylab = "Standard Error",
main = "Funnel Plot of Standard Error by Log Effect Size")
# Perform Egger's test of funnel plot asymmetry
eggers_test <- regtest(meta_results, model = "lm")
print(eggers_test)
##
## Regression Test for Funnel Plot Asymmetry
##
## Model: weighted regression with multiplicative dispersion
## Predictor: standard error
##
## Test for Funnel Plot Asymmetry: t = -1.2046, df = 140, p = 0.2304
## Limit Estimate (as sei -> 0): b = 0.7530 (CI: 0.6044, 0.9017)
# Load necessary libraries
library(metafor)
library(readxl)
library(dplyr)
library(ggplot2)
# Set the file path to your Excel file
file_path <- "C:/Users/avira/OneDrive/Documents/Research/SRMA Project/SRMA 5 (vaccine)/SRMA 5.1 analysis/DES SRMA 5.xlsx"
# Read the data from Excel
data <- read_excel(file_path)
# Prepare the data for plotting
meta_data <- data %>%
select(`I1`, `I2`, `Overall_effect_size`) %>%
mutate(
log_effect_size = as.numeric(`Overall_effect_size`),
variance = 1 / (`I1` + 1) + 1 / (`I2` + 1), # Calculate variance as inverse of the sum of inverse incidences
weight = 1 / variance # Weight is the inverse of variance
)
# Run a random-effects meta-analysis to get an overall effect size
meta_results <- rma(yi = log_effect_size, vi = variance, data = meta_data, method = "REML")
# Print the summary of the meta-analysis results
meta_summary <- summary(meta_results)
print(meta_summary)
##
## Random-Effects Model (k = 142; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -162.2725 324.5450 328.5450 334.4425 328.6319
##
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.1342)
## 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 = 141) = 5.2313, p-val = 1.0000
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.6654 0.0982 6.7752 <.0001 0.4729 0.8579 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Create a Rainforest plot using ggplot2
rainforest_plot <- ggplot(meta_data, aes(x = log_effect_size, y = 1, size = weight, fill = weight)) +
geom_point(shape = 21, color = "black", alpha = 0.6) +
scale_fill_gradient(low = "skyblue", high = "blue4") +
scale_size(range = c(2, 12), guide = "none") + # Corrected line to remove size legend
geom_vline(xintercept = meta_results$b, linetype = "dashed", color = "firebrick", size = 1) +
labs(x = "Log Effect Size", y = "", fill = "Study Weight") +
theme_minimal() +
theme(
legend.position = "right",
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
plot.title = element_text(hjust = 0.5)
)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# Print the Rainforest plot to the console
print(rainforest_plot)
# Load required libraries
library(metafor)
library(readxl)
library(dplyr)
library(grid) # Ensure grid is loaded for grid.newpage()
# Set the file path
file_path <- "C:/Users/avira/OneDrive/Documents/Research/SRMA Project/SRMA 5 (vaccine)/SRMA 5.1 analysis/DES SRMA 5.xlsx"
# Read the data from Excel file
data <- read_excel(file_path)
# Prepare the data for meta-analysis
meta_data <- data %>%
select(`PUBLICATION YEAR`, `AUTHOR`, `I1`, `I2`, `Overall_effect_size`) %>%
mutate(
log_effect_size = as.numeric(`Overall_effect_size`),
variance = 1 / (`I1` + 1) + 1 / (`I2` + 1),
pub_year = as.numeric(`PUBLICATION YEAR`),
label = paste(`AUTHOR`, `PUBLICATION YEAR`, sep = ", ")
) %>%
arrange(pub_year)
# Check the range of publication years to determine the number of plots
years <- sort(unique(meta_data$pub_year))
# Create a directory for output if it doesn't exist
output_dir <- dirname(file_path)
if (!dir.exists(output_dir)) {
dir.create(output_dir, recursive = TRUE)
}
# Initiate a PDF device to save multiple pages
pdf_filename <- file.path(output_dir, "multi_page_forest_plots.pdf")
pdf(pdf_filename, height=11, width=8.5)
# Create a forest plot for each year
for (year in years) {
# Subset data for the year
yearly_data <- subset(meta_data, pub_year == year)
# Skip if there is no data for the year
if (nrow(yearly_data) == 0) next
# Run the meta-analysis for the subset
meta_results <- rma(yi = log_effect_size, vi = variance, data = yearly_data, method = "REML")
# Print the summary to the console
print(paste("Year:", year))
print(summary(meta_results))
# Generate the forest plot for the current subset with labels
forest(meta_results,
slab = yearly_data$label,
main = paste("Forest Plot for Year", year),
xlab = "Log Effect Size",
alim = c(-3, 3),
at = seq(-3, 3, 0.5)) # Adjust the axis tick marks if needed
# Add a new page in the PDF if not the last plot
if (year != tail(years, n = 1)) {
grid.newpage()
}
}
## [1] "Year: 2015"
##
## Random-Effects Model (k = 13; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -12.9651 25.9303 29.9303 30.9001 31.2636
##
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.3305)
## 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.7517, p-val = 1.0000
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.7394 0.2723 2.7153 0.0066 0.2057 1.2731 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Year: 2016"
##
## Random-Effects Model (k = 11; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -12.4183 24.8367 28.8367 29.4418 30.5510
##
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.8083)
## 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 = 10) = 0.4112, p-val = 1.0000
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.6905 0.4064 1.6992 0.0893 -0.1060 1.4870 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Year: 2017"
##
## Random-Effects Model (k = 16; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -17.7539 35.5078 39.5078 40.9239 40.5078
##
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.5405)
## 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 = 15) = 0.8268, p-val = 1.0000
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.5080 0.3101 1.6384 0.1013 -0.0997 1.1157
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Year: 2018"
##
## Random-Effects Model (k = 12; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -12.1201 24.2402 28.2402 29.0360 29.7402
##
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.4839)
## 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 = 11) = 0.5848, p-val = 1.0000
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.6580 0.3195 2.0594 0.0395 0.0318 1.2842 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Year: 2019"
##
## Random-Effects Model (k = 10; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -10.4715 20.9429 24.9429 25.3374 26.9429
##
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.7174)
## 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 = 9) = 0.1776, p-val = 1.0000
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.6738 0.3938 1.7111 0.0871 -0.0980 1.4455 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Year: 2020"
##
## Random-Effects Model (k = 10; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -10.7359 21.4717 25.4717 25.8662 27.4717
##
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.7204)
## 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 = 9) = 0.3169, p-val = 1.0000
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.6734 0.3976 1.6937 0.0903 -0.1059 1.4527 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Year: 2021"
##
## Random-Effects Model (k = 31; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -33.6903 67.3807 71.3807 74.1831 71.8251
##
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.2720)
## 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 = 30) = 0.9180, p-val = 1.0000
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.7197 0.2021 3.5618 0.0004 0.3237 1.1157 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Year: 2022"
##
## Random-Effects Model (k = 28; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -31.4647 62.9294 66.9294 69.5211 67.4294
##
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.3414)
## 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 = 27) = 0.5929, p-val = 1.0000
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.6413 0.2272 2.8230 0.0048 0.1960 1.0865 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Year: 2023"
##
## Random-Effects Model (k = 11; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -11.6841 23.3682 27.3682 27.9734 29.0825
##
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.6461)
## 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 = 10) = 0.2088, p-val = 1.0000
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.6098 0.3705 1.6457 0.0998 -0.1164 1.3360 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Turn off the PDF device
dev.off()
## png
## 2
# Output the location of the saved PDF
cat("The forest plots have been saved to:", pdf_filename, "\n")
## The forest plots have been saved to: C:/Users/avira/OneDrive/Documents/Research/SRMA Project/SRMA 5 (vaccine)/SRMA 5.1 analysis/multi_page_forest_plots.pdf
# Install and load necessary packages
if (!require("sf")) install.packages("sf")
## Loading required package: sf
## Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE
if (!require("ggplot2")) install.packages("ggplot2")
if (!require("rnaturalearth")) install.packages("rnaturalearth")
## Loading required package: rnaturalearth
if (!require("stringdist")) install.packages("stringdist")
## Loading required package: stringdist
##
## Attaching package: 'stringdist'
## The following object is masked from 'package:tidyr':
##
## extract
if (!require("readxl")) install.packages("readxl")
library(sf)
library(ggplot2)
library(rnaturalearth)
library(stringdist)
library(readxl)
# Set file path to read the Excel file
file_path <- "C:/Users/avira/OneDrive/Documents/Research/SRMA Project/SRMA 5 (vaccine)/SRMA 5.1 analysis/DES SRMA 5.xlsx"
data <- readxl::read_excel(file_path)
# Load world countries
world <- ne_countries(scale = "medium", returnclass = "sf")
# Create a manual mapping based on approximate matching
data$matched_country <- sapply(data$Place, function(x) {
names <- world$name
index <- which.min(stringdist(tolower(names), tolower(x)))
names[index]
})
# Merge data with world map based on the matched country names
world_data <- merge(world, data, by.x = "name", by.y = "matched_country", all.x = TRUE)
# Plot using ggplot2 with color gradient based on 'Overall_effect_size'
ggplot(data = world_data) +
geom_sf(aes(fill = Overall_effect_size), color = "white") + # Ensure 'Overall_effect_size' is your intended variable
scale_fill_viridis_c(option = "C", na.value = "grey", guide = "colourbar") +
labs(title = "World Map with Overall Effect Size Variations", fill = "Overall Effect Size") +
theme_minimal() +
theme(legend.position = "bottom")
