better monotonic doc; allow for by variables

R/monotonic.R

@@ -82,6 +82,35 @@
 #'                  points = 0.5)
 #'
 #' plot(mod2, type = 'smooth', realisations = TRUE)
+#'
+#' # 'by' terms that produce a different smooth for each level of the 'by'
+#' # factor are also allowed
+#' set.seed(123123)
+#' x <- runif(80) * 4 - 1
+#' x <- sort(x)
+#'
+#' # Two different monotonic smooths, one for each factor level
+#' f <- exp(4 * x) / (1 + exp(4 * x))
+#' f2 <- exp(3.5 * x) / (1 + exp(3 * x))
+#' fac <- c(rep('a', 80), rep('b', 80))
+#' y <- c(f + rnorm(80) * 0.1,
+#'        f2 + rnorm(80) * 0.2)
+#' plot(x, y[1:80])
+#' plot(x, y[81:160])
+#'
+#' # Gather all data into a data.frame, including the factor 'by' variable
+#' mod_data <- data.frame(y, x, fac = as.factor(fac))
+#' mod_data$time <- 1:NROW(mod_data)
+#'
+#' # Fit a model with different smooths per factor level
+#' mod <- mvgam(y ~ s(x, bs = 'moi', by = fac, k = 8),
+#'              data = mod_data,
+#'              family = gaussian())
+#'
+#' # Visualise the different monotonic functions
+#' plot_predictions(mod, condition = c('x', 'fac', 'fac'),
+#'                  points = 0.5)
+#' plot(mod, type = 'smooth', realisations = TRUE)
 #' }
 smooth.construct.moi.smooth.spec <- function(object, data, knots){
 
@@ -351,7 +380,7 @@ add_mono_model_file = function(model_file,
                                    '; // monotonic basis coefficient indices\n'))
     mono_idx_data <- list(coef_indices)
     names(mono_idx_data) <- paste0('b_idx_',
-                                   mono_names_clean)
+                                   mono_names_clean[covariate])
     model_data <- append(model_data, mono_idx_data)
   }
 

R/mvgam_formulae.R

@@ -13,11 +13,15 @@
 #' \cr
 #' \cr
 #' The formulae supplied to \code{\link{mvgam}} are exactly like those supplied to
-#' \code{\link{glm}} except that smooth terms, \code{\link[mgcv]{s}},
+#' \code{\link{glm}} except that smooth terms,
+#' \code{\link[mgcv]{s}},
 #' \code{\link[mgcv]{te}},
 #' \code{\link[mgcv]{ti}} and
 #' \code{\link[mgcv]{t2}},
-#' time-varying effects using \code{\link{dynamic}}, as well as
+#' time-varying effects using \code{\link{dynamic}},
+#' monotonically increasing (using `s(x, bs = 'moi')`) or
+#' or decreasing splines (using `s(x, bs = 'mod')`; see \code{\link{monotonic}} for
+#' details), as well as
 #' Gaussian Process functions using \code{\link[brms]{gp}},
 #' can be added to the right hand side (and \code{.} is not supported in \code{mvgam} formulae).
 #' \cr
@@ -31,7 +35,8 @@
 #' \code{\link[mgcv]{jagam}},
 #' \code{\link[mgcv]{gam}},
 #' \code{\link[mgcv]{s}},
-#' \code{\link[stats]{formula}}
+#' \code{\link[stats]{formula}},
+#' \code{\link{monotonic}}
 #' @author Nicholas J Clark
 #' @name mvgam_formulae
 NULL

R/stan_utils.R

