---
title: Geostatistics practical
---


Data from the Million Deaths Study in India can be found as `india.rds` in the `data` folder of the course web site.

```{r theData}
library('terra')
x = readRDS('../data/india.rds')
values(x)[100000+1:2,] 
```

- `srs_id`: geographic identifier of sampling unit (village or census area)
- `yearCut` and `year`, year of survey, the former is a factor the latter numeric
- `ageCut` and `age`, age group, the former is a factor
- `sex` like it says
- `pop` population of age-sex group in sampling unit.  It's an estimate (and not always an integer)
- `Tuburculosis`, `Venemous.snakes` number of deaths in the relevant age-sex group
- `urban_rural`: area type of the sampling unit
- `languageGroup`:broad category of dominant language in the sampling unit


```{r amap}
xUnique = x[!duplicated(x$srs_id),]
india = geodata::gadm("IND", level=0, path=normalizePath("~/Downloads"))
indiaT = project(india, crs(xUnique))
plot(indiaT)
points(xUnique, col='red', cex=0.5)
```

```{r covariate}
indiaElev = geodata::elevation_30s('IND', path=normalizePath("~/Downloads"))
plot(indiaElev)
```

Write down the model!
$$
\begin{aligned}
Y_{itj} \sim & \text{Poisson}(\lambda_{itj}) \\
\end{aligned}
$$

```{r glgm}
xSub = x[which(x$pop > 0 & x$age >= 40),]
xSub$logPop = log(xSub$pop)
res = geostatsp::glgm(
	Tuberculosis ~ sex + ageCut + yearCut + elevation + 
		offset(logPop) +
		f(srs_id, model='iid', prior='pc.prec', param=c(1, 0.5)),
	covariates = list(elevation = indiaElev),
	data = xSub,
	family='poisson',
	grid = 50,
	prior = list(sd = 1, range = 500*1000)
	)
```

```{r glgmRes}
knitr::kable(res$parameters$summary[,3:5], digits=2)
```


```{r plotFit}
toPlot = mask(res$raster[['predict.0.5quant']], indiaT)
toPlot = exp(toPlot)*100000
thecol = mapmisc::colourScale(toPlot, col='Reds', breaks=6, dec=-2, style='equal')
mapmisc::map.new(indiaT)
plot(toPlot, col=thecol$col, breaks=thecol$breaks, legend=FALSE, add=TRUE)
plot(indiaT, add=TRUE)
mapmisc::legendBreaks("bottomleft", thecol, bty='n')
```

caluculate observed and expected

```{r obsexp}
xHavePop = x[which(x$pop>0),]
xHavePop$logPop = log(xHavePop$pop)
xHavePop$ageRelevel = relevel(xHavePop$ageCut, '[60,65)')
theRates = glm(Tuberculosis ~ ageRelevel*sex + offset(logPop),family='poisson', data=values(xHavePop))
xHavePop$pred = predict.glm(theRates, values(xHavePop), type='response')

sum(xHavePop$pred)
sum(xHavePop$Tuberculosis)

byvars = c("yearCut","srs_id",'urban_rural',"year",'languageGroup')
xAggDf = aggregate(
	values(xHavePop)[,c('pred','pop','Tuberculosis','Venomous.snakes')], 
	values(xHavePop)[,byvars],
	FUN='sum')
xUnique = xHavePop[!duplicated(values(xHavePop)[, byvars]), 
	byvars]	
xAgg = merge(xUnique, xAggDf)

values(xAgg)[24060+1:4,]
```

```{r glgm2}
xAgg$logExpected = log(xAgg$pred)
res2 = geostatsp::glgm(
	Tuberculosis ~  yearCut + elevation + 
		offset(logExpected) +
		f(srs_id, model='iid', prior='pc.prec', param=c(1, 0.5)),
	covariates = list(elevation = indiaElev/1000),
	data = xAgg,
	family='poisson',
	grid = geostatsp::squareRaster(indiaT,150), buffer=200*1000,
	prior = list(sd = 1, range = 500*1000)
	)
knitr::kable(res2$parameters$summary[,3:5], digits=3)
mapmisc::map.new(indiaT)
plot(res2$raster$random.0.5quant, add=TRUE)
plot(indiaT, add=TRUE)
```


Your Task: do something interesting with either tuberculosis or snake bites.