---
title: infections in peru
---

- Data from [peru](https://datosabiertos.gob.pe/dataset/vigilancia-epidemiol%C3%B3gica-de-enfermedad-diarreica-aguda-eda/resource/d7233b43-865c-4ba3-b126)
- counts of reportable infections, hopsitalizations and deaths
- data `peru.RData` on course web site, data folder
- need population!
- I can't find peruc ensus data, so we'll aggregate gridded data from [worldpop](https://hub.worldpop.org/geodata/summary?id=50047)
- census data at https://datosabiertos.gob.pe/dataset/población-peru

```{r peruPop}
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)
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)
```

Aggregate:

```{r agg}
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))
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)
sum(values(pop), na.rm=TRUE)
```


we'll do all hospitalizations, years 2012-2015

```{r peruForBym}
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
```{r 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)
```

```{r findNeighbours}
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(peruForModel2[which(Nneighbours==0),],  col='green')
peruForModel3 = peruForModel2[Nneighbours > 0, ]
peruForModel3$y[325] = NA
```

```{r computePopDens}
peruForModel3$area = terra::expanse(peruForModel3)/1000000
peruForModel3$popdens = peruForModel3$pop/peruForModel3$area
peruForModel3$logPopDens = round(log(peruForModel3$popdens), 2)
```

```{r bym, eval=TRUE}
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
	)
```

```{r summary}
knitr::kable(result$parameters$summary[,3:5], digits=3)
100*1000*exp(result$parameters$summary[1,grep("quant", names(result$parameters$summary))])
```



```{r data}
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')]
```

```{r popdensvar}
head(result$inla$summary.random$logPopDens)
matplot(result$inla$summary.random$logPopDens[,'ID'],
result$inla$summary.random$logPopDens[,4:6], type='l',
col='black', lty=c(2,1,2))
```

```{r plotSomeResults}
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)	
```



# Spatio-temporal!!!


Here's some code which will help you get started

```{r dataForSt}

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]
stuff$inla$mode$theta
```

- Fit a better spatio-temporal additive model
- then fit a full spatio-temporal interaction model!

```{r stbym, eval=FALSE}

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")
```
```{r stRes}
stres = readRDS("../data/stres.rds")
stres$parameters$summary[,3:5]
stres$parameters$sdSpaceTime$summarySpaceTime[,3:5]
```

```{r plotStbymsetup}
theW = stres$inla$summary.random$regionST
dimnames(theW)[-1]

# 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)
```

```{r plotStBym, fig.cap='W', fig.subcap=colnames(theWmean), fig.height=3, fig.width=7}
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')
```
