Skip to contents

About mcp formulas and models

Source: vignettes/formulas.Rmd
formulas.Rmd

mcp takes a list of formulas, and defines the change point as the point on the x-axis where the data shifts from being generated by one formula to the next. So with N formulas, you have N - 1 change points. A list with just one formula thus correspond to normal regression with 0 change points.

The formulas are called “segments” because they divide the (observed) x-axis into N segments.

Formula format for segments

The three formula parts are called response, cp, and predictor. The general format is response ~ cp ~ predictor, except for the first segment, which has no change point and therefore uses response ~ predictor. Here cp names a formula part; it is not a literal variable name. Distributional terms such as sigma() live in the predictor part.

  • The response is just the name of a column in your data.
  • The cp part can be 1 for a population-level change point or 1 + (1|group) for group-level change-point deviations.
  • The predictor part can contain 0 (no intercept), 1 (an intercept), columns from your data whose effects you want to model, and group-level effects, such as (1|group) and (factor||group).
  • Distributional predictors, such as those inside sigma(), have their own article, but follow the same predictor rules, including group-level effects.

Predictor group-level effects use familiar lme4/brms syntax. For example, y ~ 1 + (1|id) gives each id a deviation from the population intercept. Use || for several independent group coefficients: (1 + x||id) adds an intercept and numeric slope, while (factor||id) adds an intercept and treatment-coded factor contrasts. (0 + factor||id) instead adds one independent coefficient for each factor level. Each coefficient has its own population-level SD.

A group-level term carries into later segments until another term for the same distributional parameter and grouping factor replaces the whole block; use (0|id) to turn it off. Correlated multi-coefficient | terms, such as (1 + x|id), and group-level terms inside ar() or ma() are not yet supported.

Transformations in formulas are evaluated on the original predictor values, as in lm() and glm(). The change-point predictor is the exception: bare par_x and polynomial bases such as I(par_x^2) use distance from the segment onset so that they can describe joined segment shapes.

For convenience, you can omit the response and cp parts in segment 2+. The previous response and an intercept-only population-level change point are then assumed. When you call summary(fit), it will show the explicit representation. Let us see this in action for this model where we predict score as a function of time in three segments, i.e., with two change points:

library(mcp)
model = list(
  score ~ 1,  # intercept
  score ~ 1 ~ 0 + time,  # joined slope
  ~ time  # disjoined slope. "score ~ 1 ~ 1 + time" is implicit here.
)

# Interpret, but do not sample.
fit = mcp(model, data = data.frame(score = 1:10, time = 1:10), sample = FALSE)
summary(fit)
## Family: gaussian(link = 'identity')
## Segments:
##   1: score ~ 1
##   2: score ~ 1 ~ 0 + time
##   3: score ~ 1 ~ time
## 
## No samples. Nothing to summarise.

Notice how it added the response and cp parts to the last segment?

mcp is heavily inspired by brms which again is inspired by lme4::lmer. Here is a bit of history on that.

Parameter names

mcp automatically assigns names to the parameters in the format type_i where i is the segment number. Specifically:

  • Intercept_i is the intercept in the ith segment.
  • year_i is the slope in on the data column year in the ith segment. x_i is the slope on the data column x in the ith segment. The slope takes name after the data it is regressed on.
  • cp_i is the ith change point. Notice that cp_i is specified in segment i + 1. cp_1 occurs when there are two segments, and cp_2 when there are three segments, etc. OBS: future versions may start at cp_2..
  • cp_i_group is the varying deviations from cp_i. See varying change points in mcp.
  • cp_i_sd is the population-level standard deviation of the varying effects.
  • sigma_* are variance parameters about which you can read more here). With an explicit sigma() formula these are log-SD coefficients. If no sigma() is supplied, the implicit sigma_1 is instead the common response-scale SD initiated in the first segment.
  • arj_i are autocorrelation coefficients of order j for segment i (read more here).

These parameter names are saved in fit$pars. Let us specify a somewhat complex model to show off some parameter names:

model = list(
  # Intercept_1
  score ~ 1,
  
  # cp_1, cp_1_sd, cp_1_id, year_2
  1 + (1|id) ~ 0 + year,
  
  # cp_2, cp_2_sd, cp_2_condition, Intercept_3, year_3
  1 + (1|condition) ~ 1 + year
)

