begin support for ggplot plotting

Author: Nicholas Clark

DESCRIPTION

@@ -52,6 +52,7 @@ Suggests:
     rjags,
     coda,
     runjags,
+    tibble (>= 3.0.0),
     usethis,
     testthat
 Additional_repositories: https://mc-stan.org/r-packages/

NAMESPACE

@@ -14,6 +14,8 @@ S3method(as_draws_matrix,mvgam)
 S3method(as_draws_rvars,mvgam)
 S3method(coef,mvgam)
 S3method(conditional_effects,mvgam)
+S3method(eval_smooth,mod.smooth)
+S3method(eval_smooth,moi.smooth)
 S3method(find_predictors,mvgam)
 S3method(find_predictors,mvgam_prefit)
 S3method(fitted,mvgam)
@@ -130,6 +132,7 @@ importFrom(brms,pstudent_t)
 importFrom(brms,qstudent_t)
 importFrom(brms,rbeta_binomial)
 importFrom(brms,rstudent_t)
+importFrom(brms,stancode)
 importFrom(brms,student)
 importFrom(ggplot2,scale_colour_discrete)
 importFrom(ggplot2,scale_fill_discrete)
@@ -154,6 +157,14 @@ importFrom(graphics,polygon)
 importFrom(graphics,rect)
 importFrom(graphics,rug)
 importFrom(graphics,title)
+importFrom(gratia,by_level)
+importFrom(gratia,by_variable)
+importFrom(gratia,eval_smooth)
+importFrom(gratia,is_factor_by_smooth)
+importFrom(gratia,smooth_dim)
+importFrom(gratia,smooth_label)
+importFrom(gratia,smooths)
+importFrom(gratia,spline_values)
 importFrom(insight,find_predictors)
 importFrom(insight,get_data)
 importFrom(insight,get_predictors)
@@ -271,6 +282,7 @@ importFrom(stats,ts)
 importFrom(stats,update)
 importFrom(stats,update.formula)
 importFrom(stats,var)
+importFrom(tibble,add_column)
 importFrom(utils,getFromNamespace)
 importFrom(utils,head)
 importFrom(utils,lsf.str)

R/monotonic.R

@@ -412,3 +412,223 @@ add_mono_model_file = function(model_file,
   return(list(model_file = model_file,
               model_data = model_data))
 }
