spatialfusion: short demo"

In spatialfusion: Multivariate Analysis of Spatial Data Using a Unifying Spatial Fusion Framework

knitr::opts_chunk$set(echo = TRUE, cache = TRUE, fig.align = "center")

This brief demo provides code and output for fitting spatial fusion models using R package spatialfusion. The first section analyze the built-in synthetic dataset with INLA implementation while the second section analyze a simulated dataset with Stan implementation. The method argument in fusionData() function decides on which implementation to use.

1. Spatial fusion modelling with INLA on built-in synthetic data

Load libraries

library(spatialfusion)
library(tmap, quietly = T)
library(sp, quietly = T)

Load and view built-in synthetic data

summary(dataGeo)
summary(dataLattice)
summary(dataPP)

Plot data

tm_shape(dataLattice) + tm_polygons(col = "white") +
  tm_shape(dataGeo) + tm_dots(size = 0.1) +
  tm_add_legend(type = "symbol", shape = 16, size = 0.3, col = "black", label = "geostatistical") +
  tm_shape(dataPP) + tm_symbols(col = "red", shape = 4, size = 0.02) +  
  tm_add_legend(type = "symbol", shape = 4, size = 0.2, col = "red", label = "point pattern") +
  tm_layout(main.title = "dataGeo, dataPP", main.title.size = 1, 
            frame = F, fontface = 2, legend.outside = T)

tm_shape(dataLattice) +
  tm_fill(col = "mr", style = "fixed", breaks = c(0, 0.5, 2, 10, 30),
          title = "Mortality rate \n(per thousand)", legend.reverse = T) + tm_borders() +
  tm_layout(main.title = "dataLattice", main.title.size = 1, frame = F, fontface = 2,
            legend.position = c(0.77,0.8), legend.text.size = 0.5, legend.title.size = 0.5)

Data preparation

dat <- fusionData(geo.data = dataGeo, geo.formula = lungfunction ~ covariate,
                  lattice.data = dataLattice,
                  lattice.formula = mortality ~ covariate + log(pop),
                  pp.data = dataPP, distributions = c("normal", "poisson"),
                  method = "INLA")

dat

Fit a spatial fusion model

mod <- fusion(data = dat, n.latent = 1, bans = matrix(c(0,0,0), ncol = 1),
              pp.offset = 400, prior.range = c(0.1, 0.5),
              prior.sigma = c(1, 0.5),  mesh.locs = dat$locs_point,
              mesh.max.edge = c(0.05, 0.5))

Inspect the fit

mod_fit <- fitted(mod, type = "link")
par(mfrow = c(1,2))
plot(dataGeo$lungfunction, mod_fit$point1,
     xlab = "Observed", ylab = "Fitted")
abline(0,1)
plot(log(dataLattice$mortality), mod_fit$area1,
     xlab = "Observed", ylab = "Fitted")
abline(0,1)

Check parameter estimates

summary(mod, digits = 3)

Diagnostic plots

plot(mod, interactive = FALSE)

plot(mod, posterior = FALSE)

Predict latent surface

pred.locs <- spsample(dataDomain, 20000, type = "regular")
mod.pred <- predict(mod, pred.locs)
mod.pred.plot <- SpatialPointsDataFrame(coords = pred.locs, data = as.data.frame(mod.pred))
tm_shape(mod.pred.plot) +
  tm_symbols(col = "latent.s11", shape = 15, size = 0.05, style = "cont",
             midpoint = NA, legend.col.reverse = T, palette = "Oranges",
             title.col = "Latent process") +
  tm_shape(dataLattice) + tm_borders() +
  tm_layout(frame = FALSE, legend.outside = TRUE)

2. Spatial fusion modelling with Stan on simulated data

Simulate data

dat <- fusionSimulate(n.point = 200, n.area = 30, n.grid = 5, n.pred = 100,
                      psill = 1.5, phi = 1, nugget = 0, tau.sq = 0.2,
                      dimension = 10, domain = NULL, point.beta = list(rbind(1,5)),
                      area.beta = list(rbind(1, 1.5)), nvar.pp = 1,
                      distributions = c("normal","poisson"),
                      design.mat = matrix(c(2, 0.5, 1), ncol = 1), 
                      pp.offset = 0.5, seed = 1)

geo.data <- SpatialPointsDataFrame(coords = dat$mrf[dat$sample.ind, c("x","y")],
                                   data = data.frame(cov.point = dat$dat$X_point[,2], 
                                                     outcome = dat$dat$Y_point[[1]]),
                                   proj4string = CRS("+proj=longlat +ellps=WGS84"))

lattice.data <- SpatialPolygonsDataFrame(dat$poly,
                                         data = data.frame(outcome = dat$dat$Y_area[[1]], 
                                                           cov.area = dat$dat$X_area[,2]))

pp.data <- dat$data$lgcp.coords[[1]]

lattice.data@proj4string <- pp.data@proj4string <- CRS("+proj=longlat +ellps=WGS84")

Plot data

tm_shape(lattice.data) + tm_polygons(col = "white") + 
  tm_shape(geo.data) + tm_dots(size = 0.1) +
  tm_add_legend(type = "symbol", shape = 16, size = 0.3, col = "black", label = "geostatistical") +
  tm_shape(pp.data) + tm_symbols(col = "red", shape = 4, size = 0.02) +
  tm_add_legend(type = "symbol", shape = 4, size = 0.2, col = "red", label = "point pattern") +
  tm_layout(main.title = "geo.data, pp.data", main.title.size = 1, 
            frame = F, fontface = 2, legend.outside = T)

tm_shape(lattice.data) + 
  tm_fill(col="outcome", style="fixed", breaks=c(0, 10, 50, 200, 1200),
          title = "Outcome", legend.reverse = T) + tm_borders() +
  tm_layout(main.title="lattice.data", main.title.size = 1, frame = F, fontface = 2,
            legend.position = c(0.77,0.8), legend.text.size = 0.5, legend.title.size = 0.5)

Data preparation

dat2 <- fusionData(geo.data = geo.data, geo.formula = outcome ~ cov.point,
                   lattice.data = lattice.data, lattice.formula = outcome ~ cov.area,
                   pp.data = pp.data, distributions = c("normal","poisson"), 
                   method = "Stan")
dat2

Fit a spatial fusion model

mod <- fusion(data = dat2, n.latent = 1, bans = matrix(c(0,0,0), ncol = 1),
              pp.offset = 0.5, prior.phi = list(distr = "normal", pars = c(1, 1)))

Inspect the fit

mod_fit <- fitted(mod, type = "link")
par(mfrow = c(1,2))
plot(dat$data$Y_point[[1]], mod_fit$point1, xlab = "Observed", ylab = "Fitted")
abline(0,1)
plot(log(dat$data$Y_area[[1]]), mod_fit$area1, xlab = "Observed", ylab = "Fitted")
abline(0,1)

Check parameter estimates

summary(mod, digits = 2) 

Diagnostic plots

plot(mod, interactive = FALSE)

plot(mod, posterior = FALSE) 

Predict latent surface and compare with simulated truth

mod.pred <- predict(mod, new.locs = dat$pred.loc, type = "summary")
par(mfrow = c(1,1))
plot(dat$mrf[dat$pred.ind, c("sim1")], mod.pred$latent1, 
     xlab = "Truth", ylab = "Predicted")
abline(0,1)

spatialfusion documentation built on Aug. 23, 2022, 1:05 a.m.