Calculates local and integrated variance and correlation measures as introduced by Yuan et al. (2017).
Arguments
- joint
A joint
predictionof two latent model components.- prediction1
A
predictionof the first component.- prediction2
A
predictionof the second component.- samplers
An
sforSpatialPolygonobject describing the area for which to compute the cumulative variance measure.- mesh
The fmesher::fm_mesh_2d for which the prediction was performed (required for cumulative Vmeasure).
References
Y. Yuan, F. E. Bachl, F. Lindgren, D. L. Brochers, J. B. Illian, S. T. Buckland, H. Rue, T. Gerrodette. 2017. Point process models for spatio-temporal distance sampling data from a large-scale survey of blue whales. https://arxiv.org/abs/1604.06013
Examples
# \donttest{
if (bru_safe_inla() &&
require("ggplot2", quietly = TRUE) &&
require("patchwork", quietly = TRUE) &&
require("sn", quietly = TRUE) &&
require("terra", quietly = TRUE) &&
require("sf", quietly = TRUE) &&
require("RColorBrewer", quietly = TRUE) &&
require("dplyr", quietly = TRUE) &&
require("magrittr", quietly = TRUE)) {
# Load Gorilla data
gorillas <- gorillas_sf
gorillas$gcov <- gorillas_sf_gcov()
# Use RColorBrewer
# Fit a model with two components:
# 1) A spatial smooth SPDE
# 2) A spatial covariate effect (vegetation)
pcmatern <- INLA::inla.spde2.pcmatern(
gorillas$mesh,
prior.sigma = c(0.1, 0.01),
prior.range = c(0.01, 0.01)
)
cmp <- geometry ~ 0 +
vegetation(gorillas$gcov$vegetation, model = "factor_contrast") +
spde(geometry, model = pcmatern) -
Intercept(1)
fit <- lgcp(
cmp,
gorillas$nests,
samplers = gorillas$boundary,
domain = list(geometry = gorillas$mesh),
options = list(control.inla = list(int.strategy = "eb"))
)
# Predict SPDE and vegetation at a grid covering the domain of interest
pred_loc <- fm_pixels(
gorillas$mesh,
mask = gorillas$boundary,
dims = c(200, 200),
format = "sf"
)
pred <- predict(
fit,
pred_loc,
~ list(
joint = spde + vegetation,
field = spde,
veg = vegetation
)
)
pred_collect <-
rbind(
pred$joint %>% dplyr::mutate(component = "joint"),
pred$field %>% dplyr::mutate(component = "field"),
pred$veg %>% dplyr::mutate(component = "veg")
) %>%
dplyr::mutate(var = sd^2)
# Plot component mean
ggplot(pred_collect) +
geom_tile(aes(geometry = geometry, fill = mean), stat = "sf_coordinates") +
coord_equal() +
facet_wrap(~component, nrow = 1) +
theme(legend.position = "bottom")
# Plot component variance
ggplot(pred_collect) +
geom_tile(aes(geometry = geometry, fill = var), stat = "sf_coordinates") +
coord_equal() +
facet_wrap(~component, nrow = 1) +
theme(legend.position = "bottom")
# Calculate variance and correlation measure
vm <- devel.cvmeasure(pred$joint, pred$field, pred$veg)
# Compute nominal relative variance contributions; note that these can be
# greater than 100%!
vm <- dplyr::mutate(
vm,
var1_rel = if_else(var1 <= 0, NA, var1 / var.joint),
var2_rel = if_else(var2 <= 0, NA, var2 / var.joint)
)
lprange <- range(vm$var.joint, vm$var1, vm$var2)
vm <- tidyr::pivot_longer(vm,
cols = c(var.joint, var1, var2, cov, cor, var1_rel, var2_rel),
names_to = "component",
values_to = "value"
)
# Variance contribution of the components
csc <- scale_fill_gradientn(
colours = brewer.pal(9, "YlOrRd"),
limits = lprange
)
vm_ <- dplyr::filter(vm, component %in% c("var.joint", "var1", "var2"))
ggplot(vm_) +
geom_tile(aes(geometry = geometry, fill = value),
stat = "sf_coordinates"
) +
csc +
coord_equal() +
facet_wrap(~component, nrow = 1) +
theme(legend.position = "bottom")
# Relative variance contribution of the components
# When bo
vm_ <- dplyr::filter(vm, component %in% c("var1_rel", "var2_rel"))
ggplot(vm_) +
geom_tile(aes(geometry = geometry, fill = value),
stat = "sf_coordinates"
) +
scale_fill_gradientn(
colours = brewer.pal(9, "YlOrRd"),
trans = "log"
) +
coord_equal() +
facet_wrap(~component, nrow = 1)
# Where both relative contributions are larger than 1, the posterior
# correlations are strongly negative.
# Covariance and correlation of field and vegetation
vm_cov <- dplyr::filter(vm, component %in% "cov")
vm_cor <- dplyr::filter(vm, component %in% "cor")
(ggplot(vm_cov) +
geom_tile(aes(geometry = geometry, fill = value),
stat = "sf_coordinates"
) +
scale_fill_gradientn(
colours = rev(brewer.pal(9, "RdBu")),
limits = c(-1, 1) * max(abs(vm_cov$value), na.rm = TRUE)
) +
coord_equal() +
theme(legend.position = "bottom") +
ggtitle("Covariances") |
ggplot(vm_cor) +
geom_tile(aes(geometry = geometry, fill = value),
stat = "sf_coordinates"
) +
scale_fill_gradientn(
colours = rev(brewer.pal(9, "RdBu")),
limits = c(-1, 1) * max(abs(vm_cor$value), na.rm = TRUE)
) +
coord_equal() +
theme(legend.position = "bottom") +
ggtitle("Correlations")
)
# Variance and correlation integrated over space
vrt <- fm_vertices(gorillas$mesh, format = "sf")
pred_vrt <- predict(
fit,
vrt,
~ list(
joint = spde + vegetation,
field = spde,
veg = vegetation
)
)
vm.int <- devel.cvmeasure(
pred_vrt$joint,
pred_vrt$field,
pred_vrt$veg,
samplers = fm_int(gorillas$mesh, gorillas$boundary),
mesh = gorillas$mesh
)
vm.int
}
#>
#> Attaching package: ‘sn’
#> The following object is masked from ‘package:stats’:
#>
#> sd
#> terra 1.8.42
#>
#> Attaching package: ‘terra’
#> The following object is masked from ‘package:patchwork’:
#>
#> area
#>
#> Attaching package: ‘dplyr’
#> The following objects are masked from ‘package:terra’:
#>
#> intersect, union
#> The following objects are masked from ‘package:stats’:
#>
#> filter, lag
#> The following objects are masked from ‘package:base’:
#>
#> intersect, setdiff, setequal, union
#>
#> Attaching package: ‘magrittr’
#> The following objects are masked from ‘package:terra’:
#>
#> extract, inset
#> var.joint var1 var2 cor
#> 1 0.528323 0.7418273 0.2899805 -0.5427759
# }