Given a model (a list of segment formulas), mcp
infers the posterior
distributions of the parameters of each segment as well as the change points
between segments. See more details and worked examples on the mcp website.
All segments must regress on the same xvariable. Change
points are forced to be ordered using truncation of the priors. You can run
fit = mcp(model, sample=FALSE)
to avoid sampling and the need for
data if you just want to get the priors (fit$prior
), the JAGS code
fit$jags_code
, or the R function to simulate data (fit$simulate
).
mcp( model, data = NULL, prior = list(), family = gaussian(), par_x = NULL, sample = "post", cores = 1, chains = 3, iter = 3000, adapt = 1500, inits = NULL, jags_code = NULL )
model  A list of formulas  one for each segment. The first formula
has the format


data  Data.frame or tibble in long format. 
prior  Named list. Names are parameter names (

family  One of 
par_x  String (default: NULL). Only relevant if no segments contains slope (no hint at what x is). Set this, e.g., par_x = "time". 
sample  One of

cores  Positive integer or "all". Number of cores.

chains  Positive integer. Number of chains to run. 
iter  Positive integer. Number of postwarmup samples to draw. The number
of iterations per chain is 
adapt  Positive integer. Also sometimes called "burnin", this is the number of samples used to reach convergence. Set lower for greater speed. Set higher if the chains haven't converged yet or look at tips, tricks, and debugging. 
inits  A list if initial values for the parameters. This can be useful
if a model fails to converge. Read more in 
jags_code  String. Pass JAGS code to 
An mcpfit
object.
Notes on priors:
Order restriction is automatically applied to cp_\* parameters using
truncation (e.g., T(cp_1, )
) so that they are in the correct order on the
xaxis UNLESS you do it yourself. The one exception is for dunif
distributions where you have to do it as above.
In addition to the model parameters, MINX
(minimum xvalue), MAXX
(maximum xvalue), SDX
(etc...), MINY
, MAXY
, and SDY
are also available when you set priors. They are used to set uninformative
default priors.
Use SD when you specify priors for dt, dlogis, etc. JAGS uses precision
but mcp
converts to precision under the hood via the sd_to_prec()
function. So you will see SDs in fit$prior
but precision ($1/SD^2)
in fit$jags_code
# \donttest{ # Define the segments using formulas. A change point is estimated between each formula. model = list( response ~ 1, # Plateau in the first segment (int_1) ~ 0 + time, # Joined slope (time_2) at cp_1 ~ 1 + time # Disjoined slope (int_3, time_3) at cp_2 ) # Fit it. The `ex_demo` dataset is included in mcp. Sample the prior too. # options(mc.cores = 3) # Uncomment to speed up sampling ex_fit = mcp(model, data = ex_demo, sample = "both") # See parameter estimates summary(ex_fit) # Visual inspection of the results plot(ex_fit) # Visualization of model fit/predictions plot_pars(ex_fit) # Parameter distributions pp_check(ex_fit) # Prior/Posterior predictive checks # Test a hypothesis hypothesis(ex_fit, "cp_1 > 10") # Make predictions fitted(ex_fit) predict(ex_fit) predict(ex_fit, newdata = data.frame(time = c(55.545, 80, 132))) # Compare to a oneinterceptonly model (no change points) with default prior model_null = list(response ~ 1) fit_null = mcp(model_null, data = ex_demo, par_x = "time") # fit another model here ex_fit$loo = loo(ex_fit) fit_null$loo = loo(fit_null) loo::loo_compare(ex_fit$loo, fit_null$loo) # Inspect the prior. Useful for prior predictive checks. summary(ex_fit, prior = TRUE) plot(ex_fit, prior = TRUE) # Show all priors. Default priors are added where you don't provide any print(ex_fit$prior) # Set priors and rerun prior = list( int_1 = 15, time_2 = "dt(0, 2, 1) T(0, )", # tdist slope. Truncated to positive. cp_2 = "dunif(cp_1, 80)", # change point to segment 2 > cp_1 and < 80. int_3 = "int_1" # Shared intercept between segment 1 and 3 ) fit3 = mcp(model, data = ex_demo, prior = prior) # Show the JAGS model cat(ex_fit$jags_code) # }