Movatterモバイル変換

Title:

Regression with Multiple Change Points

Version:

0.3.4

Date:

2024-03-14

URL:

https://lindeloev.github.io/mcp/

BugReports:

https://github.com/lindeloev/mcp/issues

Description:

Flexible and informed regression with Multiple Change Points. 'mcp' can infer change points in means, variances, autocorrelation structure, and any combination of these, as well as the parameters of the segments in between. All parameters are estimated with uncertainty and prediction intervals are supported - also near the change points. 'mcp' supports hypothesis testing via Savage-Dickey density ratio, posterior contrasts, and cross-validation. 'mcp' is described in Lindeløv (submitted) <doi:10.31219/osf.io/fzqxv> and generalizes the approach described in Carlin, Gelfand, & Smith (1992) <doi:10.2307/2347570> and Stephens (1994) <doi:10.2307/2986119>.

License:

GPL-2

Encoding:

UTF-8

Language:

en-US

LazyData:

true

RoxygenNote:

7.3.1

Depends:

R (≥ 3.5.0)

Imports:

parallel, future (≥ 1.16), future.apply (≥ 1.4), rjags (≥4.9), coda (≥ 0.19.3), loo (≥ 2.1.0), bayesplot (≥ 1.7.0),tidybayes (≥ 3.0.0), dplyr (≥ 1.1.1), magrittr (≥ 1.5),tidyr (≥ 1.0.0), tidyselect (≥ 0.2.5), tibble (≥ 2.1.3),stringr (≥ 1.4.0), ggplot2 (≥ 3.2.1), patchwork (≥ 1.0.0),stats, rlang (≥ 0.4.1)

Suggests:

hexbin, testthat (≥ 3.1.0), purrr (≥ 0.3.4), knitr,rmarkdown

NeedsCompilation:

Packaged:

2024-03-17 19:39:34 UTC; jonas

Author:

Jonas Kristoffer Lindeløv

[aut, cre]

Maintainer:

Jonas Kristoffer Lindeløv <jonas@lindeloev.dk>

Repository:

CRAN

Date/Publication:

2024-03-17 20:10:02 UTC

mcp: Regression with Multiple Change Points

Description

logo

Flexible and informed regression with Multiple Change Points. 'mcp' can infer change points in means, variances, autocorrelation structure, and any combination of these, as well as the parameters of the segments in between. All parameters are estimated with uncertainty and prediction intervals are supported - also near the change points. 'mcp' supports hypothesis testing via Savage-Dickey density ratio, posterior contrasts, and cross-validation. 'mcp' is described in Lindeløv (submitted)doi:10.31219/osf.io/fzqxv and generalizes the approach described in Carlin, Gelfand, & Smith (1992)doi:10.2307/2347570 and Stephens (1994)doi:10.2307/2986119.

Author(s)

Maintainer: Jonas Kristoffer Lindeløvjonas@lindeloev.dk (ORCID)

Bernoulli family for mcp

Description

Bernoulli family for mcp

Usage

bernoulli(link = "logit")

Arguments

link

Link function.

Checks if all terms are in the data

Description

Checks if all terms are in the data

Usage

check_terms_in_data(form, data, i, n_terms = NULL)

Arguments

form

