theUrl = 'https://data.worldpop.org/GIS/Population/Global_2000_2020_1km_UNadj/2015/PER/per_ppp_2015_1km_Aggregated_UNadj.tif'


theFile = file.path(normalizePath("~/Downloads"), basename(theUrl))
if(!file.exists(theFile)) download.file(theUrl, theFile)
library(terra)
## terra 1.8.54
pop = rast(theFile)
pCol = mapmisc::colourScale(breaks=c(0,1,5, 10,100, 1000,10000,40000), style='fixed')
plot(pop, col=pCol$col, breaks=pCol$breaks)
plot of chunk peruPop
plot of chunk peruPop

Aggregate:

peruUrl= "http://pbrown.ca/teaching/brazil2025/data/peru.RData"
peruFile = file.path(normalizePath("~/Downloads"), basename(peruUrl))
if(!file.exists(peruFile)) download.file(peruUrl, peruFile)
(load(peruFile))
## [1] "peruWrap"    "peruHospDef"
peru = unwrap(peruWrap)
peruT = project(peru, crs(pop))
popAgg = extract(pop, peruT, fun='sum', na.rm=TRUE)
peru$pop = popAgg$per_ppp_2015
sum(peru$pop,na.rm=TRUE)
## [1] 30278996
sum(values(pop), na.rm=TRUE)
## [1] 30470739

we’ll do all hospitalizations, years 2012-2015

peruY1 = peruHospDef[peruHospDef$ano %in% 2012:2016,]
peruY = aggregate(peruY1[, grep("hosp|defu", names(peruY1))], peruY1[,'name', drop=FALSE], sum, na.rm=TRUE)
peruForModel = merge(peru, peruY, by='name')
peruForModel$y = peruForModel$hospitalizados_men5 + peruForModel$hospitalizados_may5 
peruForModel$rawRate = peruForModel$y / (5*peruForModel$pop)

multiply pop by 5, because we have 5 years data.

A plot

peruCol = mapmisc::colourScale(peruForModel$rawRate, breaks=9, digits=1, col='Reds',transform='sqrt',style='equal')
plot(peruForModel, col=peruCol$plot, lwd=0.2)
mapmisc::legendBreaks('bottomleft', peruCol)
plot of chunk plot
plot of chunk plot
peruForModel$logPop = log(5*peruForModel$pop)
peruForModel2 = peruForModel[
        which(peruForModel$pop > 0), 
    ]
adjMat1 = adjacent(peruForModel2, symmetric=TRUE)
adjMat = Matrix::sparseMatrix(i=adjMat1[,1], j=adjMat1[,2], symmetric=TRUE)
Nneighbours = apply(adjMat, 2, sum)
plot(peruForModel2, col='red', lty=0)
plot of chunk findNeighbours
plot of chunk findNeighbours
plot(peruForModel2[which(Nneighbours==0),],  col='green')
plot of chunk findNeighbours
plot of chunk findNeighbours
peruForModel3 = peruForModel2[Nneighbours > 0, ]
peruForModel3$y[325] = NA
peruForModel3$area = terra::expanse(peruForModel3)/1000000
peruForModel3$popdens = peruForModel3$pop/peruForModel3$area
peruForModel3$logPopDens = round(log(peruForModel3$popdens), 2)
result = diseasemapping::bym(
    y ~ offset(logPop) + 
        f(logPopDens, model='rw2', prior = 'pc.prec', param=c(0.1, 0.5)),
    data=peruForModel3,
    family = "poisson",
    prior = list(
        sd = 0.5, # prior medians
        propSpatial = 0.5), verbose=TRUE
    )
knitr::kable(result$parameters$summary[,3:5], digits=3)
0.025quant 0.5quant 0.975quant
(Intercept) -10.703 -10.551 -10.406
sd 1.873 2.018 2.177
propSpatial 0.532 0.632 0.729 
100*1000*exp(result$parameters$summary[1,grep("quant", names(result$parameters$summary))])
##             0.025quant 0.5quant 0.975quant
## (Intercept)   2.247267 2.616112   3.026264
cbind(
    values(peruForModel3),
    f=result$data$random.0.5quant)[
    order(result$data$random.0.5quant, decreasing=TRUE)[c(1:6, 
        seq(to=length(peruForModel3), len=6))],c('y','pop', 'f')]
##        y       pop         f
## 429   67  2687.386  6.162024
## 1515 164  8191.152  5.843529
## 241  102  4816.844  5.168600
## 973  294 25701.339  5.113633
## 1437  62  8704.869  5.080198
## 829  773 51276.419  5.043504
## 447    0 16525.565 -2.775547
## 1770   0 19893.445 -2.787422
## 557    0 39435.054 -2.893428
## 1658   0  1149.717 -3.129291
## 1659   0  3553.362 -3.156465
## 1752   0 29887.964 -3.391562
head(result$inla$summary.random$logPopDens)
##      ID       mean        sd 0.025quant   0.5quant 0.975quant       mode
## 1 -2.18 -0.1397170 0.7011565  -1.499184 -0.1478326  1.2635881 -0.1483604
## 2 -1.77 -0.2112027 0.5702219  -1.314355 -0.2177816  0.9280450 -0.2181187
## 3 -1.67 -0.2290720 0.5415962  -1.276575 -0.2352513  0.8524269 -0.2355526
## 4 -1.51 -0.2581742 0.4983793  -1.221913 -0.2637127  0.7361520 -0.2639628
## 5 -1.38 -0.2824173 0.4655872  -1.182776 -0.2874476  0.6457992 -0.2876609
## 6 -1.37 -0.2843055 0.4631497  -1.179958 -0.2892976  0.6389993 -0.2895082
##         kld
## 1 0.1462345
## 2 0.3956052
## 3 0.4834661
## 4 0.6523616
## 5 0.8194070
## 6 0.8334976
matplot(result$inla$summary.random$logPopDens[,'ID'],
result$inla$summary.random$logPopDens[,4:6], type='l',
col='black', lty=c(2,1,2))
plot of chunk popdensvar
plot of chunk popdensvar
theCol = mapmisc::colourScale(
    result$data$random.0.5quant,
    breaks=6, digits=1, col='Reds')
