14-bootstrap.Rmd

# Bootstrap {#boot}

This chapter deals with bootstrap.

The students are expected to acquire the following knowledge:

- How to use bootstrap to generate coverage intervals.

<style>
.fold-btn { 
  float: right; 
  margin: 5px 5px 0 0;
}
.fold { 
  border: 1px solid black;
  min-height: 40px;
}
</style>

<script type="text/javascript">
$(document).ready(function() {
  $folds = $(".fold");
  $folds.wrapInner("<div class=\"fold-blck\">"); // wrap a div container around content
  $folds.prepend("<button class=\"fold-btn\">Unfold</button>");  // add a button
  $(".fold-blck").toggle();  // fold all blocks
  $(".fold-btn").on("click", function() {  // add onClick event
    $(this).text($(this).text() === "Fold" ? "Unfold" : "Fold");  // if the text equals "Fold", change it to "Unfold"or else to "Fold" 
    $(this).next(".fold-blck").toggle("linear");  // "swing" is the default easing function. This can be further customized in its speed or the overall animation itself.
  })
});
</script>

```{r, echo = FALSE, warning = FALSE, message = FALSE}
togs <- T
library(ggplot2)
library(dplyr)
library(reshape2)
library(tidyr)
# togs <- FALSE
```


```{exercise}

<span style="color:blue">Ideally, a $1-\alpha$ CI would have $1-\alpha$ coverage. That is, say a 95% CI should, in the long run, contain the true value of the parameter 95% of the time. In practice, it is impossible to assess the coverage of our CI method, because we rarely know the true parameter. In simulation, however, we can. Let's assess the coverage of bootstrap percentile intervals.</span>


a. <span style="color:blue">Pick a univariate distribution with readily available mean and one that you can easily sample from.</span>

b. <span style="color:blue">Draw $n = 30$ random samples from the chosen distribution and use the bootstrap (with large enough m) and percentile CI method to construct 95% CI. Repeat the process many times and count how many times the CI contains the true mean. That is, compute the actual coverage probability (don't forget to include the standard error of the coverage probability!). What can you observe?</span>

c. <span style="color:blue">Try one or two different distributions. What can you observe?</span>

d. <span style="color:blue">Repeat (b) and (c) using BCa intervals (R package boot). How does the coverage compare to percentile intervals?</span>

e. <span style="color:blue">As (d) but using intervals based on asymptotic normality (+/- 1.96 SE).</span>
  
f. <span style="color:blue">How do results from (b), (d), and (e) change if we increase the sample size to n = 200? What about n = 5?</span>

```
<div class = "fold">
```{r, echo = togs, message = FALSE, eval = togs, warning=FALSE}
library(boot)
set.seed(0)
nit   <- 1000  # Repeat the process "many times"
alpha <- 0.05  # CI parameter
nboot <- 100   # m parameter for bootstrap ("large enough m")
# f: change this to 200 or 5.
nsample <- 30  # n = 30 random samples from the chosen distribution. Comment out BCa code if it breaks.
covers     <- matrix(nrow = nit, ncol = 3)
covers_BCa <- matrix(nrow = nit, ncol = 3)
covers_asymp_norm <- matrix(nrow = nit, ncol = 3)

isin <- function (x, lower, upper) {
  (x > lower) & (x < upper)
}

for (j in 1:nit) {  # Repeating many times
  # a: pick a univariate distribution - standard normal
  x1 <- rnorm(nsample)
  
  # c: one or two different distributions - beta and poisson
  x2 <- rbeta(nsample, 1, 2)
  x3 <- rpois(nsample, 5)
  
  X1 <- matrix(data = NA, nrow = nsample, ncol = nboot)
  X2 <- matrix(data = NA, nrow = nsample, ncol = nboot)
  X3 <- matrix(data = NA, nrow = nsample, ncol = nboot)
  for (i in 1:nboot) {
    X1[ ,i] <- sample(x1, nsample, replace = T)
    X2[ ,i] <- sample(x2, nsample, T)
    X3[ ,i] <- sample(x3, nsample, T)
  }
  X1_func <- apply(X1, 2, mean)
  X2_func <- apply(X2, 2, mean)
  X3_func <- apply(X3, 2, mean)
  X1_quant <- quantile(X1_func, probs = c(alpha / 2, 1 - alpha / 2))
  X2_quant <- quantile(X2_func, probs = c(alpha / 2, 1 - alpha / 2))
  X3_quant <- quantile(X3_func, probs = c(alpha / 2, 1 - alpha / 2))
  covers[j,1] <- (0 > X1_quant[1]) & (0 < X1_quant[2])
  covers[j,2] <- ((1 / 3) > X2_quant[1]) & ((1 / 3) < X2_quant[2])
  covers[j,3] <- (5 > X3_quant[1]) & (5 < X3_quant[2])

  mf     <- function (x, i) return(mean(x[i]))
  bootX1 <- boot(x1, statistic = mf, R = nboot)
  bootX2 <- boot(x2, statistic = mf, R = nboot)
  bootX3 <- boot(x3, statistic = mf, R = nboot)

  X1_quant_BCa <- boot.ci(bootX1, type = "bca")$bca
  X2_quant_BCa <- boot.ci(bootX2, type = "bca")$bca
  X3_quant_BCa <- boot.ci(bootX3, type = "bca")$bca
  
  covers_BCa[j,1] <- (0 > X1_quant_BCa[4]) & (0 < X1_quant_BCa[5])
  covers_BCa[j,2] <- ((1 / 3) > X2_quant_BCa[4]) & ((1 / 3) < X2_quant_BCa[5])
  covers_BCa[j,3] <- (5 > X3_quant_BCa[4]) & (5 < X3_quant_BCa[5])
  
  # e: estimate mean and standard error
  # sample mean:
  x1_bar <- mean(x1)
  x2_bar <- mean(x2)
  x3_bar <- mean(x3)
  
  # standard error (of the sample mean) estimate: sample standard deviation / sqrt(n)
  x1_bar_SE <- sd(x1) / sqrt(nsample)
  x2_bar_SE <- sd(x2) / sqrt(nsample)
  x3_bar_SE <- sd(x3) / sqrt(nsample)
  
  covers_asymp_norm[j,1] <- isin(0, x1_bar - 1.96 * x1_bar_SE, x1_bar + 1.96 * x1_bar_SE)
  covers_asymp_norm[j,2] <- isin(1/3, x2_bar - 1.96 * x2_bar_SE, x2_bar + 1.96 * x2_bar_SE)
  covers_asymp_norm[j,3] <- isin(5, x3_bar - 1.96 * x3_bar_SE, x3_bar + 1.96 * x3_bar_SE)

}
apply(covers, 2, mean)
apply(covers, 2, sd) / sqrt(nit)

apply(covers_BCa, 2, mean)
apply(covers_BCa, 2, sd) / sqrt(nit)

apply(covers_asymp_norm, 2, mean)
apply(covers_asymp_norm, 2, sd) / sqrt(nit)

```
</div>