Formula or character (tilde will be prefixed if it isn't already)

data

A data.frame or tibble

i

The segment number (integer)

n_terms

Int >= 1. Number of expected terms. Will raise error if it doesn't match.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Compute information criteria for model comparison

Description

Takes anmcpfit as input and computes information criteria using loo orWAIC. Compare models usingloo_compare andloo_model_weights.more inloo.

Usage

criterion(fit, criterion = "loo", ...)## S3 method for class 'mcpfit'loo(x, ...)## S3 method for class 'mcpfit'waic(x, ...)

Arguments

fit

Anmcpfit object.

criterion

One of"loo" (callsloo) or"waic" (callswaic).

...

Currently ignored

x

Anmcpfit object.

Value

aloo orpsis_loo object.

Functions

loo(mcpfit): Computes loo on mcpfit objects
waic(mcpfit): Computes WAIC on mcpfit objects

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Examples

# Define two models and sample them# options(mc.cores = 3)  # Speed up samplingex = mcp_example("intercepts")  # Get some simulated data.model1 = list(y ~ 1 + x, ~ 1)model2 = list(y ~ 1 + x)  # Without a change pointfit1 = mcp(model1, ex$data)fit2 = mcp(model2, ex$data)# Compute LOO for each and compare (works for waic(fit) too)fit1$loo = loo(fit1)fit2$loo = loo(fit2)loo::loo_compare(fit1$loo, fit2$loo)

Cumulative pasting of character columns

Description

Cumulative pasting of character columns

Usage

cumpaste(x, .sep = " ")

Arguments

x

A column

.sep

A character to append between pastes

Value

string.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk but Inspired byhttps://stackoverflow.com/questions/24862046/cumulatively-paste-concatenate-values-grouped-by-another-variable

Example`mcpfit` for examples

Description

This was generated usingmcp_examples("demo", sample = TRUE).

Usage

demo_fit

Format

Anmcpfit object.

Exponential family for mcp

Description

Exponential family for mcp

Usage

exponential(link = "identity")

Arguments

link

Link function (Character).

Expected Values from the Posterior Predictive Distribution

Description

Expected Values from the Posterior Predictive Distribution

Usage

## S3 method for class 'mcpfit'fitted(  object,  newdata = NULL,  summary = TRUE,  probs = TRUE,  rate = TRUE,  prior = FALSE,  which_y = "ct",  varying = TRUE,  arma = TRUE,  nsamples = NULL,  samples_format = "tidy",  scale = "response",  ...)

Arguments

object

Anmcpfit object.

newdata

Atibble or adata.frame containing predictors in the model. IfNULL (default),the original data is used.

summary

Summarise at each x-value

probs

Vector of quantiles. Only in effect whensummary == TRUE.

rate

Boolean. For binomial models, plot on raw data (rate = FALSE) orresponse divided by number of trials (rate = TRUE). If FALSE, linearinterpolation on trial number is used to infer trials at a particular x.

prior

TRUE/FALSE. Plot using prior samples? Useful formcp(..., sample = "both")

which_y

What to plot on the y-axis. One of

"ct": The central tendency which is often the mean after applying thelink function.
"sigma": The variance
"ar1","ar2", etc. depending on which order of the autoregressiveeffects you want to plot.

varying

One of:

TRUE All varying effects (fit$pars$varying).
FALSE No varying effects (c()).
Character vector: Only include specified varying parameters - seefit$pars$varying.

arma

Whether to include autoregressive effects.

TRUE Compute autoregressive residuals. Requires the response variable innewdata.
FALSE Disregard the autoregressive effects. Forfamily = gaussian(),predict() just usesigma for residuals.

nsamples

Integer orNULL. Number of samples to return/summarise.If there are varying effects, this is the number of samples from each varying group.NULL means "all". Ignored if both areFALSE. More samples trade speed for accuracy.

samples_format

One of "tidy" or "matrix". Controls the output format whensummary == FALSE.See more under "value"

scale

One of

"response": return on the observed scale, i.e., after applying the inverse link function.
"linear": return on the parameter scale (where the linear trends are modelled).

...

Currently unused

Value

Ifsummary = TRUE: Atibble with the posterior mean for each row innewdata,Ifnewdata isNULL, the data infit$data is used.
Ifsummary = FALSE andsamples_format = "tidy": Atidybayestibble with all the posteriorsamples (Ns) evaluated at each row innewdata (Nn), i.e., with⁠Ns x Nn⁠ rows. If there arevarying effects, the returned data is expanded with the relevant levels for each row.
The return columns are:
- Predictors fromnewdata.
- Sample descriptors: ".chain", ".iter", ".draw" (see thetidybayes package for more), and "data_row" (newdata rownumber)
- Sample values: one column for each parameter in the model.
- The estimate. Either "predict" or "fitted", i.e., the name of thetype argument.
Ifsummary = FALSE andsamples_format = "matrix": AnN_draws Xnrows(newdata) matrix with fitted/predictedvalues (depending ontype). This format is used bybrms and it's useful asyrep in⁠bayesplot::ppc_*⁠ functions.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Examples

fitted(demo_fit)fitted(demo_fit, probs = c(0.1, 0.5, 0.9))  # With median and 80% credible interval.fitted(demo_fit, summary = FALSE)  # Samples instead of summary.fitted(demo_fit,       newdata = data.frame(time = c(-5, 20, 300)),  # New data       probs = c(0.025, 0.5, 0.975))

Format code with one or multiple terms

Description

Take a value like "a + b" and(1) replace it with NA if na_col == NA.(2) Change to "(a + b)" if there is a "+"(3) Return itself otherwise, e.g., "a" –> "a".

Usage

format_code(col, na_col)

Arguments

col

A column

na_col

If this column is NA, return NA

Value

string

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Density geom for`plot.mcpfit()`

Description

Density geom forplot.mcpfit()

Usage

geom_cp_density(fit, facet_by, limits_y)

Arguments

fit

Anmcpfit object

facet_by

NULL or a a string, like⁠plot.mcpfit(..., facet_by = "id").⁠

Value

Aggplot2::stat_density geom representing the change point densities.

Return a geom_line representing the quantiles

Description

Called byplot.mcpfit.

Usage

geom_quantiles(samples, quantiles, xvar, yvar, facet_by, ...)

Arguments

samples

A tidybayes tibble

quantiles

Vector of quantiles (0.0 to 1.0)

xvar

An rlang::sym() with the name of the x-col insamples

yvar

An rlang::sym() with the name of the response col insamples

facet_by

String. Name of a varying group.

...

Arguments passed to geom_line

Value

Aggplot2::geom_line object.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Call`get_formula_str` for multiple ytypes and paste strings

Description

Currently used to differentiate between the JAGS model (use all) and thefit$simulate model (do not include arma).

Usage

get_all_formulas(ST, prior, par_x, ytypes = c("ct", "sigma", "arma"))

Arguments

ST

Tibble. Returned byget_segment_table.

prior

Named list. Names are parameter names (cp_i,int_i,xvar_i,'sigma“) and the values are either

A JAGS distribution (e.g.,int_1 = "dnorm(0, 1) T(0,)") indicating aconventional prior distribution. Uninformative priors based on dataproperties are used where priors are not specified. This ensures goodparameter estimations, but it is a questionable for hypothesis testing.mcp uses SD (not precision) for dnorm, dt, dlogis, etc. Seedetails. Change points are forced to be ordered through the priors usingtruncation, except for uniform priors where the lower bound should begreater than the previous change point,dunif(cp_1, MAXX).
A numerical value (e.g.,int_1 = -2.1) indicating a fixed value.
A model parameter name (e.g.,int_2 = "int_1"), indicating that this parameter is shared -typically between segments. If two varying effects are shared this way,they will need to have the same grouping variable.
A scaled Dirichlet prior is supported for change points if they are all set to⁠cp_i = "dirichlet(N)⁠ whereN is the alpha for this change point andN = 1 is most often used. This prior is less informative about thelocation of the change points than the default uniform prior, but itsamples less efficiently, so you will often need to setiter higher.It is recommended for hypothesis testing and for the estimation of morethan 5 change points.Read more.

par_x

String (default: NULL). Only relevant if no segments containsslope (no hint at what x is). Set this, e.g., par_x = "time".

ytypes

A character vector of ytypes to including in model building

Value

A string with JAGS code.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Gets code for ARMA terms, resulting in a "resid_"

Description

Developer note: Ensuring that this can be used in both simulate() and JAGSgot quite messy with a lot of if-statements. It works but some refactoringmay be good in the future.

Usage

get_ar_code(ar_order, family, is_R, xvar, yvar = NA)

Arguments

ar_order

Positive integer. The order of ARMA

family

An mcpfamily object

is_R

Bool. Is this R code (TRUE) or JAGS code (FALSE)?

Value

String with JAGS code for AR.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Extracts the order from ARMA parameter name(s)

Description

If several names are provided (vector), it returns the maximum. Ifpars_armais an empty string, it returns0.

Usage

get_arma_order(pars_arma)

Arguments

pars_arma

Character vector

Value

integer

Compute the density at a specific point.

Description

Used inhypothesis

Usage

get_density(samples, LHS, value)

Arguments

samples

An mcmc.list

LHS

Expression to compute posterior

value

What value to evaluate the density at

Value

A float

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Get a list of x-coordinates to evaluate fit$simulate at

Description

Solves two problems: if setting the number of points too high, thefunction becomes slow. If setting it too low, the posterior at large intercept-changes at change points look discrete, because they are evaluated at veryfew x in that interval.

Usage

get_eval_at(fit, facet_by)

Arguments

fit

An mcpfit object.

facet_by

String. Name of a varying group.

Details

This function makes a vector of x-values with large spacing in general,but finer resolution at change points.

Value

A vector of x-values to evaluate at.

Build an R formula (as string) given a segment table (ST)

Description

You will need to replace PAR_X for whatever your x-axis observation columnis called. In JAGS typicallyx[i_]. In R justx.

Usage

get_formula_str(ST, par_x, ytype = "ct", init = FALSE)

Arguments

ST

Tibble. Returned byget_segment_table.

par_x

String (default: NULL). Only relevant if no segments containsslope (no hint at what x is). Set this, e.g., par_x = "time".

ytype

One of "ct" (central tendency), "sigma", "ar1" (or another order), or "ma1" (or another order)

init

TRUE/FALSE. Set to TRUE for the first call. Adds segment-relativeX-codings and verbose commenting of one formula

Value

A string with JAGS code.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Adds helper variables for use in`run_jags`

Description

Returns the relevant data columns as a list and add elements with uniquevarying group levels.

Usage

get_jags_data(data, ST, jags_code, sample)

Arguments

data

A tibble

ST

A segment table (tibble), returned byget_segment_table.

jags_code

A string. JAGS model, usually returned bymake_jagscode().

sample

One of

"post": Sample the posterior.
"prior": Sample only the prior. Plots, summaries, etc. willuse the prior. This is useful for prior predictive checks.
"both": Sample both prior and posterior. Plots, summaries, etc.will default to using the posterior. The prior only has effect when doingSavage-Dickey density ratios inhypothesis.
"none" orFALSE: Do not sample. Returns an mcpfitobject without sample. This is useful if you only want to checkprior strings (fit$prior), the JAGS model (fit$jags_code), etc.

Make JAGS code for Multiple Change Point model

Description

Make JAGS code for Multiple Change Point model

Usage

get_jagscode(prior, ST, formula_str, arma_order, family, sample)

Arguments

prior

Named list. Names are parameter names (cp_i,int_i,xvar_i,'sigma“) and the values are either

A JAGS distribution (e.g.,int_1 = "dnorm(0, 1) T(0,)") indicating aconventional prior distribution. Uninformative priors based on dataproperties are used where priors are not specified. This ensures goodparameter estimations, but it is a questionable for hypothesis testing.mcp uses SD (not precision) for dnorm, dt, dlogis, etc. Seedetails. Change points are forced to be ordered through the priors usingtruncation, except for uniform priors where the lower bound should begreater than the previous change point,dunif(cp_1, MAXX).
A numerical value (e.g.,int_1 = -2.1) indicating a fixed value.
A model parameter name (e.g.,int_2 = "int_1"), indicating that this parameter is shared -typically between segments. If two varying effects are shared this way,they will need to have the same grouping variable.
A scaled Dirichlet prior is supported for change points if they are all set to⁠cp_i = "dirichlet(N)⁠ whereN is the alpha for this change point andN = 1 is most often used. This prior is less informative about thelocation of the change points than the default uniform prior, but itsamples less efficiently, so you will often need to setiter higher.It is recommended for hypothesis testing and for the estimation of morethan 5 change points.Read more.

ST

Segment table. Returned byget_segment_table().

formula_str

String. The formula string returned bybuild_formula_str.

arma_order

Positive integer. The autoregressive order.

family

One ofgaussian(),binomial(),bernoulli(), orpoission().Only default link functions are currently supported.

sample

One of

"post": Sample the posterior.
"prior": Sample only the prior. Plots, summaries, etc. willuse the prior. This is useful for prior predictive checks.
"both": Sample both prior and posterior. Plots, summaries, etc.will default to using the posterior. The prior only has effect when doingSavage-Dickey density ratios inhypothesis.
"none" orFALSE: Do not sample. Returns an mcpfitobject without sample. This is useful if you only want to checkprior strings (fit$prior), the JAGS model (fit$jags_code), etc.

Value

String. A JAGS model.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

pp_check for loo statistics

Description

pp_check for loo statistics

Usage

get_ppc_plot(fit, type, y, yrep, nsamples, draws = NULL, ...)

Arguments

type

One ofbayesplot::available_ppc("grouped", invert = TRUE) %>% stringr::str_remove("ppc_")

y

Response vector

yrep

S X N matrix of predicted responses

nsamples

Number of draws. Note that you may want to use all data for summary geoms.e.g.,pp_check(fit, type = "ribbon", nsamples = NULL).

draws

(required for loo-type plots) Indices of draws to use.

...

Arguments passed tobayesplot::ppc_type(y, yrep, ...)

Value

Aggplot2 object returned by⁠tidybayes::ppc_*(y, yrep, ...)⁠.

A string

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Get priors for all parameters in a segment table.

Description

Starts by finding all default priors. Then replace them with user priors.User priors for change points are truncated appropriately using'truncate_prior_cp“, if not done manually by the user already.

Usage

get_prior(ST, family, prior = list())

Arguments

ST

Tibble. A segment table returned byget_segment_table.

family

One ofgaussian(),binomial(),bernoulli(), orpoission().Only default link functions are currently supported.

prior

A list of user-defined priors. Will overwrite the relevantdefault priors.

Value

A named list of strings. The names correspond to the parameter namesand the strings are the JAGS code for the prior (before converting SD toprecision).

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Get JAGS code for a prior

Description

Get JAGS code for a prior

Usage

get_prior_str(prior, i, varying_group = NULL)

Arguments

prior

Named list. Names are parameter names (cp_i,int_i,xvar_i,'sigma“) and the values are either

A JAGS distribution (e.g.,int_1 = "dnorm(0, 1) T(0,)") indicating aconventional prior distribution. Uninformative priors based on dataproperties are used where priors are not specified. This ensures goodparameter estimations, but it is a questionable for hypothesis testing.mcp uses SD (not precision) for dnorm, dt, dlogis, etc. Seedetails. Change points are forced to be ordered through the priors usingtruncation, except for uniform priors where the lower bound should begreater than the previous change point,dunif(cp_1, MAXX).
A numerical value (e.g.,int_1 = -2.1) indicating a fixed value.
A model parameter name (e.g.,int_2 = "int_1"), indicating that this parameter is shared -typically between segments. If two varying effects are shared this way,they will need to have the same grouping variable.
A scaled Dirichlet prior is supported for change points if they are all set to⁠cp_i = "dirichlet(N)⁠ whereN is the alpha for this change point andN = 1 is most often used. This prior is less informative about thelocation of the change points than the default uniform prior, but itsamples less efficiently, so you will often need to setiter higher.It is recommended for hypothesis testing and for the estimation of morethan 5 change points.Read more.

i

The index inprior to get code for

varying_group

String or NULL. Null indicates a population-level prior. String indicates a varying-effects prior (one for each grouplevel).

Value

A string

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Expand samples with quantiles

Description

TO DO: implement usingfitted() andpredict() but avoid double-computing the samples? E.g.:⁠get_quantiles2 = function(fit, quantiles, facet_by = NULL) {⁠⁠fitted(fit, probs = c(0.1, 0.5, 0.9), newdata = data.frame(x = c(11, 50, 100))) %>%⁠⁠tidyr::pivot_longer(tidyselect::starts_with("Q")) %>%⁠dplyr::mutate(quantile = stringr::str_remove(name, "Q") %>% as.numeric() / 100)⁠}⁠

Usage

get_quantiles(samples, quantiles, xvar, yvar, facet_by = NULL)

Arguments

samples

A tidybayes tibble

quantiles

Vector of quantiles (0.0 to 1.0)

xvar

An rlang::sym() with the name of the x-col insamples

yvar

An rlang::sym() with the name of the response col insamples

facet_by

String. Name of a varying group.

Value

A tidybayes long format tibble with the column "quantile"

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Build a table describing a list of segments

Description

Used internally for most mcp functions.

Usage

get_segment_table(model, data = NULL, family = gaussian(), par_x = NULL)

Arguments

model

A list of formulas - one for each segment. The first formulahas the formatresponse ~ predictors while the following formulashave the formatresponse ~ changepoint ~ predictors. The responseand change points can be omitted (changepoint ~ predictors assumes sameresponse.~ predictors assumes an intercept-only change point). Thefollowing can be modeled:

Regular formulas: e.g.,~ 1 + x).Read more.
Extended formulas:, e.g.,~ I(x^2) + exp(x) + sin(x).Read more.
Variance: e.g.,~sigma(1) for a simple variance change or~sigma(rel(1) + I(x^2))) for more advanced variance structures.Read more
Autoregression: e.g.,~ar(1) for a simple onset/change in AR(1) or⁠ar(2, 0 + x⁠) for an AR(2) increasing byx.Read more

data

Data.frame or tibble in long format.

family

One ofgaussian(),binomial(),bernoulli(), orpoission().Only default link functions are currently supported.

par_x

String (default: NULL). Only relevant if no segments containsslope (no hint at what x is). Set this, e.g., par_x = "time".

Value

A tibble with one row describing each segment and the corresponding code.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Examples

model = list(  y ~ 1 + x,  1 + (1|id) ~ 1)get_segment_table(model)

Turn formula_str into a proper R function

Description

Turn formula_str into a proper R function

Usage

get_simulate(formula_str, pars, nsegments, family)

Arguments

formula_str

string. Returned byget_formula.

pars

List of user-provided parameters, in the format of fit$pars.

nsegments

Positive integer. Number of segments, typicallynrow(ST).

family

One ofgaussian(),binomial(),bernoulli(), orpoission().Only default link functions are currently supported.

Value

A string with R code for the fit$simulate() function corresponding to the model.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Internal function for summary.mcpfit, fixef.mcpfit, and ranef.mcpfit

Description

Internal function for summary.mcpfit, fixef.mcpfit, and ranef.mcpfit

Usage

get_summary(fit, width, varying = FALSE, prior = FALSE)

Arguments

fit

Anmcpfit' object.

width

Float. The width of the highest posterior density interval(between 0 and 1).

varying

Boolean. Get results for varying (TRUE) or population (FALSE)?

prior

TRUE/FALSE. Summarise prior instead of posterior?

Value

A data.frame with summaries for each model parameter.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Get formula inside a wrapper

Description

Get formula inside a wrapper

Usage

get_term_content(term)

Arguments

term

E.g., "ct(1 + x)", "sigma(0 + rel(x) + I(x^2))", etc.

Value

char formula with the content inside the brackets.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Test hypotheses on mcp objects.

Description

Returns posterior probabilities and Bayes Factors for flexible hypotheses involvingmodel parameters. The documentation for the argumenthypotheses belowshows examples of how to specify hypotheses, andread worked examples on the mcp website.For directional hypotheses,⁠hypothesis`` executes the hypothesis string in a ⁠tidybayes“ environment and summerises the proportion of samples wherethe expression evaluates to TRUE. For equals-hypothesis, a Savage-Dickeyratio is computed. Savage-Dickey requires a prior too, so remembermcp(..., sample = "both"). This function is heavily inspired by the'hypothesis' function from the 'brms' package.

Usage

hypothesis(fit, hypotheses, width = 0.95, digits = 3)

Arguments

fit

Anmcpfit object.

hypotheses

String representation of a logical test involving model parameters.Takes R code that evaluates to TRUE or FALSE in a vectorized way.

Directional hypotheses are specified using <, >, <=, or >=.hypothesisreturns the posterior probability and odds in favor of the stated hypothesis.The odds can be interpreted as a Bayes Factor. For example:

"cp_1 > 30": the first change point is above 30.
"int_1 > int_2": the intercept is greater in segment 1 than 2.
"x_2 - x_1 <= 3": the difference between slope 1 and 2 is lessthan or equal to 3.
"int_1 > -2 & int_1 < 2": int_1 is between -2 and 2 (an interval hypothesis). This can be useful as a Region Of Practical Equivalence test (ROPE).
"cp_1^2 < 30 | (log(x_1) + log(x_2)) > 5": be creative.
"`cp_1_id[1]` > `cp_1_id[2]`": id1 is greater than id2, as estimatedthrough the varying-by-"id" change point in segment 1. Note that``required for varying effects.

Hypotheses can also test equality using the equal sign (=). This runs aSavage-Dickey test, i.e., the proportion by which the probability densityhas increased from the prior to the posterior at a given value. Therefore,it requiresmcp(sample = "both"). There are two requirements:First, there can only be one equal sign, so don't use and (&) or or (|).Second, the point to test has to be on the right, and the variables on the left.

"cp_1 = 30": is the first change point at 30? Or to be more precise:by what factor has the credence in cp_1 = 30 risen/fallen whenconditioning on the data, relative to the prior credence?
"int_1 + int_2 = 0": Is the sum of two intercepts zero?
"`cp_1_id[John]`/`cp_1_id[Erin]` = 2": is the varying changepoint for John (which is relative to 'cp_1“) double that of Erin?

width

Float. The width of the highest posterior density interval(between 0 and 1).

digits

a non-null value for digits specifies the minimum number ofsignificant digits to be printed in values. The default, NULL, usesgetOption("digits"). (For the interpretation for complex numbers see signif.)Non-integer values will be rounded down, and only values greater than orequal to 1 and no greater than 22 are accepted.

Value

A data.frame with a row per hypothesis and the following columns:

hypothesis is the hypothesis; often re-arranged to test against zero.
mean is the posterior mean of the left-hand side of the hypothesis.
lower is the lower bound of the (two-sided) highest-density interval of widthwidth.
upper is the upper bound of ditto.
p Posterior probability.For "=" (Savage-Dickey), it is the BF converted to p.For directional hypotheses, it is the proportion of samples that returns TRUE.
BF Bayes Factor in favor of the hypothesis.For "=" it is the Savage-Dickey density ratio.For directional hypotheses, it is p converted to odds.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Inverse logit function

Description

Inverse logit function

Usage

ilogit(eta)

Arguments

eta

A vector of logits

Value

A vector with same length aseta

Checks if argument is an`mcpfit` object

Description

Checks if argument is anmcpfit object

Usage

is.mcpfit(x)

Arguments

x

AnR object.

Logit function

Description

Logit function

Usage

logit(mu)

Arguments

mu

A vector of probabilities (0.0 to 1.0)

Value

A vector with same length asmu

Internal function to get samples.

Description

Returns posterior samples, if available. If not, then prior samples. If not,then throw an informative error. This is useful for summary and plotting, thatworks on both.

Usage

mcmclist_samples(fit, prior = FALSE, message = TRUE, error = TRUE)

Arguments

fit

Anmcpfit object

prior

TRUE/FALSE. Summarise prior instead of posterior?

message

TRUE: gives a message if returning prior samples. FALSE = no message

error

TRUE: err if there are no samples. FALSE: return NULL

Fit Multiple Linear Segments And Their Change Points

Description

Given a model (a list of segment formulas),mcp infers the posteriordistributions of the parameters of each segment as well as the change pointsbetween segments.See more details and worked examples on the mcp website.All segments must regress on the same x-variable. Changepoints are forced to be ordered using truncation of the priors. You can runfit = mcp(model, sample=FALSE) to avoid sampling and the need fordata if you just want to get the priors (fit$prior), the JAGS codefit$jags_code, or the R function to simulate data (fit$simulate).

Usage

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)

Arguments

model

Regular formulas: e.g.,~ 1 + x).Read more.
Extended formulas:, e.g.,~ I(x^2) + exp(x) + sin(x).Read more.
Variance: e.g.,~sigma(1) for a simple variance change or~sigma(rel(1) + I(x^2))) for more advanced variance structures.Read more
Autoregression: e.g.,~ar(1) for a simple onset/change in AR(1) or⁠ar(2, 0 + x⁠) for an AR(2) increasing byx.Read more