+
+#' Evaluation of a monotonically increasing function
+#' These evaluation functions are needed so that gratia::draw methods work with mvgam
+#' monotonic smooths
+#' @importFrom gratia eval_smooth smooths by_variable smooth_label spline_values is_factor_by_smooth by_level smooth_dim
+#' @rdname monotonic
+#' @export
+eval_smooth.moi.smooth = function(smooth,
+                                  model,
+                                  n = 100,
+                                  n_3d = NULL,
+                                  n_4d = NULL,
+                                  data = NULL,
+                                  unconditional = FALSE,
+                                  overall_uncertainty = TRUE,
+                                  dist = NULL,
+                                  ...) {
+  model$cmX <- model$coefficients
+
+  ## deal with data if supplied
+  data <- process_user_data_for_eval(
+    data = data, model = model,
+    n = n, n_3d = n_3d, n_4d = n_4d,
+    id = which_smooth(
+      model,
+      smooth_label(smooth)
+    )
+  )
+
+  by_var <- gratia::by_variable(smooth) # even if not a by as we want NA later
+  if (by_var == "NA") {
+    by_var <- NA_character_
+  }
+
+  ## values of spline at data
+  eval_sm <- gratia::spline_values(smooth,
+                                   data = data,
+                                   unconditional = unconditional,
+                                   model = model,
+                                   overall_uncertainty = overall_uncertainty
+  )
+
+  ## add on info regarding by variable
+  eval_sm <- add_by_var_column(eval_sm, by_var = by_var)
+  ## add on spline type info
+  eval_sm <- add_smooth_type_column(eval_sm, sm_type = "Mono inc P spline")
+
+  # set some values to NA if too far from the data
+  if (gratia::smooth_dim(smooth) == 2L && (!is.null(dist) && dist > 0)) {
+    eval_sm <- too_far_to_na(smooth,
+                             input = eval_sm,
+                             reference = model[["model"]],
+                             cols = c(".estimate", ".se"),
+                             dist = dist
+    )
+  }
+  ## return
+  eval_sm
+}
+
+#' Evaluation of a monotonically decreasing function
+#' @rdname monotonic
+#' @param model an object of class `"gam"`
+#' @param smooth a smooth object of class `"moi.smooth"` or `"mod.smooth"`
+#' @param n numeric; the number of points over the range of the covariate at
+#'   which to evaluate the smooth.
+#' @param n_3d,n_4d numeric; the number of points over the range of last
+#'   covariate in a 3D or 4D smooth. The default is `NULL` which achieves the
+#'   standard behaviour of using `n` points over the range of all covariate,
+#'   resulting in `n^d` evaluation points, where `d` is the dimension of the
+#'   smooth. For `d > 2` this can result in very many evaluation points and slow
+#'   performance. For smooths of `d > 4`, the value of `n_4d` will be used for
+#'   all dimensions `> 4`, unless this is `NULL`, in which case the default
+#'   behaviour (using `n` for all dimensions) will be observed.
+#' @param data a data frame of covariate values at which to evaluate the
+#'   smooth.
+#' @param unconditional logical; should confidence intervals include the
+#'   uncertainty due to smoothness selection? If `TRUE`, the corrected Bayesian
+#'   covariance matrix will be used.
+#' @param overall_uncertainty logical; should the uncertainty in the model
+#'  constant term be included in the standard error of the evaluate values of
+#'  the smooth?
+#' @param dist numeric; if greater than 0, this is used to determine when
+#'   a location is too far from data to be plotted when plotting 2-D smooths.
+#'   The data are scaled into the unit square before deciding what to exclude,
+#'   and `dist` is a distance within the unit square. See
+#'   [mgcv::exclude.too.far()] for further details.
+#' @param ... ignored.
+#' @export
+eval_smooth.mod.smooth = function(smooth,
+                                  model,
+                                  n = 100,
+                                  n_3d = NULL,
+                                  n_4d = NULL,
+                                  data = NULL,
+                                  unconditional = FALSE,
+                                  overall_uncertainty = TRUE,
+                                  dist = NULL,
+                                  ...) {
+  model$cmX <- model$coefficients
+
+  ## deal with data if supplied
+  data <- process_user_data_for_eval(
+    data = data, model = model,
+    n = n, n_3d = n_3d, n_4d = n_4d,
+    id = which_smooth(
+      model,
+      smooth_label(smooth)
+    )
+  )
+
+  by_var <- gratia::by_variable(smooth) # even if not a by as we want NA later
+  if (by_var == "NA") {
+    by_var <- NA_character_
+  }
+
+  ## values of spline at data
+  eval_sm <- gratia::spline_values(smooth,
+                                   data = data,
+                                   unconditional = unconditional,
+                                   model = model,
+                                   overall_uncertainty = overall_uncertainty
+  )
+
+  ## add on info regarding by variable
+  eval_sm <- add_by_var_column(eval_sm, by_var = by_var)
+  ## add on spline type info
+  eval_sm <- add_smooth_type_column(eval_sm, sm_type = "Mono dec P spline")
+
+  # set some values to NA if too far from the data
+  if (gratia::smooth_dim(smooth) == 2L && (!is.null(dist) && dist > 0)) {
+    eval_sm <- too_far_to_na(smooth,
+                             input = eval_sm,
+                             reference = model[["model"]],
+                             cols = c(".estimate", ".se"),
+                             dist = dist
+    )
+  }
+  ## return
+  eval_sm
+}
+
+#' Utility functions; full credit goes to Gavin Simpson, the developer and
+#' maintainer of the gratia package
+#' @noRd
+`is.gamm` <- function(object) {
+  inherits(object, "gamm")
+}
+
+#' @noRd
+`is.gamm4` <- function(object) {
+  is.list(object) & (!is.null(object[["gam"]]))
+}
+
+#' @noRd
+`is.gam` <- function(object) {
+  inherits(object, "gam")
+}
+
+#' @noRd
+`is.bam` <- function(object) {
+  inherits(object, "bam")
+}
+
+#' @noRd
+`which_smooth` <- function(object, term) {
+  if (is.gamm(object) || is.gamm4(object)) {
+    object <- object[["gam"]]
+  }
+  smooths <- smooths(object)
+  which(term == smooths)
+}
+
+#' @noRd
+`process_user_data_for_eval` <- function(
+    data, model, n, n_3d, n_4d, id,
+    var_order = NULL) {
+  if (is.null(data)) {
+    data <- smooth_data(
+      model = model,
+      n = n,
+      n_3d = n_3d,
+      n_4d = n_4d,
+      id = id,
+      var_order = var_order
+    )
+  } else {
+    smooth <- get_smooths_by_id(model, id)[[1L]]
+    vars <- smooth_variable(smooth)
+    by_var <- by_variable(smooth)
+    if (!identical(by_var, "NA")) {
+      vars <- append(vars, by_var)
+    }
+    ## if this is a by variable, filter the by variable for the required
+    ## level now
+    if (gratia::is_factor_by_smooth(smooth)) {
+      data <- data %>% filter(.data[[by_var]] == gratia::by_level(smooth))
+    }
+  }
+  data
+}
+
+#' @importFrom tibble add_column
+#' @noRd
+`add_by_var_column` <- function(object, by_var, n = NULL) {
+  if (is.null(n)) {
+    n <- NROW(object)
+  }
+  insight::check_if_installed("tibble")
+  tibble::add_column(object, .by = rep(by_var, times = n), .after = 1L)
+}
+
+#' @noRd
+`add_smooth_type_column` <- function(object, sm_type, n = NULL) {
+  if (is.null(n)) {
+    n <- NROW(object)
+  }
+  insight::check_if_installed("tibble")
+  tibble::add_column(object, .type = rep(sm_type, times = n), .after = 1L)
+}

R/stan_utils.R

@@ -34,6 +34,7 @@ print.mvgammodel = function(x, ...){
 }
 
 #' @export
+#' @importFrom brms stancode
 #' @rdname code
 stancode.mvgam_prefit = function(object){