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 x-variable. 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 post-warmup 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
x-axis 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 x-value), MAXX
(maximum x-value), 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 one-intercept-only 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 re-run
prior = list(
int_1 = 15,
time_2 = "dt(0, 2, 1) T(0, )", # t-dist 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.