knitr::opts_chunk$set( collapse = TRUE, comment = "#>", dev = "png", dev.args = list(type = "cairo-png"), fig.width = 7, fig.height = 5, cache = TRUE )
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.
Load libraries
library(inlabru) library(INLA) library(ggplot2) bru_safe_sp(force = TRUE)
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")
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")
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)
summary(jfit)
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")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.