if (FALSE) {
# 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
plot(x, y)
# A standard TRPS smooth doesn't capture monotonicity
library(mgcv)
mod_data <- data.frame(y = y, x = x)
mod <- gam(
y ~ s(x, k = 16),
data = mod_data,
family = gaussian()
)
library(marginaleffects)
plot_predictions(
mod,
by = 'x',
newdata = data.frame(
x = seq(min(x) - 0.5, max(x) + 0.5, length.out = 100)
),
points = 0.5
)
# Using the 'moi' basis in mvgam rectifies this
mod_data$time <- 1:NROW(mod_data)
mod2 <- mvgam(
y ~ s(x, bs = 'moi', k = 18),
data = mod_data,
family = gaussian(),
chains = 2,
silent = 2
)
plot_predictions(
mod2,
by = 'x',
newdata = data.frame(
x = seq(min(x) - 0.5, max(x) + 0.5, length.out = 100)
),
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
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(),
chains = 2,
silent = 2
)
# Visualise the different monotonic functions
plot_predictions(
mod,
condition = c('x', 'fac', 'fac'),
points = 0.5
)
plot(mod, type = 'smooth', realisations = TRUE)
# First derivatives (on the link scale) should never be
# negative for either factor level
(derivs <- slopes(
mod,
variables = 'x',
by = c('x', 'fac'),
type = 'link'
))
all(derivs$estimate > 0)
}
Run the code above in your browser using DataLab