LGCPs - Multiple Likelihoods
Fabian E. Bachl
Generated on 2024-04-26
Source:vignettes/articles/2d_lgcp_multilikelihood.Rmd
2d_lgcp_multilikelihood.RmdIntroduction
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.
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, alpha = 0) +
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 / 2Setting 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")
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$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.10.1.9004
#> INLA version: 24.04.25-1
#> Components:
#> Common: main = spde(coordinates), group = exchangeable(1L), replicate = iid(1L)
#> Difference: main = spde(coordinates), group = exchangeable(1L), replicate = iid(1L)
#> Intercept: main = linear(1), group = exchangeable(1L), replicate = iid(1L)
#> 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.653, Running = 22.2, Post = 0.141, Total = 23
#> Fixed effects:
#> mean sd 0.025quant 0.5quant 0.975quant mode kld
#> Intercept -0.345 1.335 -2.961 -0.345 2.272 -0.345 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.939 0.577 2.019 2.868 4.278 2.703
#> Stdev for Common 2.049 0.336 1.495 2.014 2.813 1.928
#> Range for Difference 1.968 1.958 0.372 1.397 7.132 0.800
#> Stdev for Difference 0.157 0.078 0.046 0.142 0.345 0.113
#>
#> Deviance Information Criterion (DIC) ...............: 591.68
#> Deviance Information Criterion (DIC, saturated) ....: NA
#> Effective number of parameters .....................: -772.88
#>
#> Watanabe-Akaike information criterion (WAIC) ...: 1630.27
#> Effective number of parameters .................: 105.59
#>
#> Marginal log-Likelihood: -1204.58
#> 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.9004
#> INLA version: 24.04.25-1
#> Components:
#> Common: main = spde(coordinates), group = exchangeable(1L), replicate = iid(1L)
#> Difference: main = spde(coordinates), group = exchangeable(1L), replicate = iid(1L)
#> Intercept: main = linear(1), group = exchangeable(1L), replicate = iid(1L)
#> 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.514, Running = 9.46, Post = 0.118, Total = 10.1
#> Fixed effects:
#> mean sd 0.025quant 0.5quant 0.975quant mode kld
#> Intercept -0.354 1.365 -3.03 -0.354 2.322 -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.935 0.573 2.000 2.871 4.246 2.732
#> Stdev for Common 2.047 0.332 1.487 2.015 2.791 1.944
#> Range for Difference 2.044 2.798 0.140 1.197 9.148 0.373
#> Stdev for Difference 0.158 0.093 0.036 0.139 0.387 0.096
#>
#> Deviance Information Criterion (DIC) ...............: 591.16
#> Deviance Information Criterion (DIC, saturated) ....: NA
#> Effective number of parameters .....................: -773.19
#>
#> Watanabe-Akaike information criterion (WAIC) ...: 1629.67
#> Effective number of parameters .................: 105.28
#>
#> Marginal log-Likelihood: -1204.05
#> 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 has to be
numeric, so we convert the group variable to a 0/1 variable,
group_major <- group == "major".
fml.joint <- coordinates + group_major ~ Intercept + Common + (group_major - 0.5) * Difference
gorillas$nests$group_major <- gorillas$nests$group == "major"
lik_joint <- like("cp",
formula = fml.joint,
data = gorillas$nests,
samplers = gorillas$boundary,
domain = list(
coordinates = gorillas$mesh,
group_major = c(0, 1)
)
)
jfit_joint <- bru(cmp, lik_joint,
options = list(
control.inla = list(int.strategy = "eb"),
bru_max_iter = 2
)
)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$mesh,
mask = gorillas$boundary,
col = exp(jfit_joint$summary.random$Common$mean -
jfit$summary.random$Common$mean)
)
pl.minor <- ggplot() +
gg(gorillas$mesh,
mask = gorillas$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")