# Interpret, but do not sample.
fit = mcp(model, data = data.frame(score = 1:10, year = 1:10, id = 1:10, condition = rep(c("A", "B"), 5)), sample = FALSE)
str(fit$pars, vec.len = 99)  # Compact display
## List of 10
##  $ x         : chr "year"
##  $ y         : chr "score"
##  $ cp        : chr [1:2] "cp_1" "cp_2"
##  $ fixed     : chr [1:4] "Intercept_1" "year_2" "Intercept_3" "year_3"
##  $ population: chr [1:9] "cp_1" "cp_1_sd" "cp_2" "cp_2_sd" "Intercept_1" "year_2" "Intercept_3" "year_3" "sigma_1"
##  $ varying   : chr [1:2] "cp_1_id" "cp_2_condition"
##  $ sigma     : chr "sigma_1"
##  $ arma      : chr(0) 
##  $ trials    : NULL
##  $ weights   : NULL
##  - attr(*, "class")= chr [1:2] "mcplist" "list"

Modeling intercept change points

A change point is simply like an ifelse statement or multiplying with indicators (0s and 1s):

# Model parameters
x = 1:20
cp_1 = 12
Intercept_1 = 5
Intercept_2 = 10

# Ifelse version
y_ifelse = ifelse(x <= cp_1, yes = Intercept_1, no = Intercept_2)

# Indicator equivalent using dummy helpers
cp_0 = -Inf
cp_2 = Inf
y_indicator = (x > cp_0) * (x <= cp_1) * Intercept_1 +  # Between cp_0 and cp_1
              (x > cp_1) * (x <= cp_2) * Intercept_2  # Between cp_1 and cp_2

# Show it
par(mfrow = c(1,2))
plot(x, y_ifelse, main = "ifelse(x <= cp_1)")
plot(x, y_indicator, main = "(x > cp_1) * Intercept_2")

The magic of (Bayesian) MCMC sampling is that it can actually infer the change point from this simple formulation. We let mcp write the JAGS code for this simple two-plateaus model and see how it uses the indicator formulation of change points:

model = list(y ~ 1, ~ 1)
fit = mcp(model, data = data.frame(y = 1:10, x = 1:10), sample = FALSE, par_x = "x")
fit$jags_code
## model {
##   # mcp helper values
##   cp_0 = CONST1_
##   cp_2 = CONST2_
## 
##   # Priors for population-level effects
##   cp_1 ~ dunif(CONST1_, CONST2_)  # Within the observed change-point span
##   Intercept_1 ~ dt(5.5, 1/(3.7)^2, 3)   # Robustly centered mean intercept with a minimum scale of 2.5
##   Intercept_2 ~ dt(5.5, 1/(3.7)^2, 3)   # Robustly centered mean intercept with a minimum scale of 2.5
##   sigma_1 ~ dt(0, 1/(3.7)^2, 3) T(0,)  # Positive residual SD calibrated on the response scale
## 
##   # Model and likelihood
##   for (i_ in 1:length(x)) {
##     # par_x local to each segment
##     x_local_1_[i_] = min(x[i_], cp_1)
##     x_local_2_[i_] = min(x[i_], cp_2) - cp_1
##     
##     # Formula for mu
##     link_mu_[i_] =
##     
##       # Segment 1: y ~ 1
##       (x[i_] >= cp_0) * (x[i_] < cp_1) * inprod(rhs_matrix_[i_, c(1)], c(Intercept_1)) * 1 + 
##     
##       # Segment 2: y ~ 1 ~ 1
##       (x[i_] >= cp_1) * inprod(rhs_matrix_[i_, c(2)], c(Intercept_2)) * 1
##     
##     # Formula for sigma
##     link_sigma_[i_] =
##     
##       # Segment 1: y ~ 1
##       (x[i_] >= cp_0) * inprod(rhs_matrix_[i_, c(3)], c(sigma_1)) * 1
## 
##     # Likelihood and log-density for family = gaussian()
##     mu_[i_] = link_mu_[i_]
##     sigma_[i_] = max(1e-03, link_sigma_[i_])
##     y[i_] ~ dnorm(mu_[i_], 1 / sigma_[i_]^2)  # SD as precision
##   }
## }

Look at the section called # Fitted value which is the (automatically generated) model that was discussed above. Some unnecessary stuff is added to segment 1 just because it makes the code easier to generate. (x[i_] >= cp_0 when cp_0 is the smallest value of x is, of course, always true).

Modeling slope change points