```{exercise}

<span style="color:blue">
  You are given a sample of independent observations from a process of interest:
  
  | Index | 1   | 2   | 3   | 4   | 5   | 6   | 7   | 8   |
  | :-- | --: | --:	| --:	| --: | --: | --: | --: | --: |
  | X   | 7   | 2   | 4   | 6   | 4   | 5   | 9   | 10  |
  
</span>


a. <span style="color:blue">Compute the plug-in estimate of mean and 95% symmetric CI based on asymptotic normality. Use the plug-in estimate of SE.</span>

b. <span style="color:blue">Same as (a), but use the unbiased estimate of SE.</span>

c. <span style="color:blue">Apply nonparametric bootstrap with 1000 bootstrap replications and estimate the 95% CI for the mean with percentile-based CI.</span>

```
<div class = "fold">
```{r, echo = togs, message = FALSE, eval = togs, warning=FALSE}
# a
x <- c(7, 2, 4, 6, 4, 5, 9, 10)
n <- length(x)
mu <- mean(x)

SE <- sqrt(mean((x - mu)^2)) / sqrt(n)
SE
z <- qnorm(1 - 0.05 / 2)
c(mu - z * SE, mu + z * SE)

# b
SE <- sd(x) / sqrt(n)
SE
c(mu - z * SE, mu + z * SE)

# c
set.seed(0)

m  <- 1000
T_mean <- function(x) {mean(x)}

est_boot <- array(NA, m)
for (i in 1:m) {
  x_boot <- x[sample(1:n, n, rep = T)]
  est_boot[i] <- T_mean(x_boot)
}

quantile(est_boot, p = c(0.025, 0.975))
```
</div> 


