diff --git a/R/monotonic.R b/R/monotonic.R
index b1573a96..5d95634a 100644
--- a/R/monotonic.R
+++ b/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)
   }
 
diff --git a/R/mvgam_formulae.R b/R/mvgam_formulae.R
index d9889046..a16ff2f0 100644
--- a/R/mvgam_formulae.R
+++ b/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
diff --git a/R/stan_utils.R b/R/stan_utils.R
index a005ff87..60698f76 100644
--- a/R/stan_utils.R
+++ b/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],
diff --git a/man/monotonic.Rd b/man/monotonic.Rd
index 90668ffc..f65b39e7 100644
--- a/man/monotonic.Rd
+++ b/man/monotonic.Rd
@@ -97,6 +97,35 @@ plot_predictions(mod2,
                  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)
 }
 }
 \references{
diff --git a/man/mvgam_formulae.Rd b/man/mvgam_formulae.Rd
index 929b7794..de3d1019 100644
--- a/man/mvgam_formulae.Rd
+++ b/man/mvgam_formulae.Rd
@@ -21,11 +21,15 @@ are non-identifiable.
 \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 \code{s(x, bs = 'moi')}) or
+or decreasing splines (using \code{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
@@ -41,7 +45,8 @@ extensive documentation for the \code{mgcv} package.
 \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
diff --git a/src/mvgam.dll b/src/mvgam.dll
index fa4f2bc4..85f8f1ed 100644
Binary files a/src/mvgam.dll and b/src/mvgam.dll differ
diff --git a/tests/testthat/Rplots.pdf b/tests/testthat/Rplots.pdf
index c63ab130..1b6f1331 100644
Binary files a/tests/testthat/Rplots.pdf and b/tests/testthat/Rplots.pdf differ
diff --git a/tests/testthat/test-monotonic.R b/tests/testthat/test-monotonic.R
new file mode 100644
index 00000000..a38a3a46
--- /dev/null
+++ b/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)))
+})