Add method bru_response_size and ensure scalar predictors are expande…

…d to matching size when appropriate

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: inlabru
 Type: Package
 Title: Bayesian Latent Gaussian Modelling using INLA and Extensions
-Version: 2.11.1.9012
+Version: 2.11.1.9013
 Authors@R: c(
   person("Finn", "Lindgren", email = "finn.lindgren@gmail.com",
          role = c("aut", "cre", "cph"),

NAMESPACE

@@ -35,6 +35,10 @@ S3method(bru_mapper,fm_mesh_1d)
 S3method(bru_mapper,fm_mesh_2d)
 S3method(bru_mapper,inla.mesh)
 S3method(bru_mapper,inla.mesh.1d)
+S3method(bru_response_size,bru)
+S3method(bru_response_size,bru_info)
+S3method(bru_response_size,bru_like)
+S3method(bru_response_size,bru_like_list)
 S3method(bru_timings,bru)
 S3method(bru_used,bru)
 S3method(bru_used,bru_like)
@@ -303,6 +307,7 @@ export(bru_options_get)
 export(bru_options_reset)
 export(bru_options_set)
 export(bru_rerun)
+export(bru_response_size)
 export(bru_safe_inla)
 export(bru_safe_sp)
 export(bru_standardise_names)

NEWS.md

@@ -21,6 +21,8 @@
   warnings for ambiguous cases (version `2.11.1.9011`)
 * Add `[` and `]` to disallowed character set in `bru_standardise_names()`
   (version `2.11.1.9012`)
+* Add `bru_response_size` method for extracting the response size for each
+  observation `like()` object (version `2.11.1.9013`)
 
 ## Data set updates
 

R/bru.inference.R

@@ -1322,6 +1322,7 @@ extended_bind_rows <- function(...) {
 #'   `data`. Defaults to the calling environment.
 #'
 #' @return A likelihood configuration which can be used to parameterise [bru()].
+#' @seealso [bru_response_size()], [bru_used()], [component()], [component_eval()]
 #'
 #' @example inst/examples/like.R
 
@@ -1723,6 +1724,54 @@ like <- function(formula = . ~ ., family = "gaussian", data = NULL,
   lh
 }
 
+
+#' @title Response size queries
+#'
+#' @description
+#' Extract the number of response values from `bru` and related objects.
+#' @return An `integer` vector.
+#' @param object An object from which to extract response size(s).
+#' @seealso [like()]
+#' @export
+#' @examples
+#' bru_response_size(like(y ~ 1, data = data.frame(y = rnorm(10)), family = "gaussian"))
+bru_response_size <- function(object) {
+  UseMethod("bru_response_size")
+}
+
+#' @describeIn bru_response_size Extract the number of observations from a `bru_like` object.
+#' @export
+bru_response_size.bru_like <- function(object) {
+  resp <- object[["response_data"]][[object[["response"]]]]
+  if (inherits(resp, "inla.surv")) {
+    # This special case handling is needed because inla.surv objects are
+    # lists. When the INLA package changes representation of `inla.surv`
+    # to a tibble or data.frame, this special case won't be needed anymore.
+    NROW(resp[[1]])
+  } else {
+    NROW(resp)
+  }
+}
+
+#' @describeIn bru_response_size Extract the number of observations from a `bru_like_list` object.
+#' @export
+bru_response_size.bru_like_list <- function(object) {
+  vapply(object, bru_response_size, 1L)
+}
+
+#' @describeIn bru_response_size Extract the number of observations from a `bru_info` object.
+#' @export
+bru_response_size.bru_info <- function(object) {
+  bru_response_size(object[["lhoods"]])
+}
+
+#' @describeIn bru_response_size Extract the number of observations from a `bru` object.
+#' @export
+bru_response_size.bru <- function(object) {
+  bru_response_size(object[["bru_info"]])
+}
+
+
 # Placeholder for future modularised version of like()
 ## @describeIn bru_like
 ## Alias for `like()`, with `obs` standing for "observation model"
@@ -2792,7 +2841,8 @@ nonlin_predictor <- function(param, state) {
               param[["lhoods"]][[lh_idx]],
               param[["model"]][["effects"]]
             ),
-            format = "matrix"
+            format = "matrix",
+            n_pred = bru_response_size(param[["lhoods"]][[lh_idx]])
           )
         )
       }
@@ -3762,6 +3812,30 @@ iinla <- function(model, lhoods, initial = NULL, options) {
     stk.data <- INLA::inla.stack.data(stk, .response.name = "BRU.response")
   }
   inla.options$control.predictor$A <- INLA::inla.stack.A(stk)
