Full model: municipality effects

Load packages, read data and source custom scripts

rm(list = ls())
library(bamlss)

#> Loading required package: coda

#> Loading required package: colorspace

#> Loading required package: mgcv

#> Loading required package: nlme

#> This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.

#> 
#> Attaching package: 'bamlss'

#> The following object is masked from 'package:mgcv':
#> 
#>     smooth.construct

library(gamlss.dist)

#> Loading required package: MASS

library(sf)

#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

library(ggplot2)

path_proj <- day2day::git_path()
path_data <- file.path(path_proj, "data")
path_processed <- file.path(path_data, "processed")
path_modelled <- file.path(path_data, "modelled")

source(file.path(path_proj, "src", "51-bamlss.R"))

bwdata_file <- file.path(path_processed, "bwdata_41_model.fst")
model_file <- file.path(path_modelled, "pre-10-full-re.rds")
form_file <- gsub("(\\.rds)$", "-form\\1", model_file)
model_file_burned <- gsub("(\\.rds)$", "-burned\\1", model_file)

bwdata_model <- fst::read_fst(bwdata_file)
form <- readRDS(form_file)
model <- readRDS(model_file_burned)

ama <- st_read(file.path(path_processed, "ama_05_62.gpkg"))

#> Reading layer `ama_05_62' from data source `/home/rstudio/Documents/Projects/data-amazonia/70-studies/birthweight/data/processed/ama_05_62.gpkg' using driver `GPKG'
#> Simple feature collection with 62 features and 373 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: 1573261 ymin: 8909613 xmax: 3542869 ymax: 10248330
#> projected CRS:  unnamed

Municipality effects

muni_data <- tibble::tibble(res_muni = unique(bwdata_model$res_muni))
muni_data[c("lower", "mean", "upper")] <-
    predict(model, muni_data, term = "s(res_muni)", FUN = function(x) c95(x),
            intercept = FALSE)

ama <- ama %>%
    dplyr::select(mun_code, inc_study, mun_name) %>%
    dplyr::left_join(muni_data, by = c("mun_code" = "res_muni"))

ggplot(ama) +
    geom_sf(aes(fill = mean)) +
    scale_fill_fermenter(palette = "RdBu", n.breaks = 7)

ama %>%
    dplyr::filter(inc_study) %>%
    dplyr::arrange(mean) %>%
    dplyr::mutate(mun_name = factor(mun_name, mun_name)) %>%
    ggplot() +
    geom_errorbarh(aes(y = mun_name, xmin = lower, xmax = upper)) +
    geom_vline(xintercept = 1, linetype = 2) +
    geom_point(aes(y = mun_name, x = mean), col = 2) +
    theme(axis.title.y = element_blank(),
          axis.title.x = element_blank())

Testing Spatial Correlation

ama2 <- ama %>% dplyr::filter(inc_study)

listw <- spdep::nb2listw(spdep::poly2nb(as(ama2, "Spatial")))

#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on
#> GRS80 ellipsoid in CRS definition

MC <- spdep::moran.mc(ama2$mean, listw, nsim = 5000, alternative = "greater")
MC

#> 
#>  Monte-Carlo simulation of Moran I
#> 
#> data:  ama2$mean 
#> weights: listw  
#> number of simulations + 1: 5001 
#> 
#> statistic = 0.1592, observed rank = 4796, p-value = 0.04099
#> alternative hypothesis: greater

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>    user  system elapsed 
#>   6.452   0.249   6.733