Full model with t-distribution: 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, "bw-10-full-re-t.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("mu_lower", "mu_mean", "mu_upper")] <-
    predict(model, muni_data, "mu", "s(res_muni)", FUN = c95, intercept = FALSE)
muni_data[, c("sigma_lower", "sigma_mean", "sigma_upper")] <-
    predict(model, muni_data, "sigma", "s(res_muni)", FUN = function(x) c95(exp(x)),
            intercept = FALSE)

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

Visualize effects on mu

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

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

Visualize effects on sigma

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

ama %>%
    dplyr::filter(inc_study) %>%
    dplyr::arrange(sigma_mean) %>%
    dplyr::mutate(mun_name = factor(mun_name, mun_name)) %>%
    ggplot() +
    geom_errorbarh(aes(y = mun_name, xmin = sigma_lower, xmax = sigma_upper)) +
    geom_vline(xintercept = 1, linetype = 2) +
    geom_point(aes(y = mun_name, x = sigma_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_mu <- spdep::moran.mc(ama2$mu_mean, listw, nsim = 5000, alternative = "greater")
MC_mu

#> 
#>  Monte-Carlo simulation of Moran I
#> 
#> data:  ama2$mu_mean 
#> weights: listw  
#> number of simulations + 1: 5001 
#> 
#> statistic = 0.027519, observed rank = 3482, p-value = 0.3037
#> alternative hypothesis: greater

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

#> 
#>  Monte-Carlo simulation of Moran I
#> 
#> data:  ama2$mu_mean 
#> weights: listw  
#> number of simulations + 1: 5001 
#> 
#> statistic = 0.027519, observed rank = 3568, p-value = 0.2865
#> alternative hypothesis: greater

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>    user  system elapsed 
#>   8.562   0.360   9.122