+
+  N_response <- sum(bru_response_size(lhoods))
+  N_predictor <- NROW(inla.options$control.predictor$A)
+  if (N_response != N_predictor) {
+    msg <- paste0(
+      "The total number of response values (N=",
+      N_response,
+      ") and predictor values (N=",
+      N_predictor,
+      ") do not match.\n",
+      "  This is likely due to a mistake in the component or predictor constructions."
+    )
+    if ((N_response > 1L) && (N_predictor == 1L)) {
+      msg <- paste(msg, "\nPerhaps you only have components with scalar inputs?")
+    }
+    bru_log_message(
+      msg,
+      verbose = options$bru_verbose,
+      verbose_store = options$bru_verbose_store,
+      verbosity = 1
+    )
+    stop(msg)
+  }
+
   bru_log_message(
     "iinla: Model initialisation completed",
     verbose = options$bru_verbose,
@@ -4085,6 +4159,29 @@ iinla <- function(model, lhoods, initial = NULL, options) {
           )
         }
         inla.options$control.predictor$A <- INLA::inla.stack.A(stk)
+
+        N_response <- sum(bru_response_size(lhoods))
+        N_predictor <- NROW(inla.options$control.predictor$A)
+        if (N_response != N_predictor) {
+          msg <- paste0(
+            "The total number of response values (N=",
+            N_response,
+            ") and predictor values (N=",
+            N_predictor,
+            ") do not match.\n",
+            "  This is likely due to a mistake in the component or predictor constructions."
+          )
+          if ((N_response > 1L) && (N_predictor == 1L)) {
+            msg <- paste(msg, "\nPerhaps you only have components with scalar inputs?")
+          }
+          bru_log_message(
+            msg,
+            verbose = options$bru_verbose,
+            verbose_store = options$bru_verbose_store,
+            verbosity = 1
+          )
+          stop(msg)
+        }
       }
       # Store the state
       states <- c(states, list(state))

R/model.R

