Skip to contents

Extract block-averaged GCPO scores from one or more fitted bru models

Source: R/gcpo.R
bru_block_gcpo.Rd

After fitting a model with control.gcpo = list(enable = TRUE), this function reads the raw per-observation GCPO vector from fit$gcpo$gcpo and returns one score per block or observation group.

Two cases are handled automatically:

Block-based (lgcp / "cp" family)

When BRU_block is present in response_data (populated automatically by the "cp" family machinery from the samplers structure), INLA repeats the same GCPO value for every observation in a leave-out block. The mean within each block collapses the repeated values to one score per block.

Groups-based (all other families)

When BRU_block is absent, INLA builds groups internally (via num.level.sets or friends) or from a user-supplied groups list. One score per observation is returned directly from fit$gcpo$gcpo.

Scores are returned on the probability scale, consistent with fit$cpo$cpo. For multi-likelihood models, cumulative row offsets into the stacked INLA response vector are handled automatically.

Usage

bru_block_gcpo(fit, ...)

Arguments

fit

A fitted object of class bru, or a named list of such objects, or several such objects passed via \dots. Each must have been fitted with control.gcpo = list(enable = TRUE) so that fit$gcpo$gcpo is non-NULL.

...

Additional named fitted bru objects.

Value

For a single fit, a list with two elements:

blocks

Named list with one element per likelihood, each a vector of block labels (for "cp" models: unique BRU_block values; for other models: observation indices).

gcpo

GCPO scores on the probability scale, consistent with fit$cpo$cpo. A numeric vector for single-likelihood models; a named list of numeric vectors for multi-likelihood models.

For multiple fits, a named list of such objects, one per fit.

See also

bru(), bru_obs(), bru_obs_control_gcpo(), bru_gcpo_table()

Examples

# \donttest{
if (bru_safe_inla()) {
  # Block-based example (lgcp)
  cvpart <- cv_hex(
    gorillas_sf$boundary,
    cellsize = 0.5,
    n_group = 3,
    resolution = c(95, 80)
  )
  cvpart$block_ID <- seq_len(nrow(cvpart))
  cvpart$group <- NULL

  nests <- gorillas_sf$nests
  a <- sf::st_intersects(nests, cvpart)
  nests$.block <- unlist(a)

  fit1 <- lgcp(
    geometry ~ Intercept(1),
    data = nests,
    samplers = cvpart,
    domain = list(geometry = gorillas_sf$mesh),
    control.gcpo = list(enable = TRUE, type.cv = "joint")
  )
  result <- bru_block_gcpo(fit1)
  result$blocks
  result$gcpo
  sum(log(result$gcpo))

  # Multiple models
  pcmatern <- INLA::inla.spde2.pcmatern(
    gorillas_sf$mesh,
    prior.sigma = c(1, 0.01),
    prior.range = c(0.1, 0.01)
  )
  fit2 <- lgcp(
    geometry ~ Intercept(1) + field(geometry, model = pcmatern),
    data = nests,
    samplers = cvpart,
    domain = list(geometry = gorillas_sf$mesh),
    control.gcpo = list(enable = TRUE, type.cv = "joint")
  )
  result <- bru_block_gcpo(list(Model1 = fit1, Model2 = fit2))
  str(result)
}
#> List of 2
#>  $ Model1:List of 2
#>   ..$ blocks:List of 1
#>   .. ..$ lhood1: int [1:114] 39 28 33 40 62 68 51 32 29 73 ...
#>   ..$ gcpo  : num [1:114] 4.36e+40 3.91e+58 3.83e+49 1.38e+30 5.28e+89 ...
#>  $ Model2:List of 2
#>   ..$ blocks:List of 1
#>   .. ..$ lhood1: int [1:114] 39 28 33 40 62 68 51 32 29 73 ...
#>   ..$ gcpo  : num [1:114] 5.05e+48 1.83e+75 6.13e+61 1.37e+34 1.88e+125 ...
# }