```{exercise}

<span style="color:blue">
  We are given a sample of 10 independent paired (bivariate) observations:
  
  | Index   | 1    | 2     | 3    | 4    | 5     | 6     | 7     | 8     | 9     | 10 |
  | :-- | --:  | --:	 | --:	| --:	 | --:   | --:   | --:   | --:   | --:   | --:   |
  | X   | 1.26 | -0.33 | 1.33 | 1.27 | 0.41  | -1.54 | -0.93 | -0.29 | -0.01 |  2.40 |
  | Y   | 2.64 | 0.33  | 0.48 | 0.06 | -0.88 | -2.14 | -2.21 | 0.95  | 0.83  | 1.45  |

</span>


a. <span style="color:blue">Compute Pearson correlation between X and Y.</span>

b. <span style="color:blue">Use the cor.test() from R to estimate a 95% CI for the estimate from (a).</span>

c. <span style="color:blue">Apply nonparametric bootstrap with 1000 bootstrap replications and estimate the 95% CI for the Pearson correlation with percentile-based CI.</span>
  
d. <span style="color:blue">Compare CI from (b) and (c). Are they similar?</span>
  
e. <span style="color:blue">How would the bootstrap estimation of CI change if we were interested in Spearman or Kendall correlation instead?</span>

```
<div class = "fold">
```{r, echo = togs, message = FALSE, eval = togs, warning=FALSE}
x <- c(1.26, -0.33,  1.33,  1.27,  0.41, -1.54, -0.93, -0.29, -0.01,  2.40)
y <- c(2.64,  0.33,  0.48,  0.06, -0.88, -2.14, -2.21,  0.95,  0.83,  1.45)

# a
cor(x, y)

# b
res <- cor.test(x, y)
res$conf.int[1:2]

# c
set.seed(0)
m  <- 1000
n  <- length(x) 
T_cor <- function(x, y) {cor(x, y)}

est_boot <- array(NA, m)
for (i in 1:m) {
  idx <- sample(1:n, n, rep = T) # !!! important to use same indices to keep dependency between x and y
  est_boot[i] <- T_cor(x[idx], y[idx])
}

quantile(est_boot, p = c(0.025, 0.975))

# d
# Yes, but the bootstrap CI is more narrow.

# e
# We just use the functions for Kendall/Spearman coefficients instead:
T_kendall <- function(x, y) {cor(x, y, method = "kendall")}
T_spearman <- function(x, y) {cor(x, y, method = "spearman")}

# Put this in a function that returns the CI
bootstrap_95_ci <- function(x, y, t, m = 1000) {
  n <- length(x)
  est_boot <- array(NA, m)
  for (i in 1:m) {
    idx <- sample(1:n, n, rep = T) # !!! important to use same indices to keep dependency between x and y
    est_boot[i] <- t(x[idx], y[idx])
  }
  quantile(est_boot, p = c(0.025, 0.975))
}

bootstrap_95_ci(x, y, T_kendall)
bootstrap_95_ci(x, y, T_spearman)

```
</div> 



```{exercise}

<span style="color:blue">
  In this problem we will illustrate the use of the nonparametric bootstrap for estimating CIs of regression model coefficients.
</span>


a. <span style="color:blue">Load the longley dataset from base R with data(longley).</span>

b. <span style="color:blue">Use lm() to apply linear regression using "Employed" as the target (dependent) variable and all other variables as the predictors (independent). Using lm() results, print the estimated regression coefficients and standard errors. Estimate 95% CI for the coefficients using +/- 1.96 * SE.</span>

c. <span style="color:blue">Use nonparametric bootstrap with 100 replications to estimate the SE of the coefficients from (b). Compare the SE from (c) with those from (b).</span>

```
<div class = "fold">
```{r, echo = togs, message = FALSE, eval = togs, warning=FALSE}
# a
data(longley)

# b
res <- lm(Employed ~ . , longley)
tmp <- data.frame(summary(res)$coefficients[,1:2])
tmp$LB <- tmp[,1] - 1.96 * tmp[,2]
tmp$UB <- tmp[,1] + 1.96 * tmp[,2]
tmp

# c
set.seed(0)
m <- 100
n <- nrow(longley)
T_coef <- function(x) {
  lm(Employed ~ . , x)$coefficients
}

est_boot <- array(NA, c(m, ncol(longley)))
for (i in 1:m) {
  idx <- sample(1:n, n, rep = T)
  est_boot[i,] <- T_coef(longley[idx,])
}

SE <- apply(est_boot, 2, sd)
SE

# Show the standard errors around coefficients
library(ggplot2)
library(reshape2)
df <- data.frame(index = 1:7, bootstrap_SE = SE, lm_SE = tmp$Std..Error)
melted_df <- melt(df[2:nrow(df), ], id.vars = "index")  # Ignore bias which has a really large magnitude
ggplot(melted_df, aes(x = index, y = value, fill = variable)) +
  geom_bar(stat="identity", position="dodge") +
  xlab("Coefficient") +
  ylab("Standard error") # + scale_y_continuous(trans = "log") # If you want to also plot bias
  
```
</div> 