@@ -154,6 +154,8 @@ print.summary_bru_model <- function(x, ...) {
 #' computations. Default NULL, leaves it up to INLA.
 #' When seed != 0, overridden to "1:1"
 #' @param used A [bru_used()] object, or NULL (default)
+#' @param n_pred integer. If provided, scalar predictor results are expanded to
+#' vectors of length `n_pred`.
 #' @param \dots Additional arguments passed on to `inla.posterior.sample`
 #' @details * `evaluate_model` is a wrapper to evaluate model state, A-matrices,
 #' effects, and predictor, all in one call.
@@ -168,6 +170,7 @@ evaluate_model <- function(model,
                            predictor = NULL,
                            format = NULL,
                            used = NULL,
+                           n_pred = NULL,
                            ...) {
   used <- bru_used(used, labels = names(model$effects))
 
@@ -208,7 +211,8 @@ evaluate_model <- function(model,
     effects = effects,
     predictor = predictor,
     format = format,
-    used = used
+    used = used,
+    n_pred = n_pred
   )
 
   values
@@ -415,6 +419,8 @@ evaluate_effect_multi_state.component_list <- function(components,
 #' expression for a state.
 #'
 #' Default: "auto"
+#' @param n_pred integer. If provided, scalar predictor results are expanded to
+#' vectors of length `n_pred`.
 #' @details For each component, e.g. "name", the state values are available as
 #'   `name_latent`, and arbitrary evaluation can be done with `name_eval(...)`,
 #'   see [component_eval()].
@@ -427,7 +433,8 @@ evaluate_predictor <- function(model,
                                effects,
                                predictor,
                                used = NULL,
-                               format = "auto") {
+                               format = "auto",
+                               n_pred = NULL) {
   stopifnot(inherits(model, "bru_model"))
   format <- match.arg(format, c("auto", "matrix", "list"))
   pred.envir <- environment(predictor)
@@ -600,6 +607,9 @@ evaluate_predictor <- function(model,
     }
 
     result_ <- eval(predictor, envir = envir, enclos = enclos)
+    if (!is.null(n_pred) && is.numeric(result_) && length(result_) == 1) {
+      result_ <- rep(result_, n_pred)
+    }
     if (k == 1) {
       if (identical(format, "auto")) {
         if ((is.vector(result_) && !is.list(result_)) ||

R/nlinla.R

@@ -34,6 +34,9 @@ bru_compute_linearisation <- function(...) {
 #' @param eps The finite difference step size
 #' @param options A `bru_options` object. The log verbosity options
 #' are used.
+#' @param n_pred The length of the predictor expression. If not `NULL`, scalar
+#' predictor evaluations are expanded to vectors of length `n_pred`.
+#'
 #' @export
 #' @rdname bru_compute_linearisation
 bru_compute_linearisation.component <- function(cmp,
@@ -49,6 +52,7 @@ bru_compute_linearisation.component <- function(cmp,
                                                 allow_latent,
                                                 allow_combine,
                                                 eps,
+                                                n_pred = NULL,
                                                 ...,
                                                 options = NULL) {
   label <- cmp[["label"]]
@@ -59,6 +63,16 @@ bru_compute_linearisation.component <- function(cmp,
     verbosity = 4
   )
 
+  if (cmp[["main"]][["type"]] %in% c("offset", "const")) {
+    # Zero-column matrix, since a const/offset has no latent state variables.
+    return(Matrix::sparseMatrix(
+      i = c(),
+      j = c(),
+      x = c(1),
+      dims = c(NROW(pred0), 0)
+    ))
+  }
+
   if (is.null(comp_simple)) {
     A <- NULL
     assume_rowwise <- FALSE
@@ -70,29 +84,25 @@ bru_compute_linearisation.component <- function(cmp,
     )
 
     assume_rowwise <- !allow_latent && !allow_combine && is.data.frame(data)
-  }
 
-  if (assume_rowwise) {
-    if (NROW(A) == 1) {
-      if (NROW(pred0) > 1) {
+    if (assume_rowwise) {
+      if (!is.null(n_pred) && (NROW(pred0) != n_pred)) {
+        stop(
+          "Number of rows (",
+          NROW(pred0),
+          ") in the predictor for component '",
+          label,
+          "' does not match the length implied by the response data (",
+          n_pred,
+          ")."
+        )
+      }
+      if (NROW(A) == 1L) {
         A <- Matrix::kronecker(rep(1, NROW(pred0)), A)
-      } else if (is.data.frame(data)) {
-        A <- Matrix::kronecker(rep(1, NROW(data)), A)
-        pred0 <- rep(pred0, NROW(A))
       }
     }
   }
 
-  if (cmp[["main"]][["type"]] %in% c("offset", "const")) {
-    # Zero-column matrix, since a const/offset has no latent state variables.
-    return(Matrix::sparseMatrix(
-      i = c(),
-      j = c(),
-      x = c(1),
-      dims = c(NROW(pred0), 0)
-    ))
-  }
-
   triplets <- list(
     i = integer(0),
     j = integer(0),
@@ -198,7 +208,13 @@ bru_compute_linearisation.component <- function(cmp,
           },
         predictor = lhood_expr,
         used = used,
-        format = "matrix"
+        format = "matrix",
+        n_pred =
+          if (assume_rowwise) {
+            length(row_subset)
+          } else {
+            n_pred
+          }
       )
       # Store sparse triplet information
       if (symmetric_diffs) {
@@ -248,6 +264,18 @@ bru_compute_linearisation.component <- function(cmp,
     x = triplets$x,
     dims = c(NROW(pred0), NROW(state[[label]]))
   )
+  if (NROW(B) != NROW(pred0)) {
+    stop(
+      "Jacobian matrix for component '",
+      label,
+      "' has ",
+      NROW(B),
+      " rows, but expected ",
+      NROW(pred0),
+      " rows based on the predictor length."
+    )
+  }
+  B
 }
 
 #' @param lhood A `bru_like` object
@@ -271,6 +299,7 @@ bru_compute_linearisation.bru_like <- function(lhood,
   )
 
   lhood_expr <- bru_like_expr(lhood, model[["effects"]])
+  n_pred <- bru_response_size(lhood)
 
   pred0 <- evaluate_predictor(
     model,
@@ -279,22 +308,9 @@ bru_compute_linearisation.bru_like <- function(lhood,
     effects = list(effects),
     predictor = lhood_expr,
     used = used,
-    format = "matrix"
+    format = "matrix",
+    n_pred = n_pred
   )
-  if (lhood[["linear"]]) {
-    # If linear, can check if the predictor is a scalar or vector,
-    # and possibly expand to full size
-    if (length(pred0) == 1) {
-      if (is.data.frame(data) ||
-        inherits(data, c(
-          "SpatialPointsDataFrame",
-          "SpatialPolygonsDataFrame",
-          "SpatialLinesDataFrame"
-        ))) {
-        pred0 <- rep(pred0, NROW(data))
-      }
-    }
-  }
 
   # Compute derivatives for each non-const/offset component
   B <- list()
@@ -336,9 +352,13 @@ bru_compute_linearisation.bru_like <- function(lhood,
             allow_latent = label %in% used[["latent"]],
             allow_combine = lhood[["allow_combine"]],
             eps = eps,
+            n_pred = n_pred,
             ...
           )
       }
+      if ((NROW(offset) == 1L) && (NROW(B[[label]]) > 1L)) {
+        offset <- matrix(offset, NROW(B[[label]]), 1)
+      }
       offset <- offset - B[[label]] %*% state[[label]]
     }
   }