@@ -2613,7 +2613,8 @@ add_trend_predictors = function(trend_formula,
       model_data <- append(model_data, idx_data)
 
       idx_lines <- c(grep('int trend_idx', trend_model_file),
-                     grep('// gp basis coefficient indices', trend_model_file))
+                     grep('// gp basis coefficient indices', trend_model_file),
+                     grep('// monotonic basis coefficient indices', trend_model_file))
       model_file[min(grep('data {', model_file, fixed = TRUE))] <-
         paste0('data {\n',
                paste(trend_model_file[idx_lines],

tests/testthat/test-monotonic.R

@@ -0,0 +1,96 @@
+context("monotonic")
+
+# Simulate data from a monotonically increasing function
+set.seed(123123)
+x <- runif(80) * 4 - 1
+x <- sort(x)
+f <- exp(4 * x) / (1 + exp(4 * x))
+y <- f + rnorm(80) * 0.1
+mod_data <- data.frame(y = y, x = x, z = rnorm(80),
+                       time = 1:80)
+
+test_that("k must be an even integer for s(bs = 'moi')", {
+  expect_error(mvgam(y ~ s(x, bs = 'moi', k = 11),
+                data = mod_data,
+                family = gaussian()),
+               "Argument 'k(bs = 'moi')'  must be an even integer",
+               fixed = TRUE)
+
+  expect_error(mvgam(y ~ s(x, bs = 'moi', k = 1),
+                     data = mod_data,
+                     family = gaussian()),
+               "Basis dimension is too small",
+               fixed = TRUE)
+})
+
+test_that("monotonic only works for one dimensional smooths", {
+  expect_error(mvgam(y ~ s(x, z, bs = 'moi', k = 10),
+                     data = mod_data,
+                     family = gaussian()),
+               "Monotonic basis only handles 1D smooths",
+               fixed = TRUE)
+})
+
+test_that("monotonic for observation models working properly", {
+  mod <- mvgam(y ~ z + s(x, bs = 'moi', k = 18),
+               data = mod_data,
+               family = gaussian(),
+               run_model = FALSE)
+
+  # Monotonic indices should be in the model_data
+  expect_true("b_idx_s_x_" %in% names(mod$model_data))
+
+  # The smooth should be an MOI class
+  expect_true(inherits(mod$mgcv_model$smooth[[1]], 'moi.smooth'))
+
+  # The coefficients should be fixed to be non-negative
+  expect_true(any(grepl('b[b_idx_s_x_] = abs(b_raw[b_idx_s_x_]) * 1;',
+                        mod$model_file, fixed = TRUE)))
+
+  # Repeat a check for decreasing functions
+  mod <- mvgam(y ~ z + s(x, bs = 'mod', k = 18),
+               data = mod_data,
+               family = gaussian(),
+               run_model = FALSE)
+
+  # The smooth should be an MOD class
+  expect_true(inherits(mod$mgcv_model$smooth[[1]], 'mod.smooth'))
+
+  # The coefficients should be fixed to be non-positive
+  expect_true(any(grepl('b[b_idx_s_x_] = abs(b_raw[b_idx_s_x_]) * -1;',
+                        mod$model_file, fixed = TRUE)))
+})
+
+test_that("monotonic for process models working properly", {
+  mod <- mvgam(y ~ 0,
+               trend_formula = ~ z + s(x, bs = 'moi', k = 18),
+               trend_model = RW(),
+               data = mod_data,
+               family = gaussian(),
+               run_model = FALSE)
+
+  # Monotonic indices should be in the model_data
+  expect_true("b_trend_idx_s_x_" %in% names(mod$model_data))
+
+  # The smooth should be an MOI class
+  expect_true(inherits(mod$trend_mgcv_model$smooth[[1]], 'moi.smooth'))
+
+  # The coefficients should be fixed to be non-negative
+  expect_true(any(grepl('b_trend[b_trend_idx_s_x_] = abs(b_raw_trend[b_trend_idx_s_x_]) * 1;',
+                        mod$model_file, fixed = TRUE)))
+
+  # And for decreasing
+  mod <- mvgam(y ~ 0,
+               trend_formula = ~ z + s(x, bs = 'mod', k = 18),
+               trend_model = RW(),
+               data = mod_data,
+               family = gaussian(),
+               run_model = FALSE)
+
+  # The smooth should be an MOD class
+  expect_true(inherits(mod$trend_mgcv_model$smooth[[1]], 'mod.smooth'))
+
+  # The coefficients should be fixed to be non-positive
+  expect_true(any(grepl('b_trend[b_trend_idx_s_x_] = abs(b_raw_trend[b_trend_idx_s_x_]) * -1;',
+                        mod$model_file, fixed = TRUE)))
+})