```{exercise}

<span style="color:blue">
  This exercise shows a shortcoming of the bootstrap method when using the plug in estimator for the maximum.
</span>


a. <span style="color:blue">Compute the 95% bootstrap CI for the maximum of a standard normal distribution.</span>

b. <span style="color:blue">Compute the 95% bootstrap CI for the maximum of a binomial distribution with n = 15 and p = 0.2.</span>

c. <span style="color:blue">Repeat (b) using p = 0.9. Why is the result different?</span>
  
```
<div class = "fold">
```{r, echo = togs, message = FALSE, eval = togs, warning=FALSE}
# bootstrap CI for maximum

alpha <- 0.05
T_max <- function(x) {max(x)}  # Equal to T_max = max
bootstrap <- function(x, t, m = 1000) {
  n <- length(x)
  values <- rep(0, m)
  for (i in 1:m) {
    values[i] <- t(sample(x, n, replace = T))
  }
  quantile(values, probs = c(alpha / 2, 1 - alpha / 2))
}

# a
# Meaningless, as the normal distribution can yield arbitrarily large values.
x <- rnorm(100)
bootstrap(x, T_max)


# b
x <- rbinom(100, size = 15, prob = 0.2) # min = 0, max = 15
bootstrap(x, T_max)

# c
x <- rbinom(100, size = 15, prob = 0.9) # min = 0, max = 15
bootstrap(x, T_max)

# Observation: to estimate the maximum, we need sufficient probability mass near the maximum value the distribution can yield.
# Using bootstrap is pointless when there is too little mass near the true maximum.
# In general, bootstrap will fail when estimating the CI for the maximum.
```
</div> 