data

Data.frame or tibble in long format.

prior

Named list. Names are parameter names (cp_i,int_i,xvar_i,'sigma“) and the values are either

A JAGS distribution (e.g.,int_1 = "dnorm(0, 1) T(0,)") indicating aconventional prior distribution. Uninformative priors based on dataproperties are used where priors are not specified. This ensures goodparameter estimations, but it is a questionable for hypothesis testing.mcp uses SD (not precision) for dnorm, dt, dlogis, etc. Seedetails. Change points are forced to be ordered through the priors usingtruncation, except for uniform priors where the lower bound should begreater than the previous change point,dunif(cp_1, MAXX).
A numerical value (e.g.,int_1 = -2.1) indicating a fixed value.
A model parameter name (e.g.,int_2 = "int_1"), indicating that this parameter is shared -typically between segments. If two varying effects are shared this way,they will need to have the same grouping variable.
A scaled Dirichlet prior is supported for change points if they are all set to⁠cp_i = "dirichlet(N)⁠ whereN is the alpha for this change point andN = 1 is most often used. This prior is less informative about thelocation of the change points than the default uniform prior, but itsamples less efficiently, so you will often need to setiter higher.It is recommended for hypothesis testing and for the estimation of morethan 5 change points.Read more.

family

One ofgaussian(),binomial(),bernoulli(), orpoission().Only default link functions are currently supported.

par_x

String (default: NULL). Only relevant if no segments containsslope (no hint at what x is). Set this, e.g., par_x = "time".

sample

One of

"post": Sample the posterior.
"prior": Sample only the prior. Plots, summaries, etc. willuse the prior. This is useful for prior predictive checks.
"both": Sample both prior and posterior. Plots, summaries, etc.will default to using the posterior. The prior only has effect when doingSavage-Dickey density ratios inhypothesis.
"none" orFALSE: Do not sample. Returns an mcpfitobject without sample. This is useful if you only want to checkprior strings (fit$prior), the JAGS model (fit$jags_code), etc.

cores

Positive integer or "all". Number of cores.

1: serial sampling.options(mc.cores = 3) will dominatecores = 1but not larger values ofcores.
⁠>1⁠: parallel sampling on this number of cores. Ideally setchainsto the same value. Note:cores > 1 takes a few extra seconds the firsttime it's called but subsequent calls will start sampling immediately.
"all": use all cores but one and setschains to the same value. This isa convenient way to maximally use your computer's power.

chains

Positive integer. Number of chains to run.

iter

Positive integer. Number of post-warmup draws from each chain.The total number of draws isiter * chains.

adapt

Positive integer. Also sometimes called "burnin", this is thenumber of samples used to reach convergence. Set lower for greater speed.Set higher if the chains haven't converged yet or look attips, tricks, and debugging.

inits

A list if initial values for the parameters. This can be usefulif a model fails to converge. Read more injags.model.Defaults toNULL, i.e., no inits.

jags_code

String. Pass JAGS code tomcp to use directly. This is useful ifyou want to tweak the code infit$jags_code and run it within themcpframework.

Details

Notes on priors:

Order restriction is automatically applied to cp_\* parameters usingtruncation (e.g.,T(cp_1, )) so that they are in the correct order on thex-axis UNLESS you do it yourself. The one exception is for dunifdistributions 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, andSDYare also available when you set priors. They are used to set uninformativedefault priors.
Use SD when you specify priors for dt, dlogis, etc. JAGS uses precisionbutmcp converts to precision under the hood via the sd_to_prec()function. So you will see SDs infit$prior but precision ($1/SD^2)infit$jags_code

Value

Anmcpfit object.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Examples

# 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 samplingex = mcp_example("demo")  # Simulated data exampledemo_fit = mcp(model, data = ex$data, sample = "both")# See parameter estimatessummary(demo_fit)# Visual inspection of the resultsplot(demo_fit)  # Visualization of model fit/predictionsplot_pars(demo_fit)  # Parameter distributionspp_check(demo_fit)  # Prior/Posterior predictive checks# Test a hypothesishypothesis(demo_fit, "cp_1 > 10")# Make predictionsfitted(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 priormodel_null = list(response ~ 1)fit_null = mcp(model_null, data = ex$data, par_x = "time")  # fit another model heredemo_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 anyprint(demo_fit$prior)# Set priors and re-runprior = 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 modeldemo_fit$jags_code

Get example models and data

Description

Get example models and data

Usage

mcp_example(name, sample = FALSE)

Arguments

name

Name of the example. One of:

"demo": Two change points between intercepts and joined/disjoined slopes.
"ar": One change point in autoregressive residuals.
"binomial": Binomial with two change points. Much like"demo" on a logit scale.
"intercepts": An intercept-only change point.
rel_prior: Relative parameterization and informative priors.
"quadratic": A change point to a quadratic segment.
"trigonometric": Trigonometric/seasonal data and model.
"varying": Varying / hierarchical change points.
"variance": A change in variance, including a variance slope.

sample

TRUE (runfit = mcp(model, data, ...)) or FALSE.

Value

List with

model: A list of formulas
data: The simulated data
simulated: The parameters used for simulating the data.
fit: anmcpfit ifsample = TRUE,
call: the code to run the above.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Examples

ex = mcp_example("demo")plot(ex$data)  # Plot dataprint(ex$simulated)  # See true parameters used to simulateprint(ex$call)  # See how the data was simulated# Fit the model. Either...fit = mcp(ex$model, ex$data)plot(fit)ex_with_fit = mcp_example("demo", sample = TRUE)plot(ex_with_fit$fit)

Add A family object to store link functions between R and JAGS.

Description

This will make more sense once more link functions / families are added.

Usage

mcpfamily(family)

Arguments

family

A family object, e.g.,binomial(link = "identity").

Class`mcpfit` of models fitted with themcp package

Description

Models fitted with themcp function are represented asanmcpfit object which contains the user input (model, data, family),derived model characteristics (prior, parameter names, and jags code), andthe fit (prior and/or posterior mcmc samples).

Details

Seemethods(class = "mcpfit") for an overview of available methods.

User-provided information (seemcp for more details):

Slots

model: A list of formulas, making up the model.Provided by user. Seemcp for more details.
data: A data frame.Provided by user. Seemcp for more details.
family: Anmcpfamily object.Provided by user. Seemcp for more details.
prior: A named list.Provided by user. Seemcp for more details.
mcmc_post: Anmcmc.list object with posterior samples.
mcmc_prior: Anmcmc.list object with prior samples.
mcmc_loglik: Anmcmc.list object with samples of log-likelihood.
pars: A list of character vectors of model parameter names.
jags_code: A string with jags code. Usecat(fit$jags_code) to show it.
simulate: A method to simulate and predict data.
.other: Information that is used internally by mcp.

Negative binomial for mcp

Description

Parameterized asmu (mean; poisson lambda) andsize (a shape parameter),so you can dornbinom(10, mu = 10, size = 1). Read more in the doc forrnbinom,

Usage

negbinomial(link = "log")

Arguments

link

Link function (Character).

Inverse probit function

Description

Inverse probit function

Usage

phi(eta)

Arguments

eta

A vector of probits

Value

A vector with same length asmu

Plot full fits

Description

Plot prior or posterior model draws on top of data. Useplot_pars toplot individual parameter estimates.

Usage

## S3 method for class 'mcpfit'plot(  x,  facet_by = NULL,  lines = 25,  geom_data = "point",  cp_dens = TRUE,  q_fit = FALSE,  q_predict = FALSE,  rate = TRUE,  prior = FALSE,  which_y = "ct",  arma = TRUE,  nsamples = 2000,  scale = "response",  ...)

Arguments

x

Anmcpfit object

facet_by

String. Name of a varying group.

lines

Positive integer orFALSE. Number of lines (posteriordraws). FALSE orlines = 0 plots no lines. Note that lines always plotfitted values - not predicted. For prediction intervals, see theq_predict argument.

geom_data

String. One of "point", "line" (good for time-series),or FALSE (don not plot).

cp_dens

TRUE/FALSE. Plot posterior densities of the change point(s)?Currently does not respectfacet_by. This will be added in the future.

q_fit

Whether to plot quantiles of the posterior (fitted value).

TRUE Add 2.5% and 97.5% quantiles. Corresponds toq_fit = c(0.025, 0.975).
FALSE No quantiles
A vector of quantiles. For example,quantiles = 0.5plots the median andquantiles = c(0.2, 0.8) plots the 20% and 80%quantiles.

q_predict

Same asq_fit, but for the prediction interval.

rate

prior

TRUE/FALSE. Plot using prior samples? Useful formcp(..., sample = "both")

which_y

What to plot on the y-axis. One of

"ct": The central tendency which is often the mean after applying thelink function.
"sigma": The variance
"ar1","ar2", etc. depending on which order of the autoregressiveeffects you want to plot.

arma

Whether to include autoregressive effects.

TRUE Compute autoregressive residuals. Requires the response variable innewdata.
FALSE Disregard the autoregressive effects. Forfamily = gaussian(),predict() just usesigma for residuals.

nsamples

scale

One of

"response": return on the observed scale, i.e., after applying the inverse link function.
"linear": return on the parameter scale (where the linear trends are modelled).

...

Currently ignored.

Details

plot() usesfit$simulate() on posterior samples. These represent the(joint) posterior distribution.

Value

Aggplot2 object.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Examples

# Typical usage. demo_fit is an mcpfit object.plot(demo_fit)plot(demo_fit, prior = TRUE)  # The priorplot(demo_fit, lines = 0, q_fit = TRUE)  # 95% HDI without linesplot(demo_fit, q_predict = c(0.1, 0.9))  # 80% prediction intervalplot(demo_fit, which_y = "sigma", lines = 100)  # The variance parameter on y# Show a panel for each varying effect# plot(fit, facet_by = "my_column")# Customize plots using regular ggplot2library(ggplot2)plot(demo_fit) + theme_bw(15) + ggtitle("Great plot!")

Plot individual parameters

Description

Plot many types of plots of parameter estimates. See examples for typical usecases.

Usage

plot_pars(  fit,  pars = "population",  regex_pars = character(0),  type = "combo",  ncol = 1,  prior = FALSE)

Arguments

fit

Anmcpfit object.

pars

Character vector. One of:

Vector of parameter names.
"population" plots all population parameters.
"varying" plots all varying effects. To plot a particular varyingeffect, useregex_pars = "^name".

regex_pars

Vector of regular expressions. This will typically just bethe beginning of the parameter name(s), i.e., "^cp_" plots all changepoints, "^my_varying" plots all levels of a particular varying effect, and"^cp_|^my_varying" plots both.

type

String or vector of strings. Calls⁠bayesplot::mcmc_>>type<<()⁠.Common calls are "combo", "trace", and "dens_overlay". Current options include'acf', 'acf_bar', 'areas', 'areas_ridges', 'combo', 'dens', 'dens_chains','dens_overlay', 'hist', 'intervals', 'rank_hist', 'rank_overlay', 'trace','trace_highlight', and 'violin".

ncol

Number of columns in plot. This is useful when you have manyparameters and only one plottype.

prior

TRUE/FALSE. Plot using prior samples? Useful formcp(..., sample = "both")

Details

For othertype, it callsbayesplot::mcmc_type(). Use thesedirectly onfit$mcmc_post orfit$mcmc_prior if you want finercontrol of plotting, e.g.,bayesplot::mcmc_dens(fit$mcmc_post). Thereare also a number of useful plots in thecoda package, i.e.,coda::gelman.plot(fit$mcmc_post) andcoda::crosscorr.plot(fit$mcmc_post)

In any case, if you see a few erratic lines or parameter estimates, this isa sign that you may want to increase argument 'adapt' and 'iter' inmcp.

Value

Aggplot2 object.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Examples

# Typical usage. demo_fit is an mcpfit object.plot_pars(demo_fit)## Not run: # More optionsplot_pars(demo_fit, regex_pars = "^cp_")  # Plot only change pointsplot_pars(demo_fit, pars = c("int_3", "time_3"))  # Plot these parametersplot_pars(demo_fit, type = c("trace", "violin"))  # Combine plots# Some plots only take pairs. hex is good to assess identifiabilityplot_pars(demo_fit, type = "hex", pars = c("cp_1", "time_2"))# Visualize the priors:plot_pars(demo_fit, prior = TRUE)# Useful for varying effects:# plot_pars(my_fit, pars = "varying", ncol = 3)  # plot all varying effects# plot_pars(my_fit, regex_pars = "my_varying", ncol = 3)  # plot all levels of a particular varying# Customize multi-column ggplots using "*" instead of "+" (patchwork)library(ggplot2)plot_pars(demo_fit, type = c("trace", "dens_overlay")) * theme_bw(10)## End(Not run)

Posterior Predictive Checks For Mcpfit Objects

Description

Plot posterior (default) or prior (prior = TRUE) predictive checks. This is convenience wrapperaround the⁠bayesplot::ppc_*()⁠ methods.

Usage

pp_check(  object,  type = "dens_overlay",  facet_by = NULL,  newdata = NULL,  prior = FALSE,  varying = TRUE,  arma = TRUE,  nsamples = 100,  ...)

Arguments

object

Anmcpfit object.

type

One ofbayesplot::available_ppc("grouped", invert = TRUE) %>% stringr::str_remove("ppc_")

facet_by

Name of a column in data modeled as varying effect(s).

newdata

Atibble or adata.frame containing predictors in the model. IfNULL (default),the original data is used.

prior

TRUE/FALSE. Plot using prior samples? Useful formcp(..., sample = "both")

varying

One of:

TRUE All varying effects (fit$pars$varying).
FALSE No varying effects (c()).
Character vector: Only include specified varying parameters - seefit$pars$varying.

arma

Whether to include autoregressive effects.

TRUE Compute autoregressive residuals. Requires the response variable innewdata.
FALSE Disregard the autoregressive effects. Forfamily = gaussian(),predict() just usesigma for residuals.

nsamples

Number of draws. Note that you may want to use all data for summary geoms.e.g.,pp_check(fit, type = "ribbon", nsamples = NULL).

...

Further arguments passed tobayesplot::ppc_type(y, yrep, ...)

Value

Aggplot2 object for single plots. Enriched bypatchwork for faceted plots.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Examples

pp_check(demo_fit)pp_check(demo_fit, type = "ecdf_overlay")#pp_check(some_varying_fit, type = "loo_intervals", facet_by = "id")

Fits and predictions from samples and newdata

Description

Fits and predictions from samples and newdata

Usage

pp_eval(  object,  newdata = NULL,  summary = TRUE,  type = "fitted",  probs = TRUE,  rate = TRUE,  prior = FALSE,  which_y = "ct",  varying = TRUE,  arma = TRUE,  nsamples = NULL,  samples_format = "tidy",  scale = "response",  ...)

Arguments

object

Anmcpfit object.

newdata

Atibble or adata.frame containing predictors in the model. IfNULL (default),the original data is used.

summary

Summarise at each x-value

type

One of:

"fitted": return fitted values. See alsofitted()
"predict": return predicted values, using random dispersion around the central tendency(e.g.,y_predict = rnorm(N, y_fitted, sigma_fitted) forfamily = gaussian()).See alsopredict().
"residuals": same as "predict" but the observed y-values are subtracted. See alsoresiduals()

probs

Vector of quantiles. Only in effect whensummary == TRUE.

rate

prior

TRUE/FALSE. Plot using prior samples? Useful formcp(..., sample = "both")

which_y

What to plot on the y-axis. One of

"ct": The central tendency which is often the mean after applying thelink function.
"sigma": The variance
"ar1","ar2", etc. depending on which order of the autoregressiveeffects you want to plot.

varying

One of:

TRUE All varying effects (fit$pars$varying).
FALSE No varying effects (c()).
Character vector: Only include specified varying parameters - seefit$pars$varying.

arma

Whether to include autoregressive effects.

TRUE Compute autoregressive residuals. Requires the response variable innewdata.
FALSE Disregard the autoregressive effects. Forfamily = gaussian(),predict() just usesigma for residuals.

nsamples

samples_format

One of "tidy" or "matrix". Controls the output format whensummary == FALSE.See more under "value"

scale

One of

"response": return on the observed scale, i.e., after applying the inverse link function.
"linear": return on the parameter scale (where the linear trends are modelled).

...

Currently unused

Value

Ifsummary = TRUE: Atibble with the posterior mean for each row innewdata,Ifnewdata isNULL, the data infit$data is used.
Ifsummary = FALSE andsamples_format = "tidy": Atidybayestibble with all the posteriorsamples (Ns) evaluated at each row innewdata (Nn), i.e., with⁠Ns x Nn⁠ rows. If there arevarying effects, the returned data is expanded with the relevant levels for each row.
The return columns are:
- Predictors fromnewdata.
- Sample descriptors: ".chain", ".iter", ".draw" (see thetidybayes package for more), and "data_row" (newdata rownumber)
- Sample values: one column for each parameter in the model.
- The estimate. Either "predict" or "fitted", i.e., the name of thetype argument.
Ifsummary = FALSE andsamples_format = "matrix": AnN_draws Xnrows(newdata) matrix with fitted/predictedvalues (depending ontype). This format is used bybrms and it's useful asyrep in⁠bayesplot::ppc_*⁠ functions.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Samples from the Posterior Predictive Distribution

Description

Samples from the Posterior Predictive Distribution

Usage

## S3 method for class 'mcpfit'predict(  object,  newdata = NULL,  summary = TRUE,  probs = TRUE,  rate = TRUE,  prior = FALSE,  which_y = "ct",  varying = TRUE,  arma = TRUE,  nsamples = NULL,  samples_format = "tidy",  ...)

Arguments

object

Anmcpfit object.

newdata

Atibble or adata.frame containing predictors in the model. IfNULL (default),the original data is used.

summary

Summarise at each x-value

probs

Vector of quantiles. Only in effect whensummary == TRUE.

rate

prior

TRUE/FALSE. Plot using prior samples? Useful formcp(..., sample = "both")

which_y

What to plot on the y-axis. One of

"ct": The central tendency which is often the mean after applying thelink function.
"sigma": The variance
"ar1","ar2", etc. depending on which order of the autoregressiveeffects you want to plot.

varying

One of:

TRUE All varying effects (fit$pars$varying).
FALSE No varying effects (c()).
Character vector: Only include specified varying parameters - seefit$pars$varying.

arma

Whether to include autoregressive effects.

TRUE Compute autoregressive residuals. Requires the response variable innewdata.
FALSE Disregard the autoregressive effects. Forfamily = gaussian(),predict() just usesigma for residuals.

nsamples

samples_format

One of "tidy" or "matrix". Controls the output format whensummary == FALSE.See more under "value"

...

Currently unused

Value

Ifsummary = TRUE: Atibble with the posterior mean for each row innewdata,Ifnewdata isNULL, the data infit$data is used.
Ifsummary = FALSE andsamples_format = "tidy": Atidybayestibble with all the posteriorsamples (Ns) evaluated at each row innewdata (Nn), i.e., with⁠Ns x Nn⁠ rows. If there arevarying effects, the returned data is expanded with the relevant levels for each row.
The return columns are:
- Predictors fromnewdata.
- Sample descriptors: ".chain", ".iter", ".draw" (see thetidybayes package for more), and "data_row" (newdata rownumber)
- Sample values: one column for each parameter in the model.
- The estimate. Either "predict" or "fitted", i.e., the name of thetype argument.
Ifsummary = FALSE andsamples_format = "matrix": AnN_draws Xnrows(newdata) matrix with fitted/predictedvalues (depending ontype). This format is used bybrms and it's useful asyrep in⁠bayesplot::ppc_*⁠ functions.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Examples

predict(demo_fit)  # Evaluate at each demo_fit$datapredict(demo_fit, probs = c(0.1, 0.5, 0.9))  # With median and 80% credible interval.predict(demo_fit, summary = FALSE)  # Samples instead of summary.predict(  demo_fit,  newdata = data.frame(time = c(-5, 20, 300)),  # Evaluate  probs = c(0.025, 0.5, 0.975))

Print mcplist

Description

Shows a list in a more condensed format usingstr(list).

Usage

## S3 method for class 'mcplist'print(x, ...)

Arguments

x

Anmcpfit object.

...

Currently ignored

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Nice printing texts

Description

Useful forprint(fit$jags_code),print(mcp_demo$call), etc.

Usage

## S3 method for class 'mcptext'print(x, ...)

Arguments

x

Character, often with newlines.

...

Currently ignored.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Examples

mytext = "line1 = 2\n line2 = 'horse'"class(mytext) = "mcptext"print(mytext)

Probit function

Description

Probit function

Usage

probit(mu)

Arguments

mu

A vector of probabilities (0.0 to 1.0)

Value

A vector with same length asmu

Recover the levels of varying effects in mcmc.list

Description

Jags uses 1, 2, 3, ..., etc. for indexing of varying effects.This function adds back the original levels, whether numeric or string

Usage

recover_levels(samples, data, mcmc_col, data_col)

Arguments

samples

An mcmc.list with varying columns starting inmcmc_col.

data

A tibble or data.frame with the cols indata_col.

mcmc_col

A vector of strings.

data_col

A vector of strings. Has to be same length asmcmc_col.'

Remove varying or population terms from a formula

Description

WARNING: removes response side from the formula

Usage

remove_terms(form, remove)

Arguments

form

A formula

remove

Either "varying" or "population". These are removed.

Value

A formula

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Compute Residuals From Mcpfit Objects

Description

Equivalent tofitted(fit, ...) - fit$data[, fit$data$yvar] (orfitted(fit, ...) - newdata[, fit$data$yvar]),but with fixed arguments forfitted:⁠rate = FALSE, which_y = 'ct', samples_format = 'tidy'⁠.

Usage

## S3 method for class 'mcpfit'residuals(  object,  newdata = NULL,  summary = TRUE,  probs = TRUE,  prior = FALSE,  varying = TRUE,  arma = TRUE,  nsamples = NULL,  ...)

Arguments

object

Anmcpfit object.

newdata

Atibble or adata.frame containing predictors in the model. IfNULL (default),the original data is used.

summary

Summarise at each x-value

probs

Vector of quantiles. Only in effect whensummary == TRUE.

prior

TRUE/FALSE. Plot using prior samples? Useful formcp(..., sample = "both")

varying

One of:

TRUE All varying effects (fit$pars$varying).
FALSE No varying effects (c()).
Character vector: Only include specified varying parameters - seefit$pars$varying.

arma

Whether to include autoregressive effects.

TRUE Compute autoregressive residuals. Requires the response variable innewdata.
FALSE Disregard the autoregressive effects. Forfamily = gaussian(),predict() just usesigma for residuals.

nsamples

...

Currently unused

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Examples

residuals(demo_fit)residuals(demo_fit, probs = c(0.1, 0.5, 0.9))  # With median and 80% credible interval.residuals(demo_fit, summary = FALSE)  # Samples instead of summary.

Run parallel MCMC sampling using JAGS.

Description

Run parallel MCMC sampling using JAGS.

Usage

run_jags(  data,  jags_code,  pars,  ST,  cores,  sample,  n.chains,  n.iter,  n.adapt,  inits)

Arguments

data

Data.frame or tibble in long format.

jags_code

A string. JAGS model, usually returned bymake_jagscode().

pars

Character vector of parameters to save/monitor.

ST

A segment table (tibble), returned byget_segment_table.Only really used when the model contains varying effects.

cores

Positive integer or "all". Number of cores.

1: serial sampling.options(mc.cores = 3) will dominatecores = 1but not larger values ofcores.
⁠>1⁠: parallel sampling on this number of cores. Ideally setchainsto the same value. Note:cores > 1 takes a few extra seconds the firsttime it's called but subsequent calls will start sampling immediately.
"all": use all cores but one and setschains to the same value. This isa convenient way to maximally use your computer's power.

sample

One of

"post": Sample the posterior.
"prior": Sample only the prior. Plots, summaries, etc. willuse the prior. This is useful for prior predictive checks.
"both": Sample both prior and posterior. Plots, summaries, etc.will default to using the posterior. The prior only has effect when doingSavage-Dickey density ratios inhypothesis.
"none" orFALSE: Do not sample. Returns an mcpfitobject without sample. This is useful if you only want to checkprior strings (fit$prior), the JAGS model (fit$jags_code), etc.

n.chains

the number of parallel chains for the model

n.iter

number of iterations to monitor

n.adapt

the number of iterations for adaptation. Seeadapt for details. Ifn.adapt = 0 then no adaptation takes place.

inits

A list if initial values for the parameters. This can be usefulif a model fails to converge. Read more injags.model.Defaults toNULL, i.e., no inits.

Value

'mcmc.list“

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Transform a prior from SD to precision.

Description

JAGS uses precision rather than SD. This function convertsdnorm(4.2, 1.3) intodnorm(4.2, 1/1.3^2). It allows users to specifypriors using SD and then it's transformed for the JAGS code. It works for thefollowing distributions: dnorm|dt|dcauchy|ddexp|dlogis|dlnorm. In all ofthese,tau/sd is the second parameter.

Usage

sd_to_prec(prior_str)

Arguments

prior_str

String. A JAGS prior. Can be truncated, e.g.⁠dt(3, 2, 1) T(my_var, )⁠.

Value

A string

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Summarise mcpfit objects

Description

Summarise parameter estimates and model diagnostics.

Usage

## S3 method for class 'mcpfit'summary(object, width = 0.95, digits = 2, prior = FALSE, ...)fixef(object, width = 0.95, prior = FALSE, ...)ranef(object, width = 0.95, prior = FALSE, ...)## S3 method for class 'mcpfit'print(x, ...)

Arguments

object

Anmcpfit object.

width

Float. The width of the highest posterior density interval(between 0 and 1).

digits

prior

TRUE/FALSE. Summarise prior instead of posterior?

...

Currently ignored

x

Anmcpfit object.

Value

A data frame with parameter estimates and MCMC diagnostics.OBS: The change point distributions are often not unimodal and symmetric sothe intervals can be deceiving Plot them usingplot_pars(fit).

mean is the posterior mean
lower is the lower quantile of the highest-density interval (HDI) given inwidth.
upper is the upper quantile.
Rhat is the Gelman-Rubin convergence diagnostic which is often taken tobe acceptable if < 1.1. It is computed usinggelman.diag.
n.eff is the effective sample size computed usingeffectiveSize.Low effective sample sizes are also obvious as poor mixing in trace plots(seeplot_pars(fit)). Read how to deal with such problemshere
ts_err is the time-series error, taking autoregressive correlationinto account. It is computed usingspectrum0.ar.

For simulated data, the summary contains two additional columns so that itis easy to inspect whether the model can recover the parameters. Runsimulation and summary multiple times to get a sense of the robustness.

sim is the value used to generate the data.
match is"OK" ifsim is contained in the HDI interval (lower toupper).

Functions

fixef(): Get population-level ("fixed") effects of anmcpfit object.
ranef(): Get varying ("random") effects of anmcpfit object.
print(mcpfit): Print the posterior summary of anmcpfit object.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Examples

# Typical usagesummary(demo_fit)summary(demo_fit, width = 0.8, digits = 4)  # Set HDI width# Get the results as a data frameresults = summary(demo_fit)# Varying (random) effects# ranef(my_fit)# Summarise priorsummary(demo_fit, prior = TRUE)

Get tidy samples with or without varying effects

Description

Returns in a format useful forfit$simulate() with population parameters in wide formatand varying effects in long format (the number of rows will bensamples * n_levels_in_varying).

Usage

tidy_samples(  fit,  population = TRUE,  varying = TRUE,  absolute = FALSE,  prior = FALSE,  nsamples = NULL)

Arguments

fit

Anmcpfit object

population

TRUE All population effects. Same asfit$pars$population.
- FALSE No population effects. Same asc().
- Character vector: Only include specified population parameters - seefit$pars$population.

varying

One of:

TRUE All varying effects (fit$pars$varying).
FALSE No varying effects (c()).
Character vector: Only include specified varying parameters - seefit$pars$varying.

absolute

TRUE Returns the absolute location of all varying change points.
- FALSE Just returns the varying effects.
- Character vector: Only do absolute transform for these varying parameters - seefit$pars$varying.
OBS: This currently only applies to varying change points. It is not implemented forrel() regressors yet.

prior

TRUE/FALSE. Summarise prior instead of posterior?

nsamples

Value

tibble of posterior draws intidybayes format.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Convert from tidy to matrix

Description

Converts from the output oftidy_samples() orpp_eval(fit, samples_format = "tidy")to anN_draws Xnrows(newdata) matrix with fitted/predicted values. This format isused ybrms and it's useful asyrep in⁠bayesplot::ppc_*⁠ functions.

Usage

tidy_to_matrix(samples, returnvar)

Arguments

samples

Samples in tidy format

returnvar

Anrlang::sym() object.

Value

AnN_draws Xnrows(newdata) matrix.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Takes any formula-like input and returns a formula

Description

Takes any formula-like input and returns a formula

Usage

to_formula(form)

Arguments

form

Formula or character (with or without initial tilde/"~")

Value

A formula

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Unpack arma order and formula

Description

Unpack arma order and formula

Usage

unpack_arma(form_str_in)

Arguments

form_str_in

A character. These are allowed: "ar(number)" or "ar(number, formula)"

Value

A list with $order and $form_str (e.g., "ar(formula)"). The formula is ar(1) or ma(1) if no formula is given

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Takes a cp formula (as a string) and returns its properties

Description

Takes a cp formula (as a string) and returns its properties

Usage

unpack_cp(form_cp, i)

Arguments

form_cp

Segment formula as string.

i

The segment number (integer)

Value

A one-row tibble with columns:

cp_int: bool. Whether there is an intercept change in the change point.
cp_in_rel: bool. Is this intercept change relative?
cp_ran_int: bool or NA. Is there a random intercept on the change point?
cp_group_col: char or NA. Which column in data define the random intercept?

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Get the intercept of a formula

Description

Get the intercept of a formula

Usage

unpack_int(form, i, ttype)

Arguments

form

A formula

i

Segment number (integer)

ttype

The term type. One of "ct" (central tendency), "sigma" (variance),or "ar" (autoregressive)

Value

A one-row tibble describing the intercept.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Unpack right-hand side

Description

This is a pretty big function. It includes unpacking sigma(), ar(), etc.

Usage

unpack_rhs(form_rhs, i, family, data, last_segment)

Arguments

form_rhs

A character representation of a formula

i

The segment number (integer)

family

An mcpfamily object returned bymcpfamily().

data

A data.frame or tibble

last_segment

The last row in the segment table, made inget_segment_table()

Value

A one-row tibble with three columns for each ofct.sigma,ar, andma:

⁠_int⁠: NA or a one-row tibble describing the intercept.
⁠_slope⁠: NA or a tibble with a row for each slope term.
⁠_code⁠: NA or a char with the JAGS/R code to implement the slope.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Unpack the slope of a formula

Description

Makes A list of terms and applies unpack_slope_term() to each of them. Then builds the code for this segment's slope

Usage

unpack_slope(form, i, ttype, last_slope)

Arguments

form

A formula

i

Segment number (integer)

ttype

The term type. One of "ct" (central tendency), "sigma" (variance),or "ar" (autoregressive)

last_slope

The element in the slope column for this ttype in the previoussegment. I.e., probably what this function returned last time it was called!

Value

A "one-row" list with code (char) and a tibble of slopes.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Unpacks a single term

Description

Returns a row forunpack_slope().

Usage

unpack_slope_term(term, i, last_slope, ttype = "")

Arguments

term

A character, e.g., "x", "I(x^2)", or "log(x)".

i

Segment number (integer)

last_slope

The element in the slope column for this ttype in the previoussegment. I.e., probably what this function returned last time it was called!

ttype

The term type. One of "ct" (central tendency), "sigma" (variance),or "ar" (autoregressive)

Value

A "one-row" list describing a slope term.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Takes a formula and returns a string representation of y, cp, and rhs

Description

Takes a formula and returns a string representation of y, cp, and rhs

Usage

unpack_tildes(segment, i)

Arguments

segment

A formula

i

The segment number (integer)

Value

A one-row tibble with columns:

form: String. The full formula for this segment.
form_y: String. The expression for y (without tilde)
form_cp: String. The formula for the change point
form_rhs: String. The formula for the right-hand side

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Get relevant info about varying parameters

Description

Returns parameters, data columns, and implicated segments given parameter name(s) or column(s).

Usage

unpack_varying(fit, pars = NULL, cols = NULL)

Arguments

pars

NULL/FALSE for nothing.TRUE for all. A vector of varying parameter names for specifics.

cols

NULL/FALSE for nothing.TRUE for all. A vector of varying column names for specifics. Usually provided via "facet_by" argument in other functions.

Details

Returns a list with

Value

A list. See details.

Slots

pars: Character vector of parameter names.NULL if empty.
cols: Character vector of data column names.NULL if empty.
indices: Logical vector of segments in the segment table that contains the varying effect

Unpack varying effects

Description

Unpack varying effects

Usage

unpack_varying_term(term, i)

Arguments

term

A character, e.g., "x", "I(x^2)", or "log(x)".

i

Segment number (integer)

Value

A "one-row" list describing a varying intercept.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Unpacks y variable name

Description

Unpacks y variable name

Usage

unpack_y(form_y, i, family)

Arguments

form_y

Character representation of formula

i

The segment number (integer)

family

An mcpfamily object returned bymcpfamily().

Value

A one-row tibble with the columns

y: string. The y variable name.
trials: string. The trials variable name.
weights: string. The weights variable name.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Add loo if not already present

Description

Add loo if not already present

Usage

with_loo(fit, save_psis = FALSE, info = NULL)

Arguments

fit

An mcpfit object

save_psis

Logical. See documentation of loo::loo

info

Optional message if adding loo

Value

An mcpfit object with loo.

Author(s)

Jonas Kristoffer Lindeløvjonas@lindeloev.dk

Movatterモバイル変換

mcp: Regression with Multiple Change Points

Description

Author(s)

See Also

Bernoulli family for mcp

Description

Usage

Arguments

Checks if all terms are in the data

Description

Usage

Arguments

Author(s)

Compute information criteria for model comparison

Description

Usage

Arguments

Value

Functions

Author(s)

See Also

Examples

Cumulative pasting of character columns

Description

Usage

Arguments

Value

Author(s)

Examplemcpfit for examples

Description

Usage

Format

Exponential family for mcp

Description

Usage

Arguments

Expected Values from the Posterior Predictive Distribution

Description

Usage

Arguments

Value

Author(s)

See Also

Examples

Format code with one or multiple terms

Description

Usage

Arguments

Value

Author(s)

Density geom forplot.mcpfit()

Description

Usage

Arguments

Value

Return a geom_line representing the quantiles

Description

Usage

Arguments

Value

Author(s)

Callget_formula_str for multiple ytypes and paste strings

Description

Usage

Arguments

Value

Author(s)

Gets code for ARMA terms, resulting in a "resid_"

Description

Usage

Arguments

Value

Author(s)

Extracts the order from ARMA parameter name(s)

Description

Usage

Arguments

Value

Compute the density at a specific point.

Example`mcpfit` for examples

Density geom for`plot.mcpfit()`

Call`get_formula_str` for multiple ytypes and paste strings

Adds helper variables for use in`run_jags`