Create amazonas penalty for GMRF

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(dplyr)

#> 
#> Attaching package: 'dplyr'

#> The following object is masked from 'package:MASS':
#> 
#>     select

#> The following object is masked from 'package:bamlss':
#> 
#>     n

#> The following object is masked from 'package:nlme':
#> 
#>     collapse

#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag

#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

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")

bwdata_model <- fst::read_fst(file.path(path_processed, "bwdata_41_model.fst"))
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

Markov random field penalty

Penalty is constructed based on the neighbors.

ama <- ama %>%
    dplyr::filter(inc_study) %>%
    dplyr::select(mun_code)
ama <- ama[match(levels(bwdata_model$res_muni), ama$mun_code), ]
K <- neighbormatrix(ama)
colnames(K) <- ama$mun_code
rownames(K) <- ama$mun_code

saveRDS(K, file = file.path(path_processed, "ama_10_penalty.rds"))

Visualize neighbors

ama <- mutate(ama, neighbors = diag(K))
plot(ama["neighbors"])

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>    user  system elapsed 
#>  15.519   0.653  16.151