```{exercise, name = "Practical - and fictional - coverage interval comparison"}
In this exercise, we investigate how different kinds of CI's behave as we vary the number of measurements.

The story behind the data: it's 2025 and we've discovered that Slovenia has rich deposits of a rare mineral called Moustachium, which can be used to accelerate moustache growth. This mineral is highly sought, so the government has decided to contract two different companies to provide information on where to begin mining. Both companies investigated mining sites in each statistical region and gave their best estimate of the average Moustachium concentration in tonnes per square kilometer. The Data Science team has been called to estimate the uncertainty in these estimates and help avoid mining in the wrong region.

Generate synthetic data with the script below:

    set.seed(0)
    
    library(comprehenr)
    
    regions <- c("pomurska", "podravska", "koroska", "savinjska", "zasavska", "posavska", "JV Slovenija", "primorsko-notranjska", "osrednjeslovenska", "gorenjska", "goriska", "obalno-kraska")
    region_rates <- seq(1.3, 2.3, length.out=length(regions))
    region_rates <- region_rates[sample.int(length(regions), length(regions))]

    make_dataset <- function(n_contractors) {
        measurements <- matrix(nrow=length(regions), ncol=n_contractors)
        for (i in 1:length(regions)) {
            measurements[i,] <- rgamma(n_contractors, 5.0, region_rates[i])
        }

        df <- data.frame(measurements)
        row.names(df) <- regions
        names(df) <- to_vec(for(i in 1:n_contractors) paste("Contractor", i))
        return(df)
    }

    set.seed(0)
    df_2025 <- make_dataset(2)

    set.seed(0)
    df_2027 <- make_dataset(10)

    set.seed(0)
    df_2028 <- make_dataset(100)

    set.seed(0)
    df_2029 <- make_dataset(1000)

    saveRDS(df_2025, file="moustachium_2025.Rda")
    saveRDS(df_2027, file="moustachium_2027.Rda")
    saveRDS(df_2028, file="moustachium_2028.Rda")
    saveRDS(df_2029, file="moustachium_2029.Rda")

  a. <span style="color:blue">Estimate the average concentration for different regions.</span>

  b. <span style="color:blue">Estimate the average concentration uncertainty using 95% CI's (asymptotic normality with biased and unbiased standard error, standard bootstrap CI, bootstrap percentile CI).</span>

  c. <span style="color:blue">Visualize uncertainties with a histogram and discuss the best location to start mining.</span>
  
  d. <span style="color:blue">The year is 2027 and the government has decided to contract 10 companies. Rerun the code with new measurements and discuss how CI's change.</span>

  e. <span style="color:blue">Technological advancements in robotics have enabled site surveys on a massive scale. Repeat the last point for 100 surveyor robots in 2028 and 1000 surveyor robots in 2029.</span>
```
<div class = "fold">
```{r, echo = togs, message = FALSE, eval = togs, warning=FALSE}
library(ggplot2)
library(dplyr)
library(data.table)

set.seed(0)

input_dataset_path = "moustachium_2025.Rda"  # Change this for points d and e
output_plot_path = "moustachium_2025.pdf"  # Change this for points d and e

df <- readRDS(input_dataset_path)  # Data comes from here
n_contractors <- ncol(df)
results_df <- data.frame(region=row.names(df))  # Store CI bounds here

# 1. average concentration for different mining sites
results_df$average_concetration <- rowMeans(df)

# CI for the mean based on asymptotic normality (biased SE estimate)
biased_SE <- sqrt(apply(df, 1, function(vec) {sum((vec - mean(vec))^2) / length(vec)}) / n_contractors)
results_df$ci95.an.biased_var.low <- results_df$average_concetration - 1.96 * biased_SE
results_df$ci95.an.biased_var.high <- results_df$average_concetration + 1.96 * biased_SE

# CI for the mean based on asymptotic normality (unbiased SE estimate)
unbiased_SE <- sqrt(apply(df, 1, var) / n_contractors)
results_df$ci95.an.unbiased_var.low <- results_df$average_concetration - 1.96 * unbiased_SE
results_df$ci95.an.unbiased_var.high <- results_df$average_concetration + 1.96 * unbiased_SE

# Standard bootstrap CI with 1000 samples
bootstrap_variance <- function(data, n_samples) {
  # n_samples is m in pseudocode
  output <- numeric(n_samples)
  for (i in 1:n_samples) {
    index <- sample(1:length(data), length(data), rep = TRUE)
    resampled_data <- data[index]
    output[i] <- mean(resampled_data)
  }
  return(var(output))
}

bootstrap_1000_sd <- sqrt(apply(df, 1, function(vec){bootstrap_variance(vec, 1000)}))
results_df$ci95.bootstrap.standard.1000.low <- results_df$average_concetration - 1.96 * bootstrap_1000_sd
results_df$ci95.bootstrap.standard.1000.high <- results_df$average_concetration + 1.96 * bootstrap_1000_sd

# Bootstrap percentile CI with 1000 samples

bootstrap_quantile <- function(data, functional, n_samples, probs) {
    # n_samples is m in pseudocode
    output <- numeric(n_samples)
    for (i in 1:n_samples) {
        index <- sample(1:length(data), length(data), rep = TRUE)
        resampled_data <- data[index]
        output[i] <- functional(resampled_data)
    }
    return(quantile(output, probs=probs))
}

results_df$ci95.bootstrap.percentile.1000.low <- apply(df, 1, function(vec){bootstrap_quantile(vec, mean, 1000, 0.025)})
results_df$ci95.bootstrap.percentile.1000.high <- apply(df, 1, function(vec){bootstrap_quantile(vec, mean, 1000, 0.975)})

results_df

# Visualization: we use a bar chart with uncertainty bands

plot_moustachium_per_region <- function(region_names, average, ci_low, ci_high) {
    df_visualization <- data.frame(region=region_names, average=average, low=ci_low, high=ci_high)
    ggplot(df_visualization, aes(x=region, y=average)) + geom_bar(stat="identity")
}

mask <- endsWith(colnames(results_df), "low")
mask[c(1, 2)] <- T
results_df_low <- results_df[, mask]
colnames(results_df_low) <- gsub('.low','', colnames(results_df_low))

mask <- endsWith(colnames(results_df), "high")
mask[c(1, 2)] <- T
results_df_high <- results_df[, mask]
colnames(results_df_high) <- gsub('.high','', colnames(results_df_high))

long_results_df_low <- melt(setDT(results_df_low), id.vars=c("region", "average_concetration"))
names(long_results_df_low) <- c("region", "average_concentration", "variable", "low")

long_results_df_high <- melt(setDT(results_df_high), id.vars=c("region", "average_concetration"))
names(long_results_df_high) <- c("region", "average_concentration", "variable", "high")

long_results_df <- merge(long_results_df_low, long_results_df_high, by=c("region", "variable", "average_concentration"), all=T)

moustachium_plot <- ggplot(long_results_df, aes(x=region, y=average_concentration)) +
    geom_bar(stat="identity", position="dodge", alpha=0.2) +
    geom_errorbar(aes(ymin=low, ymax=high, color=variable), width=0.2, position=position_dodge(0.9)) +
    scale_x_discrete(guide = guide_axis(angle = 45)) +
    ylim(-1, 8)

# ggsave(plot=moustachium_plot, width=12, height=8, dpi=300, filename=output_plot_path)
moustachium_plot

# Visualization: we can also use a map. Circle size denotes concentration in region, low transparency denotes high uncertainty.

library(maps)

map_data_slo <- map_data('world')[map_data('world')$region == "Slovenia",]

map_df <- long_results_df[long_results_df$variable == "ci95.an.biased_var", ]

# VERY approximate longitudes and latitudes for different regions.
map_df$long <- rep(0, nrow(map_df))
map_df$lat <- rep(0, nrow(map_df))

map_df[map_df$region == "gorenjska"]$long <- 14.2
map_df[map_df$region == "gorenjska"]$lat <- 46.3

map_df[map_df$region == "goriska"]$long <- 13.85
map_df[map_df$region == "goriska"]$lat <- 46.0

map_df[map_df$region == "obalno-kraska"]$long <- 13.9
map_df[map_df$region == "obalno-kraska"]$lat <- 45.65

map_df[map_df$region == "osrednjeslovenska"]$long <- 14.5
map_df[map_df$region == "osrednjeslovenska"]$lat <- 46.

map_df[map_df$region == "primorsko-notranjska"]$long <- 14.3
map_df[map_df$region == "primorsko-notranjska"]$lat <- 45.7

map_df[map_df$region == "zasavska"]$long <- 15
map_df[map_df$region == "zasavska"]$lat <- 46.1

map_df[map_df$region == "savinjska"]$long <- 15.2
map_df[map_df$region == "savinjska"]$lat <- 46.25

map_df[map_df$region == "posavska"]$long <- 15.4
map_df[map_df$region == "posavska"]$lat <- 46

map_df[map_df$region == "koroska"]$long <- 15.1
map_df[map_df$region == "koroska"]$lat <- 46.5

map_df[map_df$region == "podravska"]$long <- 15.7
map_df[map_df$region == "podravska"]$lat <- 46.45

map_df[map_df$region == "pomurska"]$long <- 16.2
map_df[map_df$region == "pomurska"]$lat <- 46.65

map_df[map_df$region == "JV Slovenija"]$long <- 15.
map_df[map_df$region == "JV Slovenija"]$lat <- 45.7

map_df$ci_size <- (map_df$high - map_df$low)
map_df$ci_y <- map_df$lat - 0.05
map_df$ci_label <- sprintf("(%.2f, %.2f)", map_df$low, map_df$high)
map_df$avg_label <- sprintf("%.2f", map_df$average_concentration)

country_plot <- ggplot() +
    # First layer: worldwide map
    geom_polygon(data = map_data("world"),
                 aes(x=long, y=lat, group = group),
                 color = '#9c9c9c', fill = '#f3f3f3') +
    # Second layer: Country map
    geom_polygon(
        data = map_data_slo,
        aes(x=long, y=lat, group = group),
        color='darkgreen', 
        fill='green', 
        alpha=0.2
    ) +
    geom_point(data=map_df, aes(x=long, y=lat, fill=region, size=average_concentration, alpha=ci_size), color="black", pch=21) +
    geom_text(data=map_df, aes(x=long, y=ci_y, label=ci_label), size=3) +
    geom_text(data=map_df, aes(x=long, y=lat, label=avg_label), size=3) +
    scale_size_continuous(range = c(3, 12), trans = "exp") +
    scale_alpha_continuous(range = c(0.15, 0.75), trans = "reverse") +
    ggtitle("Estimated average Moustachium concentration with 95% CI") +
    coord_cartesian(xlim=c(13.2, 16.7), ylim=c(45.4, 47.))

# ggsave(plot=country_plot, width=18, height=12, dpi=300, filename="country.pdf")
country_plot
```
</div>