# data
library(readr) # reading CSV files
library(here) # constructing file paths
library(ggeffects) # extracting model predictions
library(tibble) # creating tibbles (data frames)
library(dplyr) # data manipulation
# plots
library(scales) # scale transformations and labels
library(ggplot2) # data visualization
library(patchwork) # combining plots09 - Results from generalized linear models (Figure 4)
09_V_model_Fig4.qmd
Model predictions and figures in the main text
Loading packages, data and functions
Note
The scales package needs the version 1.3 @420 to work
Reading model results for counting, DR, DC and native and NBT turnover
# Data from 03_C_data_preparation.qmd
df_country_complete6 <- readr::read_csv(here::here("data", "processed", "df_country_complete6.csv"))
# Data from 04_D_model_NBTs.qmd
mod_counting_NBT <- readRDS(here::here("output",
"models",
"model_res_counting.rds")) # NBT total countings
mod_DC <- readRDS(here::here("output",
"models",
"model_res_dc.rds"))
mod_DR <- readRDS(here::here("output",
"models",
"model_res_dr.rds"))
mod_turnover_native <- readRDS(here::here("output",
"models",
"model_res_turnover_native.rds"))
mod_turnover_nbt <- readRDS(here::here("output",
"models",
"model_res_turnover_nbt.rds"))
# function
sqrt_trans <- scales::trans_new("sqrt_with_zero",
transform = function(x) ifelse(x == 0, 0, sqrt(x)),
inverse = function(x) x^2)Total NBT counting
# Extracting predicted values
pred_mod_counting_gdp <-
ggeffects::ggpredict(mod_counting_NBT, "gdp.pad[-0.88:4 by=.05]")
pred_mod_counting_richness <-
ggeffects::ggpredict(mod_counting_NBT, "native.richness.pad[-0.4:8.6 by=.05]")
# scaling back response variables
source(here::here("R", "functions", "function_scale_back.R"))
## GDP
scale_gdp <- scale(df_country_complete6$e_migdppc)
pred_mod_counting_gdp2 <- scale_back(data = pred_mod_counting_gdp, scaled = scale_gdp)
## Native Richness
scale_native_richness <- scale(df_country_complete6$native_richness)
pred_mod_counting_richness2 <- scale_back(data = pred_mod_counting_richness, scaled = scale_native_richness)
# plotting figure 4a
a <-
tibble(
x= df_country_complete6$e_migdppc,
y = df_country_complete6$total_country_museum,
country = df_country_complete6$country,
region = df_country_complete6$region
) |>
ggplot(aes(x=x,y=y))+
geom_ribbon(data = pred_mod_counting_gdp2,
aes(x = x_original,
y = predicted,
ymin = conf.low,
ymax = conf.high),
alpha = 0.2)+
geom_line(data = pred_mod_counting_gdp2,
aes(x = x_original,
y = predicted))+
geom_point(aes(fill = region),
shape = 21)+
scale_x_continuous(expand = expansion(mult = c(0,0)),
breaks = seq(0,40000,10000),
labels = scales::label_number(scale_cut = scales::cut_short_scale()
))+
scale_y_continuous(trans = sqrt_trans,
expand = expansion(mult = c(0,0)),
breaks = c(0,100,1000,5000),
labels = scales::label_number(scale_cut = scales::cut_short_scale())
)+
guides(
y = guide_axis_logticks()
)+
scale_fill_manual(
values = c(
"Europe & Central Asia" = "#E64B35FF",
"East Asia & Pacific" = "#4DBBD5FF",
"North America" = "#3C5488FF",
"South Asia" = "#B09C85FF",
"Latin America & Caribbean" = "#00A087FF",
"Sub-Saharan Africa" = "#F39B7FFF",
"Middle East & North Africa" = "#8491B4FF"
)
)+
labs(x = "Gross Domestic Product - GDP",
y = "Number of NBT") +
theme_classic()+
theme(
axis.title = element_text(face = "bold"),
axis.text = element_text(color = "black"),
axis.line = element_blank(),
axis.ticks.x = element_blank(),
panel.grid.major = element_line(linetype = "dashed"),
legend.position = "none"
)+
coord_cartesian(clip = "off",
xlim = c(0,40000))
# plotting figure 4b
b <- tibble(
x= df_country_complete6$native_richness,
y = df_country_complete6$total_country_museum,
country = df_country_complete6$country,
region = df_country_complete6$region
) |>
ggplot(aes(x=x,y=y))+
geom_ribbon(data = pred_mod_counting_richness2,
aes(x = x_original,
y = predicted,
ymin = conf.low,
ymax = ifelse(conf.high >= 20000, 20000, conf.high)),
alpha = 0.2)+
geom_line(data = pred_mod_counting_richness2,
aes(x = x_original,
y = predicted))+
annotate(
"segment",
x = 2317.744, xend = Inf,
y = Inf, yend = Inf,
linetype = "dotted"
)+
geom_point(aes(fill = region),
shape = 21)+
scale_x_continuous(expand = expansion(mult = c(0,0)),
breaks = seq(0,3000,1000),
labels = scales::label_number(scale_cut = scales::cut_short_scale())
)+
scale_y_continuous(trans = sqrt_trans,
expand = expansion(mult = c(0,0)),
breaks = c(0,1000,10000,100000,300000),
labels = scales::label_number(scale_cut = scales::cut_short_scale()))+
guides(
y = guide_axis_logticks()
)+
scale_fill_manual(
values = c(
"Europe & Central Asia" = "#E64B35FF",
"East Asia & Pacific" = "#4DBBD5FF",
"North America" = "#3C5488FF",
"South Asia" = "#B09C85FF",
"Latin America & Caribbean" = "#00A087FF",
"Sub-Saharan Africa" = "#F39B7FFF",
"Middle East & North Africa" = "#8491B4FF"
)
)+
labs(x = "Native richness",
y = "Number of NBT") +
theme_classic()+
theme(
axis.title = element_text(face = "bold"),
axis.text = element_text(color = "black"),
axis.line = element_blank(),
axis.ticks.x = element_blank(),
panel.grid.major = element_line(linetype = "dashed"),
legend.position = "none"
)+
coord_cartesian(clip = "off",
xlim = c(-100,NA),
ylim = c(NA, 20000))Domestic Contribution
# predicted values
pred_mod_dc <-
ggeffects::ggpredict(mod_DC, "gdp.pad[-0.88:4 by=.05]")
pred_mod_dc_gdp2 <- scale_back(data = pred_mod_dc, scaled = scale_gdp)
# plotting figure 4c
c <- tibble(
x = df_country_complete6$e_migdppc,
y = df_country_complete6$prop_DC,
country = df_country_complete6$country,
region = df_country_complete6$region
) |>
ggplot(aes(x=x,y=y))+
geom_ribbon(data = pred_mod_dc_gdp2,
aes(x = x_original,
y = predicted,
ymin = conf.low,
ymax = conf.high),
alpha = 0.2)+
geom_line(data = pred_mod_dc_gdp2,
aes(x = x_original,
y = predicted))+
geom_point(aes(fill = region),
shape = 21)+
scale_x_continuous(expand = expansion(mult = c(0,0)),
breaks = seq(0,40000,10000),
labels = scales::label_number(scale_cut = scales::cut_short_scale())
) +
scale_y_continuous(expand = expansion(mult = c(0)),
labels = scales::label_number(scale_cut = scales::cut_short_scale()))+
scale_fill_manual(
values = c(
"Europe & Central Asia" = "#E64B35FF",
"East Asia & Pacific" = "#4DBBD5FF",
"North America" = "#3C5488FF",
"South Asia" = "#B09C85FF",
"Latin America & Caribbean" = "#00A087FF",
"Sub-Saharan Africa" = "#F39B7FFF",
"Middle East & North Africa" = "#8491B4FF"
)
)+
labs(x = "Gross Domestic Product - GDP",
y = "Domestic Contribution") +
theme_classic()+
theme(
axis.title = element_text(face = "bold"),
axis.text = element_text(color = "black"),
axis.line = element_blank(),
axis.ticks = element_blank(),
panel.grid.major = element_line(linetype = "dashed",
lineend = "round"),
legend.position = "none"
)+
coord_cartesian(clip = "off",
xlim = c(0,40000),
ylim = c(0,1))Domestic Retention
# predicted values
pred_mod_dr <-
ggeffects::ggpredict(mod_DR, "gdp.pad[-0.88:4 by=.05]")
pred_mod_dr_gdp2 <- scale_back(data = pred_mod_dr, scaled = scale_gdp)
# plotting figure 4d
d <- tibble(
x = df_country_complete6$e_migdppc,
y = df_country_complete6$prop_DR,
country = df_country_complete6$country,
region = df_country_complete6$region
) |>
ggplot(aes(x=x,y=y))+
geom_ribbon(data = pred_mod_dr_gdp2,
aes(x = x_original,
y = predicted,
ymin = conf.low,
ymax = conf.high),
alpha = 0.2)+
geom_line(data = pred_mod_dr_gdp2,
aes(x = x_original,
y = predicted))+
geom_point(aes(fill = region),
shape = 21)+
scale_x_continuous(expand = expansion(mult = c(0,0)),
breaks = seq(0,40000,10000),
labels = scales::label_number(scale_cut = scales::cut_short_scale())
) +
scale_y_continuous(expand = expansion(mult = c(0)),
labels = scales::label_number(scale_cut = scales::cut_short_scale()))+
scale_fill_manual(
values = c(
"Europe & Central Asia" = "#E64B35FF",
"East Asia & Pacific" = "#4DBBD5FF",
"North America" = "#3C5488FF",
"South Asia" = "#B09C85FF",
"Latin America & Caribbean" = "#00A087FF",
"Sub-Saharan Africa" = "#F39B7FFF",
"Middle East & North Africa" = "#8491B4FF"
)
)+
labs(x = "Gross Domestic Product - GDP",
y = "Domestic Retention") +
theme_classic()+
theme(
axis.title = element_text(face = "bold"),
axis.text = element_text(color = "black"),
axis.line = element_blank(),
axis.ticks = element_blank(),
panel.grid.major = element_line(linetype = "dashed",
lineend = "round"),
legend.position = "none"
)+
coord_cartesian(clip = "off",
xlim = c(0,40000),
ylim = c(0,1))Native turnover
pred_mod_beta_native <-
ggeffects::ggpredict(mod_turnover_native, "gdp.pad[-0.88:4 by=.05]")
pred_mod_beta_native2 <- scale_back(data = pred_mod_beta_native, scaled = scale_gdp)
# plotting figure 4e
e <- tibble(
x = df_country_complete6$e_migdppc,
y = df_country_complete6$native.beta.model,
country = df_country_complete6$country,
region = df_country_complete6$region
) |>
ggplot(aes(x=x,y=y))+
geom_ribbon(data = pred_mod_beta_native2,
aes(x = x_original,
y = predicted,
ymin = conf.low,
ymax = conf.high),
alpha = 0.2)+
geom_line(data = pred_mod_beta_native2,
aes(x = x_original,
y = predicted))+
geom_point(aes(fill = region),
shape = 21)+
scale_x_continuous(expand = expansion(mult = c(0,0)),
breaks = seq(0,40000,10000),
labels = scales::label_number(scale_cut = scales::cut_short_scale())
)+
scale_y_continuous(expand = expansion(mult = c(0,0)),
labels = scales::label_number(scale_cut = scales::cut_short_scale()))+
scale_fill_manual(
values = c(
"Europe & Central Asia" = "#E64B35FF",
"East Asia & Pacific" = "#4DBBD5FF",
"North America" = "#3C5488FF",
"South Asia" = "#B09C85FF",
"Latin America & Caribbean" = "#00A087FF",
"Sub-Saharan Africa" = "#F39B7FFF",
"Middle East & North Africa" = "#8491B4FF"
)
)+
labs(x = "Gross Domestic Product - GDP",
y = "Native turnover") +
theme_classic()+
theme(
axis.title = element_text(face = "bold"),
axis.text = element_text(color = "black"),
axis.line = element_blank(),
axis.ticks = element_blank(),
panel.grid.major = element_line(linetype = "dashed",
lineend = "round"),
legend.position = "none"
)+
coord_cartesian(clip = "off",
xlim = c(0,40000),
ylim = c(0,1))NBT turnover
# predictions
pred_mod_beta_type <-
ggeffects::ggpredict(mod_turnover_nbt, "gdp.pad[-0.88:4 by=.05]")
pred_mod_beta_type2 <- scale_back(data = pred_mod_beta_type, scaled = scale_gdp)
# plotting figure 4f
f <- tibble(
x = df_country_complete6$e_migdppc,
y = df_country_complete6$type.beta.model,
country = df_country_complete6$country,
region = df_country_complete6$region
) |>
ggplot(aes(x=x,y=y))+
geom_ribbon(data = pred_mod_beta_type2,
aes(x = x_original,
y = predicted,
ymin = conf.low,
ymax = conf.high),
alpha = 0.2)+
geom_line(data = pred_mod_beta_type2,
aes(x = x_original,
y = predicted))+
geom_point(aes(fill = region),
shape = 21)+
scale_x_continuous(expand = expansion(mult = c(0,0)),
breaks = seq(0,40000,10000),
labels = scales::label_number(scale_cut = scales::cut_short_scale())) +
scale_y_continuous(expand = expansion(mult = c(0)),
labels = scales::label_number(scale_cut = scales::cut_short_scale()))+
scale_fill_manual(
values = c(
"Europe & Central Asia" = "#E64B35FF",
"East Asia & Pacific" = "#4DBBD5FF",
"North America" = "#3C5488FF",
"South Asia" = "#B09C85FF",
"Latin America & Caribbean" = "#00A087FF",
"Sub-Saharan Africa" = "#F39B7FFF",
"Middle East & North Africa" = "#8491B4FF"
)
)+
labs(x = "Gross Domestic Product - GDP",
y = "NBT turnover") +
theme_classic()+
theme(
axis.title = element_text(face = "bold"),
axis.text = element_text(color = "black"),
axis.line = element_blank(),
axis.ticks = element_blank(),
panel.grid.major = element_line(linetype = "dashed",
lineend = "round"),
legend.position = "none"
)+
coord_cartesian(clip = "off",
xlim = c(0,40000),
ylim = c(0,1))Joining all plots
fig_model <-
a+b+c+d+e+f+
patchwork::plot_annotation(tag_levels = "a")+
patchwork::plot_layout(nrow = 3)&
theme(
plot.margin = margin(5.5,8,5.5,5.5,"pt")
)
fig_model
# saving model figure
ggsave(here::here("output", "figures", "Fig4_models.png"),
width = 6, height = 8,
dpi = 600, plot = fig_model)