2d_lgcp_spatiotemporal.Rmd

---
title: "LGCPs - An example in space and time"
author: "Fabian E. Bachl and Finn Lindgren"
date: "Generated on `r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{LGCPs - An example in space and time}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  dev = "png",
  dev.args = list(type = "cairo-png"),
  fig.width = 7,
  fig.height = 5
)
```

Introduction
----------------
For this vignette we are going to be working with a dataset obtained from 
the `R` package `MRSea`. We will set up a LGCP with a spatio-temporal SPDE model to estimate species distribution.


Setting things up
----------------

Load libraries
```{r results="hide",warning=FALSE,message=FALSE}
library(inlabru)
library(INLA)
library(ggplot2)
```

Get the data
-----------------------------------
Load the dataset, that has coordinates in UTM in kilometres:

```{r }
data(mrsea, package = "inlabru")
```

The points (representing animals) and the sampling regions of this dataset are associated with a season. Let's have a look at the observed points and sampling regions for all seasons:

```{r results="hide",warning=FALSE,message=FALSE}
ggplot() +
  gg(mrsea$mesh) +
  gg(mrsea$boundary) +
  gg(mrsea$samplers) +
  gg(mrsea$points, size = 0.5) +
  facet_wrap(~season) +
  ggtitle("MRSea observation seasons")
```

Integration points
-----------------------------------

The `inlabru` point process model knows how to construct the numerical integration
scheme for the LGCP likelihood.  We can also call the internal functions directly to see
what the integration scheme will look like.  Because our model will take time (season)
into account we have to construct the integration points for the LGCP accordingly.
Using the `fm_int()` we can specify the product domain space. Note that omitting
this step would simply aggregate all sampling regions over time.

```{r results="hide",warning=FALSE,message=FALSE,echo=TRUE}
ips <- fm_int(
  domain = list(coordinates = mrsea$mesh, season = 1:4),
  samplers = mrsea$samplers
)
```
Plot the integration points:
```{r results="hide",warning=FALSE,message=FALSE,echo=TRUE}
ggplot() +
  gg(mrsea$mesh) +
  gg(ips, aes(size = weight)) +
  scale_size_area(max_size = 1) +
  facet_wrap(~season)
```

Fitting the model
-----------------------------------

Fit an LGCP model to the locations of the animals. In this example we will employ a spatio-temporal SPDE. Note how the `group` and `ngroup` parameters are employed to let the SPDE model know about the name of the time dimension (season) and the total number of distinct points in time. 

```{r results="hide",warning=FALSE,message=FALSE,echo=TRUE}
matern <- inla.spde2.pcmatern(mrsea$mesh,
  prior.sigma = c(0.1, 0.01),
  prior.range = c(10, 0.01)
)

cmp <- coordinates + season ~ Intercept(1) +
  mySmooth(
    coordinates,
    model = matern,
    group = season,
    ngroup = 4
  )

fit <- lgcp(cmp,
  data = mrsea$points,
  samplers = mrsea$samplers,
  domain = list(
    coordinates = mrsea$mesh,
    season = seq_len(4)
  )
)
```

Predict and plot the intensity for all seasons:

```{r results="hide",warning=FALSE,message=FALSE,echo=TRUE}
ppxl <- fm_pixels(mrsea$mesh, mask = mrsea$boundary, format = "sp")
ppxl_all <- fm_cprod(ppxl, data.frame(season = seq_len(4)))

lambda1 <- predict(
  fit,
  ppxl_all,
  ~ data.frame(season = season, lambda = exp(mySmooth + Intercept))
)
```
```{r results="hide",warning=FALSE,message=FALSE,echo=TRUE}
pl1 <- ggplot() +
  gg(as(lambda1, "SpatialPixelsDataFrame"), aes(fill = q0.5)) +
  gg(mrsea$points, size = 0.3) +
  facet_wrap(~season) +
  coord_equal()
pl1
```