feature issue #924

DESCRIPTION

@@ -2,8 +2,8 @@ Package: brms
 Encoding: UTF-8
 Type: Package
 Title: Bayesian Regression Models using 'Stan'
-Version: 2.13.2
-Date: 2020-06-07
+Version: 2.13.3
+Date: 2020-06-09
 Authors@R: person("Paul-Christian", "Bürkner", email = "paul.buerkner@gmail.com",
   role = c("aut", "cre"))
 Depends: 

NEWS.md

@@ -1,5 +1,10 @@
 # brms 2.13.0++
 
+### New Features
+
+* Fix shape parameters across multiple monotonic terms via argument
+`id` in function `mo` to ensure conditionally monotonic effects. (#924)
+
 ### Bug Fixes
 
 * Fix generated Stan Code of models with improper global priors and

R/data-predictor.R

@@ -346,22 +346,29 @@ data_sp <- function(bterms, data, prior = brmsprior(), basis = NULL) {
     }
     out[[paste0("Jmo", p)]] <- Jmo
     # prepare prior concentration of simplex parameters
-    simo_coef <- get_simo_labels(spef)
-    for (i in seq_along(simo_coef)) {
-      simo_prior <- subset2(prior, 
-        class = "simo", coef = simo_coef[i], ls = px
-      )
-      simo_prior <- simo_prior$prior
-      if (isTRUE(nzchar(simo_prior))) {
-        simo_prior <- eval_dirichlet(simo_prior)
-        if (length(simo_prior) != Jmo[i]) {
-          stop2("Invalid Dirichlet prior for the simplex of coefficient '",
-                simo_coef[i], "'. Expected input of length ", Jmo[i], ".")
+    simo_coef <- get_simo_labels(spef, use_id = TRUE)
+    ids <- unlist(spef$ids_mo)
+    for (j in seq_along(simo_coef)) {
+      # index of first ID appearance
+      j_id <- match(ids[j], ids)
+      if (is.na(ids[j]) || j_id == j) {
+        # only evaluate priors without ID or first appearance of the ID
+        # all other parameters will be copied over in the Stan code
+        simo_prior <- subset2(prior, 
+          class = "simo", coef = simo_coef[j], ls = px
+        )
+        simo_prior <- simo_prior$prior
+        if (isTRUE(nzchar(simo_prior))) {
+          simo_prior <- eval_dirichlet(simo_prior)
+          if (length(simo_prior) != Jmo[j]) {
+            stop2("Invalid Dirichlet prior for the simplex of coefficient '",
+                  simo_coef[j], "'. Expected input of length ", Jmo[j], ".")
+          }
+        } else {
+          simo_prior <- rep(1, Jmo[j])
         }
-      } else {
-        simo_prior <- rep(1, Jmo[i])
+        out[[paste0("con_simo", p, "_", j)]] <- as.array(simo_prior)
       }
-      out[[paste0("con_simo", p, "_", i)]] <- as.array(simo_prior)
     }
   }
   out

R/formula-sp.R

@@ -104,16 +104,23 @@ mi <- function(x) {
 #' evaluate its arguments -- it exists purely to help set up a model.
 #' 
 #' @param x An integer variable or an ordered factor to be modeled as monotonic.
+#' @param id Optional character string. All monotonic terms
+#'  with the same \code{id} within one formula  will be modeled as
+#'  having the same simplex (shape) parameter vector. If all monotonic terms
+#'  of the same predictor have the same \code{id}, the resulting
+#'  predictions will be conditionally monotonic for all values of
+#'  interacting covariates (Bürkner & Charpentier, 2020).
 #'  
-#' @details For detailed documentation see \code{help(brmsformula)}
-#'   as well as \code{vignette("brms_monotonic")}.
+#' @details See Bürkner and Charpentier (2020) for the underlying theory. For
+#'   detailed documentation of the formula syntax used for monotonic terms, 
+#'   see \code{help(brmsformula)} as well as \code{vignette("brms_monotonic")}.
 #' 
 #' @seealso \code{\link{brmsformula}}
 #' 
 #' @references 
-#' Bürkner P. C. & Charpentier, E. (in review). Monotonic Effects: A Principled 
-#' Approach for Including Ordinal Predictors in Regression Models. PsyArXiv 
-#' preprint.
+#' Bürkner P. C. & Charpentier E. (2020). Modeling Monotonic Effects of Ordinal
+#' Predictors in Regression Models. British Journal of Mathematical and 
+#' Statistical Psychology. doi:10.1111/bmsp.12195
 #'   
 #' @examples   
 #' \dontrun{
@@ -127,27 +134,29 @@ mi <- function(x) {
 #' 
 #' # fit a simple monotonic model
 #' fit1 <- brm(ls ~ mo(income), data = dat)
-#' 
-#' # summarise the model
 #' summary(fit1)
 #' plot(fit1, N = 6)
 #' plot(conditional_effects(fit1), points = TRUE)
 #' 
 #' # model interaction with other variables
 #' dat$x <- sample(c("a", "b", "c"), 100, TRUE)
 #' fit2 <- brm(ls ~ mo(income)*x, data = dat)
-#' 
-#' # summarise the model
 #' summary(fit2)
 #' plot(conditional_effects(fit2), points = TRUE)
+#' 
+#' # ensure conditional monotonicity
+#' fit3 <- brm(ls ~ mo(income, id = "i")*x, data = dat)
+#' summary(fit3)
+#' plot(conditional_effects(fit3), points = TRUE)
 #' } 
 #'  
 #' @export
-mo <- function(x) {
+mo <- function(x, id = NA) {
   # use 'term' for consistency with other special terms
   term <- deparse(substitute(x))
+  id <- as_one_character(id, allow_na = TRUE)
   label <- deparse(match.call())
-  out <- nlist(term, label)
+  out <- nlist(term, id, label)
   class(out) <- c("mo_term", "sp_term")
   out
 }
@@ -299,7 +308,7 @@ tidy_spef <- function(x, data) {
   out <- data.frame(term = rm_wsp(colnames(mm)), stringsAsFactors = FALSE)
   out$coef <- rename(out$term)
   calls_cols <- paste0("calls_", all_sp_types())
-  for (col in c(calls_cols, "joint_call", "vars_mi", "Imo")) {
+  for (col in c(calls_cols, "joint_call", "vars_mi", "ids_mo", "Imo")) {
     out[[col]] <- vector("list", nrow(out))
   }
   kmo <- 0
@@ -311,13 +320,15 @@ tidy_spef <- function(x, data) {
       out$calls_mo[[i]] <- terms_split[[i]][take_mo]
       nmo <- length(out$calls_mo[[i]])
       out$Imo[[i]] <- (kmo + 1):(kmo + nmo)
+      out$ids_mo[[i]] <- rep(NA, nmo)
       kmo <- kmo + nmo
       for (j in seq_along(out$calls_mo[[i]])) {
         mo_term <- out$calls_mo[[i]][[j]]
         mo_match <- get_matches_expr(regex_sp("mo"), mo_term)
         if (length(mo_match) > 1L || nchar(mo_match) < nchar(mo_term)) {
           stop2("The monotonic term '",  mo_term, "' is invalid.")
         }
+        out$ids_mo[[i]][[j]] <- eval2(mo_term)[["id"]]
       }
     }
     # prepare me terms
@@ -349,9 +360,19 @@ tidy_spef <- function(x, data) {
 
 # extract names of monotonic simplex parameters 
 # @param spef output of tidy_spef
-get_simo_labels <- function(spef) {
-  fun <- function(i) paste0(spef$coef[i], seq_along(spef$Imo[[i]]))
-  ulapply(which(lengths(spef$Imo) > 0), fun)
+# @param use_id use the 'id' argument to construct simo labels?
+# @return a character vector of length nrow(spef)
+get_simo_labels <- function(spef, use_id = FALSE) {
+  out <- named_list(spef$term)
+  I <- which(lengths(spef$Imo) > 0)
+  for (i in I) {
+    # use the ID as label if specified
+    out[[i]] <- ifelse(
+      use_id & !is.na(spef$ids_mo[[i]]), spef$ids_mo[[i]],
+      paste0(spef$coef[i], seq_along(spef$Imo[[i]]))
+    )
+  }
+  unlist(out)
 }
 
 # standard errors of variables with missing values

R/priors.R

@@ -730,7 +730,7 @@ prior_sp <- function(bterms, data, ...) {
     prior <- prior + brmsprior(
       class = "b", coef = c("", spef$coef), ls = px
     )
-    simo_coef <- get_simo_labels(spef)
+    simo_coef <- get_simo_labels(spef, use_id = TRUE)
     if (length(simo_coef)) {
       prior <- prior + brmsprior(
         class = "simo", coef = simo_coef, ls = px

R/stan-predictor.R

@@ -920,25 +920,44 @@ stan_sp <- function(bterms, data, prior, stanvars, ranef, meef, ...) {
   }
 
   # include special Stan code for monotonic effects
-  I <- unlist(spef$Imo)
-  if (length(I)) {
-    I <- seq_len(max(I))
+  which_Imo <- which(lengths(spef$Imo) > 0)
+  if (any(which_Imo)) {
     str_add(out$data) <- glue(
       "  int<lower=1> Imo{p};  // number of monotonic variables\n",
-      "  int<lower=1> Jmo{p}[Imo{p}];  // length of simplexes\n",
-      "  // monotonic variables\n",
-      cglue("  int Xmo{p}_{I}[N{resp}];\n"),
-      "  // prior concentration of monotonic simplexes\n",
-      cglue("  vector[Jmo{p}[{I}]] con_simo{p}_{I};\n")
-    )
-    str_add(out$par) <- glue(
-      "  // simplexes of monotonic effects\n",
-      cglue("  simplex[Jmo{p}[{I}]] simo{p}_{I};\n")
-    ) 
-    str_add(out$prior) <- cglue(
-      "  target += dirichlet_lpdf(simo{p}_{I} | con_simo{p}_{I});\n"
-    )
+      "  int<lower=1> Jmo{p}[Imo{p}];  // length of simplexes\n"
+    )
+    ids <- unlist(spef$ids_mo)
+    for (i in which_Imo) {
+      for (k in seq_along(spef$Imo[[i]])) {
+        j <- spef$Imo[[i]][[k]]
+        id <- spef$ids_mo[[i]][[k]]
+        # index of first ID appearance
+        j_id <- match(id, ids)
+        str_add(out$data) <- glue(
+          "  int Xmo{p}_{j}[N{resp}];  // monotonic variable\n"
+        )
+        if (is.na(id) || j_id == j) {
+          # no ID or first appearance of the ID
+          str_add(out$data) <- glue(
+            "  vector[Jmo{p}[{j}]] con_simo{p}_{j};", 
+            "  // prior concentration of monotonic simplex\n"
+          )
+          str_add(out$par) <- glue(
+            "  simplex[Jmo{p}[{j}]] simo{p}_{j};  // monotonic simplex\n"
+          )
+          str_add(out$prior) <- glue(
+            "  target += dirichlet_lpdf(simo{p}_{j} | con_simo{p}_{j});\n"
+          )
+        } else {
+          # use the simplex shared across all terms of the same ID
+          str_add(out$tpar_def) <- glue(
+            "  simplex[Jmo{p}[{j}]] simo{p}_{j} = simo{p}_{j_id};\n"
+          )
+        }
+      }
+    }
   }
+
   stan_special_priors <- stan_special_prior_local(
     prior, class = "bsp", ncoef = nrow(spef), 
     px = px, center_X = FALSE

tests/testthat/tests.make_stancode.R

@@ -904,6 +904,16 @@ test_that("monotonic effects appear in the Stan code", {
   scode <- make_stancode(y ~ mo(x1):y, dat)
   expect_match2(scode, "mu[n] += (bsp[1]) * mo(simo_1, Xmo_1[n]) * Csp_1[n];")
 
+  # test issue #924 (conditional monotonicity)
+  prior <- c(prior(dirichlet(c(1,0.5,2)), simo, coef = "v"),
+             prior(dirichlet(c(1,0.5,2)), simo, coef = "w"))
+  scode <- make_stancode(y ~ y*mo(x1, id = "v")*mo(x2, id = "w"), 
+                         dat, prior = prior)
+  expect_match2(scode, "target += dirichlet_lpdf(simo_1 | con_simo_1);")
+  expect_match2(scode, "target += dirichlet_lpdf(simo_2 | con_simo_2);")
+  expect_match2(scode, "simplex[Jmo[6]] simo_6 = simo_2;")
+  expect_match2(scode, "simplex[Jmo[7]] simo_7 = simo_1;")
+
   expect_error(
     make_stancode(y ~ mo(x1) + (mo(x2) | x2), dat),
     "Special group-level terms require"

tests/testthat/tests.make_standata.R

@@ -367,6 +367,15 @@ test_that("make_standata correctly prepares data for monotonic effects", {
   expect_equal(sdata$con_simo_1, as.array(c(1, 0.5, 2)))
   expect_equal(sdata$con_simo_3, as.array(c(1, 0.5, 2)))
 
+  # test issue #924 (conditional monotonicity)
+  prior <- c(prior(dirichlet(c(1, 0.5, 2)), simo, coef = "v"),
+             prior(dirichlet(c(1,3)), simo, coef = "w"))
+  sdata <- make_standata(y ~ y*mo(x1, id = "v")*mo(x2, id = "w"), 
+                         data, prior = prior)
+  expect_equal(sdata$con_simo_1, as.array(c(1, 0.5, 2)))
+  expect_equal(sdata$con_simo_2, as.array(c(1, 3)))
+  expect_true(!"sdata$con_simo_3" %in% names(sdata))
+
   expect_error(
     make_standata(y ~ mo(z), data = data),
     "Monotonic predictors must be integers or ordered factors"