We can use the same principle to model change points on slopes. However, we have to “take off” where the previous slope left us on the y-axis. That is, we have to regard whatever y-value the previous segment ended with as a kind of intercept-at-x=0 in the frame of the new segment. The intercept of segment 2 is cp_1 * slope_1 and the slope in segment 2 is x * (slope_2 - cp_1).

# Model parameters
x = 1:20
cp_1 = 12
slope_1 = 2
slope_2 = -1

# Ifelse version
y_ifelse = ifelse(x <= cp_1, 
            yes = slope_1 * x,
            no = cp_1 * slope_1 + slope_2 * (x - cp_1))

# Indicator version. pmin() is a vectorized min()
cp_0 = -Inf
y_indicator = (x > cp_0) * slope_1 * pmin(x, cp_1) + 
              (x > cp_1) * slope_2 * (x - cp_1)

# Show it
par(mfrow = c(1,2))
plot(x, y_ifelse, main = "ifelse(x <= cp_1)")
plot(x, y_indicator, main = "(x > cp_1) * Intercept_2")

Let us see this in action:

model = list(y ~ 0 + x,
                ~ 0 + x)
fit = mcp(model, data = data.frame(y = 1:10, x = 1:10), sample = FALSE)
fit$jags_code
## model {
##   # mcp helper values
##   cp_0 = CONST1_
##   cp_2 = CONST2_
## 
##   # Priors for population-level effects
##   cp_1 ~ dunif(CONST1_, CONST2_)  # Within the observed change-point span
##   x_1 ~ dt(0, 1/(0.8222222)^2, 3)   # Regularizing mean coefficient scaled to a reference predictor change
##   x_2 ~ dt(0, 1/(0.8222222)^2, 3)   # Regularizing mean coefficient scaled to a reference predictor change
##   sigma_1 ~ dt(0, 1/(3.7)^2, 3) T(0,)  # Positive residual SD calibrated on the response scale
## 
##   # Model and likelihood
##   for (i_ in 1:length(x)) {
##     # par_x local to each segment
##     x_local_1_[i_] = min(x[i_], cp_1)
##     x_local_2_[i_] = min(x[i_], cp_2) - cp_1
##     
##     # Formula for mu
##     link_mu_[i_] =
##     
##       # Segment 1: y ~ 0 + x
##       (x[i_] >= cp_0) * inprod(rhs_matrix_[i_, c(1)], c(x_1)) * x_local_1_[i_] + 
##     
##       # Segment 2: y ~ 1 ~ 0 + x
##       (x[i_] >= cp_1) * inprod(rhs_matrix_[i_, c(2)], c(x_2)) * x_local_2_[i_]
##     
##     # Formula for sigma
##     link_sigma_[i_] =
##     
##       # Segment 1: y ~ 0 + x
##       (x[i_] >= cp_0) * inprod(rhs_matrix_[i_, c(3)], c(sigma_1)) * 1
## 
##     # Likelihood and log-density for family = gaussian()
##     mu_[i_] = link_mu_[i_]
##     sigma_[i_] = max(1e-03, link_sigma_[i_])
##     y[i_] ~ dnorm(mu_[i_], 1 / sigma_[i_]^2)  # SD as precision
##   }
## }

Again, look at the #Fitted value to see the indicator-version in action. And again, mcp adds something about cp_0 = -Inf and cp_2 = Inf, just for internal convenience.

You will find the exact same formula for y_ = ... if you do print(fit$simulate), though this function contains a whole lot of other stuff too.

Modeling relative slopes and intercepts

Relative parameters and changepoints can be modeled by specifying prior relationships or parameter offsets. Relative slopes are easy: just replace x_2 with x_1 + x_2. You could do the same if all segments are intercept-only. However, if the previous segment had a slope, we want the intercept to be relative to where that “ended”. The mcp solution is to only “turn off” the “hanging intercept” from that slope’s ending (pmin(x, cp_i)) when the model encounters an absolute intercept. An indicator does this.

# Model parameters
x = 1:20
cp_1 = 12
Intercept_1 = 5
Intercept_2 = 3  # let's model this as relative

# Indicator version.
cp_0 = -Inf
y_indicator = (x > cp_0) * Intercept_1 +  # Note: no (x < cp_1)
              (x > cp_1) * (Intercept_2)

# Plot it
plot(x, y_indicator, main = "Relative intercept")

You can look at fit$jags_code and fit$simulate() to see this in action.