##### META-ANALYSIS FOR LITERATURE REVIEW #####
## Load required packages
library(metafor) # Package for meta-analysis
library(readxl) # Package for reading Excel files
# Fix font issue
windowsFonts("Arial" = windowsFont("Arial"))
## Read in excel file with soil+nutrients data
SoilNutrients <- read_xlsx("Grant et al 2022. Literature Review and Meta-Analysis.xlsx", sheet = "SOIL")
############################# NITRATE #################################
# Response Ratio for Nitrate/NO3
RR_effect_sizes_No3 <- escalc("ROM",
m1i = Mean_seabird_NO3,
n1i = Seabird_sample_NO3,
sd1i = SD_seabird_NO3,
m2i = Mean_control_NO3,
n2i = Control_sample_NO3,
sd2i = SD_control_NO3,
data = SoilNutrients)
# Create unique ID number for each paper
RR_effect_sizes_No3$ID <- seq.int(nrow(RR_effect_sizes_No3))
##### Run the random effects model
REM_NO3 <- rma.mv(yi, vi,
random = ~ 1 | ID,
slab = paste(Author, Year, sep = ", "),
data = RR_effect_sizes_No3)
# View results
REM_NO3
# Create forest plot of nitrate for all studies
tiff(file = 'Nitrate_Plot.tiff', width = 25, height = 29, units = "cm", res = 500)
forest_NO3 <- forest(RR_effect_sizes_No3$yi, RR_effect_sizes_No3$vi,
annotate = TRUE, slab = REM_NO3$slab,
xlab = "ln(Repsonse Ratio) - Nitrate",
cex = 1,
cex.lab = 1.25,
ylim = c(-1, 28),
xlim = c(-7, 10),
alim = c(-3, 7),
at = c(-3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7),
)
addpoly(REM_NO3, row= -1, cex=1.25, mlab="Summary", fonts = "Arial")
abline(h=0)
text(-5.1, 27, "1st Author and Year", cex = 1.25, font = 2)
text(8.7, 27, "LRR [95% CI]", cex = 1.25, font = 2)
dev.off()
#### Run a meta-regression with island/mainland as a predictor
mix_effect_island_NO3 <-
rma.mv(
yi, # this is the column from our dataframe that contains calculated effect sizes
vi, # calculated variances
mods = ~ Island_Mainland - 1, # Specify which moderators the model should take into account
random = ~ 1 | ID,
method = "REML", #Model fitting method
digits = 4, # round calculations to 4 digits
data = RR_effect_sizes_No3)
# View results
mix_effect_island_NO3
# Create vector with island and mainland sample sizes
IslandMainland_N_NO3 <- c(23, 2)
# Create vector with p-values for island and mainland
IslandMainland_P_NO3 <- c("<0.0001", 0.0696)
# Forest plot for island v mainland
tiff(file = 'Nitrate_Island_Plot.tiff', width = 20, height = 8, units = "cm", res = 500)
forest_IM <- forest(mix_effect_island_NO3$b,
ci.lb = mix_effect_island_NO3$ci.lb,
ci.ub = mix_effect_island_NO3$ci.ub,
annotate = TRUE,
xlab = "ln(Response Ratio) - Nitrate",
slab = c("Island", "Mainland"),
xlim = c(-7, 8),
cex = 1,
digits = c(2L, 2L),
ylim = c(-1, 5),
top = 1)
text(-3.5, rev(seq(2:1)), IslandMainland_N_NO3, cex = 1)
text(-2, rev(seq(2:1)), IslandMainland_P_NO3, cex = 1)
text(-2, 4, "P-value", cex = 1.25, font = 2)
text(-3.5, 4, "n", cex = 1.25, font = 2)
text(-5.8, 4, "Sub Group", cex = 1.25, font = 2)
text(6.5, 4, "LRR [95% CI]", cex = 1.25, font = 2)
abline(h=0)
abline(h=3)
dev.off()
#### Run a meta-regression with seabird order as a predictor
mix_effect_order_NO3 <-
rma.mv(
yi, # this is the column from our dataframe that contains calculated effect sizes
vi, # calculated variances
mods = ~ Order - 1, # Specify which moderators the model should take into account
random = ~ 1 | ID,
method = "REML", #Model fitting method
digits = 4, # round calculations to 4 digits
data = RR_effect_sizes_No3)
# View results
mix_effect_order_NO3
# Create vector with order sample sizes
Order_N_NO3 <- c(4, 2, 9, 6)
# Create vector with p-values for order
Order_P_NO3 <- c(0.0544, 0.0981, "<0.0001", "<0.0001")
# Forest plot for order
tiff(file = 'Nitrate_Order_Plot.tiff', width = 20, height = 8.5, units = "cm", res = 500)
forest(mix_effect_order_NO3$b,
ci.lb = mix_effect_order_NO3$ci.lb,
ci.ub = mix_effect_order_NO3$ci.ub,
annotate = TRUE,
xlab = "ln(Response Ratio) - Nitrate",
slab = c("Charadriiformes", "Procellariiformes", "Sphenisciformes", "Suliformes"),
cex = 1,
alim = c(-1, 5),
at = c(-1, 0, 1, 2, 3, 4, 5),
xlim = c(-7, 8),
ylim = c(-1.5, 5),
top = 1)
text(-3.5, rev(seq(4:1)), Order_N_NO3, cex = 1)
text(-2, rev(seq(4:1)), Order_P_NO3, cex = 1)
addpoly(REM_NO3, row= -1, cex=1, mlab="Full Model")
text(-3.5, -1, 25, cex = 1)
text(-2, -1, "<0.0001", cex = 1)
abline(h=0)
dev.off()
#### Run a meta-regression with Zone as a predictor
mix_effect_zone_NO3 <-
rma.mv(
yi, # this is the column from our dataframe that contains calculated effect sizes
vi, # calculated variances
mods = ~ Zone - 1, # Specify which moderators the model should take into account
random = ~ 1 | ID,
method = "REML", #Model fitting method
digits = 4, # round calculations to 4 digits
data = RR_effect_sizes_No3)
# View results
mix_effect_zone_NO3
# Create vector with zone sample sizes
Zone_N_NO3 <- c(9, 14, 2)
# Create vector with p-values for zone
Zone_P_NO3 <- c(0.0043, "<0.0001", 0.0113)
# Create forest plot of zone
tiff(file = 'Nitrate_Zone_Plot.tiff', width = 20, height = 7.5, units = "cm", res = 500)
forest_Z_NO3 <- forest(mix_effect_zone_NO3$b,
ci.lb = mix_effect_zone_NO3$ci.lb,
ci.ub = mix_effect_zone_NO3$ci.ub,
annotate = TRUE,
xlab = "ln(Response Ratio) - Nitrate",
slab = c("Polar", "Temperate", "Tropical"),
cex = 1,
alim = c(-1, 5),
at = c(-1, 0, 1, 2, 3, 4, 5),
xlim = c(-7, 8),
ylim = c(-1.5, 4),
top = 1)
text(-3.5, rev(seq(3:1)), Zone_N_NO3, cex = 1)
text(-2, rev(seq(3:1)), Zone_P_NO3, cex = 1)
abline(h=0)
dev.off()
################################ AMMONIUM ######################################
# Response Ratio for Ammonium/NH4
RR_effect_sizes_NH4 <- escalc("ROM",
m1i = Mean_seabird_NH4,
n1i = Seabird_sample_NH4,
sd1i = SD_seabird_NH4,
m2i = Mean_control_NH4,
n2i = Control_sample_NH4,
sd2i = SD_control_NH4,
data = SoilNutrients)
# Create unique ID number for each paper
RR_effect_sizes_NH4$ID <- seq.int(nrow(RR_effect_sizes_NH4))
##### Run the random effects model
REM_NH4 <- rma.mv(yi, vi,
random = ~ 1 | ID,
slab = paste(Author, Year, sep = ", "),
data = RR_effect_sizes_NH4)
# View results
REM_NH4
# Create forest plot of nitrate for all studies
tiff(file = 'Ammonium_Plot.tiff', width = 25, height = 29, units = "cm", res = 500)
forest_NH4 <- forest(RR_effect_sizes_NH4$yi, RR_effect_sizes_NH4$vi,
annotate = TRUE, slab = REM_NH4$slab,
xlab = "ln(Repsonse Ratio) - Ammonium",
cex = 1,
cex.lab = 1.25,
ylim = c(-1, 28),
xlim = c(-7, 10),
alim = c(-3, 7),
at = c(-3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7),
)
addpoly(REM_NH4, row= -1, cex=1.25, mlab="Summary", fonts = "Arial")
abline(h=0)
text(-5.1, 27, "1st Author and Year", cex = 1.25, font = 2)
text(8.7, 27, "LRR [95% CI]", cex = 1.25, font = 2)
dev.off()
#### Run a meta-regression with island/mainland as a predictor
mix_effect_island_NH4 <-
rma.mv(
yi, # this is the column from our dataframe that contains calculated effect sizes
vi, # calculated variances
mods = ~ Island_Mainland - 1, # Specify which moderators the model should take into account
random = ~ 1 | ID,
method = "REML", #Model fitting method
digits = 4, # round calculations to 4 digits
data = RR_effect_sizes_NH4)
# View results
mix_effect_island_NH4
# Create vector with island and mainland sample sizes
IslandMainland_N_NH4 <- c(23, 2)
# Create vector with p-values for island and mainland
IslandMainland_P_NH4 <- c("<0.0001", 0.0257)
# Forest plot for island v mainland
tiff(file = 'Ammonium_Island_Plot.tiff', width = 20, height = 8, units = "cm", res = 500)
forest_IM_NH4 <- forest(mix_effect_island_NH4$b,
ci.lb = mix_effect_island_NH4$ci.lb,
ci.ub = mix_effect_island_NH4$ci.ub,
annotate = TRUE,
xlab = "ln(Response Ratio) - Ammonium",
slab = c("Island", "Mainland"),
xlim = c(-7, 8),
alim = c(-1, 6),
cex = 1,
digits = c(2L, 2L),
ylim = c(-1, 5),
top = 1)
text(-3.5, rev(seq(2:1)), IslandMainland_N_NH4, cex = 1)
text(-2, rev(seq(2:1)), IslandMainland_P_NH4, cex = 1)
text(-2, 4, "P-value", cex = 1.25, font = 2)
text(-3.5, 4, "n", cex = 1.25, font = 2)
text(-5.8, 4, "Sub Group", cex = 1.25, font = 2)
text(6.5, 4, "LRR [95% CI]", cex = 1.25, font = 2)
abline(h=0)
abline(h=3)
dev.off()
#### Run a meta-regression with seabird order as a predictor
mix_effect_order_NH4 <-
rma.mv(
yi, # this is the column from our dataframe that contains calculated effect sizes
vi, # calculated variances
mods = ~ Order - 1, # Specify which moderators the model should take into account
random = ~ 1 | ID,
method = "REML", #Model fitting method
digits = 4, # round calculations to 4 digits
data = RR_effect_sizes_NH4)
# View results
mix_effect_order_NH4
# Create vector with order sample sizes
Order_N_NH4 <- c(4, 3, 9, 5)
# Create vector with p-values for order
Order_P_NH4 <- c(0.1137, "<0.0001", "<0.0001", "0.0002")
# Forest plot for order
tiff(file = 'Ammonium_Order_Plot.tiff', width = 20, height = 8.5, units = "cm", res = 500)
forest(mix_effect_order_NH4$b,
ci.lb = mix_effect_order_NH4$ci.lb,
ci.ub = mix_effect_order_NH4$ci.ub,
annotate = TRUE,
xlab = "ln(Response Ratio) - Ammonium",
slab = c("Charadriiformes", "Procellariiformes", "Sphenisciformes", "Suliformes"),
cex = 1,
alim = c(-1, 6),
at = c(-1, 0, 1, 2, 3, 4, 5, 6),
xlim = c(-7, 8),
ylim = c(-1.5, 5),
top = 1)
text(-3.5, rev(seq(4:1)), Order_N_NH4, cex = 1)
text(-2, rev(seq(4:1)), Order_P_NH4, cex = 1)
addpoly(REM_NH4, row= -1, cex=1, mlab="Full Model")
text(-3.5, -1, 25, cex = 1)
text(-2, -1, "<0.0001", cex = 1)
abline(h=0)
dev.off()
#### Run a meta-regression with Zone as a predictor
mix_effect_zone_NH4 <-
rma.mv(
yi, # this is the column from our dataframe that contains calculated effect sizes
vi, # calculated variances
mods = ~ Zone - 1, # Specify which moderators the model should take into account
random = ~ 1 | ID,
method = "REML", #Model fitting method
digits = 4, # round calculations to 4 digits
data = RR_effect_sizes_NH4)
# View results
mix_effect_zone_NH4
# Create vector with zone sample sizes
Zone_N_NH4 <- c(9, 14, 2)
# Create vector with p-values for zone
Zone_P_NH4 <- c("<0.0001", "<0.0001", 0.0646)
# Create forest plot of zone
tiff(file = 'Ammmonium_Zone_Plot.tiff', width = 20, height = 7.5, units = "cm", res = 500)
forest_Z_NH4 <- forest(mix_effect_zone_NH4$b,
ci.lb = mix_effect_zone_NH4$ci.lb,
ci.ub = mix_effect_zone_NH4$ci.ub,
annotate = TRUE,
xlab = "ln(Response Ratio) - Ammonium",
slab = c("Polar", "Temperate", "Tropical"),
cex = 1,
alim = c(-1, 6),
at = c(-1, 0, 1, 2, 3, 4, 5, 6),
xlim = c(-7, 8),
ylim = c(-1.5, 4),
top = 1)
text(-3.5, rev(seq(3:1)), Zone_N_NH4, cex = 1)
text(-2, rev(seq(3:1)), Zone_P_NH4, cex = 1)
abline(h=0)
dev.off()
################################ NITROGEN ######################################
# Response Ratio for Nitrogen/N
RR_effect_sizes_N <- escalc("ROM",
m1i = Mean_seabird_N,
n1i = Seabird_sample_N,
sd1i = SD_seabird_N,
m2i = Mean_control_N,
n2i = Control_sample_N,
sd2i = SD_control_N,
data = SoilNutrients)
# Create unique ID number for each paper
RR_effect_sizes_N$ID <- seq.int(nrow(RR_effect_sizes_N))
##### Run the random effects model
REM_N <- rma.mv(yi, vi,
random = ~ 1 | ID,
slab = paste(Author, Year, sep = ", "),
data = RR_effect_sizes_N)
# View results
REM_N
# Create forest plot of nitrogen for all studies
tiff(file = 'Nitrogen_Plot.tiff', width = 25, height = 29, units = "cm", res = 500)
forest_N <- forest(RR_effect_sizes_N$yi, RR_effect_sizes_N$vi,
annotate = TRUE, slab = REM_N$slab,
xlab = "ln(Repsonse Ratio) - Nitrogen",
cex = 1,
cex.lab = 1.25,
ylim = c(-1, 28),
xlim = c(-7, 10),
alim = c(-3, 7),
at = c(-3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7),
)
addpoly(REM_N, row= -1, cex=1.25, mlab="Summary", fonts = "Arial")
abline(h=0)
text(-5.1, 27, "1st Author and Year", cex = 1.25, font = 2)
text(8.7, 27, "LRR [95% CI]", cex = 1.25, font = 2)
dev.off()
#### Run a meta-regression with island/mainland as a predictor
mix_effect_island_N <-
rma.mv(
yi, # this is the column from our dataframe that contains calculated effect sizes
vi, # calculated variances
mods = ~ Island_Mainland - 1, # Specify which moderators the model should take into account
random = ~ 1 | ID,
method = "REML", #Model fitting method
digits = 4, # round calculations to 4 digits
data = RR_effect_sizes_N)
# View results
mix_effect_island_N
# Create vector with island and mainland sample sizes
IslandMainland_N_N <- c(18, 7)
# Create vector with p-values for island and mainland
IslandMainland_P_N <- c("<0.0001", 0.2972)
# Forest plot for island v mainland
tiff(file = 'Nitrogen_Island_Plot.tiff', width = 20, height = 8, units = "cm", res = 500)
forest_IM_N <- forest(mix_effect_island_N$b,
ci.lb = mix_effect_island_N$ci.lb,
ci.ub = mix_effect_island_N$ci.ub,
annotate = TRUE,
xlab = "ln(Response Ratio) - Nitrogen",
slab = c("Island", "Mainland"),
xlim = c(-7, 8),
alim = c(-1, 5),
at = c(-1, 0, 1, 2, 3, 4, 5),
cex = 1,
digits = c(2L, 2L),
ylim = c(-1, 5),
top = 1)
text(-3.5, rev(seq(2:1)), IslandMainland_N_N, cex = 1)
text(-2, rev(seq(2:1)), IslandMainland_P_N, cex = 1)
text(-2, 4, "P-value", cex = 1.25, font = 2)
text(-3.5, 4, "n", cex = 1.25, font = 2)
text(-5.8, 4, "Sub Group", cex = 1.25, font = 2)
text(6.5, 4, "LRR [95% CI]", cex = 1.25, font = 2)
abline(h=0)
abline(h=3)
dev.off()
#### Run a meta-regression with seabird order as a predictor
mix_effect_order_N <-
rma.mv(
yi, # this is the column from our dataframe that contains calculated effect sizes
vi, # calculated variances
mods = ~ Order - 1, # Specify which moderators the model should take into account
random = ~ 1 | ID,
method = "REML", #Model fitting method
digits = 4, # round calculations to 4 digits
data = RR_effect_sizes_N)
# View results
mix_effect_order_N
# Create vector with order sample sizes
Order_N_N <- c(2, 5, 6, 10)
# Create vector with p-values for order
Order_P_N <- c(0.1855, 0.1177, "<0.0001", 0.4589)
# Forest plot for order
tiff(file = 'Nitrogen_Order_Plot.tiff', width = 20, height = 8.5, units = "cm", res = 500)
forest(mix_effect_order_N$b,
ci.lb = mix_effect_order_N$ci.lb,
ci.ub = mix_effect_order_N$ci.ub,
annotate = TRUE,
xlab = "ln(Response Ratio) - Nitrogen",
slab = c("Charadriiformes", "Procellariiformes", "Sphenisciformes", "Suliformes"),
cex = 1,
alim = c(-1, 5),
at = c(-1, 0, 1, 2, 3, 4, 5),
xlim = c(-7, 8),
ylim = c(-1.5, 5),
top = 1)
text(-3.5, rev(seq(4:1)), Order_N_N, cex = 1)
text(-2, rev(seq(4:1)), Order_P_N, cex = 1)
addpoly(REM_N, row= -1, cex=1, mlab="Full Model")
text(-3.5, -1, 25, cex = 1)
text(-2, -1, "<0.0001", cex = 1)
abline(h=0)
dev.off()
#### Run a meta-regression with Zone as a predictor
mix_effect_zone_N <-
rma.mv(
yi, # this is the column from our dataframe that contains calculated effect sizes
vi, # calculated variances
mods = ~ Zone - 1, # Specify which moderators the model should take into account
random = ~ 1 | ID,
method = "REML", #Model fitting method
digits = 4, # round calculations to 4 digits
data = RR_effect_sizes_N)
# View results
mix_effect_zone_N
# Create vector with zone sample sizes
Zone_N_N <- c(6, 16, 3)
# Create vector with p-values for zone
Zone_P_N <- c("<0.0001", 0.0937, 0.0254)
# Create forest plot of zone
tiff(file = 'Nitrogen_Zone_Plot.tiff', width = 20, height = 7.5, units = "cm", res = 500)
forest_Z_N <- forest(mix_effect_zone_N$b,
ci.lb = mix_effect_zone_N$ci.lb,
ci.ub = mix_effect_zone_N$ci.ub,
annotate = TRUE,
xlab = "ln(Response Ratio) - Nitrogen",
slab = c("Polar", "Temperate", "Tropical"),
cex = 1,
alim = c(-1, 5),
at = c(-1, 0, 1, 2, 3, 4, 5),
xlim = c(-7, 8),
ylim = c(-1.5, 4),
top = 1)
text(-3.5, rev(seq(3:1)), Zone_N_N, cex = 1)
text(-2, rev(seq(3:1)), Zone_P_N, cex = 1)
abline(h=0)
dev.off()
################################ PHOSPHOROUS ######################################
# Response Ratio for Phosphorous/P
RR_effect_sizes_P <- escalc("ROM",
m1i = Mean_seabird_P,
n1i = Seabird_sample_P,
sd1i = SD_seabird_P,
m2i = Mean_control_P,
n2i = Control_sample_P,
sd2i = SD_control_P,
data = SoilNutrients)
# Create unique ID number for each paper
RR_effect_sizes_P$ID <- seq.int(nrow(RR_effect_sizes_P))
##### Run the random effects model
REM_P <- rma.mv(yi, vi,
random = ~ 1 | ID,
slab = paste(Author, Year, sep = ", "),
data = RR_effect_sizes_P)
# View results
REM_P
#Back transform model (Exp) for easier interpretation
predict(REM_P, transf=exp, digits = 2)
# Create forest plot of phosphorous for all studies
tiff(file = 'Phosphorous_Plot.tiff', width = 25, height = 21, units = "cm", res = 500)
forest_P <- forest(RR_effect_sizes_P$yi, RR_effect_sizes_P$vi,
annotate = TRUE, slab = REM_P$slab,
xlab = "ln(Repsonse Ratio) - Phosphorous",
cex = 1,
cex.lab = 1.25,
ylim = c(-1, 20),
xlim = c(-7, 10),
alim = c(-3, 7),
at = c(-3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7),
)
addpoly(REM_P, row= -1, cex=1.25, mlab="Summary", fonts = "Arial")
abline(h=0)
text(-5.1, 19, "1st Author and Year", cex = 1.25, font = 2)
text(8.7, 19, "LRR [95% CI]", cex = 1.25, font = 2)
dev.off()
#### Run a meta-regression with island/mainland as a predictor
mix_effect_island_P <-
rma.mv(
yi, # this is the column from our dataframe that contains calculated effect sizes
vi, # calculated variances
mods = ~ Island_Mainland - 1, # Specify which moderators the model should take into account
random = ~ 1 | ID,
method = "REML", #Model fitting method
digits = 4, # round calculations to 4 digits
data = RR_effect_sizes_P)
# View results
mix_effect_island_P
# Create vector with island and mainland sample sizes
IslandMainland_N_P <- c(14, 3)
# Create vector with p-values for island and mainland
IslandMainland_P_P <- c("<0.0001", 0.0251)
# Forest plot for island v mainland
tiff(file = 'Phosphorous_Island_Plot.tiff', width = 20, height = 8, units = "cm", res = 500)
forest_IM_P <- forest(mix_effect_island_P$b,
ci.lb = mix_effect_island_P$ci.lb,
ci.ub = mix_effect_island_P$ci.ub,
annotate = TRUE,
xlab = "ln(Response Ratio) - Phosphorous",
slab = c("Island", "Mainland"),
xlim = c(-7, 8),
cex = 1,
alim = c(-1, 5),
at = c(-1, 0, 1, 2, 3, 4, 5),
digits = c(2L, 2L),
ylim = c(-1, 5),
top = 1)
text(-3.5, rev(seq(2:1)), IslandMainland_N_P, cex = 1)
text(-2, rev(seq(2:1)), IslandMainland_P_P, cex = 1)
text(-2, 4, "P-value", cex = 1.25, font = 2)
text(-3.5, 4, "n", cex = 1.25, font = 2)
text(-5.8, 4, "Sub Group", cex = 1.25, font = 2)
text(6.5, 4, "LRR [95% CI]", cex = 1.25, font = 2)
abline(h=0)
abline(h=3)
dev.off()
#### Run a meta-regression with seabird order as a predictor
mix_effect_order_P <-
rma.mv(
yi, # this is the column from our dataframe that contains calculated effect sizes
vi, # calculated variances
mods = ~ Order - 1, # Specify which moderators the model should take into account
random = ~ 1 | ID,
method = "REML", #Model fitting method
digits = 4, # round calculations to 4 digits
data = RR_effect_sizes_P)
# View results
mix_effect_order_P
# Create vector with order sample sizes
Order_N_P <- c(4, 3, 4, 5)
# Create vector with p-values for order
Order_P_P <- c("0.0090", "0.0049", "0.0001", "0.0001")
# Forest plot for order
tiff(file = 'Phosphorous_Order_Plot.tiff', width = 20, height = 8.5, units = "cm", res = 500)
forest(mix_effect_order_P$b,
ci.lb = mix_effect_order_P$ci.lb,
ci.ub = mix_effect_order_P$ci.ub,
annotate = TRUE,
xlab = "ln(Response Ratio) - Phosphorous",
slab = c("Charadriiformes", "Procellariiformes", "Sphenisciformes", "Suliformes"),
cex = 1,
alim = c(-1, 5),
at = c(-1, 0, 1, 2, 3, 4, 5),
xlim = c(-7, 8),
ylim = c(-1.5, 5),
top = 1)
text(-3.5, rev(seq(4:1)), Order_N_P, cex = 1)
text(-2, rev(seq(4:1)), Order_P_P, cex = 1)
addpoly(REM_P, row= -1, cex=1, mlab="Full Model")
text(-3.5, -1, 17, cex = 1)
text(-2, -1, "<0.0001", cex = 1)
abline(h=0)
dev.off()
#### Run a meta-regression with Zone as a predictor
mix_effect_zone_P <-
rma.mv(
yi, # this is the column from our dataframe that contains calculated effect sizes
vi, # calculated variances
mods = ~ Zone - 1, # Specify which moderators the model should take into account
random = ~ 1 | ID,
method = "REML", #Model fitting method
digits = 4, # round calculations to 4 digits
subset = c(3, 6, 10, 13, 19, 20, 24, 25, 26, 27, 28, 33, 37, 38, 39, 40),
data = RR_effect_sizes_P)
# View results
mix_effect_zone_P
# Create vector with zone sample sizes
Zone_N_P <- c(3, 13)
# Create vector with p-values for zone
Zone_P_P <- c(0.0052, "<0.0001")
# Create forest plot of zone
tiff(file = 'Phosphorous_Zone_Plot.tiff', width = 20, height = 6.5, units = "cm", res = 500)
forest_Z_P <- forest(mix_effect_zone_P$b,
ci.lb = mix_effect_zone_P$ci.lb,
ci.ub = mix_effect_zone_P$ci.ub,
annotate = TRUE,
xlab = "ln(Response Ratio) - Phosphorous",
slab = c("Polar", "Temperate"),
cex = 1,
alim = c(-1, 5),
at = c(-1, 0, 1, 2, 3, 4, 5),
xlim = c(-7, 8),
ylim = c(-1.5, 3),
top = 1)
text(-3.5, rev(seq(2:1)), Zone_N_P, cex = 1)
text(-2, rev(seq(2:1)), Zone_P_P, cex = 1)
abline(h=0)
dev.off()