Modelo_espacial
Delia Ortega
Preparación de los Datos
Casos de cancer de labio observados y esperados y la proporción de población que tiene actividades relacionadas con la agricultura o pesca (AFF). Para calcular los valores esperados se debe tener la población y los casos por cada estrato (grupo).
library(SpatialEpi)
data(scotland)
class(scotland)## [1] "list"names(scotland)## [1] "geo" "data" "spatial.polygon" "polygon"geo es un dataframe con los nombres y coordenadas de los centroides de cada condado y el objeto spatial.polygon es el mapa de Escocia, mientras data contiene los datos del número de casos y la covariable.
head(scotland$data)## county.names cases expected AFF
## 1 skye-lochalsh 9 1.4 0.16
## 2 banff-buchan 39 8.7 0.16
## 3 caithness 11 3.0 0.10
## 4 berwickshire 9 2.5 0.24
## 5 ross-cromarty 15 4.3 0.10
## 6 orkney 8 2.4 0.24map <- scotland$spatial.polygon
plot(map)
El mapa no contiene información sobre el sistema de coordenadas de referencia.
codes <- rgdal::make_EPSG()
codes[which(codes$code == "27700"), ]## code note
## 2475 27700 OSGB 1936 / British National Grid
## prj4
## 2475 +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs +type=crs
## prj_method
## 2475 Transverse Mercatorproj4string(map) <- "+proj=tmerc +lat_0=49 +lon_0=-2
+k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36
+units=km +no_defs"Se transforma el mapa usando la proyección WGS84 para tener los datos de longitud y latitud
map <- spTransform(map,
CRS("+proj=longlat +datum=WGS84 +no_defs"))Se realiza la selección de los datos y se asignan en un dataframe d:
d <- scotland$data[,c("county.names", "cases", "expected", "AFF")]
names(d) <- c("county", "Y", "E", "AFF")cálculo de los SIR (razón de incidencia estandarizada)
d$SIR <- d$Y / d$E
head(d)## county Y E AFF SIR
## 1 skye-lochalsh 9 1.4 0.16 6.428571
## 2 banff-buchan 39 8.7 0.16 4.482759
## 3 caithness 11 3.0 0.10 3.666667
## 4 berwickshire 9 2.5 0.24 3.600000
## 5 ross-cromarty 15 4.3 0.10 3.488372
## 6 orkney 8 2.4 0.24 3.333333Usando el objeto \(spatialpolygons\) y el dataframe \(d\) se genera otro objeto que contenga la variable creada SIR usando ID que corresponde al nombre del condado.
library(sp)
rownames(d) <- d$county
map <- SpatialPolygonsDataFrame(map, d, match.ID = TRUE)
head(map@data)## county Y E AFF SIR
## skye-lochalsh skye-lochalsh 9 1.4 0.16 6.428571
## banff-buchan banff-buchan 39 8.7 0.16 4.482759
## caithness caithness 11 3.0 0.10 3.666667
## berwickshire berwickshire 9 2.5 0.24 3.600000
## ross-cromarty ross-cromarty 15 4.3 0.10 3.488372
## orkney orkney 8 2.4 0.24 3.333333Modelo BYM
Se estima un modelo poisson con la siguiente forma:
\[\begin{align*} Y_{i} \sim \text{Poisson}(E_{i}\theta_{i}), i=1,2,...,n\\ \log(\theta_{i})=x_{i}\beta+u_{i}+v_{i}\\ \mu_{i}|\mu_{-i}\sim N\left(\overline{\mu}_{\delta_{i}},\frac{\sigma_\mu^2}{n\delta_{i}}\right) \end{align*}\]
Donde \(d_{i}=(1, d_{i1},.....,d_{ip})\) es el vector del intercepto y p covariables del area \(i\) y \(\beta=(\beta_{0}, \beta_{1},....,\beta_{p})\) es el vector de coeficientes. Por cada unidad de incremento de la covariable \(d_{j}\), \(j=1,2,...,p\), el riesgo relativo aumenta por un factor de \(\exp(\beta_{1})\). El término \(u_{i}\) representa el efecto aleatorio específico del area \(i\) que modela la dependencia espacial entre los riesgos relativos. A este componente se asocia una distribución condicional autoregresiva (CAR) que asume una ditribución normal con media igual al promedio del efecto espacial en los vecinos adyacentes \((\overline{\mu}_{\delta_{i}})\) y varianza inversamente proporcional al número de vecinos \((\frac{\sigma_\mu^2}{n\delta_{i}})\). \(v_{i}\) es el error no correlacionado que se asume independiente e identicamente distribuido con media cero y varianza \(\sigma_{v}^2\).
Matriz de vecinos Con la función \(poly2nb\) para crear la lista de vecinos basadas en las areas con fronteras contiguas, cada elemento de nb representa un areay contiene los indices de sus vecinos.
library(spdep)
library(INLA)
nb <- poly2nb(map)
head(nb)## [[1]]
## [1] 5 9 19
##
## [[2]]
## [1] 7 10
##
## [[3]]
## [1] 12
##
## [[4]]
## [1] 18 20 28
##
## [[5]]
## [1] 1 12 19
##
## [[6]]
## [1] 0La función \(nb2INLA()\) convierte la lista en un archivo con la representación de la matriz de vecinos como se requiere en R-INLA.
nb2INLA("map.adj", nb)
g <- inla.read.graph(filename = "map.adj")Usando INLA
Se incluyen dos vectores en los datos para denotar los indices de estos vectores aleatorios \(u_{i}\) y \(v_{i}\).
map$idareau <- 1:nrow(map@data)
map$idareav <- 1:nrow(map@data)el modelo:
formula <- Y ~ AFF +
f(idareau, model = "besag", graph = g, scale.model = TRUE) +
f(idareav, model = "iid")Estimación:
res <- inla(formula,
family = "poisson", data = map@data,
E = E, control.predictor = list(compute = TRUE)
)Resultados:
summary(res)##
## Call:
## c("inla(formula = formula, family = \"poisson\", data = map@data, ", "
## E = E, control.predictor = list(compute = TRUE))")
## Time used:
## Pre = 0.625, Running = 0.717, Post = 0.0628, Total = 1.41
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -0.305 0.120 -0.539 -0.306 -0.069 -0.307 0
## AFF 4.330 1.277 1.744 4.357 6.770 4.409 0
##
## Random effects:
## Name Model
## idareau Besags ICAR model
## idareav IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision for idareau 4.15 1.45 2.02 3.92 7.63 3.49
## Precision for idareav 19416.69 19411.20 1352.54 13675.08 71097.63 3702.39
##
## Expected number of effective parameters(stdev): 28.54(3.53)
## Number of equivalent replicates : 1.96
##
## Marginal log-Likelihood: -189.70
## Posterior marginals for the linear predictor and
## the fitted values are computedDistribución Posterior
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.0.5marginal <- inla.smarginal(res$marginals.fixed$AFF)
marginal <- data.frame(marginal)
ggplot(marginal, aes(x = x, y = y)) + geom_line() +
labs(x = expression(beta[1]), y = "Density") +
geom_vline(xintercept = 0, col = "black") + theme_bw()