About mcp formulas and models

Source: vignettes/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 format of each segment is generally response ~ changepoint ~ predictors + sigma(variance_predictors), except for the first segment where there is no change point (response ~ predictors + sigma(variance_predictors)). Here is a simple and a complex formula

  • The response is just the name of a column in your data.
  • The change point can be 1 for a population-level change point or 1 + (1|group) varying change points sampled around the population-level change point.
  • The predictors can be 0 (no change in intercept), 1 (change in intercept), and any column in your data which you want to model a slope on.
  • The variance predictors has it’s own article, but they follow the same rules as the predictors.

For convenience, you can omit the response and the change point in segment 2+, in which case the former response and an intercept-only (as opposed to random/varying) change point is 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, 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 the change point 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:

  • int_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). Note that for family = gaussian() and other families with an SD residual, there will always be an sigma_1, i.e., the common sigma 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(
  # int_1
  score ~ 1,

  # cp_1, cp_1_sd, cp_1_id, x_2
  1 + (1|id) ~ 0 + year,

  # cp_2, cp_2_sd, cp_2_condition, int_2, x_2
  1 + (1|condition) ~ 1 + rel(year)
)

# Intepret, but do not sample.
fit = mcp(model, sample = FALSE)
str(fit$pars, vec.len = 99)  # Compact display
## List of 9
##  $ x         : chr "year"
##  $ y         : chr "score"
##  $ trials    : NULL
##  $ weights   : NULL
##  $ varying   : chr [1:2] "cp_1_id" "cp_2_condition"
##  $ sigma     : chr "sigma_1"
##  $ arma      : chr(0) 
##  $ reg       : chr [1:8] "cp_1" "cp_1_sd" "cp_2" "cp_2_sd" "int_1" "int_3" "year_2" "year_3"
##  $ population: chr [1:9] "cp_1" "cp_1_sd" "cp_2" "cp_2_sd" "int_1" "int_3" "year_2" "year_3" "sigma_1"
##  - 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
int_1 = 5
int_2 = 10

# Ifelse version
y_ifelse = ifelse(x <= cp_1, yes = int_1, no = int_2)

# Indicator equivalent using dummy helpers
cp_0 = -Inf
cp_2 = Inf
y_indicator = (x > cp_0) * (x <= cp_1) * int_1 +  # Between cp_0 and cp_1
              (x > cp_1) * (x <= cp_2) * int_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) * int_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, sample = FALSE, par_x = "x")
fit$jags_code
## model {
## 
##   # Priors for population-level effects
##   cp_0 = MINX  # mcp helper value.
##   cp_2 = MAXX  # mcp helper value.
## 
##   cp_1 ~ dunif(MINX, MAXX)
##   int_1 ~ dt(0, 1/(3*SDY)^2, 3) 
##   int_2 ~ dt(0, 1/(3*SDY)^2, 3) 
##   sigma_1 ~ dnorm(0, 1/(SDY)^2) T(0, )
## 
## 
##   # Model and likelihood
##   for (i_ in 1:length(x)) {
##     X_1_[i_] = min(x[i_], cp_1)
##     X_2_[i_] = min(x[i_], cp_2) - cp_1
##     
##     # Fitted value
##     y_[i_] = 
##     
##       # Segment 1: y ~ 1
##       (x[i_] >= cp_0) * (x[i_] < cp_1) * int_1 + 
##     
##       # Segment 2: y ~ 1 ~ 1
##       (x[i_] >= cp_1) * int_2 
##     
##     # Fitted standard deviation
##     sigma_[i_] = max(0, 
##       (x[i_] >= cp_0) * sigma_1 )
## 
##     # Likelihood and log-density for family = gaussian()
##     y[i_] ~ dnorm((y_[i_]), 1 / sigma_[i_]^2)  # SD as precision
##     loglik_[i_] = logdensity.norm(y[i_], (y_[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) * int_2")

Let us see this in action:

model = list(y ~ 0 + x,
                ~ 0 + x)
fit = mcp(model, sample = FALSE)
fit$jags_code
## model {
## 
##   # Priors for population-level effects
##   cp_0 = MINX  # mcp helper value.
##   cp_2 = MAXX  # mcp helper value.
## 
##   cp_1 ~ dunif(MINX, MAXX)
##   x_1 ~ dt(0, 1/(SDY/(MAXX-MINX))^2, 3) 
##   x_2 ~ dt(0, 1/(SDY/(MAXX-MINX))^2, 3) 
##   sigma_1 ~ dnorm(0, 1/(SDY)^2) T(0, )
## 
## 
##   # Model and likelihood
##   for (i_ in 1:length(x)) {
##     X_1_[i_] = min(x[i_], cp_1)
##     X_2_[i_] = min(x[i_], cp_2) - cp_1
##     
##     # Fitted value
##     y_[i_] = 
##     
##       # Segment 1: y ~ 0 + x
##       (x[i_] >= cp_0) * x_1 * X_1_[i_] + 
##     
##       # Segment 2: y ~ 1 ~ 0 + x
##       (x[i_] >= cp_1) * x_2 * X_2_[i_] 
##     
##     # Fitted standard deviation
##     sigma_[i_] = max(0, 
##       (x[i_] >= cp_0) * sigma_1 )
## 
##     # Likelihood and log-density for family = gaussian()
##     y[i_] ~ dnorm((y_[i_]), 1 / sigma_[i_]^2)  # SD as precision
##     loglik_[i_] = logdensity.norm(y[i_], (y_[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

mcp allows for specifying relative intercepts and change points through rel(). 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
int_1 = 5
int_2 = 3  # let's model this as relative

# Indicator version.
cp_0 = -Inf
y_indicator = (x > cp_0) * int_1 +  # Note: no (x < cp_1)
              (x > cp_1) * (int_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.