plot(peruForModel3, col=theCol$plot, lty=0)
mapmisc::legendBreaks('left', theCol, inset=0.1)
plot of chunk plotSomeResults
plot of chunk plotSomeResults

Spatio-temporal!!!

Here’s some code which will help you get started

peruForModel3 = peruForModel3[!duplicated(peruForModel3$name), ]
peruSt = merge(
    peruHospDef,
    values(peruForModel3)[,c("name","pop","logPop","logPopDens")],
    all=FALSE)


peruNb = terra::adjacent(peruForModel3, 
  type='intersects')
attributes(peruNb)$region.id = peruForModel3$name


peruSt$y = peruSt$hospitalizados_men5 + 
    peruSt$hospitalizados_may5 

stuff = diseasemapping::bym(
    y ~ offset(logPop) + logPopDens + ano,
    # need to add f(name, model = 'iid', replicate, ...)
    data = peruSt,
    adjMat = peruNb,
    region.id = 'name',
    prior = list(sd=1, propSpatial = 0.5),
    control.mode = list(theta=c(-1, -.36), restart=TRUE),
    verbose=TRUE
    )

stuff$parameters$summary[,3:5]
##              0.025quant    0.5quant  0.975quant
## (Intercept) 47.00397470 48.29970820 49.59536749
## logPopDens   0.24407835  0.30141231  0.35931364
## ano         -0.03087001 -0.03023302 -0.02959602
## sd           1.55061880  1.64673234  1.74978153
## propSpatial  0.50285744  0.58925982  0.67135992
stuff$inla$mode$theta
## Log precision for region.indexS     Logit phi for region.indexS 
##                      -0.9969214                       0.3585528
args(diseasemapping:::stbym)

stres = diseasemapping:::stbym(
    y ~ offset(logPop) + logPopDens + 
       f(ano, model="rw2", prior='pc.prec', param = c(0.01, 0.5)),
    data = peruSt,
    adjMat = peruNb,
    region.id = 'name', time.id='ano',
    prior = list(
        spaceTime = list(
            ar = list(prior='pccor0', param = c(0.25, 0.5)),
            sd = c(0.1, 0.5),
            propSpatial = c(0.5, 0.5)
            ),
        sd=1, propSpatial = 0.5),
    verbose=TRUE
#   control.mode = list(theta=c(-1,0.5,2,-0.6, 2.5), restart=TRUE)
    )
saveRDS(stres, "../data/stres.rds")
stres = readRDS("../data/stres.rds")
stres$parameters$summary[,3:5]
##             0.025quant   0.5quant   0.975quant
## (Intercept) 89.2400662 99.8915244 110.55380254
## logPopDens   0.2688482  0.3256111   0.38309494
## ano         -0.0614316 -0.0561247  -0.05082395
## sd           1.2187576  1.3367843   1.45699581
## propSpatial  0.4180214  0.5505086   0.66552197
stres$parameters$sdSpaceTime$summarySpaceTime[,3:5]
##                                0.025quant  0.5quant 0.975quant
## Spatial frac for regionST       0.5893585 0.6336577  0.6797766
## AR coef for regionST            0.7593312 0.7774910  0.7946385
## Spatial frac for region.indexS  0.4180214 0.5505086  0.6655220
## SD for regionST                 1.6217579 1.6846216  1.7523963
## SD for region.indexS            1.2192984 1.3369373  1.4563604
theW = stres$inla$summary.random$regionST
dimnames(theW)[-1]
## $time
##  [1] "2000" "2001" "2002" "2003" "2004" "2005" "2006" "2007" "2008" "2009"
## [11] "2010" "2011" "2012" "2013" "2014" "2015" "2016" "2017" "2018" "2019"
## [21] "2020" "2021" "2022" "2023"
## 
## $quantile
## [1] "mean"       "sd"         "0.025quant" "0.5quant"   "0.975quant"
## [6] "mode"       "kld"
# 1 means the 'bym' part, 2 means posterior mean
theWmean = theW[,,"0.5quant"]
theBreaks = c(0,0.5, 1, 2, 10000)
theBreaks = sort(unique(c(theBreaks, -theBreaks)))
peruR = project(peruForModel3, mapmisc::omerc(peruForModel3, angle=60))
peruRcrop = crop(peruR, ext(-1e7, 1e7, -1e7, 4e5))
apply(as.matrix(ext(peruRcrop)),1,diff)
## [1] 1996956.4  823014.8
for(Dtime in 1:ncol(theWmean)) {
    thecol = mapmisc::colourScale(
        theWmean[,Dtime],
        breaks = theBreaks, style='fixed', col='Reds')
    mapmisc::map.new(peruRcrop)
    plot(peruR, col=thecol$plot, lwd=0.1, border='#FFFFFF20', add=TRUE)
    mtext(colnames(theWmean)[Dtime], line=-1, adj=0.9)
}
mapmisc::scaleBar(peruR, 'top', bty='n')
mapmisc::legendBreaks('topright', thecol, bty='n')

WWWWWWWWWWWWWWWWWWWWWWWW