Need help: Monotonic effects

[DominiqueMakowski]

The main remaining feature I could see in this package are Monotonic effects, used to model Likert scales as predictors.

Here is the stan code:

library(brms)

m <- brm(mpg ~ mo(gear), data = mtcars, refresh=0)

stancode(m)

// generated with brms 2.21.1
functions {
  /* compute monotonic effects
   * Args:
   *   scale: a simplex parameter
   *   i: index to sum over the simplex
   * Returns:
   *   a scalar between 0 and rows(scale)
   */
  real mo(vector scale, int i) {
    if (i == 0) {
      return 0;
    } else {
      return rows(scale) * sum(scale[1:i]);
    }
  }
}
data {
  int<lower=1> N;  // total number of observations
  vector[N] Y;  // response variable
  int<lower=1> Ksp;  // number of special effects terms
  int<lower=1> Imo;  // number of monotonic variables
  array[Imo] int<lower=1> Jmo;  // length of simplexes
  array[N] int Xmo_1;  // monotonic variable
  vector[Jmo[1]] con_simo_1;  // prior concentration of monotonic simplex
  int prior_only;  // should the likelihood be ignored?
}
transformed data {
}
parameters {
  real Intercept;  // temporary intercept for centered predictors
  simplex[Jmo[1]] simo_1;  // monotonic simplex
  vector[Ksp] bsp;  // special effects coefficients
  real<lower=0> sigma;  // dispersion parameter
}
transformed parameters {
  real lprior = 0;  // prior contributions to the log posterior
  lprior += student_t_lpdf(Intercept | 3, 19.2, 5.4);
  lprior += dirichlet_lpdf(simo_1 | con_simo_1);
  lprior += student_t_lpdf(sigma | 3, 0, 5.4)
    - 1 * student_t_lccdf(0 | 3, 0, 5.4);
}
model {
  // likelihood including constants
  if (!prior_only) {
    // initialize linear predictor term
    vector[N] mu = rep_vector(0.0, N);
    mu += Intercept;
    for (n in 1:N) {
      // add more terms to the linear predictor
      mu[n] += (bsp[1]) * mo(simo_1, Xmo_1[n]);
    }
    target += normal_lpdf(Y | mu, sigma);
  }
  // priors including constants
  target += lprior;
}
generated quantities {
  // actual population-level intercept
  real b_Intercept = Intercept;
}

[DominiqueMakowski]

Here is a start...

using RDatasets
data = dataset("datasets", "mtcars")

# Variables
y = data[!, :MPG]
x = data[!, :Gear]

# Define the monotonic function
function monotonic(scale::Vector{Float64}, i::Int)
    if i == 0
        return 0.0
    else
        return length(scale) * sum(scale[1:i])
    end
end