library(mgcv)
library(raster)
library(magrittr)
library(sp)

# Generate a dummy probability surface using Swaziland as an example
swz_elev <- raster::getData('alt', country="SWZ")
swz_prob <- swz_elev / cellStats(swz_elev, max)

# Generate some random points to act as sampling sites
candidate_sites <- coordinates(swz_elev)[which(!is.na(swz_elev[])),]
sample_pixels <- sample(1:nrow(candidate_sites), 100)
sampling_sites <- candidate_sites[sample_pixels,] %>% as.data.frame()

# Take a binomial sample of 100 at each site
sampling_sites$prob <- raster::extract(swz_prob, sampling_sites[,c("x", "y")])
sampling_sites$n_pos <- rbinom(100, 100, sampling_sites$prob )
sampling_sites$n_neg <- 100 - sampling_sites$n_pos

# Fit a GAM with bivaruate smooth on lat/lng
gam_mod <- gam(cbind(n_pos, n_neg) ~ s(x, y, bs="gp"),
               data = sampling_sites,
              family = "binomial")

# Look at predicted mean
prediction_data <- as.data.frame(candidate_sites)
predicted_prob <- predict(gam_mod, prediction_data, type="response")
predicted_prob_raster <- swz_elev
predicted_prob_raster[which(!is.na(swz_elev[]))] <- predicted_prob
plot(predicted_prob_raster)

# Simulate from the posterior
Cg <- predict(gam_mod, prediction_data, type = "lpmatrix")
sims <- rmvn(1000, mu = coef(gam_mod), V = vcov(gam_mod, unconditional = TRUE))
fits <- Cg %*% t(sims)
fits_prob <- exp(fits) / (1 + exp(fits))

#### Aggregate simulated values over districts
# Get district boundary polygons
swz_districts <- raster::getData("GADM", level=2, country="SWZ")

# Identify which district each pixel is in
which_district <- sp::over(SpatialPoints(prediction_data, proj4string = crs(swz_districts)), swz_districts)

# Calculate the district mean per draw from the posterior
district_posterior <- apply(fits_prob, 2, function(x){tapply(x, which_district$GID_2, mean)})

# Look at the posterior for the first district
hist(district_posterior[1,])