Skip to contents

Introduction

For this vignette we are going to be working with the inlabru’s ´gorillas´ dataset which was originally obtained from the R package spatstat. The data set contains two types of gorillas nests which are marked as either major or minor. We will set up a multi-likelihood model for these nests which creates two spatial LGCPs that share a common intercept but have employ different spatial smoothers.

Setting things up

Load libraries

Get the data

For the next few practicals we are going to be working with a dataset obtained from the R package spatstat, which contains the locations of 647 gorilla nests. We load the dataset like this:

data(gorillas, package = "inlabru")

Plot the nests and visualize the group membership (major/minor) by color:

ggplot() +
  gg(gorillas$mesh) +
  gg(gorillas$nests, aes(color = group)) +
  gg(gorillas$boundary) +
  coord_fixed() +
  ggtitle("Gorillas nests and group membership")

Fiting the model

First, we define all components that enter the joint model. That is, the intercept that is common to both LGCPs and the two different spatial smoothers, one for each nest group.

matern <- inla.spde2.pcmatern(gorillas$mesh,
  prior.range = c(0.1, 0.01),
  prior.sigma = c(1, 0.01)
)

cmp <- ~
  Common(coordinates, model = matern) +
    Difference(coordinates, model = matern) +
    Intercept(1)

Given these components we define the linear predictor for each of the likelihoods. (Using “.” indicates a pure additive model, and one can use include/exclude options for like() to indicate which components are actively involved in each model.)

fml.major <- coordinates ~ Intercept + Common + Difference / 2
fml.minor <- coordinates ~ Intercept + Common - Difference / 2

Setting up the Cox process likelihoods is easy in this example. Both nest types were observed within the same window:

lik_minor <- like("cp",
  formula = fml.major,
  data = gorillas$nests[gorillas$nests$group == "major", ],
  samplers = gorillas$boundary,
  domain = list(coordinates = gorillas$mesh)
)
lik_major <- like("cp",
  formula = fml.minor,
  data = gorillas$nests[gorillas$nests$group == "minor", ],
  samplers = gorillas$boundary,
  domain = list(coordinates = gorillas$mesh)
)

… which we provide to the ´bru´ function.

jfit <- bru(cmp, lik_major, lik_minor,
  options = list(
    control.inla = list(int.strategy = "eb"),
    bru_max_iter = 1
  )
)
library(patchwork)
pl.major <- ggplot() +
  gg(gorillas$mesh,
    mask = gorillas$boundary,
    col = exp(jfit$summary.random$Common$mean)
  )
pl.minor <- ggplot() +
  gg(gorillas$mesh,
    mask = gorillas$boundary,
    col = exp(jfit$summary.random$Difference$mean)
  )
(pl.major + scale_fill_continuous(trans = "log")) +
  (pl.minor + scale_fill_continuous(trans = "log")) &
  theme(legend.position = "right")

jfit0 <- jfit
jfit <- bru_rerun(jfit)
library(patchwork)
pl.major <- ggplot() +
  gg(gorillas$mesh,
    mask = gorillas$boundary,
    col = exp(jfit$summary.random$Common$mean)
  )
pl.minor <- ggplot() +
  gg(gorillas$mesh,
    mask = gorillas$boundary,
    col = exp(jfit$summary.random$Difference$mean)
  )
(pl.major + scale_fill_continuous(trans = "log")) +
  (pl.minor + scale_fill_continuous(trans = "log")) &
  theme(legend.position = "right")

summary(jfit0)
#> inlabru version: 2.7.0.9004
#> INLA version: 23.01.31
#> Components:
#> Common: main = spde(coordinates)
#> Difference: main = spde(coordinates)
#> Intercept: main = linear(1)
#> Likelihoods:
#>   Family: 'cp'
#>     Data class: 'SpatialPointsDataFrame'
#>     Predictor: coordinates ~ Intercept + Common - Difference/2
#>   Family: 'cp'
#>     Data class: 'SpatialPointsDataFrame'
#>     Predictor: coordinates ~ Intercept + Common + Difference/2
#> Time used:
#>     Pre = 1.13, Running = 116, Post = 0.118, Total = 117 
#> Fixed effects:
#>             mean    sd 0.025quant 0.5quant 0.975quant   mode kld
#> Intercept -0.349 1.351     -2.998   -0.349      2.299 -0.349   0
#> 
#> Random effects:
#>   Name     Model
#>     Common SPDE2 model
#>    Difference SPDE2 model
#> 
#> Model hyperparameters:
#>                       mean    sd 0.025quant 0.5quant 0.975quant  mode
#> Range for Common     2.989 0.722      1.870     2.89      4.686 2.680
#> Stdev for Common     2.065 0.415      1.390     2.02      3.015 1.914
#> Range for Difference 1.919 1.330      0.628     1.54      5.477 1.070
#> Stdev for Difference 0.141 0.069      0.041     0.13      0.304 0.103
#> 
#> Deviance Information Criterion (DIC) ...............: 589.25
#> Deviance Information Criterion (DIC, saturated) ....: NA
#> Effective number of parameters .....................: -774.86
#> 
#> Watanabe-Akaike information criterion (WAIC) ...: 1628.71
#> Effective number of parameters .................: 104.73
#> 
#> Marginal log-Likelihood:  -1205.61 
#>  is computed 
#> Posterior summaries for the linear predictor and the fitted values are computed
#> (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
summary(jfit)
#> inlabru version: 2.7.0.9004
#> INLA version: 23.01.31
#> Components:
#> Common: main = spde(coordinates)
#> Difference: main = spde(coordinates)
#> Intercept: main = linear(1)
#> Likelihoods:
#>   Family: 'cp'
#>     Data class: 'SpatialPointsDataFrame'
#>     Predictor: coordinates ~ Intercept + Common - Difference/2
#>   Family: 'cp'
#>     Data class: 'SpatialPointsDataFrame'
#>     Predictor: coordinates ~ Intercept + Common + Difference/2
#> Time used:
#>     Pre = 0.937, Running = 22.6, Post = 0.101, Total = 23.6 
#> Fixed effects:
#>            mean    sd 0.025quant 0.5quant 0.975quant  mode kld
#> Intercept -0.35 1.372      -3.04    -0.35      2.339 -0.35   0
#> 
#> Random effects:
#>   Name     Model
#>     Common SPDE2 model
#>    Difference SPDE2 model
#> 
#> Model hyperparameters:
#>                       mean    sd 0.025quant 0.5quant 0.975quant  mode
#> Range for Common     2.983 0.657      1.937    2.899      4.510 2.726
#> Stdev for Common     2.055 0.353      1.466    2.018      2.850 1.939
#> Range for Difference 1.662 1.086      0.501    1.372      4.523 0.984
#> Stdev for Difference 0.148 0.079      0.039    0.134      0.339 0.101
#> 
#> Deviance Information Criterion (DIC) ...............: 587.71
#> Deviance Information Criterion (DIC, saturated) ....: NA
#> Effective number of parameters .....................: -775.97
#> 
#> Watanabe-Akaike information criterion (WAIC) ...: 1627.56
#> Effective number of parameters .................: 104.15
#> 
#> Marginal log-Likelihood:  -1205.67 
#>  is computed 
#> Posterior summaries for the linear predictor and the fitted values are computed
#> (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')