mcp  R Documentation 
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 draws from each chain.
The total number of draws 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 
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
An mcpfit
object.
Jonas Kristoffer Lindeløv jonas@lindeloev.dk
get_segment_table
# 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 and sample the prior too.
# options(mc.cores = 3) # Uncomment to speed up sampling
ex = mcp_example("demo") # Simulated data example
demo_fit = mcp(model, data = ex$data, sample = "both")
# See parameter estimates
summary(demo_fit)
# Visual inspection of the results
plot(demo_fit) # Visualization of model fit/predictions
plot_pars(demo_fit) # Parameter distributions
pp_check(demo_fit) # Prior/Posterior predictive checks
# Test a hypothesis
hypothesis(demo_fit, "cp_1 > 10")
# Make predictions
fitted(demo_fit)
predict(demo_fit)
predict(demo_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$data, par_x = "time") # fit another model here
demo_fit$loo = loo(demo_fit)
fit_null$loo = loo(fit_null)
loo::loo_compare(demo_fit$loo, fit_null$loo)
# Inspect the prior. Useful for prior predictive checks.
summary(demo_fit, prior = TRUE)
plot(demo_fit, prior = TRUE)
# Show all priors. Default priors are added where you don't provide any
print(demo_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$data, prior = prior)
# Show the JAGS model
demo_fit$jags_code
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.