# LGCPs - Multiple Likelihoods

#### Fabian E. Bachl

#### Generated on 2024-05-23

Source:`vignettes/articles/2d_lgcp_multilikelihood.Rmd`

`2d_lgcp_multilikelihood.Rmd`

## Introduction

For this vignette we are going to be working with the inlabru’s
´gorillas_sf´ 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.

## 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_sf, package = "inlabru")`

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

## 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_sf$mesh,
prior.range = c(0.1, 0.01),
prior.sigma = c(1, 0.01)
)
cmp <- ~
Common(geometry, model = matern) +
Difference(geometry, 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 <- geometry ~ Intercept + Common + Difference / 2
fml.minor <- geometry ~ 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_sf$nests[gorillas_sf$nests$group == "major", ],
samplers = gorillas_sf$boundary,
domain = list(geometry = gorillas_sf$mesh)
)
lik_major <- like("cp",
formula = fml.minor,
data = gorillas_sf$nests[gorillas_sf$nests$group == "minor", ],
samplers = gorillas_sf$boundary,
domain = list(geometry = gorillas_sf$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_sf$mesh,
mask = gorillas_sf$boundary,
col = exp(jfit$summary.random$Common$mean)
)
pl.minor <- ggplot() +
gg(gorillas_sf$mesh,
mask = gorillas_sf$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")
```

## Rerunning

Rerunning with the previous estimate as starting point sometimes improves the accuracy of the posterior distribution estimation.

```
jfit0 <- jfit
jfit <- bru_rerun(jfit)
```

```
library(patchwork)
pl.major <- ggplot() +
gg(gorillas_sf$mesh,
mask = gorillas_sf$boundary,
col = exp(jfit$summary.random$Common$mean)
)
pl.minor <- ggplot() +
gg(gorillas_sf$mesh,
mask = gorillas_sf$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.10.1.9007
#> INLA version: 24.05.22
#> Components:
#> Common: main = spde(geometry), group = exchangeable(1L), replicate = iid(1L)
#> Difference: main = spde(geometry), group = exchangeable(1L), replicate = iid(1L)
#> Intercept: main = linear(1), group = exchangeable(1L), replicate = iid(1L)
#> Likelihoods:
#> Family: 'cp'
#> Data class: 'sf', 'data.frame'
#> Predictor: geometry ~ Intercept + Common - Difference/2
#> Family: 'cp'
#> Data class: 'sf', 'data.frame'
#> Predictor: geometry ~ Intercept + Common + Difference/2
#> Time used:
#> Pre = 0.667, Running = 18.3, Post = 0.186, Total = 19.2
#> Fixed effects:
#> mean sd 0.025quant 0.5quant 0.975quant mode kld
#> Intercept -0.354 1.365 -3.029 -0.354 2.32 -0.354 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.936 0.572 2.001 2.871 4.246 2.731
#> Stdev for Common 2.047 0.332 1.487 2.015 2.791 1.944
#> Range for Difference 1.950 2.274 0.192 1.267 7.878 0.513
#> Stdev for Difference 0.159 0.085 0.038 0.144 0.356 0.106
#>
#> Deviance Information Criterion (DIC) ...............: 593.25
#> Deviance Information Criterion (DIC, saturated) ....: 591.61
#> Effective number of parameters .....................: -771.62
#>
#> Watanabe-Akaike information criterion (WAIC) ...: 1631.56
#> Effective number of parameters .................: 106.35
#>
#> Marginal log-Likelihood: -1204.29
#> 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.10.1.9007
#> INLA version: 24.05.22
#> Components:
#> Common: main = spde(geometry), group = exchangeable(1L), replicate = iid(1L)
#> Difference: main = spde(geometry), group = exchangeable(1L), replicate = iid(1L)
#> Intercept: main = linear(1), group = exchangeable(1L), replicate = iid(1L)
#> Likelihoods:
#> Family: 'cp'
#> Data class: 'sf', 'data.frame'
#> Predictor: geometry ~ Intercept + Common - Difference/2
#> Family: 'cp'
#> Data class: 'sf', 'data.frame'
#> Predictor: geometry ~ Intercept + Common + Difference/2
#> Time used:
#> Pre = 0.518, Running = 7.92, Post = 0.129, Total = 8.56
#> Fixed effects:
#> mean sd 0.025quant 0.5quant 0.975quant mode kld
#> Intercept -0.355 1.367 -3.033 -0.355 2.324 -0.355 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.935 0.573 1.999 2.871 4.244 2.732
#> Stdev for Common 2.047 0.332 1.486 2.015 2.791 1.945
#> Range for Difference 2.091 2.889 0.140 1.217 9.413 0.371
#> Stdev for Difference 0.159 0.093 0.036 0.139 0.387 0.097
#>
#> Deviance Information Criterion (DIC) ...............: 591.17
#> Deviance Information Criterion (DIC, saturated) ....: 589.53
#> Effective number of parameters .....................: -773.17
#>
#> Watanabe-Akaike information criterion (WAIC) ...: 1629.64
#> Effective number of parameters .................: 105.26
#>
#> Marginal log-Likelihood: -1204.04
#> 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)')
```

## Single-likelihood version

In this particular model, we can also rewrite the problem as a single
point process over a product domain over space and `group`

.
In versions `<= 2.7.0`

, the integration domain had to be
numeric, so we convert the group variable to a 0/1 variable,
`group_major <- group == "major"`

, which is also useful in
the predictor expression:

`fml.joint <- geometry + group_major ~ Intercept + Common + (group_major - 0.5) * Difference`

```
gorillas_sf$nests$group_major <- gorillas_sf$nests$group == "major"
lik_joint <- like("cp",
formula = fml.joint,
data = gorillas_sf$nests,
samplers = gorillas_sf$boundary,
domain = list(
geometry = gorillas_sf$mesh,
group_major = c(0, 1)
)
)
```

```
jfit_joint <- bru(cmp, lik_joint,
options = list(
control.inla = list(int.strategy = "eb") # Approximate for faster vignette
)
)
```

Plotting the ratios of `exp(Common)`

and
`exp(Difference)`

between the new fit and the old confirms
that the results are the same up to small numerical differences.

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