Package {rmoriebricklayer}

Title:	Reproducible Data Capsules with Provenance and Fallback
Version:	1.0.0
Description:	For analyzing published and unpublished tables, from open-data portals and statistical agencies to a researcher's own data, so that a result can be traced to the exact data behind it and checked again afterwards. A 'capsule' records where each source came from, its Secure Hash Algorithm 256 ('SHA-256') hash, the schema it should have and an archived copy in the Internet Archive 'Wayback Machine' ('web.archive.org'); sources are found through the Comprehensive Knowledge Archive Network ('CKAN', 'ckan.org'). Run again, a capsule says whether each source still matches and, if it does not, how its distribution moved (Kolmogorov-Smirnov, chi-square, population stability index and 'Benford' first-digit screens), and when a source cannot be reached it runs on schema-driven synthetic data that is marked as such. Capsules can be signed with keyed hashes ('HMAC-SHA-256', RFC 2104) or with the post-quantum signature standards 'ML-DSA' (FIPS 204), 'SLH-DSA' (FIPS 205) and 'XMSS' (RFC 8391), implemented from the standards with no system cryptography library and checked against 'OpenSSL' 3.5 and the known-answer vectors of 'NIST' (the US National Institute of Standards and Technology). For the administrative tables these capsules usually hold, it computes period-over-period change with exact intervals, trends (Mann-Kendall with 'Theil-Sen' slopes), rates, and the bounds that rounding and suppression in a release imply; small-area rates with the empirical Bayes shrinkage of Clayton and 'Kaldor' (1987) DOI 10.2307/2532003; self-exciting event series as 'Hawkes' processes; custodial populations by stock-and-flow measures; and the burden of air pollution from dispersion to attributable deaths.
License:	AGPL-3
Copyright:	The bundled excerpts of published statistics in inst/extdata belong to their publishers; see inst/COPYRIGHTS.
Encoding:	UTF-8
Language:	en-US
Depends:	R (≥ 4.1.0)
Imports:	graphics, grDevices, methods, parallel, stats, tools, utils
Suggests:	withr, bslib (≥ 0.9.0), digest, jsonlite, knitr, later, MASS, openssl, pkgdown, rmarkdown, sf, promises, shiny, stringi, survival, survey, testthat (≥ 3.0.0)
VignetteBuilder:	knitr
Config/testthat/edition:	3
SystemRequirements:	C++17, libcurl (deb: libcurl4-openssl-dev, rpm: libcurl-devel, brew: curl), zlib (deb: zlib1g-dev, rpm: zlib-devel).
Config/roxygen2/version:	8.0.0
NeedsCompilation:	yes
Packaged:	2026-10-21 19:52:40 UTC; rootcoder
Author:	Vansh Singh Ruhela [aut, cre]
Maintainer:	Vansh Singh Ruhela <vsruhela@proton.me>
Repository:	CRAN
Date/Publication:	2026-10-21 22:00:02 UTC

Admissions implied by a population and a length of stay

Description

Works out how many admissions a given average daily population and average length of stay imply over a period. Published statistics often give the population held and the average stay but not the number admitted (or the other way round); this recovers the missing quantity from the identity that links all three.

Usage

admits(adp, alos, t = 365)

Arguments

Details

Rearranging the identity \bar{X}_{hc} = N_a \bar{X}_t / t (Lakner 1976, p.18-20) for N_a. Useful when two of the three quantities are published and the third is not. The identity holds exactly only in a steady state (admissions and stays roughly constant across the period); with a population that is growing or shrinking quickly the result is an approximation.

Errors when any argument is not numeric, is NA or Inf, is not positive, or has more than one value.

Value

A numeric scalar: the implied number of admissions over the period, \mathrm{adp} \times t / \mathrm{alos}. It is not rounded.

See Also

adp(), alos(), stkflow()

Examples

# Lakner (p.20) runs this the other way: t = 365, an average daily
# population of 25 and 1,750 admissions imply a stay of 5.2 days.
alosimplied <- 25 * 365 / 1750
round(alosimplied, 1)

# and back again
admits(25, alosimplied)

# a 90-day quarter with 40 held on average and a 12-day average stay
admits(40, 12, t = 90)

# an average stay of zero is refused
try(admits(25, 0))

Average daily population

Description

The mean number of person-days served per day over a period: a STOCK, answering how many people are held at one time. Use it to report how much confinement (or bed occupancy, shelter use, open caseload) a system used over a year or any other window, from the days each person spent in it. It is the measure capacity planning needs, because it counts beds in use rather than people who passed through.

Usage

adp(days, t = 365)

Arguments

Details

This is Lakner's \bar{X}_{hc} = \sum X_i / t (1976, p.15), where X_i is the days served by person i inside the period. The denominator is TIME, which is what makes it a stock. Compare alos(), whose denominator is people: the two are linked by \mathrm{ADP} = N \times \mathrm{ALOS} / t (Little's law), so admits() and stkflow() can move between them.

Only days served inside the period should be counted: a stay that began before the period or ends after it contributes its days within the window.

Errors when days is not numeric (text and factors are refused), is empty, has NA, Inf or a negative value, or when t is not a single positive number.

Value

A numeric scalar: person-days per day, that is the average number of people held on any one day of the period.

References

Lakner, E. (1976) A Manual of Statistical Sampling Methods for Corrections Planners. University of Illinois at Urbana-Champaign.

See Also

alos(), stkflow(), adpcnts()

Examples

# Lakner's own worked example (p.15): 3,000 inmates served 13,500
# detention days in a year.
adp(13500)

# per-person days give the same total
adp(c(10, 20, 30), t = 30)

# a leap year, with the period length taken from the dates
adp(13500, t = perdays("2024-01-01", "2024-12-31"))

# negative days are refused
try(adp(c(10, -2)))

Person-days estimated from periodic headcounts

Description

When only a count of people present on certain days is available – not a record per person – the person-days over the period are the mean of those counts scaled to the period's length. Use it with adp() or alos() when the source is a daily or weekly count sheet rather than admission and release dates.

Usage

adpcnts(cnts, t = 365)

Arguments

Details

Lakner's X_t = (\frac{1}{C}\sum N_i) t (1976, p.21). The counts need not cover every day: counting on weekdays only is the usual case, and the mean of the days counted stands in for the days not counted. That substitution is an assumption, and it fails if the days counted differ systematically from the days missed – weekday-only counting where weekend admissions are released before Monday, for instance.

Value

A numeric scalar: the estimated person-days over the period, mean(counts) * t. Despite the function's name this is a total of person-days, not an average daily population; divide by t (or pass it to adp()) for the population held on an average day.

See Also

adp()

Examples

# Lakner (p.21): counts on 255 days summing to 34,935 imply just over
# fifty thousand detention days across the year.
adpcnts(rep(34935 / 255, 255))

# weekly counts over a 31-day month, then the average daily population
pd <- adpcnts(c(120, 130, 125, 128), t = 31)
pd
adp(pd, t = 31)

# a missing count must be dropped or filled first
try(adpcnts(c(120, NA, 125)))

Agent-assisted reproducibility-bundle help

Description

Asks a language model to plan a brick-proof, reproducible data bundle, through llmask: your own key, your own endpoint or a local Ollama server, with the packages' help at hand (the model reads the help pages of the functions it proposes before it answers). Use it to get a first draft of the steps and R code for packaging an analysis: finding the portal's download URL, fetching with a Wayback fallback, recording hashes and a manifest, writing a synthetic stand-in, then packaging, signing and verifying the bundle.

Usage

agentbndl(req, model = NULL, route = NULL)

Arguments

Details

The request is sent with a fixed system prompt that describes that task and asks for a short plan and R code. The route is chosen as in llmask(); when no route is set up, the function returns a message saying what to set up (a provider key, a local Ollama or your own endpoint) instead of an error, and makes no network request. When a route is set up, the request and the prompt are sent to that provider or server. The answer is model output: check the code it proposes against the help pages before running it.

Errors: req not a single non-empty string; route not one of the names listed; route = "hosted" (the hosted tier is retired); and errors raised by llmask().

Value

Character scalar: the answer, or a message saying what to set up when no route answers.

Examples

# Needs a language model: a provider key (e.g. ANTHROPIC_API_KEY), your own
# endpoint or a local Ollama.
agentbndl("scaffold a bundle for analysis.R from the Toronto CKAN data")
agentbndl("sign and verify my capsule", route = "ollama")



# With every route switched off (no key in the environment and no saved
# settings: the settings folder is pointed at an empty one) the call returns
# a setup hint, not an error, and reaches no network -- safe to run anywhere:
keys <- c("OLLAMA_HOST", "MORIE_LLM_BASE_URL", "XDG_CONFIG_HOME",
          llmprvd()$key_env,
          llmprvd()$alt_env[nzchar(llmprvd()$alt_env)])
old <- Sys.getenv(keys, unset = NA)
Sys.setenv(OLLAMA_HOST = "off", MORIE_LLM_BASE_URL = "off", XDG_CONFIG_HOME = tempfile())
Sys.unsetenv(setdiff(keys, c("OLLAMA_HOST", "MORIE_LLM_BASE_URL", "XDG_CONFIG_HOME")))
cat(agentbndl("hello"))
for (v in keys) {
  if (is.na(old[[v]])) Sys.unsetenv(v) else do.call(Sys.setenv, as.list(old[v]))
}

# argument errors come first
try(agentbndl(""))
try(agentbndl("plan a bundle", route = "hosted"))

Directly age-standardised rates, with Fay-Feuer gamma intervals

Description

The rate a population would have with the age structure of a standard population, for one or more groups, with Fay and Feuer's gamma confidence interval.

Usage

agestd(events, pop, std, grp = NULL, per = 1e+05, conflvl = 0.95)

Arguments

Details

With d_i events and n_i people in age group i and standard weights w_i (the standard population's shares), the rate is y = \sum w_i d_i / n_i, its variance v = \sum (w_i / n_i)^2 d_i, and with w_M = \max_i w_i / n_i the limits are \frac{v}{2y} \chi^2_{2y^2/v}(\alpha/2) and \frac{v + w_M^2}{2(y + w_M)} \chi^2_{2(y + w_M)^2/(v + w_M^2)}(1 - \alpha/2) (Fay and Feuer 1997). The interval is conservative (coverage at least nominal in Fay and Feuer's simulations) and reduces to the exact Garwood interval for a single age group: suited to the few events of small communities. With no events the lower limit is 0. Rates are per per people. Age groups are matched to std by position: order the rows of each group by age, in the order of std.

The choice of standard matters: Statistics Canada's comparisons of First Nations, Metis and Inuit mortality standardise to the Indigenous population's own age structure (Kumar and Tjepkema 2019), which is younger than Canada's.

Value

A data frame with group, events, rate, lower, upper, and the attributes var (the variances, on the per-person scale) for rategap.

References

Fay MP, Feuer EJ (1997). Confidence intervals for directly standardized rates: a method based on the gamma distribution. Statistics in Medicine 16(7):791-801.
DOI 10.1002/(sici)1097-0258(19970415)16:7<791::aid-sim500>3.0.co;2-#

Kumar MB, Tjepkema M (2019). Suicide among First Nations people, Metis and Inuit (2011-2016): Findings from the 2011 Canadian Census Health and Environment Cohort (CanCHEC). Statistics Canada 99-011-X2019001.

See Also

rategap.

Examples

# made-up counts for one small community, by five age groups, to show the computation
std <- c(30000, 25000, 20000, 15000, 10000)   # a younger standard population
agestd(events = c(4, 9, 11, 6, 3), pop = c(9000, 8000, 6000, 4000,
                                                           2000),
                std = std)

Cross-fitted augmented inverse probability weighting for a binary treatment

Description

The average effect of a binary treatment (on everyone, or on the treated) from one row per person, with the outcome and the propensity fitted by learners on other folds, combined in the doubly robust (Neyman-orthogonal) score of Robins, Rotnitzky and Zhao; with survey or sampling weights if given.

Usage

aipwcf(
  data,
  y,
  trt,
  covs,
  estd = c("ATE", "ATT"),
  lrn = "lasso",
  basis = c("linear", "poly2"),
  w = NULL,
  des = NULL,
  folds = 5L,
  trim = 0.01,
  seed = 1L,
  conflvl = 0.95
)

Arguments

Details

With \hat\mu_d(X) the outcome fitted among those with treatment d and \hat e(X) the propensity, both out of fold, the score for the average effect is

\psi = \hat\mu_1 - \hat\mu_0 + \frac{D (Y - \hat\mu_1)}{\hat e} - \frac{(1 - D)(Y - \hat\mu_0)}{1 - \hat e} - \theta,

and for the effect on the treated

\psi = \frac{D (Y - \hat\mu_0) - \hat e (1 - D)(Y - \hat\mu_0)/(1 - \hat e)}{P(D = 1)} - \frac{D}{P(D = 1)}\theta.

\hat\theta solves the weighted mean of the score at zero; its variance is the weighted mean of \psi^2 over n, the influence-function variance of Chernozhukov et al. (2018). Propensities are clipped to [trim, 1 - trim]. The estimator is consistent when either the outcome or the propensity model is right, and with cross-fitting the learners can be flexible (lasso or ridge on a second-order basis, or any learner function).

The average effect equals DoubleML's interactive regression model given the same folds and learners; for the effect on the treated, DoubleML estimates P(D = 1) within each fold and this function in the whole sample, so the two differ slightly.

The propensity learner is logistic regression on the covariates (with basis = "poly2", on their squares and products too). With w the point estimate is design-weighted, but the standard error treats the weights as fixed and the sample as independent.

For a complex survey, pass the design from svydes instead: its weight gives the point estimate and its replicate weights (or strata and primary sampling units) the standard error, with the variance formula of the replicate type, the confidence interval on the design's degrees of freedom. The learners are fitted once, without the weights, which enter only the average of the scores, so each replicate estimate is the replicate-weighted average of the same cross-fitted scores: the same numbers as re-running the function on every replicate weight, from one fit instead of one per replicate. Rows with a missing value are left out as a domain (their scores count as zero in every replicate), not dropped from the design.

Value

A data frame of class rmbl_aipw: estimand, estimate, se, lower, upper, p_value; attributes n, treated, psrange, trimmed (how many propensities were clipped), variance ("influence", or the design's type) and propensity (a data frame of each row's out-of-fold propensity before clipping, its treatment and its weight, which plot() draws).

References

Robins JM, Rotnitzky A, Zhao LP (1994). Estimation of regression coefficients when some regressors are not always observed. Journal of the American Statistical Association 89(427):846-866. DOI 10.1080/01621459.1994.10476818

Chernozhukov V et al. (2018), see plndml.

See Also

plrcf, plndml, baltab.

Examples

set.seed(2)
n <- 800
d <- data.frame(age = runif(n, 18, 70), income = rlnorm(n, 10, 0.5))
d$belong <- rbinom(n, 1, plogis(-2 + 0.03 * d$age + 0.2 * log(d$income)))
d$wellbeing <- 50 + 5 * d$belong + 0.2 * d$age + 2 * log(d$income) +
  rnorm(n, 0, 8)
aipwcf(d, "wellbeing", "belong", c("age", "income"), lrn = "ols")
# with 50 bootstrap replicate weights: one fit, the replicate standard error
d$w <- runif(n, 0.5, 2)
reps <- sapply(1:50, function(b) d$w * rpois(n, 1))
colnames(reps) <- paste0("bsw", 1:50)
d <- cbind(d, reps)
des <- svydes(d, w = "w", reps = colnames(reps))
aipwcf(d, "wellbeing", "belong", c("age", "income"), lrn = "ols",
              des = des)

Average length of stay

Description

The mean number of person-days served per person: the FLOW side of the same person-days. Use it to say how long a typical admission lasts, and together with adp() to tell whether a change in the population held came from more people or from longer stays. For the full distribution of stays (median, spread, interval) use staysmry().

Usage

alos(days, n)

Arguments

Details

Lakner's \bar{X}_t = \sum X_i / N' (1976, p.16), with a caveat worth repeating (p.16-17): the period must be longer than the longest stay people actually serve, or the average is biased DOWNWARD, since the longest stays are the ones that fail to finish inside the window. For short-stay facilities a year is comfortable; for long sentences it is not, and the period has to be set from the records.

Errors when days is not numeric, is empty, or has NA, Inf or negative values, and when n is not a single positive number.

Value

A numeric scalar: days per person, \sum X_i / n.

References

Lakner, E. (1976) A Manual of Statistical Sampling Methods for Corrections Planners. University of Illinois at Urbana-Champaign.

See Also

adp(), stkflow()

Examples

# Lakner (p.17): 2,700 inmates admitted and released served 12,150
# detention days between them.
alos(12150, 2700)

# from one value per released person
stays <- c(2, 5, 1, 30, 4, 7)
alos(stays, length(stays))

# zero people cannot have an average stay
try(alos(12150, 0))

Analyse a published administrative table in one call

Description

The questions asked of a published table are always the same: what is in it, what changed between periods and how sure can one be, what the rates are if there is an exposure, whether a short series trends, and whether this release differs from the one held before. Each has a function in this package; this runs them together, with the pieces a reviewer needs and that are easy to forget on a deadline: exact intervals on every change, multiple-comparison adjustment over the groups scanned, the envelope that rounding and suppression in the release imply, and a drift screen against the prior capsule.

Usage

anltab(
  data,
  value,
  prd,
  by = NULL,
  pop = NULL,
  prior = NULL,
  units = c("count", "continuous", "percent"),
  rnd = NULL,
  rndkind = c("nearest", "random"),
  suplim = NULL,
  dir = c("neutral", "higher", "lower"),
  per = 1000,
  conflvl = 0.95,
  alpha = 0.05,
  adjust = "BH"
)

Arguments

Details

The steps, in order:

1. profile: profcols() and misssmry() on data;

2. change: yoy() between consecutive periods within each by group. For units = "count" the exact p-values are then adjusted across all comparisons with scanadj(), and, when rnd or suplim is given, yoybnd() adds the publication envelope (see pubbnd()). For "continuous" and "percent" units neither step runs: rnd and suplim are then ignored, with a warning;

3. rates (only with pop): rtcalc() per group and period, and ratechg() between periods, adjusted with scanadj();

4. trend: per group, the values summed within each period, then trndtst() (Mann-Kendall tau, p-value and Theil-Sen slope) and, for counts, cnttrnd() (Poisson rate ratio per period). A group with fewer than three periods gets NA; a test that fails gets NA rather than an error;

5. drift (only with prior): cpsdrft() of data against prior.

Nothing is written and nothing is fetched; a progress bar for the trend step is drawn on stderr unless options(morie.quiet = TRUE). Errors: data or prior not a data frame; a named column not in data; an invalid units, dir or rndkind; a conflvl outside (0, 1); and any error of the functions called (for example negative counts in yoy()).

Value

A list of class bricklayer_analysis (with a print method that summarises it in a few lines):

profile

profcols() of data.

missingness

named numeric from misssmry().

change

the yoy() table; for counts with the added columns p_value, padj and significant, and with env_* and compctlo/compctup when publication bounds were requested.

rates, ratechg

the rtcalc() and ratechg() tables when population is given, otherwise NULL.

trend

data frame with one row per group: the by columns, periods, first, last, tau, trndp, slope, rateratio, rateratiolo, rateratioup.

drift

cpsdrft() result against prior, or NULL.

meta

list of the arguments used, nrows, the sorted periods, prior (logical) and when (a time stamp).

See Also

rptanl() writes the result as a Markdown or HTML report.

Examples

path <- system.file("extdata", "otis_a01_individuals.csv",
                    package = "rmoriebricklayer")
otis <- read.csv(path)
a <- anltab(otis, value = "individuals", prd = "year",
                   by = c("table", "group"))
a
a$change[, c("table", "group", "year", "pctchg", "padj")]
a$trend[, c("group", "tau", "trndp", "rateratio")]

# counts published rounded to base 5: the intervals widen
b <- anltab(otis, value = "individuals", prd = "year",
                   by = c("table", "group"), rnd = 5)
b$change[1:3, c("group", "year", "pctlo", "compctlo")]

# rates with an exposure column
d <- data.frame(year = rep(2019:2022, 2), site = rep(c("A", "B"), each = 4),
                n = c(10, 12, 15, 9, 30, 28, 35, 40),
                pop = rep(c(1000, 2000), each = 4))
r <- anltab(d, "n", "year", by = "site", pop = "pop", per = 1000)
r$ratechg[, c("site", "year", "rate", "rateratio", "padj")]

# a column that is not in the data is an error
try(anltab(otis, value = "count", prd = "year"))

Capture files the Internet Archive cannot reach, and publish them as an archive.org item

Description

Downloads each URL from here, records its SHA-256, size and time in a manifest, timestamps the manifest with an RFC 3161 time-stamping authority, and, with archive.org keys, uploads the files, the manifest and its timestamp as one public archive.org item. The item is the mirror that fetchvrf falls back to, and the timestamp shows that these were the bytes on that date.

Usage

arccap(
  urls,
  dir,
  item = NULL,
  meta = list(),
  tsa = "https://freetsa.org/tsr",
  ak = NULL,
  sk = NULL,
  tmout = 600L
)

Arguments

Details

Some servers refuse the Internet Archive's crawler while serving everyone else: Save Page Now then reports that the target “cannot find” the page, or “blocks access” to it, although the page is there. Those URLs cannot be put in the Wayback Machine, but they can be captured by anyone who can reach them, and the capture can be published. Pass a wbsave result as urls to capture just the URLs it failed on with such a refusal.

The manifest (manifest.json) lists, for each URL, the file, its SHA-256, its size and the time it was fetched. The timestamp (manifest.tsr) is the authority's signed statement that the manifest's SHA-256 existed at its time; tsvrf checks it, as does openssl ts -verify. Each file's digest in the manifest ties the file to that time.

The upload uses the Internet Archive's S3-like API: each file is sent with a PUT, streamed from disk, to
https://s3.us.archive.org/<item>/<file>.
The first request creates the item with the metadata given. The keys are the same as for wbsave, so they are read from IA_S3_ACCESS_KEY and IA_S3_SECRET_KEY, or their _FILE twins. They travel only in the authorization header and are never written or printed. Publish only what the source's licence lets you redistribute.

Value

A data frame of class rmbl_archive_capture, one row per URL: url, file, sha256, bytes, fetched_utc, status ("captured", "uploaded" or "failed: ...") and mirror (the archive.org download address, for fetchvrf); attributes manifest, timestamp (a tsvrf result, or NULL) and item.

References

Adams C, Cain P, Pinkas D, Zuccherato R (2001). Internet X.509 Public Key Infrastructure Time-Stamp Protocol (TSP). RFC 3161. DOI 10.17487/RFC3161

Internet Archive. ias3: the Internet Archive's S3-like API. archive.org/developers/ias3.html

See Also

fetchvrf, wbsave, tsvrf.

Examples

w <- wbsave(urls)
lic <- "https://www.statcan.gc.ca/en/terms-conditions/open-licence"
meta <- list(title = "Statistics Canada tables (captures)",
             creator = "Statistics Canada", licenseurl = lic)
cap <- arccap(w, file.path(tempdir(), "capture"),
                       item = "my-statcan-tables-2026", meta = meta)
fetchvrf(cap$url[1], tempfile(), sha256 = cap$sha256[1],
               mirrors = cap$mirror[1])

Use Text As-Is, Falling Back to ASCII When It Cannot Be Represented

Description

Returns x unchanged when it is valid, well-formed text (so legitimate UTF-8 such as an accented name is preserved), and only transliterates to plain ASCII via toascii() when the text is not valid UTF-8 (an encoding error) or when force = TRUE (for ASCII-only destinations such as a package DESCRIPTION) . This lets author and supervisor names keep their accents wherever UTF-8 is supported while degrading gracefully instead of erroring where it is not.

Usage

asciifb(x, force = FALSE)

Arguments

Details

Each element is tested on its bytes with validUTF8(), so the outcome does not depend on the session locale. Valid elements (including plain ASCII and properly encoded accented text) are returned unchanged; invalid ones go through toascii(), which drops the bytes that are not UTF-8 and transliterates the rest. With force = TRUE every element goes through toascii(). A missing value is kept as NA unless force = TRUE, when it becomes "".

Value

A character vector the length of x: each element unchanged when it is valid UTF-8 (and force is not TRUE), otherwise its ASCII transliteration.

See Also

toascii(), which does the transliteration.

Examples

# By default valid UTF-8 is preserved (accents kept where supported).
asciifb("\u00c1ngela")               # "\u00c1ngela"

# force = TRUE always transliterates (for ASCII-only destinations
# such as a package DESCRIPTION).
asciifb("\u00c1ngela", force = TRUE)  # "Angela"

# Plain ASCII is returned unchanged either way.
asciifb("plain name")

# Vectorised; each element handled independently.
asciifb(c("caf\u00e9", "resume"), force = TRUE)

# Bytes that are not UTF-8 (a Latin-1 file read as UTF-8) are repaired.
bad <- rawToChar(as.raw(c(0x63, 0x61, 0x66, 0xe9)))
validUTF8(bad)
asciifb(bad)

Population attributable fraction

Description

The share of all cases in a population that is due to an exposure: the fraction that would not occur if the exposure were removed and nothing else changed. It combines how strongly the exposure raises risk (the relative risk) with how common it is (the prevalence); multiplied by the number of cases it gives the attributable cases, as polbrdn does.

Usage

attrfrac(rr, expoprev, rrlo = NULL, rrup = NULL)

Arguments

Details

Levin's formula (Levin 1953; Rothman, Greenland and Lash 2008, chapter 5):

PAF = p (RR - 1) / (1 + p (RR - 1)),

with p the proportion exposed. With everyone exposed it is 1 - 1/RR. It assumes the relative risk is causal and unconfounded and that removing the exposure changes nothing else. A relative risk below 1 gives a negative fraction (the exposure protects).

When both confidence limits of the RR are given, the fraction's limits are the formula applied to them: the formula increases with the RR, so the interval of the fraction is the image of the interval of the risk. Refusals: an rr or limit that is not a single finite positive number (logical values included), a prevalence outside [0, 1], only one of the two limits, or limits that do not enclose rr.

Value

A single number, the attributable fraction: in [0, 1) for rr >= 1, negative for a protective exposure. When the limits were given it carries an attribute ci, a named numeric vector (lower, upper).

References

Rothman, K. J., Greenland, S. and Lash, T. L. (2008). Modern Epidemiology, 3rd ed., chapter 5. Levin, M. L. (1953). The occurrence of lung cancer in man. Acta Unio Internationalis Contra Cancrum 9, 531-541.

See Also

crfpm25 for a relative risk from a concentration, polbrdn for the attributable cases.

Examples

attrfrac(1.5, 0.4)
attrfrac(2, 1)      # everyone exposed: 1 - 1/RR

# With the confidence limits of the relative risk.
paf <- attrfrac(1.08, 1, rrlo = 1.06, rrup = 1.09)
attr(paf, "ci")

# A protective exposure gives a negative fraction.
attrfrac(0.8, 0.5)

# A prevalence above 1 is refused.
try(attrfrac(1.5, 40))

Audit the categorical columns of a data frame for coding hazards

Description

Checks every categorical column of a data frame for the coding problems that silently change results: numbers that are really category codes whose labels were lost, labels still in foreign value-label form, labels that differ only by case or by an invisible space, empty strings or "NA" stored as a category, unused levels, and identifiers mistaken for categories. Run it once after every import, before any model. It also reports which level R will use as the reference category.

Usage

audcats(data, cols = NULL)

Arguments

Details

For each column the levels are levels(x) for a factor and the sorted distinct values (radix sort, as factor() would order them) otherwise. A hazard is reported when:

• the column is haven labelled (it still carries value labels and must be decoded with deccodes or declbld);

• every level looks like a number (^[0-9.]+$);

• every level starts with a code prefix such as "1. White";

• a level has leading or trailing whitespace, including non-breaking and other Unicode spaces, or two levels are equal once trimmed;

• a level is empty, or is a missing-value sentinel (NA, N/A, NAN, NULL, NONE, ., -, ?, compared case-insensitively);

• the reference (first) level is empty, a sentinel, or padded;

• two different trimmed levels agree once lower-cased;

• a factor has levels that never occur;

• there are more than 50 levels.

The function only reports; it never changes data. It errors when data is not a data frame.

Value

A data frame of class c("bricklayer_category_audit", "data.frame"), one row per audited column, with character columns column, storage (the column's class), levels (up to eight levels joined with "|"), reference (the first level) and hazards (all hazards joined with " ;; ", or ""), and integer n_levels. The attribute clean is TRUE when no hazard was found. The print method lists each column and its hazards.

See Also

deccodes, grdlvls, xfrvrf.

Examples

# codes without labels, and a case-variant city name
a <- audcats(data.frame(
  race = factor(c("1", "2", "2", "3")),
  city = c("Toronto", "toronto", "Ottawa")[c(1, 2, 3, 3)]
))
a
attr(a, "clean")

# padded, empty and sentinel labels, and an unused level
audcats(data.frame(
  g = factor(c(" A", "A", "NA", ""), levels = c("", " A", "A", "NA", "B"))
))

# a clean frame; non-categorical columns are skipped by default
audcats(data.frame(sex = c("F", "M", "F"), age = c(31, 45, 27)))

try(audcats(1:3))

Base64 encoding

Description

Encodes and decodes base64, in both the standard alphabet and the URL-safe variant. Computed in R with no dependency, so a capsule can embed binary content in a text manifest wherever this package runs.

Usage

jsb64enc(input)

jsb64dec(input)

jsb64uenc(input)

jsb64udec(input)

Arguments

Details

jsb64enc() breaks its output into 72-character lines (a newline after each complete block of 54 input bytes), matching jsonlite's encoder byte for byte; the decoder ignores line breaks, other white space and = padding, so either form round trips. Any other character outside the alphabet is an error rather than being skipped, because skipping it would decode a corrupted body into a shorter valid-looking one.

The URL-safe variant substitutes - and ⁠_⁠ for + and / and drops the = padding, which is what makes it safe in a URL path, a query string or a filename.

Base64 is an ENCODING, not encryption or a digest: it hides nothing and anyone can reverse it. Use csha256() to pin content and chmac() to authenticate it.

Errors: an encoder given anything but a raw vector, a character vector or NULL; jsb64dec() given anything but text or raw bytes, or text with a character outside the standard alphabet; jsb64udec() given anything but character text (raw bytes are refused), NA, or a character outside the URL-safe alphabet.

Value

The encoders return a character scalar ("" for empty input, NA_character_ for NULL); the decoders return a raw vector (use rawToChar() for text).

See Also

jsgzenc(), which composes this with gzip.

Examples

# Round trip through the standard alphabet.
b <- jsb64enc("hello capsule")
b
rawToChar(jsb64dec(b))

# Raw input works the same way.
jsb64enc(charToRaw("abc"))
jsb64dec(jsb64enc(charToRaw("abc")))

# Padding appears when the length is not a multiple of three.
jsb64enc("a")
jsb64enc("ab")
jsb64enc("abc")

# The URL-safe variant has no "+", "/" or "=" to escape.
jsb64uenc(as.raw(c(255, 224, 63)))
jsb64enc(as.raw(c(255, 224, 63)))
jsb64udec(
  jsb64uenc("path/safe?yes")
)

# Long input is wrapped, and the decoder ignores the breaks.
long <- jsb64enc(strrep("x", 200))
grepl("\n", long)
rawToChar(jsb64dec(long)) == strrep("x", 200)

# A character outside the alphabet is refused, not skipped.
try(jsb64dec("ab$c"))

The Goodman-Bacon decomposition of a two-way fixed-effects estimate

Description

Splits the two-way fixed-effects difference-in-differences estimate of a staggered adoption into the two-group, two-period comparisons it averages, with their weights: treated against never treated, earlier against later treated, and later against earlier treated (the comparisons that use already-treated units as controls).

Usage

bacondcmp(data, unit, time, y, trtd)

Arguments

Details

For a balanced panel with treatment that stays on and no covariates, the estimate of plndid is exactly \hat\beta = \sum_{k} s_{kU}\hat\beta_{kU} + \sum_{k \ne l} s_{kl}\hat\beta_{kl}, with weights that sum to one and grow with the groups' sizes n and the variance of their treatment over time \bar D(1 - \bar D) (Goodman-Bacon 2021, theorem 1). Units treated from the first period act as an already-treated group. A large weight on “later vs earlier” comparisons, together with effects that grow over time, is the warning sign that the two-way estimate is biased.

Value

A data frame of class rmbl_bacon: type ("treated vs never", "earlier vs later", "later vs earlier"), treated, control (the adoption periods; Inf for never), estimate, weight; attribute twfe (the weighted sum, equal to the two-way estimate).

References

Goodman-Bacon A (2021). Difference-in-differences with variation in treatment timing. Journal of Econometrics 225(2):254-277. DOI 10.1016/j.jeconom.2021.03.014

See Also

plndid, stgatt.

Examples

p <- polex()
b <- bacondcmp(p, "district", "month", "crmrate", "reformed")
b
tapply(b$weight, b$type, sum)

A covariate balance table

Description

Means of each covariate in the treated and comparison groups, their standardised mean difference, and Welch's test, optionally weighted (after matching or weighting).

Usage

baltab(data, grp, covs, w = NULL)

Arguments

Details

The standardised mean difference is (\bar x_1 - \bar x_0) / s_1 with s_1 the unweighted standard deviation in the treated group (the convention for the effect on the treated, so that it does not move when weights change the comparison group); |SMD| < 0.1 is the usual threshold for balance. With weights the means are weighted and the test uses the weighted variances and Kish's effective sample sizes.

Value

A data frame of class rmbl_balance: covariate, meantrtd, meanctl, smd, p_value.

References

Welch BL (1947). The generalization of ‘Student's’ problem when several different population variances are involved. Biometrika 34(1-2):28-35. DOI 10.1093/biomet/34.1-2.28

Examples

d <- data.frame(t = rep(0:1, each = 50),
                age = c(stats::rnorm(50, 30), stats::rnorm(50, 33)))
baltab(d, "t", "age")

How much a result depends on the open band's assumed cap

Description

Recomputes a statistic of banded data (by default the Gini coefficient) for a range of assumed caps on the open top band, and reports how far the answer moves. Every statistic computed from banded data depends on where the open band is assumed to end; this function measures that dependence, so you can report a range instead of a single figure when the dependence is large.

Usage

bandsens(bands, cnts, stat = ginicoef, caps = NULL, rule = "midpoint")

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

Arguments

Details

Let b be the smallest lower bound among the open top bands. The default caps are 12 points spaced evenly on the log scale from 1.05 b to 20 b. For each cap c, bandvals(bands, rule, openupcap = c) gives one value per band, each value is repeated cnts times, and stat is applied to the result. A statistic that errors or does not return a single number gives NA; warnings are suppressed.

The span is \max - \min over the finite values and the relative span is the span divided by the median absolute value (NA when any value is zero). The print method adds a note when the relative span exceeds 0.1.

Errors: cnts not numeric or of a different length from the bands; no open top band.

Value

A data frame of class c("rmbl_band_sensitivity", "data.frame") with numeric columns cap and value, and attributes span, relspan and base (the lower bound b of the open band).

Examples

b <- prsbands(c("1", "2 to 5", "6 to 10", "Greater than 10"))
counts <- c(1200, 430, 110, 38)

s <- bandsens(b, counts)
s
attr(s, "relspan")

# A statistic that barely moves has been measured; one that swings
# has not.
bandsens(b, counts, stat = mean, caps = c(12, 20, 50))

# nothing to vary without an open band
try(bandsens(c("1", "2 to 5"), c(10, 4)))

Representative values for banded categories

Description

Gives each band a single number that a calculation can use (a mean, a concentration index, a regression), under a rule you choose. For closed bands the rule places the value at the midpoint, the lower or upper bound, or the geometric mean of the bounds. An open band ("50+", “under 18”) has no midpoint, so its missing bound is filled from an explicit assumption, and the rows that rest on that assumption are flagged in the result instead of being absorbed into it.

Usage

bandvals(
  bands,
  rule = c("midpoint", "lower", "upper", "geometric"),
  openupcap = NULL,
  openupfac = 2,
  openloflr = NULL
)

Arguments

Details

For a band with bounds [l, u] the value is (l + u)/2 ("midpoint"), l ("lower"), u ("upper") or \sqrt{l u} ("geometric"; when l \le 0 the arithmetic midpoint is used instead, because the geometric mean is undefined there). The "lower" and "upper" rules bracket any value that could come from the bands.

An open upper band gets u = openupcap if given, else u = l \times openupfac. An open lower band gets l = openloflr if given, else \min(0, u). A cap below the band's lower bound, or a floor above its upper bound, is an error, because the value would lie outside the band. Unparsed bands give NA.

Use bandsens to see how much a result moves with the cap, and expbands to turn the values into one value per unit.

Value

The band table with two added columns: value (numeric, the representative value) and assumed (logical, TRUE for the open bands whose value rests on the cap or floor).

See Also

prsbands, bandsens, expbands

Examples

b <- prsbands(c("1", "2 to 5", "6 to 10", "Greater than 10"))
bandvals(b)

# The lower and upper rules bracket the truth.
bandvals(b, rule = "lower")$value
bandvals(b, rule = "upper", openupcap = 40)$value

# geometric placement, and an explicit floor for an open bottom band
bandvals(b, rule = "geometric", openupcap = 40)$value
bandvals(c("under 5", "5 to 9", "10+"), openloflr = 1)

# a cap inside the open band is refused
try(bandvals(b, openupcap = 5))

Benford first-digit test

Description

Compares the distribution of leading significant digits in x with Benford's law, P(d) = \log_{10}(1 + 1/d). Naturally occurring quantities that span several orders of magnitude follow it closely; figures that were rounded, truncated, capped, re-scaled, or invented typically do not. That makes it a cheap screen for a numeric column that arrived looking plausible but is not the measurement it claims to be.

It is a SCREEN, not a verdict. Columns with a narrow range, a unit floor or ceiling, or an assigned-identifier structure (postcodes, year fields, prices ending in 99) legitimately violate Benford's law. Treat a small p-value as a reason to look, never as evidence of fabrication.

Usage

benfordtst(x)

Arguments

Details

Zeros and non-finite values have no leading significant digit and are excluded; the sign is ignored.

With n usable values, O_d the count with leading digit d and E_d = n \log_{10}(1 + 1/d),

X^2 = \sum_{d=1}^{9} \frac{(O_d - E_d)^2}{E_d}

on 8 degrees of freedom. The smallest Benford share is 0.046 (digit 9), so the chi-square approximation needs about n \ge 110 for every expected count to reach 5; below that treat the p-value as rough. Digits are counted in compiled code in blocks of a million values, and a long vector shows a progress bar in an interactive session (silenced by options(morie.quiet = TRUE)). Errors: non-numeric input, or no finite non-zero value.

Value

A list of class bricklayer_benford with counts (observed digit frequencies 1–9), expected, proportion, statistic, df, p_value, logpvalue (the log of the upper tail, which keeps its information where p_value underflows to 0), and n (the number of values used). counts, expected and proportion are named numeric vectors with names "1" to "9". The print method shows the verdict and a bar chart of observed against expected shares.

References

Benford F (1938). The law of anomalous numbers. Proceedings of the American Philosophical Society 78(4), 551–572.

See Also

drftchisq(), cpsdrft().

Examples

# A quantity spanning several orders of magnitude follows the law.
set.seed(3)
benfordtst(10^stats::runif(2000, 0, 6))

# Digits drawn uniformly do not.
benfordtst(as.numeric(paste0(sample(1:9, 2000, TRUE), "000")))

# The expected proportions are the closed form.
b <- benfordtst(10^stats::runif(500, 0, 5))
all.equal(unname(b$expected) / b$n, log10(1 + 1 / (1:9)))

# Zeros carry no leading digit and are excluded from n.
benfordtst(c(0, 0, 1, 2, 3))$n

try(benfordtst(c(0, NA)))

Run the rmoriebricklayer command line

Description

Dispatches the verbs of the rmoriebricklayer launcher:

login [--token] [--email ADDRESS [--code CODE]] [--no-browser]

sign in for the curated tables at data.rmorie.com (a free MORIE data key, issued at www.rmorie.com/access/): a key you paste with --token (read from the terminal or a pipe when KEY is omitted), a code sent to --email, or the GitHub device flow

logout

forget the data key

doctor

report the language-model routes available here (own endpoint, each provider with your key, local Ollama), in the order ask tries them

models [--route ROUTE]

list the models your provider (or server) offers, the one in use marked

data list

the curated tables at data.rmorie.com

data pull db/table

download one table as CSV

app [PATH] [--port N] [--no-browser]

open the analysis dashboard (rmblapp) on this machine, with PATH, a CSV file or a capsule directory, loaded; needs the shiny package

config [setup | set KEY VALUE | unset KEY | get KEY | help]

show or change the settings of the language model (llmcfg): the route ask uses and the address, key and model of each route; setup walks through them

ask [--route ROUTE] [--model NAME] [--fresh] PROMPT...

send a prompt to the model (or the named one, on the named route) and print the reply; the same question again is answered from this machine (--fresh asks anew). aye and tellme are the same verb, and a quoted question with no verb (rmbl "how do I weight a survey?") is asked too

explain-error MESSAGE...

where an error of these packages comes from (explerr), without a model or the network

bundle REQUEST...

agentbndl from the shell

wayback FILE [--out FILE.jsonl] [--jobs N] [--within 30d] [--max N]

save every URL listed in FILE (one per line; blank lines and # comments skipped) in the Internet Archive's Wayback Machine with wbsave, results in --out (default FILE.wayback.jsonl), so a stopped run carries on; wayback --retry FILE.jsonl [--include-gone] retries the failures of an earlier run (wbrtry)

verify PATH [--json] [--fail-on LEVEL] [--public-key HEX]

check every capsule under PATH as the rmbl GitHub Action does (cpschk): each file against its JSON Schema (schemapath), vrfcps and cpsbndlvrf. The secret of an HMAC bundle is read from RMBL_KEY, which is then removed from the environment. Prints the Markdown report, or JSON with --json; exits 1 when a failed row is at or above --fail-on (error, the default, drift or warn). --base MANIFEST|DIR compares each manifest with its base copy, --manifests GLOB checks only the capsule files matching GLOB, --snapshot looks up each source's latest Wayback snapshot and --no-footprint leaves the footprint (co2fp) out of the report. --seen-log FILE adds the signatures an XMSS seen-log lists to those under PATH: an XMSS index one key signed twice fails its xmss_index row

drift PATH [--json] [--fail-on LEVEL] [--html FILE] [--app]

verify, then fetch each source again and compare it column by column (cpsdrft, Benford): a column past --psi-warn X (0.1) is a warning, a flagged one past --psi-fail X (0.25) drift. --allow-file accepts file:// sources, --cache DIR keeps each source and reads it when a fetch fails; the options of verify apply. --html FILE writes the interactive report (drfthtml) and --app opens it as a Shiny app (drftapp)

refresh PATH

pin each source's new bytes (cpsrfrsh) and print what changed (mnfdiff)

report PATH [--format md|json]

the verify report, exit 0; the options of verify apply

sign PATH [--scheme S] [--key-file F] [--overwrite]

pin each provenance on its first run and sign each capsule into capsule_bundle.json (cpsseal), then verify as verify does. The key is RMBL_KEY, removed from the environment at once, or the contents of --key-file, and is never printed: the shared secret for hmac (the default), else the JSON of a fipskg key, or its secret hex with --public-key, or the JSON of a pqckg key for xmss-sha256 (see the section XMSS keys). An existing bundle is refused without --overwrite

functions [PATTERN]

list the exported functions with their one-line titles, optionally filtered by a regular expression

describe NAME

title, usage, arguments and description of one function, from its help page

examples NAME

print the examples of one function

version

print the package version (also --version, -v)

help [VERB]

this list, or one verb's usage (also --help, -h; VERB --help shows a verb's usage without running it)

The launcher installed by instcli calls this function with the shell's arguments and exits with its status. Calling it from R is the same as typing the command, which is how the examples below work.

Usage

blcli(args = commandArgs(trailingOnly = TRUE), out = cat)

Arguments

Details

An option a verb does not take is refused with exit status 2 rather than ignored or sent as part of a prompt; everything after a bare -- is positional (for a prompt that itself starts with --) and only ask, bundle and explain-error accept it. A first argument that is plainly a question (two or more words, or ending in ?, with no path, line break or key-shaped word) and no verb is asked as ask would ask it.

Errors are not raised: every error is written to out as "PROG: message" and turned into the exit status. The status is 0 on success, 1 when the verb failed (an error, nothing found, no route), 2 for a usage error (an unknown verb or option, a missing value) and 130 when interrupted (Ctrl-C). The command name in messages is the environment variable RMBL_PROG, else rmoriebricklayer.

Side effects: during the call options(morie.progress = TRUE) is set, so downloads draw a progress bar; when RMBL_PROG is set (the launcher sets it) and options(rmoriebricklayer.dcache) is not, the data.rmorie.com manifest and tables are kept in tools::R_user_dir("rmoriebricklayer", "cache") rather than a fresh temporary folder. The verbs login, logout and config write the shared credentials and settings files; ask, models, data, bundle and wayback use the network, and drift and refresh fetch each capsule's source (refresh writes the capsule), sign writes each capsule's provenance and bundle; verify and report only read the capsule files; functions, describe, examples, explain-error, version and help only read the installed help.

Value

The exit status (integer scalar), invisibly: 0 on success, 1 on failure, 2 for a usage error, 130 when interrupted.

XMSS keys

An XMSS key (pqckg) signs each of its 2^height indices once; two signatures at one index let anyone forge a third. Its private key never changes, but its state, the indices already used, changes with every signature, and a copy of an old state signs used indices again (RFC 10033, “Hash-Based Signatures: State and Backup Management”, sections 2 and 4). rmbl sign --scheme xmss-sha256 --key-file KEY therefore treats the key file as the state:

• Before anything is signed it rewrites KEY, under a lock, through a temporary file in the same directory that is flushed to the disk and renamed over it, with next_index advanced and the state in a used field: highest (the highest index used), generation (one more at every write), writer (a digest of the host name and boot id), the last reserved block and, once split, sectors. A crash after that write loses the reserved indices and never reuses one. A KEY that cannot be rewritten (not writable, a read-only file system) is refused.

• --state-file FILE keeps the state in FILE (an absolute path to a writable regular file, or a new file in a writable directory) and only reads KEY; a key in RMBL_KEY needs it. The command line never uses the configuration directory that cpssign uses in an R session; it reads and updates the record in RMBL_XMSS_INDEX_FILE only when that is set.

• --sectors K --sector NAME (K a power of two, at most half the indices) splits the indices into K disjoint ranges and signs only in one of them, which keeps its own highest index; the range of NAME is SHA-256 of NAME and the key fingerprint, modulo K. Two names land in one sector by chance, so for real use give each runner its own --sector-index I. A sector that is used up is refused. Once split, sign only through rmbl sign.

• --reserve N marks the next N indices used in one write (at least one per capsule signed); the capsules of the run are signed from that block, and the indices it does not use are lost, as RFC 10033, section 5.9 allows.

• A stale copy is refused: when a bundle under PATH, or the --seen-log FILE that every rmbl sign appends its signatures to, shows this key signed an index above what the state says, nothing is signed. A note says when another host wrote the state last.

• When N or fewer indices are left in the range (--warn-left N, by default a tenth of the range) a warning is printed in the report and on standard error (RFC 10033, section 3.4).

A CI secret is a fixed copy of the key that every job starts from, so each job would sign the same indices again: in CI sign with ML-DSA (fipskg), or keep the XMSS state in a file that persists between jobs and give each concurrent runner its own sector. rmbl verify shows each XMSS signature's index and key fingerprint, and xmssidxseen lists indices signed twice. NIST SP 800-208 approves XMSS only when the key is made and used inside a hardware cryptographic module that never exports it; a key file is software, so this follows the software guidance of RFC 10033 and is not a validated SP 800-208 deployment.

See Also

instcli to put the command on your PATH.

Examples

blcli("version")
blcli("help")

# one verb's usage
blcli(c("help", "data"))

# capture the output instead of printing it
con <- textConnection("lines", "w")
status <- blcli(c("functions", "^drbg_"), out = function(x) cat(x, file = con))
close(con)
status
lines

# an unknown verb is a usage error: status 2, nothing run
blcli("nonsense")

# check a capsule as the rmbl GitHub Action does
d <- usecpstmpl(tempfile("capsule-"), ex = TRUE)
invisible(cpsrfrsh(d))
blcli(c("verify", d))


# these read every help page, which takes a few seconds
blcli(c("functions", "json"))
blcli(c("describe", "anltab"))
blcli(c("examples", "perdays"))

Download a file with live progress

Description

Downloads one file from the web to a local path, showing a progress bar (percent, size and rate) while a person is watching. It is the single download routine of the package and its sibling packages, and it applies the package's network safety rules: HTTPS only, no local or private addresses, a size cap, and checked redirects. Use it when a script must fetch a source file and you want a visible, safe transfer.

Usage

bldl(
  url,
  dest,
  headers = NULL,
  label = basename(dest),
  size = NULL,
  tmout = 3600,
  quiet = NULL,
  tty = NULL,
  allowfile = FALSE,
  maxbytes = 2^31
)

Arguments

Details

Progress is an animated bar, redrawn in place, when the console can redraw a line: a terminal, or the RStudio, Positron, R.app or RKWard console (never a knitr or Jupyter notebook; the options rmoriebricklayer.dynamic or cli.dynamic, or the environment variable R_CLI_DYNAMIC, override the choice). In the RStudio or Positron console the bar is drawn on stdout, elsewhere on stderr. Otherwise milestone lines go to stderr (every 10 percent, or every 50 MB when the size is unknown), and nothing is shown when quiet. By default output is shown in an interactive session or when options(morie.progress = TRUE) is set, and suppressed by options(morie.quiet = TRUE).

Network: the transfer runs in the package's compiled HTTP client. The URL is checked, its host resolved and pinned, and every redirect is checked again (no switch from https to http; credentials are not passed to another host). Plain http:// is allowed only when rmoriebricklayer.allowhttp, an option, is TRUE. A transfer slower than 64 bytes a second for 30 seconds is abandoned. A file:// URL is copied locally through the same progress display, only with allowfile = TRUE.

Errors: a URL that is not https (or allowed http/file), or that points at a loopback, link-local or private address; headers not a named character vector; tmout below 1 or above 2147483647; maxbytes below 1; a network failure; an HTTP status of 400 or more; a body larger than maxbytes. dest is overwritten. A server busy for a moment (HTTP 429, 500, 502, 503 or 504) is asked again up to three times, after the wait its Retry-After header asks for (at most a minute) or else 1, 2 and 4 seconds.

Value

dest, invisibly.

See Also

blfetch for a download with a Wayback Machine fallback.

Examples

# offline: a local file through the same routine
src <- tempfile(fileext = ".txt")
writeLines("hello", src)
dest <- tempfile(fileext = ".txt")
bldl(paste0("file://", src), dest, quiet = TRUE,
                    allowfile = TRUE)
readLines(dest)

# the safety checks run before any connection is made
try(bldl("http://example.org/table.csv", tempfile()))
try(bldl("https://127.0.0.1/table.csv", tempfile()))
try(bldl(paste0("file://", src), tempfile()))

# a real download needs the network


bldl("https://www.ontario.ca/some/open-data/table.csv",
                    "table.csv", maxbytes = 50e6)

Fetch a URL to disk with an Internet Archive fallback (C++/libcurl)

Description

The shared data-fetch foundation of the morie ecosystem. Downloads url to dest; if the live download fails (404, network), it retries the Wayback Machine snapshot – so a rotated or removed source file (CIHI, open-data portals) stays retrievable. Backed by the package's C++ kernel rmbl_fetch_with_fallback (libcurl), which rmorie and morie share through LinkingTo. The return value says which source served the file, so a report can note that it used an archived copy.

Usage

blfetch(url, dest, wb = "", tmout = 120L)

Arguments

Details

Before any network use the arguments are checked: url must be one https URL to a public host (plain http only when the option rmoriebricklayer.allowhttp is TRUE; loopback, private, link-local and metadata addresses are refused), the directory of dest must exist and be writable, and tmout must be at least one second. The compiled layer then resolves the host, pins the connection to the checked addresses and re-checks every redirect.

If the live download fails, the fallback is the wb URL when given, else a snapshot looked up through the Wayback availability API (as wbsnpnat() does). dest is overwritten. A spinner is drawn on stderr while waiting, unless options(morie.quiet = TRUE).

Errors: invalid arguments as above; and, when both sources fail, an error saying whether the live URL could not be reached or answered an HTTP status, and whether the Wayback Machine had no snapshot or its copy could not be downloaded.

Value

Invisibly, "live" or "wayback": the source that served the file. Throws on total failure.

Examples

# Inputs are validated before any network access:
try(blfetch("", tempfile()))          # empty url -> error
try(blfetch("ftp://example.org/data.csv", tempfile()))
try(blfetch("https://127.0.0.1/data.csv", tempfile()))
try(blfetch("https://example.org/data.csv",
                     file.path(tempdir(), "no-such-dir", "data.csv")))


# needs the network
dst <- tempfile(fileext = ".xlsx")

# Live download; auto-resolves a Wayback snapshot only if the live URL fails.
blfetch(
  paste0("https://www.cihi.ca/sites/default/files/document/",
         "hospital-beds-2024-2025-data-tables-en.xlsx"),
  dst)

# Pin an explicit Wayback snapshot to fall back to, and a shorter timeout.
blfetch(
  "https://example.org/rotated-file.csv", tempfile(fileext = ".csv"),
  wb = "https://web.archive.org/web/2024id_/https://example.org/rotated-file.csv",
  tmout = 60)

# The return value tells you which source served the file.
status <- blfetch("https://cloud.r-project.org/", tempfile())
status   # "live" or "wayback"

Plot styling for the package's plot methods

Description

Every plot() method in the package draws with base graphics and the same small set of style choices, which these functions expose. blpal() returns colour-blind-safe palettes: "okabe" (the Okabe-Ito categorical set), "diverging" (blue, grey, orange), "sequential" (light to dark blue) and "grey". plotopt() sets the defaults every method uses, for the session: palette, grid, cex, las and bg. Each method also takes main, col (which overrides the palette), palette and ... passed to the underlying base-graphics call, so a single plot can be changed without changing the session. Use blpal() also to colour your own plots so they match the package's.

Usage

blpal(name = c("okabe", "diverging", "sequential", "grey"), n = NULL)

plotopt(...)

Arguments

Details

The palettes: "okabe" is the eight-colour set of Okabe and Ito (2008), distinguishable under the common forms of colour-vision deficiency. "diverging" has three anchors (blue #2166AC, grey #8C8C8C, orange-red #D6604D) and is used for changes (decrease, no change, increase). The palette "sequential" has five blues from light to dark; "grey" five greys from dark to light. With n given, "okabe" and "grey" repeat their colours (rep_len()) and "diverging" and "sequential" are interpolated with grDevices::colorRampPalette().

The session options are stored in options(rmoriebricklayer.plot); the defaults are palette = "okabe", grid = TRUE, cex = 1, las = 1, bg = "white" and panel = "#EEF1F6". Setting some options keeps the others.

Text is drawn at the session cex. Titles, axis titles, row labels and keys are made smaller to fit a small device, but never below 0.6 of that size, so that they stay readable; for larger text set, for example, plotopt(cex = 1.2), and for a table with many rows use a taller device.

Errors: an unknown palette name; n not a positive whole number; an option name other than the six above, or an unnamed value; an unknown palette option.

Value

blpal(): a character vector of hexadecimal colours. plotopt(): a named list of the six options in force (after the change), invisibly when setting.

References

Okabe, M. and Ito, K. (2008). Color Universal Design (CUD): how to make figures and presentations that are friendly to colorblind people. jfly.uni-koeln.de/color

Examples

blpal("okabe", 3)
old <- plotopt(palette = "grey", grid = FALSE)
plotopt()
plotopt(palette = old$palette, grid = old$grid)

# an interpolated sequential scale, shown as swatches
cols <- blpal("sequential", 7)
barplot(rep(1, 7), col = cols, border = NA, axes = FALSE, names.arg = cols, las = 2,
        cex.names = 0.6)

# the same palette in a base plot of your own
plot(1:8, pch = 19, cex = 3, col = blpal("okabe"))

# names are checked
try(blpal("viridis"))
try(plotopt(colour = "red"))

Live progress for long-running work

Description

Makes a progress bar for a loop or a long computation: the animated bar of bldl(), redrawn in place at no more than 12 frames a second, when the console can redraw a line – a terminal, or the RStudio, Positron, R.app or RKWard console (never a knitr or Jupyter notebook; the options rmoriebricklayer.dynamic or cli.dynamic, or the environment variable R_CLI_DYNAMIC, override the choice). In the RStudio or Positron console the bar is drawn on stdout, elsewhere on stderr. Otherwise a milestone line goes to stderr at each tenth, and nothing is shown when quiet. With total = 0 the work has no count, and the bar is a spinner with the time so far and the current stage (the info text), with a line at each new stage when the console cannot redraw.

Usage

blprog(label, total, quiet = NULL, color = "blue", delay = NULL)

Arguments

Details

Progress is shown when a person is watching: an interactive session, or a command line that sets options(morie.progress = TRUE). options(morie.quiet = TRUE) or quiet = TRUE silences it.

Nothing is drawn until the work has run for delay seconds, so a fast call never flashes a bar. A bar made inside a function while a bar made by a calling function is still running stays silent, so nested work shows one bar, for the outer loop.

Updates are cheap: when quiet an update only stores the value, and when drawing, frames are throttled, so $tick() can be called on every iteration of a loop. Call $done() through on.exit() right after making the bar, so the line is cleared and the bar is released even when the work fails or is interrupted:

p <- blprog("bootstrap", B)
on.exit(p$done(), add = TRUE)
for (b in seq_len(B)) {
  ...
  p$tick(info = sprintf("resample 
}

$done() may be called more than once; only the first call acts.

Value

An object of class rmbl_progress, a list of functions: update(done, info = NULL) sets the steps done (and the text after the bar), tick(by = 1, info = NULL) adds by steps, wait(seconds) sleeps while keeping the animation moving, and done(summary = NULL) ends the bar, clears the line and, when the bar was shown, writes summary on its own line. Each returns the object invisibly.

Examples

p <- blprog("squares", 50, delay = 0)
for (i in 1:50) {
  x <- sum(seq_len(i)^2)
  p$tick(info = sprintf("step %d/50", i))
}
p$done("squares: 50 steps")

# work with no count: a spinner and the stage
f <- function() {
  p <- blprog("model", 0, delay = 0)
  on.exit(p$done(), add = TRUE)
  p$update(info = "fitting")
  p$wait(0.2)
  p$update(info = "writing")
  invisible(TRUE)
}
f()

Treated, comparison and excluded places around sites

Description

Assigns each place (a block, an address) to the treated group when it lies within a radius of a site (a store that donated, a new facility), to the comparison group when it lies within the same radius of a comparison point and farther than a buffer from every site, and excludes it otherwise, so that spillover from the sites does not reach the comparison.

Usage

bufasgn(x, y, sitex, sitey, comx, comy, radius, buffer = 3 * radius)

Arguments

Details

Distances are Euclidean, so coordinates must be projected (metres or feet, not degrees). The radius and the buffer are in the same unit. The buffer should be several times the radius: the design it follows used 330 feet and 990 feet.

Value

A data frame with the columns group ("treated", "comparison" or "excluded"), site_distance and comparison_distance.

References

Applied to blocks near donor stores: Zorro Medina, Hackett, Green and Vargas (2025), see evstudy.

Examples

set.seed(4)
x <- runif(200, 0, 5000)
y <- runif(200, 0, 5000)
a <- bufasgn(x, y, c(1000, 3500), c(1200, 3000), c(2500, 4500), c(500, 4000),
                   radius = 330)
table(a$group)

Where each loaded package came from

Description

Records, for every package named, its version and enough about its provenance to find the same one again: the library it was loaded from, the repository it was installed from, and – for a package installed from a remote – the remote's URL and commit. Store it in a manifest next to capenv() when an analysis depends on packages installed from GitHub or a local build.

Usage

capdeps(pkgs = loadedNamespaces())

Arguments

Details

Why versions alone are not enough. Two installations can report the same version and differ: one built from CRAN, one from a fork, one from a local R CMD INSTALL of a working tree with uncommitted changes. A version number identifies an intention; the repository and commit identify what was actually loaded.

The fields are read with utils::packageDescription() from each package's installed DESCRIPTION: Repository (for example "CRAN"), RemoteUrl and RemoteSha (written by remotes, pak and devtools for packages installed from a Git host) and Built. Packages are sorted and de-duplicated. A progress bar is drawn on stderr for many packages unless options(morie.quiet = TRUE). Nothing is written. Error: pkgs that is not a character vector.

Value

A data frame with one row per package and character columns package, version, library (the library directory), repository, remote_url, remote_sha and built (R version, platform and time of the build). Unavailable fields are NA rather than omitted, so a reader can see that the information was absent rather than forgotten; for a package that is not installed, every field except package is NA. An empty packages gives a data frame with no rows and the same columns.

See Also

capenv(), envdiff().

Examples

deps <- capdeps(c("stats", "utils"))
deps[, c("package", "version")]

# a package with no repository field records that fact
is.na(capdeps("stats")$repository)

# a package that is not installed is a row of NAs, not an error
capdeps("notapackage123")[, c("package", "version")]

# names must be character
try(capdeps(1))

Capture the Analysis Environment for a Manifest

Description

Records the facts a replicator needs to rebuild the session in which a result was produced: the R version, platform and operating system, a UTC timestamp, the random number generator in use and whether it was seeded, and the versions of the packages involved. Store the result in a results manifest (for example with wrmnfjs()) next to the numbers it produced, and compare two captures with envdiff().

Usage

capenv(pkgs = loadedNamespaces())

Arguments

Details

Package versions come from utils::packageVersion(); a package that is not installed gets NA rather than an error. The packages are sorted and de-duplicated. The generator kinds come from RNGkind() (R has changed its default sample() algorithm before, so a recorded seed means nothing without the kind it was fed to), and rng_seeded says whether .Random.seed exists in the global environment. Nothing is written and the random state is not touched.

Value

A list of seven elements:

r_version

character, e.g. "4.5.1".

platform

character, R.version$platform.

os

character, the system name and release.

captured_utc

character, ISO 8601 time YYYY-MM-DDThh:mm:ssZ.

rng_kind

character, the three RNGkind() values joined by commas.

rng_seeded

logical, whether .Random.seed exists.

packages

named list of version strings, sorted by name (NA for a package that is not installed). A list, not a vector, so the names survive the manifest's JSON: a named vector is written as a bare array.

Examples

# Record specific packages' versions alongside the session facts.
env <- capenv(c("stats", "utils"))
env$r_version
env$os
env$packages          # named list of versions

# Default captures every currently loaded namespace.
names(capenv())

# the generator and whether it was seeded
set.seed(42)
capenv("stats")[c("rng_kind", "rng_seeded")]

# a package that is not installed is recorded as NA, not an error
capenv(c("stats", "notapackage123"))$packages

Case flow through a court: clearance, backlog and time to clear

Description

Turns the counts a court publishes each period – cases received, cases disposed and, when it is published, cases pending – into the measures used to judge whether a court is keeping up: the clearance rate (cases disposed per case received), its rolling value over a window, the net flow into the inventory and, with the pending count, the number of periods the current inventory would take to clear at the recent pace of disposition.

Usage

caseflow(recv, disp, pend = NULL, prd = NULL, win = 12L, conflvl = 0.95)

Arguments

Details

A clearance rate above 1 means the court disposed of more cases than it received, so its inventory fell. The interval treats received and disposed as Poisson counts over the same period: given their sum, the disposed count is binomial, and the exact (Clopper-Pearson) limits for that proportion give exact limits for the ratio. The interval measures the period-to-period noise of counts this size; it does not make received and disposed the same cases, which they are not (a case disposed this month was usually received months before).

pertoclear divides the pending inventory by the mean number disposed over the last win periods (fewer at the start of the series). The published pending count need not equal the previous count plus the net flow – courts reactivate cases on bench warrants and transfer them – so penchg is reported next to net rather than derived from it.

The result carries attribute "caseflow": win, conflvl and the totals received and disposed. Errors: recv empty; a count vector that is not numeric, has a missing, negative or infinite value, or a length other than recv's; prd of another length, with missing or repeated values; win not a whole number of at least 1; conflvl outside (0, 1).

Value

A data frame of class rmbl_case_flow, one row per period, with period, received, disposed, clearance with lower and upper, net (received less disposed: positive grows the inventory), rolclr with rollo and rolup (NA until a full window is available), and either pending, penchg and pertoclear, or, without pending, cumnet (the inventory change since the first period).

References

Ostrom, B. J. and Hanson, R. A. (1999). Efficiency, Timeliness, and Quality: A New Perspective from Nine State Criminal Trial Courts. National Center for State Courts, Williamsburg, VA. (Clearance rate and time-to-disposition measures.)

Ulm, K. (1990). A simple method to calculate the confidence interval of a standardized mortality ratio. American Journal of Epidemiology, 131(2), 373-375. (Exact limits for a ratio of Poisson counts through the beta distribution.)

See Also

lifetbl, dispmix, stkflow for custody populations.

Examples

f <- read.csv(system.file("extdata", "ocj_flow_ontario.csv",
                          package = "rmoriebricklayer"))
cf <- caseflow(f$received, f$disposed, pend = f$pending, prd = f$month)
cf
plot(cf)

# without the pending count, a quarterly rolling window and 90% intervals
cq <- caseflow(f$received, f$disposed, prd = f$month, win = 3, conflvl = 0.90)
head(as.data.frame(cq)[, c("period", "clearance", "rolclr", "cumnet")])

# one value per period in every count
try(caseflow(c(100, 120, 90), c(95, 130)))

BLAKE2b digest, optionally keyed

Description

Computes the BLAKE2b hash of text or bytes as lowercase hex, with any digest length from 1 to 64 bytes. With a key it becomes a message authentication code: only someone holding the key can produce or check the digest. Use it for fingerprints and keyed checks inside your own workflow; use csha256() or chmac() where other tools must reproduce the value.

Usage

cblake2b(x, key = NULL, length = 32L)

Arguments

Details

BLAKE2b (RFC 7693): a modern cryptographic hash, faster than SHA-256 in software, that takes a key natively and produces any digest length from 1 to 64 bytes.

The native key is the interesting part. A keyed BLAKE2b IS a message authentication code, with no HMAC wrapper and so no doubled hashing – cblake2b(msg, key = k) does the job of chmac() at lower cost. Use SHA-256 or HMAC where interoperability with other tools matters; use this where it does not and speed does.

Character input is hashed element by element, as the bytes of each string; NA elements give NA. A raw vector is one message. The digest length is a parameter of BLAKE2b, so a 32-byte digest is not the first half of the 64-byte digest of the same input. Errors: length outside 1 to 64; a key that is not a single string or raw vector, or is longer than 64 bytes.

Value

A character vector of lowercase hex digests – one per element for character input, length-1 for raw input.

References

Saarinen MJ, Aumasson JP (2015). The BLAKE2 Cryptographic Hash and Message Authentication Code (MAC). RFC 7693. DOI 10.17487/RFC7693

See Also

csha256(), csha512(), chmac()

Examples

# The RFC 7693 test vector for BLAKE2b-512 of "abc".
cblake2b("abc", length = 64)

# 32 bytes by default, the same width as SHA-256.
cblake2b("abc")
nchar(cblake2b("abc"))

# Any digest length, which SHA-2 cannot do.
cblake2b("abc", length = 8)

# Keyed, so it authenticates without an HMAC construction.
cblake2b("manifest", key = "secret")
cblake2b("manifest", key = "secret") ==
  cblake2b("manifest", key = "other")

# Vectorised, and raw input hashes the same bytes.
cblake2b(c("a", "b"))
identical(cblake2b("abc"), cblake2b(charToRaw("abc")))

try(cblake2b("abc", length = 65))

Bootstrap replicate means (C backend)

Description

B resamples of x, drawn with replacement and each the same length as x, with the mean of every resample returned. The spread of these replicate means estimates the standard error of the sample mean, and their quantiles give a percentile confidence interval, without assuming the data are normal (Efron and Tibshirani, 1993). The resampling uses the core's own 64-bit Mersenne Twister seeded by seed, NOT R's RNG, so a given seed reproduces the same replicates in every binding of the core and R's own random stream is left untouched.

Usage

cbootmean(x, B = 1000L, seed = 42L)

Arguments

Details

Replicate b draws n indices uniformly from 1, \dots, n with replacement (C++ std::mt19937_64 seeded with seed, then std::uniform_int_distribution, whose algorithm is the C++ library's own, so the replicates for a seed can differ between libstdc++ and libc++) and returns the mean of the drawn values, accumulated in extended precision. The replicates are a single stream: the first B replicates for a seed are the same whatever larger B is asked for later.

Edge cases: an empty x gives B values of NaN; an NA in x makes every replicate that draws it NA, so remove missing values first. B must be a whole number of at least 1 ("B must be a single whole number >= 1"); a character, factor or list x is refused ("must be numeric"). seed must be a whole number from 0 to below 2^{53}: NA, a negative or a fractional seed is refused.

Performance: the work is B \times n draws (about 0.4 s for B = 2000, n = 20000 on a laptop). The kernel polls for a Ctrl-C after about every million draws, so a long run can be interrupted; while it runs, an interactive session (or one with options(morie.progress = TRUE)) shows a spinner on stderr from 0.75 s on, and options(morie.quiet = TRUE) silences it.

Value

A numeric vector of length B: the replicate means, in the units of x.

References

Efron B, Tibshirani RJ (1993). An Introduction to the Bootstrap. Chapman and Hall.

See Also

exprspplr() for a bootstrap of a model estimate.

Examples

set.seed(1)
x <- stats::rnorm(50, mean = 5)

reps <- cbootmean(x, B = 500, seed = 7)
length(reps)

# The replicates centre on the sample mean, and their spread estimates
# the standard error.
c(sample = mean(x), bootstrap = mean(reps))
c(bootstrapse = stats::sd(reps),
  formulase = stats::sd(x) / sqrt(length(x)))

# A percentile confidence interval for the mean.
stats::quantile(reps, c(0.025, 0.975))

# Reproducible: the same seed gives the same replicates, and R's own
# random stream is not consumed.
identical(cbootmean(x, 100, seed = 1),
          cbootmean(x, 100, seed = 1))

# A longer run continues the same stream.
identical(cbootmean(x, 10, seed = 1),
          cbootmean(x, 100, seed = 1)[1:10])

# At least one replicate is needed.
try(cbootmean(x, B = 0))

Covariance matrix of a numeric matrix (C backend)

Description

The covariance of every pair of columns of a numeric matrix or data frame, computed by the compiled core. Covariances describe how variables move together; the matrix is the input to principal components, Mahalanobis distances and multivariate drift checks. Column names are carried through to both dimensions of the result.

Usage

ccovm(x)

Arguments

Details

Entry (j, k) is \frac{1}{n - 1} \sum_i (x_{ij} - \bar x_j)(x_{ik} - \bar x_k), the n - 1 denominator of stats::cov(), with which it agrees to rounding. The diagonal holds the column variances.

Missing values are not removed: an NA in a column makes every covariance involving that column NA (pairwise or complete-case deletion is not done; use stats::cov() with ⁠use =⁠ for those). A single row gives NaN throughout. Refusals: a vector that is not a matrix ("must be a matrix or data frame") and a data frame with a non-numeric column ("every column of x must be numeric"). Logical and integer matrices are converted to double. The cost is O(n p^2) for n rows and p columns.

Value

A symmetric ncol(x) by ncol(x) numeric matrix, with the column names of x (if any) as both row and column names; entry units are the product of the two columns' units.

See Also

ccorr() for a single Pearson correlation, cmom() for one column's variance.

Examples

X <- cbind(a = c(1, 2, 3, 4), b = c(2, 4, 7, 8), c = c(5, 3, 2, 1))
ccovm(X)

# Agrees with stats::cov().
all.equal(ccovm(X), stats::cov(X))

# The diagonal is the column variances.
all.equal(diag(ccovm(X)), apply(X, 2, stats::var))

# Data frames are accepted.
ccovm(data.frame(u = 1:5, v = c(2, 1, 4, 3, 6)))

# A plain vector is not a matrix of variables.
try(ccovm(1:5))

CRC-32 checksum (C backend)

Description

Computes the CRC-32 checksum of character or raw input: the CRC of ITU V.42, zip and PNG (reflected polynomial 0xEDB88320, initial value and final XOR 0xFFFFFFFF). Use it as a fast first check that a large file or capsule member arrived intact, or to compare with the CRC a zip archive records for its members.

Usage

ccrc32(x)

Arguments

Details

A CRC is NOT a cryptographic digest: it detects accidental corruption – a truncated download, a flipped bit on disk – but anyone can construct a different input with the same value, so it must never be used where csha256() is meant. It is here because it is far cheaper than SHA-256 over gigabyte-scale capsule members, which makes it the right first pass when the question is only “did this file arrive intact”.

Character input is checked element by element, as its bytes, with NA giving NA; other atomic input is converted with as.character(). A raw vector gives one checksum. For files use crc32file().

Value

A numeric vector of unsigned 32-bit checksums (one per element for character input, NA for NA; length 1 for raw input). Returned as double rather than integer because values above 2^31 - 1 are not representable as an R integer.

See Also

csha256() for integrity against tampering rather than accident.

Examples

# The standard check value: CRC-32("123456789") is 0xCBF43926.
ccrc32("123456789")
ccrc32("123456789") == 0xCBF43926

# The empty input has checksum zero.
ccrc32("")

# Vectorised, and sensitive to a single changed byte.
ccrc32(c("brick", "brack"))

# Raw bytes work the same way.
ccrc32(charToRaw("123456789"))

Verify a certificate chain

Description

Builds the path from a leaf certificate to a trust anchor you supply, verifying each signature, each validity window, and the constraints that stop a certificate being used for something it was not issued for.

Usage

certchnvrf(
  leaf,
  trust,
  inter = list(),
  attm = Sys.time(),
  purp = NULL,
  crls = list(),
  pols = NULL,
  rvk = c("supplied", "fetch", "none")
)

Arguments

Details

attm is the point the validity windows are checked against, and it matters which one you pass. For a timestamp token the right answer is the time the token asserts: a token signed in 2020 by a certificate that expired in 2021 was validly signed, and checking it against today would reject it for no good reason. For a signature you are being asked to rely on now, pass now.

The path is built by following issuer names from the leaf through the intermediates until a certificate in trust is reached, at most ten steps (a cycle among supplied certificates cannot loop). Then every link is checked: the validity window at attm, the signature under the issuer's key (RSA PKCS#1 v1.5 with SHA-1, SHA-256, SHA-384 or SHA-512; ECDSA with SHA-256, SHA-384 or SHA-512 on P-256, P-384 and P-521), that every issuer is a CA with a path length that allows the chain, name constraints, the policy rules of RFC 5280 section 6.1, the requested purp (extended key usage of the leaf), that the anchor is self-signed, and revocation as rvk says. Each check is one row of checks; ok is TRUE only when every row passed. A failed check is a row, not an error.

Errors: trust empty; a certificate that certprs() cannot read; an unknown rvk. With revocation = "fetch" the network is used (distribution points and OCSP responders named in the certificates); otherwise nothing is fetched.

Value

A list of class bricklayer_certpath_check: ok (logical, the verdict), path (the list of certificates from the leaf to the anchor, as built), checks (a data frame with check, logical ok and detail, one row per check), attm (the time used) and policies (the outcome of policy processing). print() shows the verdict and the checks.

See Also

certprs(), tsvrf().

Examples

# A small test chain shipped as hexadecimal DER: a self-signed P-256 root
# and a leaf it issued (a throwaway test authority; trust it for nothing).
roothex <- paste0(
  "308201973082013ca003020102021452157dfd900c163c5027372ed4090b9d11748286300a06",
  "082a8648ce3d040302302931123010060355040a0c09726d626c207465737431133011060355",
  "04030c0a454332353620526f6f74301e170d3236313030323038343334315a170d3336303932",
  "393038343334315a302931123010060355040a0c09726d626c20746573743113301106035504",
  "030c0a454332353620526f6f743059301306072a8648ce3d020106082a8648ce3d0301070342",
  "0004624ad2e98f4566cd8f003a188f822875fa4760ed894507e8a4130a08fcb59280e958aa6c",
  "e6147989d155ce20d1a6d94fe109d3d2f2bc13d1394cea3520197fdda3423040300f0603551d",
  "130101ff040530030101ff300e0603551d0f0101ff040403020106301d0603551d0e04160414",
  "2ef598166f5b9e6bb65d3252cf127b2c992348ae300a06082a8648ce3d040302034900304602",
  "21009cd3d7968f65630677b3c7437cae5187021d553c756820d34b6a8a1240e7efb7022100a3",
  "c37f0a8e1f38a5263b14414e6f916962bbd36bcbaca40b07d371352de745c0")
leafhex <- paste0(
  "308201ae30820154a00302010202142fa99d48ead4c033c568837b2288f9e9bf137068300a06",
  "082a8648ce3d040302302931123010060355040a0c09726d626c207465737431133011060355",
  "04030c0a454332353620526f6f74301e170d3236313030323038343334315a170d3331313030",
  "313038343334315a302931123010060355040a0c09726d626c20746573743113301106035504",
  "030c0a4543323536204c6561663059301306072a8648ce3d020106082a8648ce3d0301070342",
  "00045f0ad14e657dedf3f6cc18530f9d905823b6e9f492ecc0eb1989c9ad8a11a859ebb7681a",
  "e5b1f579974e3c6709e824245608c448da6fe248ebf67bf727d598f4a35a305830090603551d",
  "1304023000300b0603551d0f040403020780301d0603551d0e04160414b1d22b6dd5987f0dfa",
  "0ff3a5d0140351d752e0fb301f0603551d230418301680142ef598166f5b9e6bb65d3252cf12",
  "7b2c992348ae300a06082a8648ce3d0403020348003045022100bc9d6ec452554adb8f33e4ac",
  "53d5e29620fc19528074153028ce849d98c8c15202205c841de357dbad5042daba7fb4af3724",
  "78d684f648e6ef23c68a3d764894c8fc")
hex2raw <- function(h) {
  as.raw(strtoi(substring(h, seq(1, nchar(h), 2), seq(2, nchar(h), 2)), 16L))
}
root <- certprs(hex2raw(roothex))
leaf <- certprs(hex2raw(leafhex))

# the leaf chains to the root, checked a day after it was issued
res <- certchnvrf(leaf, trust = root, attm = leaf$not_before + 86400)
res$ok
res$checks[, c("check", "ok")]

# after the leaf expired, the same chain fails on validity
certchnvrf(leaf, trust = root, attm = leaf$not_after + 86400)$ok

# a purpose the leaf was not issued for
r <- certchnvrf(leaf, trust = root, attm = leaf$not_before + 86400,
                       purp = "timeStamping")
r$checks[!r$checks$ok, ]

# a chain must end at something you trust
try(certchnvrf(leaf, trust = list()))


# with revocation looked up on the network
certchnvrf("tsa.crt", trust = "ca.crt", purp = "timeStamping",
                  rvk = "fetch")

Parse an X.509 certificate

Description

Reads the fields a verifier needs: who issued it, who it is for, when it is valid, what key it carries, and what that key is allowed to do. Use it to inspect the certificate of a timestamping authority or a signer before you decide to trust it, or as the first step of certchnvrf(). The DER decoder is the package's own; no system library or network is used.

Usage

certprs(cert)

Arguments

Details

The certificate is decoded as DER (RFC 5280 section 4.1); a file whose text contains a -----BEGIN line is read as PEM and its base64 body decoded first (only the first certificate in the file is used). Version 1 certificates, without the explicit version field, are accepted. Extensions read: basic constraints, key usage, extended key usage, subject alternative name, name constraints, certificate policies, policy mappings and policy constraints, and inhibit any-policy. The signature is not checked here: that needs the issuer's key, which certchnvrf() supplies.

Errors: cert neither a raw vector nor the path of an existing file (a directory is refused); bytes that are not well-formed DER or carry trailing bytes; a DER structure that is not a certificate or has no signature.

Value

A list of class bricklayer_certificate with these elements:

• version (integer) and serial (hex);

• issuer and subject (distinguished names as text);

• not_before and not_after (POSIXct, UTC);

• signature_algorithm (a name such as "ECDSA-SHA256") and signature_oid;

• signature and tbs (raw; the signed part);

• key (a list: type "RSA" or "EC", bits, and the key values, with curve for EC);

• is_ca and path_len;

• key_usage and extended_key_usage (character; EKU as OIDs);

• subject_rdns, issuer_rdns, san and name_constraints;

• the policy fields policies, policy_mappings and require_explicit_policy;

• the policy fields inhibit_policy_mapping and inhibit_any_policy;

• self_signed (logical: issuer equals subject);

• der (the bytes).

print() shows a summary.

See Also

certchnvrf(), tsvrf().

Examples

# A small test chain shipped as hexadecimal DER: a self-signed P-256 root
# and a leaf it issued (a throwaway test authority; trust it for nothing).
roothex <- paste0(
  "308201973082013ca003020102021452157dfd900c163c5027372ed4090b9d11748286300a06",
  "082a8648ce3d040302302931123010060355040a0c09726d626c207465737431133011060355",
  "04030c0a454332353620526f6f74301e170d3236313030323038343334315a170d3336303932",
  "393038343334315a302931123010060355040a0c09726d626c20746573743113301106035504",
  "030c0a454332353620526f6f743059301306072a8648ce3d020106082a8648ce3d0301070342",
  "0004624ad2e98f4566cd8f003a188f822875fa4760ed894507e8a4130a08fcb59280e958aa6c",
  "e6147989d155ce20d1a6d94fe109d3d2f2bc13d1394cea3520197fdda3423040300f0603551d",
  "130101ff040530030101ff300e0603551d0f0101ff040403020106301d0603551d0e04160414",
  "2ef598166f5b9e6bb65d3252cf127b2c992348ae300a06082a8648ce3d040302034900304602",
  "21009cd3d7968f65630677b3c7437cae5187021d553c756820d34b6a8a1240e7efb7022100a3",
  "c37f0a8e1f38a5263b14414e6f916962bbd36bcbaca40b07d371352de745c0")
leafhex <- paste0(
  "308201ae30820154a00302010202142fa99d48ead4c033c568837b2288f9e9bf137068300a06",
  "082a8648ce3d040302302931123010060355040a0c09726d626c207465737431133011060355",
  "04030c0a454332353620526f6f74301e170d3236313030323038343334315a170d3331313030",
  "313038343334315a302931123010060355040a0c09726d626c20746573743113301106035504",
  "030c0a4543323536204c6561663059301306072a8648ce3d020106082a8648ce3d0301070342",
  "00045f0ad14e657dedf3f6cc18530f9d905823b6e9f492ecc0eb1989c9ad8a11a859ebb7681a",
  "e5b1f579974e3c6709e824245608c448da6fe248ebf67bf727d598f4a35a305830090603551d",
  "1304023000300b0603551d0f040403020780301d0603551d0e04160414b1d22b6dd5987f0dfa",
  "0ff3a5d0140351d752e0fb301f0603551d230418301680142ef598166f5b9e6bb65d3252cf12",
  "7b2c992348ae300a06082a8648ce3d0403020348003045022100bc9d6ec452554adb8f33e4ac",
  "53d5e29620fc19528074153028ce849d98c8c15202205c841de357dbad5042daba7fb4af3724",
  "78d684f648e6ef23c68a3d764894c8fc")
hex2raw <- function(h) {
  as.raw(strtoi(substring(h, seq(1, nchar(h), 2), seq(2, nchar(h), 2)), 16L))
}
root <- certprs(hex2raw(roothex))
root
leaf <- certprs(hex2raw(leafhex))
leaf$subject
leaf$key$curve
c(leaf$not_before, leaf$not_after)

# the same certificate from a PEM file
pem <- tempfile(fileext = ".pem")
writeLines(c("-----BEGIN CERTIFICATE-----",
             strsplit(jsb64enc(hex2raw(roothex)), "\n")[[1]],
             "-----END CERTIFICATE-----"), pem)
certprs(pem)$self_signed

# bytes that are not a certificate
try(certprs(charToRaw("not a certificate")))

Regularized incomplete gamma function (C backend)

Description

The lower regularized incomplete gamma function P(a, x), which is the CDF of a Gamma distribution with shape a and unit rate. Exposed because it is the building block of the chi-square tail used by drftchisq() and benfordtst(): the p-value of a chi-square statistic q on k degrees of freedom is 1 - P(k/2, q/2).

Usage

cgamcdf(shape, x)

Arguments

Details

P(a, x) = \frac{1}{\Gamma(a)} \int_0^x t^{a-1} e^{-t} \, dt,

evaluated by the shared core (a series for small x, a continued fraction otherwise, the standard pairing), and agreeing with stats::pgamma(x, shape) to rounding. For a gamma distribution with rate r, pass x * r.

shape must be a single positive number ("shape must be a single positive number"). A negative x gives 0, a missing x gives NA, and Inf gives 1, as in stats::pgamma().

Value

A numeric vector the length of x, each entry a probability in [0, 1].

See Also

stats::pgamma(), the base R equivalent.

Examples

cgamcdf(3.5, c(0.5, 1, 4, 12))

# Identical to the unit-rate gamma CDF.
all.equal(cgamcdf(3.5, c(0.5, 1, 4, 12)),
          stats::pgamma(c(0.5, 1, 4, 12), shape = 3.5))

# Shape 1 is the exponential distribution.
all.equal(cgamcdf(1, c(0.5, 2)), stats::pexp(c(0.5, 2)))

# A chi-square tail on k degrees of freedom is 1 - P(k/2, q/2).
1 - cgamcdf(2 / 2, 5.99 / 2)      # about 0.05 on 2 df

# The shape must be positive.
try(cgamcdf(0, 1))

Envelope of a difference, a percent change or a rate under publication bounds

Description

Given the intervals from pubbnd() for the current and the previous cell, the set of differences and percent changes that the observed counts could have produced is an interval too: its ends are reached at the ends of the input intervals, because the difference and the ratio are monotone in each argument. This is interval arithmetic, not a confidence interval; it says nothing about sampling and everything about what the release withheld.

Usage

chgenv(now, prev, pop = NULL, prepop = NULL, per = 1000)

Arguments

Details

With current interval [a_l, a_u] and previous interval [b_l, b_u]: the change lies in [a_l - b_u, a_u - b_l] and the percent change in [100(a_l / b_u - 1), 100(a_u / b_l - 1)]. An end of the percent interval is NA when its denominator could be zero (an upper bound needs b_l > 0, a lower bound b_u > 0), since a change from zero has no percentage. A rate interval is [per \cdot a_l / pop, per \cdot a_u / pop]; the population is taken as exact.

Errors when now and prev have different numbers of rows. The populations are not checked; give one per row (or one value).

Value

A data frame, one row per cell, with numeric columns chglo, chgup, pctlo, pctup (NA where the previous cell could have been zero); with population also ratelo and rateup; with prepop as well preratelo and prerateup.

See Also

pubbnd(), yoybnd().

Examples

now <- pubbnd(c(45, 120), rnd = 5)
prev <- pubbnd(c(40, 100), rnd = 5)
chgenv(now, prev)

# with rates per 1,000 residents in both years
chgenv(now, prev, pop = c(10000, 20000),
                prepop = c(10000, 19000))

# a suppressed previous cell could have been zero: no upper percentage
chgenv(pubbnd(12), pubbnd("<5"))

# the two sides must line up
try(chgenv(now, prev[1, ]))

Tamper-evident chain of capsule manifests

Description

Keeps a history of runs (one manifest per run) as a hash chain: each entry records the digest of the entry before it, so the chain's integrity covers the ORDER and COMPLETENESS of the history, not merely the contents of each manifest. Deleting an entry, inserting one, or editing one breaks the links from that point onward, and chnvrf() reports the first index where the break occurs. Use it when the record of which analyses were run, and in what order, has to be auditable.

Usage

chnnew()

chnapp(chain, entry, label = NA_character_)

chnhead(chain)

chnseal(chain)

chnvrf(chain)

Arguments

Details

chnnew() starts an empty chain, chnapp() adds an entry, chnhead() returns the last link, chnvrf() checks every link and chnseal() returns one digest to sign.

This is the structure behind an append-only audit log, and it is what a per-manifest digest alone cannot give you: individually valid manifests say nothing about whether any were removed.

Each entry stores payload, the SHA-256 of the entry (of its UTF-8 bytes for a single string, otherwise of serialize(entry, version = 2, xdr = TRUE) with the four bytes that name the writer's R version set to zero), prev, the previous link (64 zeros for the first), scheme (2) and digest = SHA-256 of "v2|<prev>|<payload>|<label>", where <label> is "0" for a missing label and "1" followed by the label otherwise. The entry itself is not kept, only its digest; the label is covered, so editing it breaks the link. Zeroing the writer-version bytes makes the payload the same on every R version, for any R object (a canonical JSON encoding would cover only the objects JSON can hold); the rest of a version-2 stream depends only on the object.

Entries written by 0.5.12 and earlier carry no scheme. Their digest is the SHA-256 of the string "<prev>:<payload>" (the label not covered, and the payload of a non-string entry taken over serialized bytes that include the writer's R version). chnvrf() and chnseal() still accept them under those rules, with a warning that they were sealed with the weaker scheme. Appending to such a chain adds scheme-2 entries.

chnvrf() walks the chain from the start: an entry fails when its prev is not the previous link, its scheme is neither 2 nor absent, or its digest does not recompute (compared in constant time with cdgsteq()).

What the head covers, and what it does not

chnhead() is the STORED digest of the last entry. Deleting an entry from the MIDDLE leaves that value untouched – the stored digests do not change, only the links between them stop agreeing – so the head alone will not notice. chnvrf() will, and names the index.

Conversely, truncating from the END leaves a perfectly valid prefix that chnvrf() accepts, while the head changes.

The two failures are complementary, which is why chnseal() exists: it verifies the links AND folds the entry count and every link digest into one value (SHA-256 of the count and the links joined by "|"), so a single signature over the seal detects an edit, a deletion, a reordering, an insertion and a truncation alike. Sign the seal, not the head.

Signing is what turns tamper-EVIDENT into tamper-PROOF: without a signature an attacker who rewrites the whole chain leaves it internally consistent, because recomputing every link is cheap.

Error: a chain that does not come from chnnew() (for chnapp(), chnhead() and chnvrf()).

Value

chnnew() and chnapp() return an object of class c("bricklayer_chain", "list") with entries (a list of entries, each with label, payload, prev, scheme and digest) and genesis (64 zeros). chnvrf() returns a list of class bricklayer_chain_check with valid (logical), n (integer, entries), broken_at (integer index of the first bad entry, NA when intact) and head. chnhead() returns the head digest (64 hex characters; the genesis value for an empty chain). chnseal() returns a single digest over the verified chain's length and links, or NA if the chain does not verify.

See Also

cpssign() to sign the seal, mrkroot() for pinning the contents of one capsule rather than a history.

Examples

op <- options(rmoriebricklayer.xmssidxfile = tempfile())  # keep this example's
# one-time index record out of your configuration directory; leave the default in use
ch <- chnnew()
ch <- chnapp(ch, "manifest for run 1", label = "run-1")
ch <- chnapp(ch, "manifest for run 2", label = "run-2")
ch <- chnapp(ch, "manifest for run 3", label = "run-3")
ch

# An intact chain verifies.
chnvrf(ch)$valid

# Editing an entry breaks it, and names where.
edited <- ch
edited$entries[[2]]$digest <- csha256("something else")
chnvrf(edited)$valid
chnvrf(edited)$broken_at

# So does deleting one, which a per-manifest digest would not catch.
dropped <- ch
dropped$entries[[2]] <- NULL
chnvrf(dropped)$valid

# Sign the SEAL, which covers the links and the length together.
key <- pqckg(height = 2)
sig <- cpssign(chnseal(ch), key)
cpsvrf(chnseal(ch), sig, signpubkey(key))

# Truncating the chain still seals, but to a different value, so the
# signature no longer verifies.
truncated <- ch
truncated$entries[[3]] <- NULL
cpsvrf(chnseal(truncated), sig, signpubkey(key))

# A chain whose links disagree has no seal to present at all.
chnseal(dropped)

# any object can be an entry, such as a whole manifest list
ch2 <- chnapp(chnnew(), list(run = 1, mean = 2.5), label = "list entry")
chnhead(ch2)

# only chains are accepted
try(chnvrf(list(entries = list())))
options(op)

Split a file into fixed-size chunks

Description

Reads a file as bytes and returns it as consecutive chunks of a fixed size, ready for mrkroot(), mrklvs() and mrkprf(). Use it to pin a large data file as a Merkle root, so that a later change can be traced to the chunk that moved.

Usage

chnkfile(path, chnkbytes = 1048576L)

Arguments

Details

The file is read sequentially with readBin() in pieces of chnkbytes; every chunk has exactly that many bytes except the last, which holds the remainder. All chunks are held in memory at once, so the memory used is about the size of the file. The default 1 MiB chunk is a compromise: smaller chunks localise a change more precisely but make the tree and its proofs larger. An edit that changes the file's length shifts every later chunk boundary, so only edits that keep the length localise to one chunk.

Errors: path that does not name an existing file (or names a directory); chnkbytes that is not a positive whole number.

Value

A list of raw vectors, in file order. A zero-length file gives an empty list. Chunks are returned as bytes rather than strings because a file is bytes: an R string cannot hold a zero byte, so a character chunk could not represent an arbitrary binary file at all.

See Also

mrkroot(), sha512file()

Examples

p <- tempfile()
writeLines(rep("some capsule content", 50), p)

# One small chunk size to show the splitting.
ch <- chnkfile(p, chnkbytes = 128)
length(ch)
lengths(ch)

# The chunks pin the file as a Merkle root.
mrkroot(ch)

# They are the file's bytes, so they concatenate back to it exactly.
identical(unlist(ch), readBin(p, "raw", file.size(p)))

# Overwriting bytes in place changes exactly one leaf.
b <- readBin(p, "raw", file.size(p))
b[300] <- as.raw(0x58)
q <- tempfile()
writeBin(b, q)
which(mrklvs(ch) != mrklvs(chnkfile(q, chnkbytes = 128)))

# a chunk size must be positive
try(chnkfile(p, chnkbytes = 0))
unlink(c(p, q))

Fit a Hawkes process by maximum likelihood (analytic gradient, fast routes)

Description

Fits a self-exciting (Hawkes) point process to a series of event times, such as the times of shootings, break-ins or calls for service, and returns the estimated background rate, the branching ratio (the share of events triggered by earlier events) and the shape of the triggering kernel, with AIC, BIC and a goodness-of-fit test. Use it to ask whether events cluster in time beyond what a steady or seasonal rate explains.

Usage

chwkfit(
  times,
  hor,
  kernel = c("exponential", "weibull", "gamma", "lomax"),
  base = c("constant", "sinusoidal"),
  meth = "auto",
  eps = 1e-09,
  start = NULL
)

Arguments

Details

The model is the univariate Hawkes process

\lambda(t) = \nu(t) + \eta \sum_{t_j < t} g(t - t_j)

with a constant baseline (\nu = e^{a_0}) or a sinusoidal one (\nu(t) = \exp(a_0 + a_1 t/T + a_2 \sin(2\pi t/365.25) + a_3 \cos(2\pi t/365.25)), a trend plus a yearly cycle when times are in days, its integral by the trapezoid rule on max(64, floor(T) + 1) points), branching ratio \eta \in [0.001, 0.99] and one of four normalised kernels g: exponential (decay beta), Weibull (shape alpha, scale lambda), gamma (shape alpha, rate beta) or Lomax (shape alpha, scale c). The log-likelihood is \sum_i \log \lambda(t_i) - \int_0^T \lambda(t)\,dt. Every route uses the analytic gradient of the shared C++ core (the code morie's Python arm calls), so one likelihood evaluation per point instead of one per parameter. The optimiser is a projected BFGS inside a fixed parameter box (a_0 \in [-15, 15], the other baseline terms in [-5, 5], kernel parameters bounded as listed in the source), at most 2000 iterations.

"exact"

Ozaki's O(n) recursion for the exponential kernel; for Weibull and gamma the double sum stops where the kernel underflows to exactly 0 (the same value as the full sum); Lomax the full O(n^2) sum.

"soe"

the completely monotone kernels (Lomax; gamma with shape < 1) as a sum of exponentials (Beylkin & Monzon 2010): relative error eps on every intensity, O(n K); a gamma kernel with shape >= 1 is truncated at eps.

"truncate"

each event excites only lags with kernel tail mass above eps: an approximation for light-tailed kernels, O(n w).

"em"

the EM algorithm (Veen & Schoenberg 2008), finished by the projected BFGS on the exact likelihood from EM's point (EM's steps shrink before it reaches the maximum). It is an R loop of up to 500 iterations and much slower than the other routes (tens of seconds for a few hundred events); because it starts from the same point but moves differently, it can end at a different local maximum.

"inar"

Kirchner's (2017) INAR(p) least-squares estimator on binned counts (constant baseline only): a different, fast, approximate estimator.

"auto"

"exact" for the exponential kernel, "truncate" for Weibull, "soe" for gamma and Lomax.

For the Weibull, gamma and Lomax kernels without start, a second deterministic start at the exponential fit (the shape-1 member) is also run and the better optimum kept, because the likelihood is multimodal. The exponential kernel without start is also run from a second start (branching ratio 0.1 and a decay a quarter of the reciprocal mean gap), and the better optimum kept: on data with little clustering its likelihood can have more than one maximum, one of them with the decay on its bound. A parameter on its bound (see atbound) is still a reason to refit from another start and compare nll.

The reported nll is the exact likelihood at the estimate whatever the route, so AIC is comparable across routes. An event exactly at hor is part of the record and counted, as in morie's Python fit. The Kolmogorov-Smirnov test compares the time-rescaling residuals of chwkres() with the uniform distribution (exact p-value up to 10,000 events without ties, asymptotic otherwise).

Events at the same instant do not excite each other: the intensity sums over strictly earlier events, as the definition of a Hawkes process (a simple point process, with no two events at one time) has it. Tied times are fitted that way, with a warning: data recorded to a resolution, such as daily dates, carry ties the process itself never makes, and treating them as exact biases the fit (Filimonov & Sornette 2015), so spread them uniformly across their interval first with chwkjit().

Errors: times with NA, unsorted, or outside [0, horizon]; hor not one positive finite number; fewer than 50 events; an event rate n / horizon outside [e^{-15}, 100 e^{15}] per unit time (rescale the times); eps not in (0, 1); start of the wrong length or with a non-finite likelihood; method = "soe" with a kernel other than Lomax or gamma; method = "inar" with the sinusoidal baseline. Long fits show a progress spinner (EM a bar) on stderr unless options(morie.quiet = TRUE).

Value

A list:

theta

numeric, all estimates in the order baseline, eta, kernel.

basepar

numeric, a0 (constant) or a0..a3 (sinusoidal), on the log scale.

braratio

numeric, \eta, the expected number of events each event triggers.

kerpar

numeric, the kernel parameters (beta; alpha, lambda; alpha, beta; or alpha, c).

nll, aic, bic

numeric: the exact negative log-likelihood, 2k + 2\,nll and k \log n + 2\,nll.

n, horizon, kpar

the number of events (integer), the window end and the number of parameters (integer).

ksstat, kspva

numeric, the Kolmogorov-Smirnov test of the time-rescaling residuals against the uniform.

kernel, baseline, method

character, the choices used (method after "auto" is resolved).

eps

numeric, eps for "soe" and "truncate", NA otherwise.

converged

logical, the optimiser stopped at a maximum within the parameter box (before its iteration limit, at a finite likelihood). A maximum on the box counts: atbound names the parameters that lie on it.

atbound

character, the names of parameters that lie on the box (character(0) for an interior maximum).

A warning is given when \eta reaches its upper bound 0.99 (a near-critical, probably non-stationary fit). An estimate on the box is not an interior maximum of the likelihood: day-dated (tied) times drive a shape or decay to its wall, and a Lomax fitted to exponential data runs to its exponential limit.

References

Ozaki T (1979). Maximum likelihood estimation of Hawkes' self-exciting point processes. Ann. Inst. Statist. Math. 31, 145–155. DOI 10.1007/BF02480272

Beylkin G, Monzon L (2010). Approximation by exponential sums revisited. Appl. Comput. Harmon. Anal. 28, 131–149. DOI 10.1016/j.acha.2009.08.011

Veen A, Schoenberg FP (2008). Estimation of space-time branching process models in seismology using an EM-type algorithm. JASA 103, 614–624. DOI 10.1198/016214508000000148

Kirchner M (2017). An estimation procedure for the Hawkes process. Quant. Finance 17, 571–595. DOI 10.1080/14697688.2016.1211312

Filimonov V, Sornette D (2015). Apparent criticality and calibration issues in the Hawkes self-excited point process model: application to high-frequency financial data. Quant. Finance 15(8), 1293–1314. DOI 10.1080/14697688.2015.1032544

See Also

chwkjit(), for times recorded to a resolution; chwkres() for the residuals.

Examples

set.seed(1)
times <- sort(stats::runif(400, 0, 100))
fit <- chwkfit(times, 100, "exponential")
fit$braratio
fit$atbound # an interior maximum: the better of two starts
# one start alone can stop with the decay on its bound, at a worse nll
chwkfit(times, 100, "exponential", start = c(log(2.4), 0.4, 4))$atbound
fit$nll
chwkfit(times, 100, "lomax", meth = "soe", eps = 1e-8)$nll
chwkfit(times, 100, "exponential", base = "sinusoidal")$aic

# the fast approximate INAR estimator
chwkfit(times, 100, meth = "inar")$braratio

# daily dates: spread the ties across each day first
days <- floor(times)
jit <- chwkjit(days, hor = 100)
chwkfit(jit, 100)$converged

# too few events, or times past the horizon, are refused
try(chwkfit(times[1:20], 100))
try(chwkfit(times, 50))

Spread event times recorded to a resolution across their interval

Description

Daily (or hourly, or second) timestamps put every event of an interval at one time, but a Hawkes process never has two events at one instant, and fitting the rounded times as exact biases the estimates (Filimonov & Sornette 2015). This function spreads each recorded time uniformly over the interval it stands for, so the times can be passed to chwkfit().

Usage

chwkjit(times, res = 1, seed = 42, hor = Inf)

Arguments

Details

Each time becomes t + w * U, with w = min(resolution, horizon - t) and U uniform on [0, 1), and the result is sorted; this is the remedy Filimonov and Sornette use. The times are sorted (radix order) before the uniforms are assigned, so the result does not depend on the input order. The uniforms are the splitmix64 stream of cunif(), so the result is reproducible and, on daily dates with seed = 42, the same as the jitter rmorie and morie (Python) apply before their TPS Hawkes fits. R's random number generator is not used.

Errors: times with NA or non-finite values; res not one positive finite number; hor not one number; a time past hor.

Value

A numeric vector of the jittered times, the same length as times, sorted increasingly.

References

Filimonov V, Sornette D (2015). Apparent criticality and calibration issues in the Hawkes self-excited point process model: application to high-frequency financial data. Quant. Finance 15(8), 1293–1314. DOI 10.1080/14697688.2015.1032544

See Also

chwkfit(), cunif()

Examples

days <- c(0, 0, 0, 1, 3, 3)
chwkjit(days)
anyDuplicated(chwkjit(days))

# a horizon keeps the last day inside the window
chwkjit(c(0, 1, 2), hor = 2)

# hourly resolution, in days, and another seed
chwkjit(c(0, 0, 0), res = 1 / 24, seed = 7)

# times past the horizon are refused
try(chwkjit(c(0, 5), hor = 3))

Hawkes-process negative log-likelihood (C backend)

Description

The negative log-likelihood of a univariate self-exciting Hawkes process with constant baseline on ⁠[0, horizon]⁠, for the event times times. A Hawkes process is the natural model for arrivals that trigger further arrivals – repeat calls to a service, aftershocks, retweet cascades, revisions to an open-data release.

Usage

chwknll(
  times,
  hor,
  kernel = c("exponential", "weibull", "lomax", "gamma"),
  par
)

Arguments

Details

Four triggering kernels are available. "exponential" is memoryless and evaluates by an O(n) recursion; the other three are not, so they cost O(n^2).

Parameters are passed on the scales the kernel is defined on: par = c(a0, eta, ...) where a0 is the LOG baseline intensity ( \nu = e^{a0}) and eta the branching ratio in (0, 1) – the expected number of children per event, so the process is stationary only for eta < 1. The remaining entries are the kernel's own shape parameters: beta (exponential), ⁠alpha, lambda⁠ (Weibull: shape and scale), ⁠alpha, c⁠ (Lomax: shape and scale, g(u) = \alpha c^\alpha (u + c)^{-(\alpha + 1)}), ⁠alpha, beta⁠ (gamma: shape and rate). This is the likelihood chwkfit() maximises, on the same parameters.

The conditional intensity is \lambda(t) = \nu + \eta \sum_{t_i < t} g(t - t_i) with g a probability density on (0, \infty) (so \eta is the branching ratio), and the value returned is

-\ell = -\sum_i \log \lambda(t_i) + \int_0^T \lambda(t)\,dt,

with T = hor (Hawkes, 1971). For the exponential kernel g(u) = \beta e^{-\beta u}. With no events the value is the compensator alone, \nu T.

The feasible region, outside which the sentinel 1e12 is returned instead of a value: -20 < a_0 < 20; 10^{-6} < \eta < 0.999; and for the kernel parameters, exponential 0.05 \le \beta \le 30; Weibull shape in [0.05, 20] and scale in [0.001, 1000]; Lomax shape in [0.05, 30] and scale in [0.001, 100]; gamma shape in [0.05, 20] and rate in [0.05, 30]. A missing parameter also gives the sentinel.

Refusals (errors, not the sentinel): times with a missing value, not sorted, or outside [0, T]; a hor that is not a single positive finite number; a kernel not among the four; and a par of the wrong length for the kernel. The quadratic kernels poll for a Ctrl-C as they run, so a long evaluation can be interrupted.

Value

A length-1 numeric: the negative log-likelihood, to be MINIMISED. A parameter set outside the feasible region returns the sentinel 1e12 rather than erroring, so the value can be handed straight to stats::optim() without the optimiser walking off the domain.

References

Hawkes AG (1971). Spectra of some self-exciting and mutually exciting point processes. Biometrika 58(1), 83–90. DOI 10.1093/biomet/58.1.83

Examples

set.seed(4)
times <- sort(stats::runif(40, 0, 10))

# Exponential kernel: log-baseline -0.5, branching 0.3, decay 1.2.
chwknll(times, 10, "exponential", c(-0.5, 0.3, 1.2))

# Lower is better, so this is what an optimiser minimises.
nll <- function(p) chwknll(times, 10, "exponential", p)
fit <- stats::optim(c(-0.5, 0.3, 1.2), nll)
fit$par

# A branching ratio at or above 1 is not a stationary process, and is
# reported as the infeasible sentinel rather than a number.
chwknll(times, 10, "exponential", c(-0.5, 1.5, 1.2)) == 1e12

# The heavier-tailed kernels take two shape parameters.
chwknll(times, 10, "gamma", c(-0.5, 0.3, 2, 1.5))

# Unsorted times are an error, not the sentinel.
try(chwknll(rev(times), 10, "exponential", c(-0.5, 0.3, 1.2)))

Time-rescaling residuals of a Hawkes process

Description

Transforms the event times into residuals that are independent and uniform on (0, 1) when the fitted Hawkes model is right, by the time-rescaling theorem. Plot them, or test them against the uniform with stats::ks.test(), to check the fit of chwkfit(); departures show where the model misses clustering or a trend.

Usage

chwkres(
  times,
  hor,
  kernel = c("exponential", "weibull", "gamma", "lomax"),
  par,
  base = c("constant", "sinusoidal")
)

Arguments

Details

With the compensator \Lambda(t) = \int_0^t \lambda(s)\,ds, the residual of event i is

U_i = 1 - \exp(-(\Lambda(t_i) - \Lambda(t_{i-1})))

with t_0 = 0 (Brown et al. 2002). The rescaled gaps \Lambda(t_i) - \Lambda(t_{i-1}) are unit-rate exponential under the model, so the U_i are uniform. The model and kernels are those of chwkfit(); the shared C++ core computes the residuals in O(n) for the exponential kernel and O(n w) for the others. Tied times give tied residuals.

Errors: times with NA, unsorted or outside [0, horizon]; par of the wrong length (2 + kernel parameters for the constant baseline, 5 + kernel parameters for the sinusoidal one).

Value

A numeric vector of length length(times), values in (0, 1).

References

Brown EN, Barbieri R, Ventura V, Kass RE, Frank LM (2002). The time-rescaling theorem and its application to neural spike train data analysis. Neural Comput. 14, 325–346. DOI 10.1162/08997660252741149

Examples

set.seed(2)
times <- sort(stats::runif(200, 0, 50))
U <- chwkres(times, 50, "exponential", c(log(3), 0.2, 2))
stats::ks.test(U, "punif")$p.value

# residuals of a fitted model, and a quick uniform Q-Q check
fit <- chwkfit(times, 50, "weibull")
U2 <- chwkres(times, 50, "weibull", fit$theta)
summary(U2)
plot(stats::ppoints(length(U2)), sort(U2), xlab = "uniform quantile", ylab = "residual")
abline(0, 1)

# par must match the kernel and baseline
try(chwkres(times, 50, "weibull", c(log(3), 0.2, 2)))

Simulate event times from a fitted Hawkes process

Description

Draws a new event history from a chwkfit() result by its cluster representation: immigrant events from the baseline (a Poisson process, thinned from its upper bound for a sinusoidal baseline), then each event's offspring, a Poisson number with mean the branching ratio, at delays drawn from the fitted kernel, generation after generation, until none fall before the horizon.

Usage

chwksim(fit, hor = NULL, seed = 1L)

Arguments

Value

The sorted event times (numeric), with the attribute generation (0 for an immigrant, 1 for its children, ...).

See Also

plothwk(), chwkfit().

Examples

set.seed(1)
times <- sort(stats::runif(300, 0, 100))
fit <- chwkfit(times, 100)
sim <- chwksim(fit, seed = 2)
length(sim)
table(attr(sim, "generation"))
try(chwksim(list(a = 1)))

Trimmed inverse-probability weights (C backend)

Description

The inverse-probability-of-treatment weights 1/e for the treated and 1/(1-e) for the untreated, with the propensity score clamped into ⁠[trimlo, trimhi]⁠ FIRST. Clamping matters: an untrimmed score near 0 or 1 produces a weight large enough for one observation to dominate the entire estimate. These weights reweight a treated and an untreated group so that each resembles the whole sample on the covariates behind the propensity score, the basis of inverse-probability-weighted estimates of a treatment effect (Rosenbaum and Rubin, 1983).

Usage

cipw(treat, ps, trimlo = 0.01, trimhi = 0.99)

Arguments

Details

For unit i with treatment T_i and propensity e_i, the clamped score is \tilde e_i = \min(\max(e_i, l), u) with l = trimlo and u = trimhi, and the weight is 1 / \tilde e_i when T_i = 1 and 1 / (1 - \tilde e_i) when T_i = 0. With the defaults no weight exceeds 100.

Refusals: vectors of different lengths, a treatment coded other than 0/1 ("treat must be coded 0/1", naming the offending values), a propensity outside [0, 1], and bounds that do not satisfy 0 < l < u < 1. A missing treatment or propensity gives a missing weight for that unit only. The weights are not normalised; divide by their sum within each arm for stabilised (Hajek) estimates.

Value

A numeric vector of weights the length of treat, each at least 1 (and at most 1 / trimlo or 1 / (1 - trimhi)), NA where the treatment or propensity is missing.

References

Rosenbaum PR, Rubin DB (1983). The central role of the propensity score in observational studies for causal effects. Biometrika 70(1), 41–55. DOI 10.1093/biomet/70.1.41

See Also

cwtd() to average an outcome with these weights.

Examples

treat <- c(1, 0, 1, 0)
e <- c(0.5, 0.25, 0.02, 0.9)

cipw(treat, e)

# A balanced score gives weight 2 to either arm.
cipw(c(1, 0), c(0.5, 0.5))

# Without trimming the third observation would carry weight 50; the
# default clamp holds it to 100 at the 0.01 floor, and a looser floor
# tames it further.
cipw(treat, e, trimlo = 0.10)

# Logical treatment indicators work too.
cipw(c(TRUE, FALSE), c(0.4, 0.4))

# A treatment coded 1/2 is refused rather than read as all-control.
try(cipw(c(1, 2), c(0.5, 0.5)))

Generate a Data Citation From Provenance

Description

Builds a ready-to-paste data citation, as plain text and as a BibTeX @misc entry, from the dataset and resource blocks of a provenance record. Use it so that every report and paper built from a capsule cites the data set in the same words, with the retrieval date and licence that were actually pinned.

Usage

citecps(prov)

Arguments

Details

The fields used, with their fallbacks:

• author: dataset$publisher, else "Unknown publisher";

• title: resource$name, else dataset$package_slug, else "Untitled dataset";

• year: the first four characters of captured_at_utc, else the current year; retrieval date: its first ten characters (omitted if absent);

• link: https://doi.org/<dataset$doi> if a DOI is recorded, else resource$direct_url, else dataset$catalogue_page;

• licence: dataset$licence_short, else dataset$licence_name.

The text form is Publisher (year). Title [Data set]. link Retrieved date. Licence: x. The BibTeX key is the package slug with non-alphanumeric characters removed, followed by the year. In the BibTeX entry the licence is written in the note field only when a retrieval date is present. No character in the fields is escaped for LaTeX. Error: prov that is not a list.

Value

A list with text and bibtex, each a length-1 character string (the BibTeX entry has embedded newlines), or NULL if provenance is NULL.

Examples

prov <- list(
  captured_at_utc = "2026-06-23T04:41:40Z",
  dataset = list(publisher = "Ontario Ministry of the Solicitor General",
                 licence_short = "OGL-Ontario",
                 package_slug = "data-on-inmates-in-ontario"),
  resource = list(name = "Restrictive Confinement - Detailed Dataset",
                  direct_url = "https://data.ontario.ca/example.csv")
)
cit <- citecps(prov)

# Plain-text citation ready to paste.
cat(cit$text)

# BibTeX @misc entry for LaTeX bibliographies.
cat(cit$bibtex)

# a recorded DOI replaces the URL
prov$dataset$doi <- "10.5683/SP3/EXAMPLE"
cat(citecps(prov)$text)

# NULL provenance returns NULL (composes safely).
citecps(NULL)

# anything else that is not a list is an error
try(citecps("data-on-inmates-in-ontario"))

Tidy the column names of an ingested table

Description

Converts column names to a consistent, syntactically valid form: transliterated to ASCII, non-alphanumerics collapsed to a single separator, and duplicates disambiguated with a numeric suffix. The counterpart of janitor::clean_names(). Use it right after reading an open-data file, so the rest of the script can refer to columns without backticks.

Usage

clncolnms(
  data,
  case = c("snake", "lower_camel", "upper_camel", "screaming_snake", "none"),
  sep = "_"
)

Arguments

Details

Open-data extracts arrive with names like "Total Population (2021)" and "% change", which need backticks everywhere and break silently when a re-release renames "% change" to "% Change" or adds a trailing space. Normalising once at ingestion makes the schema stable against that.

The steps, in order: transliterate to ASCII (iconv(..., "ASCII//TRANSLIT"), keeping the original where that fails); replace % by pct, & by and and # by n as words; split camelCase at a lower-to-upper boundary; turn every run of other characters into one word break; an empty name becomes x; join the words in the chosen case; prefix x to a name starting with a digit; and give repeated names the suffixes <sep>2, <sep>3, ... in order of appearance (the first keeps its name). The transliteration depends on the platform's iconv, so an accented letter may come out slightly differently on different systems.

Record the mapping in the capsule manifest. Cleaning names changes what a downstream script must refer to, so an unrecorded cleaning is itself a reproducibility hazard – the returned object carries the original names in its "orignms" attribute for exactly that.

Value

For a data frame, the data frame with new names and an "orignms" attribute (character, the names before cleaning); for a character vector, the cleaned character vector.

Examples

clncolnms(c("Total  Population (2021)", "% change",
                     "Ville / City", "dup", "dup"))

# Applied to a data frame, with the original names retained.
df <- data.frame(`Total Pop` = 1:2, `% change` = 3:4,
                 check.names = FALSE)
cleaned <- clncolnms(df)
names(cleaned)
attr(cleaned, "orignms")

# Other cases.
clncolnms(c("Total Pop"), case = "lower_camel")
clncolnms(c("Total Pop"), case = "upper_camel")
clncolnms(c("Total Pop"), case = "screaming_snake")

# a different separator, and a name that starts with a digit
clncolnms(c("2021 count", "Rate per 1000"), sep = ".")

# an unknown case is an error
try(clncolnms("a b", case = "kebab"))

Decompose a change in a clearance rate into rate and case-mix effects

Description

An aggregate clearance rate – cases disposed per case received by a court, or offences cleared per offence reported to a police service – is a weighted mean of the rates of its parts (offence groups, courts, case types), weighted by each part's share of the denominator. It can therefore fall while the rate of every part rises, if the work shifts toward the parts that clear more slowly: Simpson's paradox. This function compares two periods (or two jurisdictions), splits the change in the aggregate exactly into the change in the parts' rates and the change in the mix, reports the rate standardised to a common mix, tests the within-part change, and says whether the comparison is a reversal.

Usage

clrdcmp(num, den, grp, cmp, base = NULL, conflvl = 0.95)

Arguments

Details

With part rates r_g and denominator shares w_g, the aggregate is R = \sum_g w_g r_g. Kitagawa's symmetric decomposition writes the change between comparisons A and B as

R_B - R_A = \sum_g \bar w_g (r_{gB} - r_{gA}) + \sum_g \bar r_g (w_{gB} - w_{gA}),

with \bar w_g and \bar r_g the averages over A and B: the first sum is the rate effect (the change the aggregate would show with the mix held at its average), the second the mix effect. The two add to the observed change with no residual. Because the shares' changes add to zero, each part's mix term is reported as (\bar r_g - \bar R)(w_{gB} - w_{gA}), centred on the average aggregate rate \bar R: the total is the same, and a part's term is negative when work moved toward a part that clears below average.

The within-part change is estimated by the Mantel-Haenszel rate ratio, treating each part's numerator as a Poisson count over its denominator, with the Greenland-Robins variance; the p-value is for the hypothesis that the ratio is 1. Each part's own rate ratio has the exact conditional (binomial) interval. A court's clearance rate can exceed 1, so it is treated as a rate rather than a proportion; for police clearance, where it is a proportion, the Poisson treatment is conservative.

The decomposition is only as good as the parts: a different grouping (courts instead of offence groups) can tell a different story, and both can be right. Demers and Rossmo found 210 reversals in fifteen years of clearance rates for Canada's fifty largest police services, four of them reversed in every crime category, and advise comparing clearance rates only within comparable mixes.

Value

A list of class rmbl_clearance_decomposition:

groups

one row per part: the counts, rates and denominator shares (weight) in both comparisons, the rate ratio with its exact interval, and the part's rateeff and mixeff.

overall

the aggregate rates, their change, its rateeff and mixeff (which add to it exactly), the comparison rate standardised to the baseline mix (staonbase) and to the average mix (staonave).

rateratio

the crude rate ratio and the Mantel-Haenszel rate ratio over the parts, each with its interval, and the p-value for the stratified ratio.

simpson

"full" when every part moved one way and the aggregate the other; "mix" when the aggregate moved against the mix-standardised change (the mix outweighed the rates); "none" otherwise; with the number of parts that rose and fell.

References

Simpson, E. H. (1951). The interpretation of interaction in contingency tables. Journal of the Royal Statistical Society, Series B, 13(2), 238-241.

Kitagawa, E. M. (1955). Components of a difference between two rates. Journal of the American Statistical Association, 50(272), 1168-1194. DOI 10.2307/2281213

Das Gupta, P. (1993). Standardization and Decomposition of Rates: A User's Manual. U.S. Bureau of the Census, Current Population Reports P23-186. Washington, DC.

Greenland, S. and Robins, J. M. (1985). Estimation of a common effect parameter from sparse follow-up data. Biometrics, 41(1), 55-68.

Demers, S. and Rossmo, D. K. (2015). Simpson's paradox in Canadian police clearance rates. Canadian Journal of Criminology and Criminal Justice, 57(3), 424-434. DOI 10.3138/cjccj.2014-E29

See Also

caseflow for the clearance rate over time.

Examples

g <- read.csv(system.file("extdata", "ocj_offence_groups.csv",
                          package = "rmoriebricklayer"))
g <- g[g$year %in% c(2014, 2022), ]
cd <- clrdcmp(g$disposed, g$received, g$offence_group, g$year)
cd
plot(cd)

# Simpson's paradox in four numbers: each part clears better in B, yet the
# aggregate falls, because B's workload moved to the part that clears less
toy <- clrdcmp(num = c(90, 20, 19, 60),
                               den = c(100, 50, 20, 140),
                               grp = c("minor", "serious", "minor", "serious"),
                               cmp = c("A", "A", "B", "B"))
toy

# every part must appear in both comparisons
try(clrdcmp(c(5, 6, 7), c(10, 10, 10), c("x", "y", "x"), c(1, 1, 2)))

Cumulative risk of a first event by age

Description

The chance that a person experiences an event (a first arrest, a first conviction, a first jail stay, a first child protection placement) by each age, from first events and the population by age, as a synthetic cohort; with mortality as a competing risk when deaths are given.

Usage

cmlrisk(age, fireve, pop, atrisk = NULL, dthrate = NULL, conflvl = 0.95)

Arguments

Details

With E_a first events at age a and N_a persons at risk (those who have not had the event), the hazard is h_a = E_a/N_a. When only the population P_a is known, the persons at risk are taken as P_a(1 - CR_{a-1}) (the synthetic-cohort assumption). Without mortality, CR_a = 1 - \prod_{x \le a}(1 - h_x). With death rates q_x, the multiple-decrement table carries survivors l_{x+1} = l_x(1 - h_x - q_x) and CR_a = \sum_{x \le a} l_x h_x / l_0, the risk of the event before death. The interval (without mortality) is Greenwood's on the log of the survival: \mathrm{Var}(\log S_a) = \sum_{x \le a} h_x/(N_x(1 - h_x)).

Counting first events needs records that link a person's events over time; an agency's counts of all events overstate first events. Events of people who do not live in the area, counted against its population, bias the risk upward.

Value

A data frame: age, hazard, atrisk, cmlrisk, lower, upper.

References

Bonczar TP (2003). Prevalence of imprisonment in the U.S. population, 1974-2001. Bureau of Justice Statistics, NCJ 197976.

Applied to arrests, convictions and jail in New York City: Hepburn P, Kohler-Hausmann I, Zorro Medina A (2019). Cumulative risks of multiple criminal justice outcomes in New York City. Demography 56(3):1161-1171. DOI 10.1007/s13524-019-00781-7

Examples

age <- 16:25
first <- c(900, 1400, 1700, 1600, 1400, 1100, 900, 700, 600, 500)
pop <- rep(100000, 10)
cmlrisk(age, first, pop)

Mean, variance, skewness and kurtosis (C backend)

Description

The four standard descriptive moments of a numeric vector – where it is centred, how spread out it is, whether it leans to one side and how heavy its tails are – in one call to the compiled core. These are the four numbers a data description or a drift check usually reports together.

Usage

cmom(x)

Arguments

Details

The mean is base R's (an extended-precision sum followed by a corrective pass); the deviations from it are then divided by their largest magnitude before any power is taken, and the scale is restored at the end. That avoids both the cancellation of the textbook “sum of squares minus square of sum” formula and overflow of fourth powers for values beyond about 10^{77}. It reads the data three times (mean, largest deviation, moments), in O(n) time. The variance uses the n - 1 denominator, matching stats::var(). The shape statistics use the sample moment definitions

g_1 = m_3 / m_2^{3/2}, \qquad g_2 = m_4 / m_2^2 - 3,

with m_k = n^{-1} \sum_i (x_i - \bar x)^k, so kurtosis is reported as EXCESS kurtosis: a normal sample sits near zero, not near three. These are the uncorrected moment ratios (type 1 in the classification of Joanes and Gill, 1998), not the bias-adjusted versions some software reports.

Edge cases are NaN, never an error: an empty vector gives four NaN; one value gives its mean and NaN for the rest; skewness needs at least 3 values and kurtosis at least 4; any NA or NaN in x makes every entry NaN. A constant vector has variance 0 and NaN shape statistics. A character, factor or list x is refused ("must be numeric").

Value

A named length-4 numeric vector: mean (units of x), variance (squared units of x, n - 1 denominator), skewness and kurtosis (both unitless; kurtosis as excess over the normal's 3).

References

Joanes DN, Gill CA (1998). Comparing measures of sample skewness and kurtosis. Journal of the Royal Statistical Society: Series D (The Statistician) 47(1), 183–189. DOI 10.1111/1467-9884.00122

See Also

csdev(), cqnt() and cmad() for robust alternatives that one extreme value cannot move.

Examples

cmom(c(2, 4, 4, 4, 5, 5, 7, 9))

# The first two entries agree with base R.
m <- cmom(1:10)
all.equal(m[["mean"]], mean(1:10))
all.equal(m[["variance"]], stats::var(1:10))

# A symmetric sample has no skew; excess kurtosis is near 0 for normal
# data and positive for a heavy-tailed sample.
cmom(c(-2, -1, 0, 1, 2))[["skewness"]]
cmom(c(rep(0, 20), -8, 8))[["kurtosis"]] > 0

# Too short to define a shape statistic: NaN rather than a guess.
cmom(c(1, 2))

Normal log-density (C backend)

Description

The logarithm of the normal (Gaussian) density at each value of x, computed directly from the formula rather than as log(dnorm(x)). Use it for log-likelihoods: summing log-densities is how a normal model is scored, and computing them directly keeps the sum finite far into the tails, where the density itself underflows to zero and its logarithm would be -Inf.

Usage

cnormlpdf(x, mean = 0, sd = 1)

Arguments

Details

Each value is

\log \phi(x; \mu, \sigma) = -\log \sigma - \tfrac12 \log(2\pi) - \frac{(x - \mu)^2}{2\sigma^2},

evaluated by the shared core, and agrees with stats::dnorm(x, mean, sd, log = TRUE) to rounding.

mean and sd are single numbers: a vector of either is refused, so this is not vectorised over the parameters the way stats::dnorm() is. A non-finite mean, and an sd that is not a positive finite number, are errors, as in cnormpdf(); a missing x gives NA. A character, factor or list x is refused ("must be numeric").

Value

A numeric vector the length of x: the natural-log density at each value (in log units per unit of x). Equivalent to stats::dnorm(x, mean, sd, log = TRUE).

See Also

cmom() for estimating mean and sd from data.

Examples

cnormlpdf(c(-1, 0, 1))

# Identical to stats::dnorm(log = TRUE).
all.equal(cnormlpdf(-2:2, 0.3, 1.7),
          stats::dnorm(-2:2, 0.3, 1.7, log = TRUE))

# Still finite where the density itself underflows to zero.
stats::dnorm(50)                  # 0
cnormlpdf(50)            # about -1251

# A log-likelihood: the sum of the log-densities of a sample.
y <- c(4.8, 5.1, 5.6, 4.9, 5.3)
sum(cnormlpdf(y, mean = mean(y), sd = stats::sd(y)))

# A standard deviation of zero has no density, and the parameters are
# single numbers.
try(cnormlpdf(1, sd = 0))
try(cnormlpdf(1, mean = c(0, 1)))
try(cnormlpdf("1"))

Normal density (C backend)

Description

Evaluates the density of the normal distribution with a given mean and standard deviation at each value of x, in compiled code. It is the same kernel that other packages reach through LinkingTo: rmoriebricklayer, and for a positive standard deviation it agrees with dnorm.

Usage

cnormpdf(x, mean = 0, sd = 1)

Arguments

Details

f(x) = \frac{1}{\sigma\sqrt{2\pi}} \exp\left(-\frac{(x - \mu)^2}{2\sigma^2}\right).

The function is vectorised over x only; mean must be a single finite number and sd a single positive number, and anything else is an error (dnorm() recycles them, and gives Inf or NaN for sd = 0 or below). NA in x gives NA; x = Inf gives 0. Character, factor or list input for x is an error.

Value

A numeric vector the length of x.

See Also

cstat, dnorm.

Examples

# Standard normal density at a few quantiles.
cnormpdf(c(-1, 0, 1))

# Peak of the standard normal is 1/sqrt(2*pi) at x = 0.
cnormpdf(0)

# Shift and scale via mean and sd.
cnormpdf(c(0, 5, 10), mean = 5, sd = 2)

# Identical to stats::dnorm().
all.equal(cnormpdf(-2:2, 0, 1), stats::dnorm(-2:2, 0, 1))

# a non-positive sd and text are refused
try(cnormpdf(0, sd = 0))
try(cnormpdf("0"))

A Community Well-Being score from the census profile

Description

Indigenous Services Canada's Community Well-Being (CWB) index (cwbscr) computed from the public 2021 Census Indigenous Population Profile, for each geography and group.

Usage

cnscwb(
  x,
  dwell = NULL,
  groups = c("First Nations", "Métis", "Inuit", "Non-Indigenous"),
  bounds = c(2650, 75000),
  resid = "Total - Residence on or off reserve",
  minpop = 250
)

Arguments

Details

The Department computes the CWB from census microdata; the public profile gives close but not the same inputs, and this function says where it departs:

income

per capita: the number of total income recipients aged 15 and over times their average total income, over the population in private households (all ages), as the CWB does.

education

the shares with at least a high school diploma (all but those with no certificate, diploma or degree) and with a bachelor's degree or higher, both among those aged 25 and over: the profile has no age group of 20 and over, which the CWB uses for the first.

labour force

the participation rate and the employment rate of participants (100 less the unemployment rate) of those aged 25 to 64: the CWB uses 20 to 64.

housing

the shares of households with one person or fewer per room and needing no major repairs, from the household profile (dwell), which splits Indigenous from non-Indigenous households only, so First Nations, Metis and Inuit people all get the pooled Indigenous households' housing: the CWB counts persons, not households. Without dwell the housing component and the score are NA.

A group of fewer than minpop people in a geography gets no score. The Department scores communities of 65 people or more; the default here is 250, because the public profile's counts are rounded to multiples of 5 and its long-form counts rest on a few dozen households at that size. Geographies are matched by DGUID, so places that share a name stay apart. The inputs are those of the 2021 profile (98-510-X2021001 to -016): other editions name their characteristics differently. Read the result as a CWB-style description of socio-economic conditions, close to but not the Department's index, and never as a measure of the wellness of First Nations, Metis or Inuit people as they define it (see cwbscr).

Value

A data frame of class rmbl_census_cwb: geo, dguid, group; the components income, education, housing, labour and the score cwb (0 to 100, NA where not scored); and the inputs incpercap (dollars), highsch25p, uni25p, population, participation, employment, nomajrep, notcro (percentages).

References

Indigenous Services Canada (2019). Report on trends in First Nations communities, 1981 to 2016 (Community Well-Being index).
sac-isc.gc.ca/eng/1345816651029/1557323327644

See Also

cwbscr.

Examples

p <- cnsprof(indiex("censuspop"))
h <- cnsprof(indiex("censushh"))
cnscwb(p, h)[, c("geo", "group", "income", "education", "housing", "labour",
                     "cwb")]

Read the 2021 Census Indigenous Population Profile

Description

Reads a table of the Indigenous Population Profile, 2021 Census of Population (Statistics Canada catalogue 98-510-X2021001 to 98-510-X2021016), as the data service mirrors it, into one row per geography, gender, age group, characteristic and group: First Nations, Metis, Inuit, the other identity groups and non-Indigenous people (or, in the household tables, Indigenous and non-Indigenous households), each in its own rows. The data service keeps these tables under census2021_ipp/.

Usage

cnsprof(x)

Arguments

Details

The profile gives each characteristic in a hierarchy shown only by indentation (“Total - Highest certificate, diploma or degree”, “High (secondary) school diploma or equivalency certificate”, ...): the level, the parent and the topic's total (rootid) are worked out from it, so a share can be taken of its parent. A topic marked “25% sample data” comes from the long-form sample (sample = "25%"); the others from the full count ("100%"). Counts are randomly rounded to multiples of 5 (of 10 below 10) by Statistics Canada, and a cell suppressed to protect confidentiality (x), too unreliable (F), not available (..) or not applicable (...) has no value here; E values are kept and marked. The global non-response rates of the geography (tnrsf, tnrlf, in percent) are kept: Statistics Canada suppresses an area whose long-form rate is 50% or more, and a high rate means a larger error than the sample alone.

Value

A data frame of class rmbl_census_profile: geo, dguid, geolvl, residence (on or off reserve, in the tables that split it), gender, age, chaid, characteristic, level, parid, rootid, sample, group, value, quality, tnrsf, tnrlf.

References

Statistics Canada (2023). Indigenous Population Profile, 2021 Census of Population. Statistics Canada Catalogue no. 98-510-X2021001. Released June 21, 2023. Published at www12.statcan.gc.ca/census-recensement/2021/dp-pd/ipp-ppa/index.cfm?Lang=E (the site answers scripts with 403); the Wayback Machine holds a copy of that address captured on 8 March 2026 (web.archive.org, snapshot 20260308015010).

Statistics Canada (2022). Guide to the Census of Population, 2021. Catalogue no. 98-304-X.

See Also

profshr, cnsstd, cnscwb, indiex.

Examples

p <- cnsprof(indiex("censuspop"))
head(p[p$characteristic == "Foster children" & p$geo == "Canada", ])

An Indigenous standard population from the census

Description

The age structure of a group in a geography of the 2021 Census Indigenous Population Profile, as the standard of agestd or stdprev.

Usage

cnsstd(
  x,
  grp = "Indigenous identity",
  geo = "Canada",
  ages = c("0 to 4 years", "5 to 9 years", "10 to 14 years", "15 to 19 years",
    "20 to 24 years", "25 to 34 years", "35 to 44 years", "45 to 54 years",
    "55 to 64 years", "65 to 74 years", "75 years and over"),
  resid = "Total - Residence on or off reserve"
)

Arguments

Details

Comparisons of First Nations, Metis and Inuit health are standardised to the Indigenous population's own, younger, age structure (Kumar and Tjepkema 2019): a standard as old as Canada's gives weight to ages where few Indigenous people are. This takes the population in private households of each age group (the profile's first characteristic) for the group and geography given. Age groups that do not overlap are chosen by default (0 to 4, 5 to 9, 10 to 14, 15 to 19, 20 to 24, 25 to 34, 35 to 44, 45 to 54, 55 to 64, 65 to 74, 75 and over); collapse them to the age groups of the rates being standardised.

Value

A data frame: age, population, weight (shares summing to 1).

References

Kumar MB, Tjepkema M (2019). Suicide among First Nations people, Metis and Inuit (2011-2016): findings from the 2011 Canadian Census Health and Environment Cohort (CanCHEC). Statistics Canada catalogue no. 99-011-X2019001.

See Also

agestd, stdprev.

Examples

p <- cnsprof(indiex("censuspop"))
cnsstd(p)

Trend in a count series, as a rate ratio per period

Description

Estimates the trend in a series of counts (arrests, admissions, incidents per year) as a rate ratio per period: 0.85 means the count falls by 15% a year. It fits a Poisson log-linear model and, when the counts vary more than a Poisson model allows, widens the interval. An offset carries the denominator when the exposure varies, so the trend is then in the rate rather than in the count.

Usage

cnttrnd(y, x = NULL, offset = NULL, conflvl = 0.95)

Arguments

Details

The model is \log \mu_t = \beta_0 + \beta_1 (x_t - \bar x) + \log o_t, with o_t the offset (1 when none), fitted by iteratively reweighted least squares (at most 50 iterations). The rate ratio is e^{\beta_1}, with a Wald interval e^{\beta_1 \pm z\, \mathrm{se}}.

The dispersion is reported because a Poisson fit assumes it is one. It is the Pearson statistic \sum (y - \mu)^2/\mu divided by n - 2. When it exceeds 1.5 the quasi-Poisson interval is returned instead: the standard error is multiplied by the square root of the dispersion and the t distribution on n - 2 degrees of freedom replaces the normal, and overdispersed is set.

Missing values are dropped. Errors: negative or fractional counts; x and y of different lengths; infinite x; every count zero; an offset of the wrong length or not positive; fewer than three periods; a conflvl outside (0, 1). When the counts are zero at one end of the series only (c(0, 0, 5)) the likelihood rises without bound and there is no finite trend: a warning says so, the rate ratio is Inf (or 0) and the interval and p-value are NA. There is a plot() method for the result.

Value

A list of class c("rmbl_count_trend", "list") with rateratio (per period), lower, upper, p_value (test of no trend), logslope (\beta_1), se (its standard error, scaled when overdispersed), fitted (fitted counts), dispersion, overdispersed (logical), n, data (list of the x and y used) and method.

References

Bilder, C. R. and Loughin, T. M. Analysis of Categorical Data with R, 2nd edn. Chapman and Hall/CRC, on the quasi-likelihood treatment of an overdispersed Poisson fit: the variance is scaled by an estimated dispersion, which widens the interval while leaving the point estimate alone. That is the behaviour reported here through dispersion and overdispersed.

See Also

trndtst, ratechg.

Examples

# A count falling by about 15% a year.
set.seed(3)
y <- stats::rpois(8, lambda = 200 * 0.85^(0:7))
fit <- cnttrnd(y)
round(fit$rateratio, 3)
c(fit$lower, fit$upper)

# With a varying denominator the trend is in the rate.
cnttrnd(c(20, 25, 30), offset = c(1000, 1500, 2500))$rateratio

# Counts that swing far more than Poisson noise get the wider interval.
wild <- cnttrnd(c(10, 50, 5, 60, 8, 70))
wild$dispersion
wild$method

try(cnttrnd(c(1.5, 2, 3)))

Energy and carbon footprint of a computation

Description

Evaluates expr and reports the electricity it used and the CO2-equivalent that electricity implies, by one of two methods, for the grid it ran on, offline. Use it to report the energy and emissions of an analysis next to its results (as journals and funders increasingly ask), or to compare two ways of running the same computation. The result keeps the value of expr, and lists every figure that was assumed rather than measured.

Usage

co2fp(
  expr,
  loc = NULL,
  meth = c("auto", "rapl", "model"),
  cpupwrw = NULL,
  usage = "process",
  loadcur = c("linear", "codecarbon"),
  cores = NULL,
  memgb = NULL,
  mempwrwpergb = 0.3725,
  pue = 1,
  rapldir = "/sys/class/powercap",
  procstat = "/proc/stat"
)

Arguments

Details

Measured (method = "rapl"): the CPU package energy counters (RAPL, microjoules, read from rapldir) before and after, with wrap-around handled through max_energy_range_uj. This is what codecarbon does on Linux; the counters cover the CPU packages only, so memory is still modelled. Available when raplok() is TRUE.

Modelled (method = "model"): the Green Algorithms formula (Lannelongue, Grealey and Inouye 2021),

E = t \,(P_{cpu}\, u + M\, P_{mem})\, PUE

with t the running time, P_{cpu} the CPU's thermal design power in watts, u its utilisation, M the memory in GB, P_{mem} = 0.3725 W/GB and PUE the data-centre overhead. When the CPU's power is unknown (cpupwrw = NULL) the fallback is codecarbon's: 50% of a constant 85 W, that is 42.5 W charged for the whole running time whatever the measured utilisation (the utilisation is still reported, but not applied).

Utilisation is measured, not assumed, in one of two ways. usage = "process" (the default) is this R process's own CPU time, children included, divided by the wall time and the core count: codecarbon's process tracking mode, which charges only what the computation itself used. usage = "machine" is the whole machine's busy share over the run, from /proc/stat on Linux: codecarbon's machine tracking mode, right when the computation is the only thing the machine is doing or when the question is what the machine drew. Where /proc/stat is absent the process figure is used and the assumptions say so. A number between 0 and 1 fixes the utilisation, as the Green Algorithms calculator asks for it. loadcur chooses how power scales with that utilisation: linearly (Green Algorithms) or along codecarbon's curves, a floor at 10% of P_{cpu} for a process and 0.1 + 0.9 u^3 for a machine, so a figure can be made comparable with either tool.

Cross-checked against codecarbon 3.3.1 on one workload with the same forced inputs: memory energy identical (5.96 W for 16 GB), CPU energy within 8%. loadcur = "codecarbon" reproduces its power-from-load curves, so what remains of the difference is the load measurement itself (codecarbon samples the load during the run; this function integrates the CPU time).

Emissions are E \times the carbon intensity of the location (co2int: a bundled table of every country at its latest published year and of sub-national zones, so no network access is needed anywhere). With location = NULL the location is detected from the environment, the time zone or the locale (dtctloc) and the detection is recorded. Every default that stood in for a measurement is listed in assumptions, so a report can say what was measured and what was assumed.

The energy figures include the PUE factor (CPU, memory and total alike). The running time has a floor of the CPU time divided by the cores, so a computation faster than the clock's resolution does not come out at zero. Short computations (well under a second) give figures dominated by timer resolution; measure a representative workload.

Every argument is checked, and the carbon intensity looked up, before expr is evaluated, so a bad argument never costs the computation. Errors: an unknown loc; pue below 1 or not a single finite number; cpupwrw or cores not a single positive number (cores also whole); memgb or mempwrwpergb negative; usage other than "process", "machine" or one number in [0, 1]; method = "rapl" where the counters are not readable. An error inside expr propagates unchanged. Nothing is written and nothing is fetched.

Value

A list of class rmbl_footprint with elements value (the result of expr), durs (numeric, wall seconds), cputms (numeric, CPU seconds of this process and its children), usage (numeric in [0, 1], the utilisation), usagemode ("process", "machine" or "fixed"), cores (integer), memgb (numeric), enekwh (a list of numerics cpu, memory and total, in kWh, PUE included), co2eg (numeric, grams of CO2-equivalent), co2int (the list returned by co2int, with its year and source), method ("rapl" or "model", the method actually used), assumptions (character vector of every default or fallback that stood in for a measurement) and citations (character vector). The print method shows a three-line summary.

References

Lannelongue, L., Grealey, J. and Inouye, M. (2021). Green Algorithms: quantifying the carbon footprint of computation. Advanced Science 8(12), 2100707. Courty, B. et al. (2024). CodeCarbon: estimate and track carbon emissions from machine learning computing. DOI 10.5281/zenodo.4658424. Lacoste, A., Luccioni, A., Schmidt, V. and Dandres, T. (2019). Quantifying the carbon emissions of machine learning. arXiv:1910.09700.

Examples

fp <- co2fp(sum(sqrt(seq_len(2e5))), loc = "CA-ON")
fp$enekwh$total
fp$co2eg
fp$co2int$year
fp$assumptions
print(fp)
# wherever you are, offline: the location is detected and recorded
co2fp(1 + 1, cpupwrw = 45, memgb = 8)$assumptions
# the whole machine's utilisation, on codecarbon's machine curve
co2fp(1 + 1, loc = "FR", usage = "machine", loadcur = "codecarbon",
                  cpupwrw = 45, memgb = 8)$usagemode
# a fixed utilisation and a data-centre PUE, with a carbon intensity given as a number
fp2 <- co2fp(sum(sqrt(seq_len(2e5))), loc = 30, usage = 0.5,
                         cpupwrw = 65, memgb = 16, pue = 1.67, meth = "model")
fp2$enekwh
fpeqv(fp2$co2eg)
# a bad argument is refused before the expression runs
try(co2fp(Sys.sleep(10), pue = 0.9))

Carbon intensity of electricity by location, offline, anywhere

Description

The grams of CO2-equivalent emitted per kilowatt-hour of electricity at a location, with the year and source of the figure. The table is bundled with the package, so it works offline in any country: every country at the latest year published by Our World in Data from Ember's yearly electricity data (the latest year each grid reports: 2025 or 2024 for nearly all, older for a few grids, Ukraine (2022) among them, named in year), sub-national zones from Electricity Maps' 2024 yearly data as redistributed in the Green Algorithms data release v3.1 (Canadian provinces and territories, US balancing authorities, Australian states, Indian and Japanese regions, Brazilian and Chilean systems, ...), and the world average. All figures are lifecycle intensities (generation plus upstream), the basis both sources publish.

Usage

co2int(loc = "WORLD")

co2inttab()

Arguments

Details

A location is matched by its code ("CA", "CA-ON", "US-CAL-CISO"), its ISO 3166-1 alpha-3 code ("CAN") or its name ("Canada"), case-insensitively. An unknown sub-national code falls back to its country and says so in note. NULL detects the location with dtctloc.

Value

co2int(): a list with g_per_kwh, location, name, year, basis, source and note (a character vector, empty unless a fallback or detection was involved). co2inttab(): the full table as a data frame with the same columns plus iso3.

Provenance and licence of the table

The country rows are Ember's yearly electricity data as published by Our World in Data (CC BY 4.0); the sub-national zones and the Canadian provinces come from the Electricity Maps yearly files redistributed in the Green Algorithms data repository (v3.1), under that repository's terms, and the default PUE from the same repository. Every row names its source and year. Values are reproduced unchanged: Ember reports three small diesel grids (Montserrat, Saint Helena, Uzbekistan) at exactly 1000 gCO2e/kWh and Turkmenistan's gas fleet at 1306, and the table does not round or cap them. Kosovo appears once, under the code XK.

References

Ember (2026). Yearly Electricity Data, via Our World in Data, ourworldindata.org/grapher/carbon-intensity-electricity. Electricity Maps (2025). 2024 Yearly Carbon Intensity Data, via Green Algorithms data v3.1. Lannelongue, L., Grealey, J. and Inouye, M. (2021). Green Algorithms: quantifying the carbon footprint of computation. Advanced Science 8(12), 2100707.

Examples

co2int("CA-ON")
co2int("Canada")$year
co2int("CA-XX")$note
co2int(120)$source
nrow(co2inttab())
# the Canadian provinces in the table, lowest intensity first
tab <- co2inttab()
ca <- tab[startsWith(tab$location, "CA-"), c("location", "name", "g_per_kwh", "year")]
head(ca[order(ca$g_per_kwh), ])
# a location the table does not know is an error, not a silent world average
try(co2int("Atlantis"))

Raw outcome trends by adoption cohort

Description

The mean outcome of each adoption cohort by time since adoption, smoothed by a moving average and scaled to 1 in the period before adoption: a check, before any model, that the cohorts moved together before the policy.

Usage

cohtrnd(data, unit, time, y, ftreat, win = c(-12, 12), ma = 3L, norm = TRUE)

Arguments

Value

A data frame of class rmbl_cohort_trends: cohort, evtm, value. plot() draws one line per cohort, with the period before adoption marked.

Examples

p <- polex()
ct <- cohtrnd(p, "district", "month", "crmrate", "refmonth",
                    win = c(-10, 10))
plot(ct)

A community-weighted index of protective and risk factors

Description

Combines published indicators into one index for each group or place, with weights the community sets: protective factors count for, risk factors against. The weights are never estimated.

Usage

commidx(data, grp, ind, est, lower, upper, w, qual = NULL, conflvl = 0.95)

Arguments

Details

Inuit Tapiriit Kanatami's National Inuit Suicide Prevention Strategy (2016) names risk factors (among them crowded housing, food insecurity, childhood adversity, mental distress) and protective factors (among them cultural continuity, social equity, family strength, healthy child development). How much each matters is a judgement for the people the index describes, so the function takes the weights as given: positive for protective, negative for risk. The index is \sum_i w_i p_i / \sum_i |w_i|, with standard error \sqrt{\sum_i w_i^2 s_i^2} / \sum_i |w_i| from the published intervals (independent estimates). It is a summary of the indicators under the weights given; with other weights it ranks places otherwise, and it is never a ranking of Nations or communities on their culture. For Inuit regions, use Inuit data only.

Value

A data frame of class rmbl_community_index: group, index, lower, upper, usewithcau; attributes weights and conflvl. plot() draws it.

References

Inuit Tapiriit Kanatami (2016). National Inuit Suicide Prevention Strategy. Ottawa. itk.ca/national-inuit-suicide-prevention-strategy

Carroll SR, et al. (2020), see holbal.

Examples

# weights a community might set for adults, 2019/2022 (table 13-10-0924-01); others
# would differ
h <- indiex("health")
h <- h[h$period == "2019/2022" & h$age == "Total, 18 years and over", ]
w <- c(2, 1, -2, -1)
names(w) <- c("Sense of belonging to local community, somewhat strong or very strong",
              "Has a regular healthcare provider",
              "Perceived life stress, most days quite a bit or extremely stressful",
              "Heavy drinking")
ci <- commidx(h, "group", "indicator", "estimate", "lower", "upper",
                      qual = "quality", w = w)
ci
# an indicator the table does not publish for the period leaves the index NA, never
# filled: food insecurity is ".." for 2019/2022
w2 <- c(w, "Household food insecure" = -2)
commidx(h, "group", "indicator", "estimate", "lower", "upper",
                qual = "quality",
                w = w2)$index
plot(ci)

Concentration of a total across units

Description

Measures how unevenly a total is spread across units: whether a few people account for most of the placements, or a few regions for most of the days. ginicoef() gives the Gini coefficient, lorenzc() the points of the Lorenz curve (for plotting, for example with plotlrnz()), and topshr() the share of the total held by the largest units (“the top 1% of people account for 20% of the days”). The three are computed in compiled code from the same sorted values.

Usage

ginicoef(x, na.rm = TRUE)

lorenzc(x, na.rm = TRUE)

topshr(x, fracs = c(0.01, 0.05, 0.1, 0.25), na.rm = TRUE)

Arguments

Details

Gini is the mean absolute difference between pairs of units over twice the mean, G = \sum_i \sum_j |x_i - x_j| / (2 n^2 \bar{x}), which is also twice the area between the Lorenz curve and the diagonal. Zero is a perfectly even spread; the maximum for n units is 1 - 1/n, not 1, so a Gini near 1 requires many units as well as an uneven spread.

A Gini computed from banded or rounded data inherits that coarsening: compare with bandsens() when the values are bands. For a tail that a few units dominate, hilltlidx() says whether the largest values follow a power law.

Errors: x not numeric, with Inf or a negative value, or with no non-missing value; ginicoef() when every value is zero (the Gini of a zero total is undefined); topshr() when fractions is not numeric, is empty, has NA or lies outside [0, 1].

Value

ginicoef(): a numeric scalar in [0, 1 - 1/n]. lorenzc(): a data frame with numeric columns population (cumulative share of units, from 0 to 1) and value (cumulative share of the total, units in ascending order), n + 1 rows including the origin. topshr(): a data frame with columns fraction (the requested fraction of units), units (how many of the largest units that fraction is) and share (the share of the total those units hold).

References

Hedderich, J. and Sachs, L. (2020). Applied Statistics: Methods Using R. Springer-Verlag, Berlin Heidelberg. Section 3.14, p. 117, gives the construction used here: the units are placed at equal intervals on the horizontal axis and the cumulated, ascendingly ordered shares of the total on the vertical one, so that the curve is the diagonal exactly when p percent of the units account for p percent of the total, and sags further the greater the concentration.

See Also

hilltlidx(), bandsens()

Examples

# Ten units holding one each: no concentration.
ginicoef(rep(1, 10))

# One unit holding everything: the maximum for ten units.
ginicoef(c(rep(0, 9), 1))
1 - 1 / 10

placements <- c(rep(1, 1200), rep(3, 430), rep(8, 110), rep(20, 38))
ginicoef(placements)

# What the most frequent few account for.
topshr(placements, c(0.01, 0.05, 0.1))

head(lorenzc(placements))

# missing values are dropped
ginicoef(c(1, NA, 3))

# a zero total has no Gini coefficient
try(ginicoef(c(0, 0, 0)))

Connection across places: does one indicator rise where another does?

Description

For two indicators published for the same places (regions, provinces, Inuit regions), their rank correlation across places, an exact permutation p-value, and an interval that carries the sampling error of the published estimates: for connection to land and wellness, community belonging and mental health, education and health.

Usage

connassoc(
  data,
  place,
  x,
  y,
  xlo = NULL,
  xup = NULL,
  ylo = NULL,
  yup = NULL,
  reps = 10000,
  seed = 1,
  conflvl = 0.95
)

Arguments

Details

Mashford-Pringle and Stewart (2019) describe the land as “an instrumental part of all aspects of Indigenous life including health and well-being”; Mashford-Pringle (2016) describes Indigenous health and education as interconnected. Across places, that connection has a testable shadow: where more people are on the land (or belong, or finish school), do more report wellness? This function measures it by Spearman's \rho. Its p-value is exact when there are at most 8 places (over all n! orderings) and a Monte Carlo estimate otherwise ((b + 1)/(R + 1) for b of R = reps random orderings at least as extreme). Each estimate carries sampling error, so rholo and rhoup are the quantiles of \rho over reps draws of both indicators from normal distributions with their published standard errors. That is a sensitivity band, not a confidence interval for a correlation across places: it shows how much the published sampling error alone moves \rho for these places, it is pulled toward 0 (noise is added to estimates that are already noisy), and it ignores the sampling of the places themselves. Quality flags are not carried: filter suppressed rows out first.

This is an association across places (an ecological one). It cannot say that connection to land makes a person well, which needs data on persons; a handful of places gives little power; and a correlation is not a measure of the relationship the authors describe. It is a description of how the published indicators move together.

Value

A list of class rmbl_connection: data (the places used), rho, p_value, exact, rholo, rhoup, n.

Acknowledgement

This estimator puts to quantitative use the work of Professor Angela Mashford-Pringle (Timiskaming First Nation) and Professor Suzanne Stewart (Yellowknives Dene First Nation), both of the University of Toronto.

See the acknowledgements of the first-nations-metis-inuit-health vignette.

References

Mashford-Pringle A, Stewart SL (2019). Akiikaa (it is the land): exploring land-based experiences with university students in Ontario. Global Health Promotion 26(3_suppl):64-72. DOI 10.1177/1757975919828722

Mashford-Pringle A (2016). Aboriginal health and Aboriginal education are interconnected. Research in the Sociology of Health Care 34:85-96. DOI 10.1108/s0275-495920160000034005

See Also

dualcont.

Examples

# across regions, 2022 (table 41-10-0067-01): where First Nations people report
# excellent or
# very good mental health, do Metis people living in the same region too?
mh <- indiex("mentalhealth")
mh <- mh[mh$rating == "Excellent or very good" & mh$region != "Canada", ]
fn <- mh[mh$group == "First Nations", ]
me <- mh[mh$group == "M\u00e9tis", ]
d <- merge(fn, me, by = "region", suffixes = c("fn", "me"))
ca <- connassoc(d, "region", "estimatefn", "estimateme", "lowerfn",
                             "upperfn",
                             "lowerme", "upperme")
ca
plot(ca)

Full pairwise correlation table

Description

Computes the correlation of every pair of numeric columns and returns it in long form – one row per pair, which is easier to sort, filter and join than a matrix – with the number of complete pairs each figure rests on. The counterpart of corrr::correlate().

Usage

cortbl(data, meth = c("spearman", "pearson"), minprs = 3L)

Arguments

Details

Non-numeric columns are ignored. For each pair, rows where either value is missing are dropped (pairwise deletion), so different pairs can rest on different rows; nprs says how many. Spearman's rank correlation (the default) is the Pearson correlation of the ranks, robust to outliers and to monotone transformations; see topcors() for why it is the default. A pair with fewer than minprs complete rows, or with a constant column (no variance), gets NA. The pairs are listed in column order (utils::combn()). A progress bar is drawn on stderr for many pairs unless options(morie.quiet = TRUE).

Errors: data not a data frame; fewer than two numeric columns; minprs not a whole number of at least 2; an invalid meth.

Value

A data frame of class c("bricklayer_cortable", "data.frame") with one row per pair of numeric columns and the columns x and y (character, the column names), correlation (numeric, in [-1, 1] or NA) and nprs (integer); the method is in attr(, "method").

See Also

topcors() for just the strongest, ccovm() for the matrix form.

Examples

set.seed(1)
df <- data.frame(a = stats::rnorm(100), b = stats::rnorm(100))
df$c <- df$a + stats::rnorm(100, sd = 0.2)

cortbl(df)

# Pearson instead of Spearman
cortbl(df, meth = "pearson")

# nprs shows how much data each figure rests on.
gappy <- df
gappy$a[1:80] <- NA
cortbl(gappy)

# A pair with too little overlap is NA, not a number from nothing.
thin <- df
thin$a[1:99] <- NA
cortbl(thin)$correlation

# one numeric column cannot be correlated
try(cortbl(data.frame(a = 1:5, b = letters[1:5])))

Attest a capsule, binding a signature to what it covers

Description

cpsatst() signs a manifest's canonical digest and records everything a verifier needs alongside it: the scheme, the public key, the context string, the pre-hash, a note and the time of signing. cpschkatst() checks the result against a manifest without being told anything further. Use them to publish a capsule together with proof that its manifest has not changed since you signed it.

Usage

cpsatst(man, key, ctx = NULL, prehash = "none", note = NULL)

cpschkatst(attest, man, keyexp = NULL, key = NULL)

Arguments

Details

What this adds over cpssign(). A bare signature leaves three things implicit – which key, which scheme, and which bytes. A verifier who has to be told those out of band cannot check anything they were not already given, which makes the signature a formality. An attestation states them, so the check is cpschkatst(attestation, manifest) and nothing else.

What is signed: the payload digest (mnfdgst() of the manifest), scheme, public, ctx (hex for a raw context), prehash, note and signed_utc is put in canonical form (mnfcanon()) and its SHA-256 is signed with cpssign(). Every field, the note and the time included, is therefore covered: editing any of them breaks the signature. With an XMSS key from pqckg() the scheme is "xmss-sha256" and the context and pre-hash are not used by the signature; signing uses up one of the key's one-time indices, as every XMSS signature does.

With a shared secret (a string or raw vector) as key the scheme is "hmac": HMAC-SHA-256 under that secret, as cpssign(scheme = "hmac"). public is then a key id, the hex of HMAC-SHA-256 of the string "public" under the secret: it names the key, so keyexp can pin it, without revealing it. The secret itself is never recorded, so checking needs it again (key in cpschkatst()). context is "" and prehash "none", the bytes morie's capsule signer writes. HMAC is symmetric: anyone who can check the signature can also make one. It is not one of the post-quantum signature schemes and carries none of their guarantees.

The check, one row each: attestation_complete (every field is present), manifest_digest (the manifest's digest equals the attested one), key_as_expected (only when keyexp is given) and signature (the signature verifies under the attested public key; for "hmac", under key, whose key id must equal public). ok is TRUE when every row passes.

What it does NOT establish. That the public key belongs to whoever you think it does: an attestation is only as good as the channel the key arrived on, which is why keyexp exists. And that the manifest is true – only that it has not changed since it was signed. cpsfals() is for the other question. No chain seal is recorded; to attest a run history, sign chnseal() directly.

Errors: man not a list; key not from fipskg() or pqckg() and not a non-empty string or raw vector; an invalid prehash or ctx (from cpssign()), or either given with an HMAC key; attestation not from cpsatst(). A signature that fails to verify is a failed row, not an error.

Value

cpsatst() a list of class bricklayer_attestation with the character fields digest, scheme, public (hex), context, prehash, note and signed_utc (ISO 8601 UTC time), and signature (a bricklayer_signature from cpssign()).

cpschkatst() a list of class bricklayer_attestation_check with ok (logical) and checks, a data frame with columns check, ok and detail, one row per check.

See Also

cpssign(), mnfdgst(), cpsfals(), chnseal().

Examples

m <- mkmnf(list(dataset = "otis", rows = 1200L),
  env = FALSE
)
key <- fipskg("ML-DSA-65")
att <- cpsatst(m, key, note = "counts as published")
att

# a verifier needs the attestation and the manifest, nothing else
res <- cpschkatst(att, m)
res$ok
res$checks

# any change to the manifest breaks it
m2 <- m
m2$meta$rows <- 1201L
cpschkatst(att, m2)$ok

# and so does presenting a different key
cpschkatst(att, m,
  keyexp = fipskg("ML-DSA-65")$public
)$ok

# the note is signed too: editing it breaks the signature row
forged <- att
forged$note <- "counts corrected"
cpschkatst(forged, m)$checks[, c("check", "ok")]

# a context string and a pre-hash
att2 <- cpsatst(m, key, ctx = "release 2026-10", prehash = "sha512")
cpschkatst(att2, m, keyexp = key)$ok

# a shared secret signs with HMAC; `public` is its key id, not the secret
att3 <- cpsatst(m, "shared-secret")
att3$public
cpschkatst(att3, m, key = "shared-secret", keyexp = att3$public)$ok
cpschkatst(att3, m, key = "another-secret")$checks[, c("check", "ok", "detail")]

# a key that is not a signing key is refused
try(cpsatst(m, 42))

Bundle a capsule into one signed artifact

Description

Records a digest of every file named, the manifest's canonical digest, and an attestation covering both, as a single JSON file. cpsbndlvrf() re-hashes the files on disk and checks everything against it.

Usage

cpsbndl(
  dir,
  man,
  key,
  files = NULL,
  ctx = NULL,
  prehash = "none",
  note = NULL,
  path = NULL
)

cpsbndlrd(path)

cpsbndlvrf(bundle, dir, man = NULL, keyexp = NULL, key = NULL, seen = NULL)

Arguments

Details

The file digests are inside the signed payload, not alongside it. That is the whole point: a list of hashes that is not itself signed can be rewritten to match whatever the files now say, and a recipient who re-hashes the files and compares them to that list learns only that the list is consistent with itself.

What it establishes: the named files have not changed since signing, the manifest has not changed, and both were signed together by the holder of that key. What it does not: that the key belongs to anyone in particular, or that the manifest's claims are true.

How it is built. Each file is hashed with SHA-256 (sha256file()); the payload is the list of path, sha256 and size_bytes, the manifest's mnfdgst() and the UTC time; cpsatst() signs a manifest whose meta is that payload, and the bundle (payload plus attestation) is written as pretty JSON with doubles at 17 digits. With files = NULL every regular file under dir is hashed, recursively and including dot-files (a .Rprofile dropped into a capsule runs code), except the bundle file itself. A progress bar is shown when options(morie.progress) allows it.

How it is checked. cpsbndlvrf() re-hashes every listed file (a path that is absolute or climbs out of dir fails its row), reports every file on disk the bundle does not list (the bundle's own file excepted), checks the attestation with cpschkatst() (signature, digest, and keyexp when given; an "hmac" signature needs the shared secret as key), and, when man is given, compares its digest with the bundled one. Every check is a row of checks; ok is TRUE only when all pass. A failed check is a row, not an error.

Errors: dir missing or not a directory; no files to bundle, or a named file not in dir; a key cpsatst() cannot sign with; cpsbndlrd() given a file that is not a bundle; cpsbndlvrf() given something other than a bundle or a path. Reading a bundle whose signature records a stateful key's state warns (see cpsatst()).

Value

cpsbndl(): the bundle, invisibly, a list of class bricklayer_bundle with files (a list of path, sha256, size_bytes), manifest_digest, bundled_utc and attestation (from cpsatst()), and the path it was written to as attribute "path".

cpsbndlrd(): the same structure read from the file.

cpsbndlvrf(): a list of class bricklayer_bundle_check with ok (logical) and checks (a data frame of check, logical ok and detail: one file:PATH row per listed file, no_unlisted_files, the ⁠attestation:⁠ rows, with a manifest manifest_digest, and with seen xmss_index). For an XMSS bundle the attestation:signature row's detail names the signature's index and the key fingerprint.

See Also

cpsatst(), vrfcps(), mnfdgst().

Examples

dir <- tempfile()
dir.create(dir)
write.csv(data.frame(x = 1:3), file.path(dir, "data.csv"),
          row.names = FALSE)
m <- mkmnf(list(dataset = "demo"), env = FALSE)
key <- fipskg("ML-DSA-44")

b <- cpsbndl(dir, m, key, note = "as published")
cpsbndlvrf(attr(b, "path"), dir, man = m)$ok

# touch a byte of the data and the bundle no longer holds
write.csv(data.frame(x = 1:4), file.path(dir, "data.csv"),
          row.names = FALSE)
res <- cpsbndlvrf(attr(b, "path"), dir, man = m)
res$ok
res$checks[!res$checks$ok, c("check", "detail")]

# read back, and checked against the key you expected to have signed
b2 <- cpsbndlrd(attr(b, "path"))
b2
cpsbndlvrf(b2, dir, keyexp = "00ff")$ok

# a file added beside the bundle is reported too
writeLines("x", file.path(dir, "extra.txt"))
cpsbndlvrf(b2, dir)$checks[, c("check", "ok")]

# a bundle signed with a shared secret (HMAC), as morie signs
dir2 <- tempfile()
dir.create(dir2)
writeLines("x,y", file.path(dir2, "d.csv"))
h <- cpsbndl(dir2, m, "shared-secret")
h$attestation$public  # the key id, not the secret
cpsbndlvrf(attr(h, "path"), dir2, key = "shared-secret",
           keyexp = h$attestation$public)$ok
# without the secret the signature row cannot pass
cpsbndlvrf(attr(h, "path"), dir2)$ok
unlink(dir2, recursive = TRUE)

# a directory that does not exist
try(cpsbndl(file.path(tempdir(), "no-such-capsule"), m, key))
unlink(dir, recursive = TRUE)

Check every capsule under a folder

Description

Runs the checks of the rmbl GitHub Action on every capsule (a folder holding data_provenance.json or capsule_bundle.json) under path: each capsule file against its JSON Schema (schemapath), then vrfcps and cpsbndlvrf. With mode = "drift" or "refresh" each source is fetched again and compared with the pin, column by column (cpsdrft and benfordtst); "refresh" then writes and pins new bytes. rmbl verify, drift and report call it. format() renders the Markdown report or the JSON, and ghann, ghcmt and ghout the GitHub pieces.

Usage

cpschk(
  path = ".",
  mode = c("verify", "drift", "refresh", "report"),
  key = "",
  pubkey = "",
  psiwarn = 0.1,
  psifail = 0.25,
  base = NULL,
  cache = NULL,
  allowfile = FALSE,
  snapshot = FALSE,
  wayback = FALSE,
  glob = "",
  fp = TRUE,
  seenlog = ""
)

## S3 method for class 'rmbl_cpschk'
format(x, fmt = c("md", "json"), failon = NULL, ...)

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

Arguments

Details

Each check is one row with a severity: error, drift, warn or note (a note never changes a result). A capsule's result is the worst severity among its failed rows, and status the worst over all capsules: ok, warn, drift or error. Each failed row adds a note, and one more from explerr where it finds what the row means or where its message is raised.

Drift: when the fetched source's SHA-256 differs from the pin and the pinned data file is a CSV, cpsdrft() runs with psithr = psifail. A column is drift when cpsdrft() flags it and its PSI is above psifail, and warn when it is flagged below that or its PSI is above psiwarn; an added or removed column is drift. New bytes whose columns all stay within the bands are a note; new bytes that cannot be compared are a warning.

Refresh: the new bytes are written as the data file and pinned in the provenance as cpsrfrsh does (refused when they fail the schema the provenance pins, valsch); with wayback = TRUE an http(s) source is first captured with wbsave. The capsule is then verified again.

Manifests: with base, each capsule's manifest.json and results/manifest.json is compared with its base copy: the base's recorded results against this copy's observed values, each within the base's tolerance (mnfrcmp), and the resources of the two (mnfdiff). A manifest only one side has is added or removed.

Cache: a fetched source is kept in the cache folder under source- and the first 16 hex digits of the SHA-256 of its URL, one file per source overwritten by each fetch, and is read only when the live fetch fails.

Side effects: "drift" and "refresh" fetch each source (the network for http(s) ones), as snapshot = TRUE asks the Wayback Machine once per capsule (wbsnpurl); "refresh" writes the capsules; a cache folder is written to. Errors: path is not a directory; a bad argument. A capsule whose checks stop is one failed row, not an error.

Value

A list of class rmbl_cpschk: status, mode, path, glob, rows (a data frame of every check: capsule, part, file, check, ok, detail, severity and line, the line of the file where known), capsules (a list per capsule: capsule, source, pinned_sha256, fetched_sha256, changed, result, checks_passed, checks_total and notes, and where made drift with its columns table, drift_summary, snapshot, snapshot_utc and manifests), base (its label), psiwarn, psifail and footprint (the co2fp() result without its value, or NULL). The cpsdrft() result of each capsule whose columns were compared, with its provenance and manifest diff, is in the attribute drifts (named by capsule), for drfthtml and drftapp. format() returns one string and print() writes the Markdown report.

See Also

blcli (the verbs), ghann, cpsrfrsh, explerr.

Examples

root <- tempfile("capsules-")
d <- usecpstmpl(file.path(root, "one"), ex = TRUE)
invisible(cpsrfrsh(d))
x <- cpschk(root, fp = FALSE)
x$status
x$rows[, c("check", "ok")]
x

# the data changed: the digest row fails, and the notes say what it means
p <- loadprov(file.path(d, "data_provenance.json"))
cat("1,2\n", file = file.path(d, p$resource$filename), append = TRUE)
x <- cpschk(root, fp = FALSE)
x$status
x$capsules[[1]]$notes
substr(format(x, fmt = "json"), 1, 120)

# drift: fetch the source again (here a local file) and compare
cpschk(root, "drift", psiwarn = 0.05, fp = FALSE)$status

try(cpschk(tempfile()))
try(cpschk(root, psiwarn = 1))

Compare a fetched data frame with the one a capsule was pinned against

Description

Checks whether the data in a freshly downloaded table still look like the data an analysis was built on, column by column. A digest only says whether the bytes changed; this says whether the distributions changed, which is what invalidates a result. Use it when re-running an analysis on a re-released extract.

Usage

cpsdrft(
  ref,
  cur,
  alpha = 0.01,
  psithr = 0.25,
  psiminn = 1000L,
  bins = 10L,
  idelvls = 100L,
  seed = 1L
)

Arguments

Details

Runs the appropriate drift test on every shared column and collects the verdicts in one report: drftks() plus drftpsi() for a numeric column, drfthomog() for a categorical one. Columns that appeared or vanished are listed separately, since no test applies to them.

The categorical test is the two-sample homogeneity test, NOT drftchisq() 's goodness-of-fit against a known distribution: the reference here is itself a finite sample, and ignoring its sampling error would report drift too readily.

This is the check that a byte-level digest cannot make. A re-released extract legitimately has a different SHA-256 while being the same data statistically; conversely a column can keep its name, type and row count while having been silently rescaled. cpsdrft() asks whether the DATA moved.

Column types: a column numeric on both sides is "numeric"; a character or factor column whose every value is an ISO date (YYYY-MM-DD) or date-time on both sides is converted and tested as "date"; a categorical column with more than idelvls distinct values, or with at least 20 distinct values that number more than a fifth of the rows, is an "identifier" and is not tested; anything else is "categorical". Missing values are dropped. For a categorical column the PSI is computed over the reference's categories with every unseen current value pooled into one extra bin.

A column is drifted when its p-value is below alpha, or when its PSI exceeds psithr and both samples have at least psiminn values. anydrft is also TRUE when a column was added or removed. Each column runs one or two tests, so with many columns some will fall below alpha by chance; see scanadj() for a multiplicity correction and drftcal() for an empirical false-alarm rate. A progress bar is shown in an interactive session for many columns.

Errors: either argument not a data frame; alpha not strictly between 0 and 1; psiminn negative; idelvls below 2.

Value

A list of class bricklayer_drift with elements:

columns

a data frame, one row per shared column, with column, type ("numeric", "date", "categorical" or "identifier"), statistic (KS D or chi-square), p_value, psi, jsdiv, unseen_share (share of current values the reference never had; categorical and identifier columns only), drifted (logical, NA when no test applies) and note (why a column was not tested, or "");

added, removed

character vectors of column names;

nref, ncur

row counts;

alpha, psithr, psiminn

the settings used;

bands

the PSI bands the plots and reports colour by: warn (0.1, or psithr when lower) and fail (psithr);

summ

one binned summary per shared column, for plot(), drfthtml() and drftapp(): missing shares, sizes, summary statistics or the ten largest categories, the PSI bins with each bin's term of the index, a histogram and first-digit counts. Never the rows;

anydrft

logical.

The print method lists the per-column verdicts, drifted columns first.

See Also

valsch() for the structural check, which this complements rather than replaces.

Examples

set.seed(5)
ref <- data.frame(
  value = stats::rnorm(300),
  size = stats::runif(300, 1, 10),
  grade = sample(c("a", "b", "c"), 300, TRUE)
)

# A fresh draw from the same process: no drift.
same <- data.frame(
  value = stats::rnorm(300),
  size = stats::runif(300, 1, 10),
  grade = sample(c("a", "b", "c"), 300, TRUE)
)
d <- cpsdrft(ref, same)
d$anydrft

# A silently rescaled column, and a new category, are both caught.
moved <- same
moved$size <- moved$size * 3
moved$grade[1:100] <- "z"
cpsdrft(ref, moved)

# The PSI is always reported as an effect size, but on a sample this
# small it is not allowed to raise the flag by itself -- binning noise
# alone would clear 0.25. Lower psiminn to override that.
cpsdrft(ref, same)$columns$psi

# Structural changes are reported rather than tested.
cpsdrft(ref, same[, c("value", "grade")])$removed

try(cpsdrft(ref, as.matrix(same)))

Run falsification controls against a statistic

Description

Recomputes stat under conditions in which its value is known in advance, and reports whether it behaved. Four controls, each answering a different way of being wrong:

Usage

cpsfals(data, stat, trt = NULL, n = 199L, subfrac = 0.8, seed = NULL)

Arguments

Details

• permutation – shuffle trt and the association it carries is destroyed, so the statistic should fall to its null distribution. Reports where the observed value sits in that distribution, and the smallest p-value the number of permutations could have produced.

• random common cause – add a column of noise. It cannot possibly matter, so a statistic that moves is reading something other than the data.

• placebo treatment – replace trt with a random draw of the same shape. The effect should vanish; if it does not, the statistic is picking up structure that has nothing to do with the exposure.

• subset stability – recompute on random subsets. Wide spread means the finding rests on particular rows, which is worth knowing before it rests on a conclusion.

What a pass means. Only that these controls did not catch anything. They are falsification tests, so they can refute and cannot confirm: passing all four is consistent with a statistic that is wrong for a reason none of them probes.

When each control passes:

• permutation – the two-sided permutation p-value (1 + \#\{|s_b| \ge |s_{obs}|\}) / (1 + B) over the B usable permutations is below 0.05 (the observed statistic is unlike its null);

• placebo – the statistic under one more permuted treatment lies within the 95th percentile of the absolute permutation values;

• random common cause – adding a column .rmbl_random_common_cause of standard normal noise changes the statistic by at most 10^{-8} \max(1, |s_{obs}|);

• subset stability – the observed value lies inside the middle 95% of the statistic over n random subsets of round(subfrac * nrow(data)) rows (at least two).

A permutation or subset on which the statistic errors or is not finite is dropped from that control. The statistic is evaluated ⁠2 n + 3⁠ times, so its cost sets the run time; a progress bar is shown when options(morie.progress) allows it. With seed the random number stream is set for the call and the caller's stream restored afterwards.

Errors: data not a data frame or with fewer than 4 rows; stat not a function, or failing or not returning one finite number on the data as given, the placebo data or the noise-augmented data; n below 9; subfrac not strictly between 0 and 1; trt not a column of data.

Value

A list of class bricklayer_falsification: observed (the statistic on the data as given), controls (a data frame with one row per control: control, logical passed (NA when skipped), numeric value (the p-value, the placebo or noise statistic, or the subset spread) and detail), permutation (the vector of permuted statistics, NULL without treatment), n, treatment, rng_kind and seed (NA when none was given). print() shows one line per control; plot() draws the permutation distribution.

See Also

cpspwr() for the positive control – whether a real effect would have been seen at all; cpsatst() for the different question of whether the record is intact rather than whether the finding survives.

Examples

set.seed(1)
d <- data.frame(x = rnorm(200))
d$y <- 0.8 * d$x + rnorm(200)
# a real association survives its controls
real <- cpsfals(d, function(z) cor(z$x, z$y),
                        trt = "x", n = 199, seed = 42)
real$controls[, c("control", "passed")]

# and a statistic that ignores the data fails the ones that can see it
fake <- cpsfals(d, function(z) 0.5, trt = "x",
                        n = 199, seed = 42)
fake$controls[fake$controls$control == "permutation", "passed"]

# without a treatment column two controls are skipped, and say so
fit <- cpsfals(d, function(z) mean(z$y), n = 49, seed = 1)
fit$controls[, c("control", "passed")]

# the statistic must return one finite number
try(cpsfals(d, function(z) c(1, 2), n = 49))

Detect an injected effect of known size

Description

Adds an effect of each given size to the data, reruns the whole detection procedure, and reports how often it was found. The result is a power curve, and the smallest size detected reliably is the smallest effect the analysis could have seen. Use it before reporting a null result, to show whether the analysis could have found an effect of the size that matters.

Usage

cpspwr(
  data,
  stat,
  inject,
  sizes = c(0, 0.2, 0.5),
  trt = NULL,
  n = 99L,
  reps = 10L,
  alpha = 0.05,
  seed = NULL
)

Arguments

Details

This is the positive control that cpsfals() lacks: its four are all negative. Those establish that the finding is not an artefact; this establishes that a real effect would not have been missed. A study that passes every falsification control and has no power against the effect it was looking for has not found that the effect is absent – it has found nothing either way, and the two are routinely reported as the same thing.

For each size (0 is always added) and each of reps repetitions, the statistic is computed on the data passed through the injector, as inject(data, size), and a permutation null is built by shuffling the trt column n times and recomputing the statistic. The two-sided permutation p-value is

p = \frac{1 + \#\{|T^*| \ge |T|\}}{1 + n_{ok}}

over the n_{ok} permutations with a finite statistic (Phipson and Smyth's correction, so p is never 0), and the effect counts as detected when p \le \alpha. A repetition whose injection or observed statistic fails (error or non-finite value) is skipped and not counted in reps. Since the data are the same in every repetition, inject should add randomness of its own (new noise, a resampled effect) when repetitions are meant to differ; otherwise only the permutations vary.

The cost is |sizes| \times reps \times (n + 1) evaluations of stat. A progress bar is drawn on stderr unless options(morie.quiet = TRUE). With seed, R's random number state is set for the run and restored afterwards; without it, the run uses and advances the session's generator.

Errors: data not a data frame of at least 4 rows; stat or inject not a function; trt missing or not a column; n below 9; reps below 1; n so small that 1/(n+1) > \alpha, so that nothing could ever be detected.

Value

A list of class c("bricklayer_power", "list"):

curve

data frame with one row per size: size, detected (integer), reps (integer, the usable repetitions) and rate (detected / reps, NA when none was usable).

alpha, n, reps, treatment

the settings used.

smadet

the smallest positive size found at a rate of at least 0.8, or NA when no size reached it.

false_positive_rate

the rate at size 0.

See Also

cpsfals() for the negative controls.

Examples

set.seed(1)
d <- data.frame(x = rnorm(120), y = rnorm(120))
# inject a linear effect of x on y
inj <- function(z, size) {
  z$y <- z$y + size * z$x
  z
}
pw <- cpspwr(d, function(z) cor(z$x, z$y), inj,
                    sizes = c(0, 0.3, 0.6), trt = "x",
                    n = 99, reps = 5, seed = 42)
pw$curve
pw

# The rate at size 0 estimates the false-positive rate. With few
# repetitions it is usually 0, because alpha is small -- five draws
# at 0.05 come up empty about three times in four -- so read it as a
# sanity check that it is not LARGE, not as an estimate of alpha.
pw$curve[pw$curve$size == 0, "rate"]

# an injection with fresh noise each time, and effects too small to be
# detected in four repetitions: the smallest detected size is NA
inj2 <- function(z, size) {
  z$y <- z$y + size * z$x + stats::rnorm(nrow(z), sd = 0.5)
  z
}
cpspwr(d, function(z) cor(z$x, z$y), inj2, sizes = c(0.1, 0.2),
              trt = "x", n = 49, reps = 4, seed = 7)$smadet

# too few permutations to ever reach alpha = 0.05
try(cpspwr(d, function(z) cor(z$x, z$y), inj, trt = "x", n = 9))

Fetch a capsule's source again and pin the new bytes

Description

Downloads the source named in a capsule's data_provenance.json (resource$direct_url, else source$url; an https:// URL, or a local file inside the capsule), and when its bytes differ from the pin writes them as the capsule's data file and pins them in the provenance, in the form the analysis.R of usecpstmpl writes on its first run: source$sha256, and resource$filename, sha256, size_bytes and retrieved_utc (row_count_data_rows too when it was recorded). vrfcps then verifies the capsule against the new data. From the shell: rmbl refresh PATH.

Usage

cpsrfrsh(dir, allowfile = FALSE)

Arguments

Details

Nothing is written when the bytes are unchanged. A download goes through fdl, with the provenance's wayback$csv_snapshot as its fallback (no other Wayback copy is looked up). A source is fetched over https:// only. A local path (relative to the capsule) or a file:// URL is read only when it resolves, symbolic links followed, to a file inside the capsule, inside the allowfile folder, or in the package's extdata; never under /proc, /sys or /dev. A gzip source may inflate to at most the option rmoriebricklayer.maxcsvbytes (1 GiB by default) when it is read. A capsule bundle is not signed again: that needs your key (cpsbndl). Side effects: writes the data file and data_provenance.json in dir. Errors: dir is not a directory or holds no provenance; the provenance names no source, or a TODO one, or a file:// one without allowfile; a URL that is not https://; a local source outside those folders; the source cannot be fetched, or inflates past the limit; resource$filename is not a plain relative path inside dir; data_provenance.json is a link out of dir; the new data fail the schema the provenance pins (valsch).

Value

The mnfdiff of the provenance before and after: all "unchanged" when the source still matches its pin.

See Also

mnfdiff, vrfcps, usecpstmpl.

Examples

d <- usecpstmpl(tempfile("capsule-"), ex = TRUE)
cpsrfrsh(d)            # the first pin: the data file is added
vrfcps(d)$ok
cpsrfrsh(d)$change     # the source has not changed since

try(cpsrfrsh(tempfile()))

Assess a capsule in one call

Description

Runs the checks this package provides over one data frame and collects the findings into a single report: the structural schema check, the distributional comparison against a reference, the missingness picture, the multivariate outliers, a Benford screen on the wide numeric columns, and the integrity digests. Use it as the first look at a capsule or a fresh extract, before any analysis, and keep the Markdown version (rptmd()) with the outputs.

Usage

cpsrpt(
  data,
  ref = NULL,
  schema = NULL,
  rules = NULL,
  chunks = NULL,
  sig = NULL,
  key = NULL,
  alpha = 0.01,
  maxrowsoutl = 20000L
)

Arguments

Details

Nothing here is new arithmetic. The value is that the answers arrive together and ordered by severity, because the failure mode this is built against is a person running one check, seeing it pass, and concluding the data is fine.

Value

A list of class c("bricklayer_report", "list"): findings (a data frame of character columns severity, check ("schema", "rule", "drift", "missing", "shape", "outliers", "benford" or "signature"), subject (the column, rule or row count concerned) and detail), verdict (the worst severity: "fatal", "warn", "note", or "clean" when there is no finding), profile (from profcols()), missingness (from misssmry()), drift (the cpsdrft() result, or NULL without a reference), digest (a list with data_digest, and with chunks also merkle_root, merkle_scheme (the mrkroot() hashing scheme, 2) and n_chunks, and, with a signature and key, the logical signature_valid), and nrows / n_cols. print() shows the verdict and the findings; summary() gives the verdict and the count of each severity.

What severity means

• fatal – nothing downstream can run as it stands: a schema check or a rule of severity "fatal" failed (for example a required column is missing), or a column of the reference is absent from the data.

• warn – something moved or is badly incomplete: a schema check, or a rule of severity "warning", failed (a value left its pinned range, a type changed), a column drifted from the reference, a column is entirely missing or more than half missing, or a signature did not verify (or came with no key).

• note – worth a look but not necessarily wrong: a column not in the reference, a column whose drift could not be tested, a contiguous run of missing values (one outage rather than scattered gaps), a constant column, values beyond the Tukey fences, multivariate outliers or collinear columns, a Benford departure.

How each check runs

The drift comparison is cpsdrft() at level alpha. The missing run is the longest run from missruns() of at least max(2, nrow(data) %/% 20) rows. The outlier scan is mahaoutl() on the numeric columns that are not constant or entirely missing, run only when there are at least two of them, at most maxrowsoutl rows and more complete rows than columns. The Benford screen (benfordtst()) runs on a numeric column with at least 200 non-zero finite values whose largest absolute value is more than 1,000 times its smallest, and reports p below 0.01. The digest is dgstobj() of the data; a signature is checked with cpsvrf() against that digest.

Findings are sorted by severity, then check, then subject. Nothing is written and nothing is fetched. Errors when data is not a data frame or has no columns, and when schema or rules are malformed (from valsch() and valrules()).

A clean report is not proof the data is correct. It means these particular checks found nothing, and every one of them has a stated blind spot – see cpsdrft() on statistical power and benfordtst() on why a departure is a screen rather than a verdict.

See Also

cpsdrft(), profcols(), valsch(), rptmd() to write it out.

Examples

set.seed(1)
ref <- data.frame(
  id = 1:200,
  score = stats::runif(200, 0, 10),
  grade = sample(c("a", "b", "c"), 200, TRUE),
  stringsAsFactors = FALSE
)

# A fresh extract from the same process: nothing to report.
cur <- ref
cur$score <- stats::runif(200, 0, 10)
cpsrpt(cur, ref = ref, schema = infsch(ref))

# One rescaled column and a new category: both surface, worst first.
bad <- cur
bad$score <- bad$score * 5
bad$grade[1:80] <- "z"
r <- cpsrpt(bad, ref = ref, schema = infsch(ref))
r
r$verdict
r$findings[, c("severity", "check", "subject")]

# A missing column is fatal, because nothing downstream can run.
cpsrpt(bad[, c("id", "score")], ref = ref,
               schema = infsch(ref))$verdict

# validation rules join the findings, and chunks are pinned with a Merkle root
r2 <- cpsrpt(bad, rules = list(ruleinset("grade", c("a", "b", "c"))),
                     chunks = c("part one", "part two"))
summary(r2)
r2$digest$merkle_root

# a data frame is required
try(cpsrpt(1:10))

Pin and sign every capsule under a folder

Description

Seals each capsule (a folder holding data_provenance.json) under path, so that a program in any language only has to write the data file and a plain provenance naming it ({"resource": {"filename": "data.csv"}}): on its first run the provenance is pinned to the data file's bytes as cpsrfrsh pins it (resource$sha256, size_bytes, retrieved_utc, the row count when one is recorded), then cpsbndl signs every file of the capsule into capsule_bundle.json, over its manifest.json when it has one. rmbl sign calls it and then verifies (cpschk).

Usage

cpsseal(path = ".", key, overwrite = FALSE, glob = "")

Arguments

Details

A provenance already pinned (one with resource$sha256) is left as it is. Side effects: writes data_provenance.json (first run) and capsule_bundle.json in each capsule. Errors, before anything is written: no provenance under path; a capsule that already has a bundle, unless overwrite; then, per capsule, a resource$filename that is not a plain relative path to a file in the capsule, and the refusals of cpsbndl() (a key it cannot sign with).

Value

A data frame, one row per capsule: capsule, pinned (whether the provenance was pinned now), files (the number signed) and public (the HMAC key id or the public key; never the secret).

See Also

cpsbndl, cpschk, blcli (rmbl sign).

Examples

root <- tempfile("capsules-")
d <- file.path(root, "one")
dir.create(d, recursive = TRUE)
utils::write.csv(data.frame(x = 1:3), file.path(d, "data.csv"),
                 row.names = FALSE)
writeLines('{"resource": {"filename": "data.csv"}}',
           file.path(d, "data_provenance.json"))
s <- cpsseal(root, "a shared secret")
s
cpschk(root, key = "a shared secret", fp = FALSE)$status

# a second run needs overwrite, and a post-quantum key signs too
try(cpsseal(root, "a shared secret"))
cpsseal(root, fipskg("ML-DSA-44"), overwrite = TRUE)$pinned

Sign a capsule manifest

Description

Authenticates msg – normally a manifest digest, or the whole manifest text – so a verifier can tell that it came from the holder of the key and has not been altered since. A digest alone shows that data were not corrupted; a signature also shows who vouched for them, because only the key holder can make it. Verify with cpsvrf().

Usage

cpssign(
  msg,
  key,
  scheme = NULL,
  ctx = NULL,
  determ = FALSE,
  prehash = "none",
  masked = getOption("rmoriebricklayer.masked", TRUE)
)

Arguments

Details

Three kinds of key are accepted: a shared secret (HMAC), an XMSS key from pqckg(), or a standardised ML-DSA or SLH-DSA key from fipskg(). For new work the standardised keys are the simplest: they sign any number of messages with no state to keep.

With a standardised key the signature is FIPS 204 ML-DSA.Sign or FIPS 205 slh_sign over the message with its context (hedged by default, deterministic on request). ML-DSA signs in milliseconds; SLH-DSA from a few hundred milliseconds (the f sets) to several seconds (the s sets at 192 and 256 bits); an XMSS signature rebuilds the key's tree, about 3 s at the default height 10. Slow signatures show a spinner on stderr after 0.75 s when a person is watching (an interactive session or options(morie.progress = TRUE); options(morie.quiet = TRUE) silences it), and Ctrl-C stops them.

Refusals: a message that is not one string or a raw vector; a key of the wrong kind for the scheme; a context over 255 bytes; an unknown pre-hash; an exhausted XMSS key, an index the used-index record shows as used, a key whose record cannot be read or written, or an XMSS key the record does not know (see xmssadv() and xmssregkey()).

With scheme = "hmac" the key is a shared secret string and the result is an HMAC-SHA-256 tag. Symmetric, so anyone who can verify can also sign.

With scheme = "xmss" the key is a pqckg() object and the result is a post-quantum one-time signature under the key's Merkle root. Asymmetric: a verifier holding only the public root cannot forge.

Value

A list of class bricklayer_signature: scheme, signature, and for XMSS also auth, index, root, height, randomizer (the per-signature R of RFC 8391, which a verifier needs and which therefore travels with the signature) and wire – the signature in the RFC's own byte order, index || R || WOTS || auth, hex-encoded. wire is byte-identical to what the XMSS reference implementation produces from the same key material, so it can be handed to another implementation as bytes.

The key must be advanced after every signature

An XMSS signature consumes a leaf, and subsequent signing must use the advanced key – reusing an index breaks the scheme. The signature object carries no key material (it is what verifiers are handed); the advanced key is xmssadv(key, signature), and that is what to save. The key is an R value, so cpssign() cannot advance the caller's copy; instead it records, in the session and in the index file described under xmssadv(), the highest index each key has signed with, and refuses that index or an earlier one. Passing an exhausted key is an error, not a silent wrap-around.

Security

The implementation is hand-written. Its standardised schemes are checked byte for byte against OpenSSL 3.5 and NIST known-answer vectors in the test suite. Side channels are checked, not assumed: inst/ctcheck runs every operation that touches a secret under valgrind memcheck with the secret marked undefined (the ctgrind method), so a branch or a memory address that depends on it is a reported error, and afterwards scans the dead stack for copies of the secret. Both checks run in CI on every change, with GCC and with Clang. Timing is also measured on x86-64 and arm64 hardware (inst/dudect), power leakage is assessed in simulation under the value and the transition models (inst/tvla), and ML-KEM decapsulation and ML-DSA signing are first-order masked by default. These are checks of this code on those compilers and under those leakage models; the README's security section lists what each check does and does not cover.

See Also

cpsvrf(), chmac()

Examples

op <- options(rmoriebricklayer.xmssidxfile = tempfile())  # keep this example's
# one-time index record out of your configuration directory; leave the default in use
# Symmetric: one shared secret.
sig <- cpssign("sha256:abc123", key = "shared-secret",
                    scheme = "hmac")
cpsvrf("sha256:abc123", sig, "shared-secret")
cpsvrf("sha256:TAMPERED", sig, "shared-secret")

# Post-quantum: the verifier needs only the public root.
key <- pqckg(height = 2)
s1 <- cpssign("manifest-1", key)
cpsvrf("manifest-1", s1, signpubkey(key))

# Advance the key past the signature before signing again.
key <- xmssadv(key, s1)
key$next_index
s2 <- cpssign("manifest-2", key)
cpsvrf("manifest-2", s2, signpubkey(key))

# A signature does not transfer to another message.
cpsvrf("manifest-1", s2, signpubkey(key))

# A standardised key: no state, any number of signatures.
fk <- fipskg("ML-DSA-44")
s3 <- cpssign("manifest-3", fk, ctx = "release", determ = TRUE)
cpsvrf("manifest-3", s3, fipspubkey(fk), ctx = "release")

# An XMSS signature needs an XMSS key.
try(cpssign("manifest-4", "shared-secret", scheme = "xmss"))
options(op)

Verify a capsule manifest signature

Description

Checks sig against msg. For "hmac" the comparison is constant-time (verified under valgrind memcheck in CI: inst/ctcheck, cases digest-equal and hmac-sha256). For XMSS the Winternitz chains are walked to their ends and the authentication path replayed to the Merkle root; the digest is bound to both the leaf index and the root, so a signature cannot be replayed at another index or under another key.

Usage

cpsvrf(msg, sig, key, ctx = NULL, prehash = NULL)

Arguments

Details

Returns FALSE rather than erroring on a malformed or truncated signature: a verifier must treat unparseable input as “not verified”, never as an exception to be caught and ignored.

TRUE means the signature is valid for msg under key (and, for the standardised schemes, ctx and the pre-hash). Verification needs only public material, except for HMAC, where the verifier holds the same secret as the signer and could therefore also sign. HMAC tags are compared in constant time.

Only two inputs are errors rather than FALSE: a sig that is not a bricklayer_signature at all, and, for an XMSS signature, a key that is not a key object (a string, say) or a msg that is not one string or a raw vector. A key or signature whose class, scheme or fields were edited returns FALSE. An XMSS signature object written by 0.5.9 to 0.5.12 (one that carries a key_state field) still verifies, with a warning that it exposed the signer's private key.

Value

TRUE or FALSE (a length-1 logical, never NA).

See Also

cpssign(), signpubkey(), fipspubkey().

Examples

op <- options(rmoriebricklayer.xmssidxfile = tempfile())  # keep this example's
# one-time index record out of your configuration directory; leave the default in use
key <- pqckg(height = 2)
sig <- cpssign("pinned-manifest", key)
pub <- signpubkey(key)

cpsvrf("pinned-manifest", sig, pub)

# Every way of being wrong returns FALSE.
cpsvrf("edited-manifest", sig, pub)              # message changed
bad <- sig; bad$signature <- paste0("ff", substring(bad$signature, 3))
cpsvrf("pinned-manifest", bad, pub)              # signature edited
cpsvrf("pinned-manifest", sig,
               signpubkey(pqckg(height = 2)))  # foreign key

# A truncated signature is not verified, and does not error.
trunc <- sig; trunc$signature <- substring(sig$signature, 1, 64)
cpsvrf("pinned-manifest", trunc, pub)

# HMAC: the verifier holds the same secret.
tag <- cpssign("pinned-manifest", "shared-secret", scheme = "hmac")
cpsvrf("pinned-manifest", tag, "shared-secret")

# Something that is not a signature at all is an error, not FALSE.
try(cpsvrf("pinned-manifest", list(signature = "00"), pub))
options(op)

Association in a contingency table

Description

Cramer's V, with the bias correction of Bergsma (2013) available, and the small-expected-count condition reported rather than assumed away. Use it to say how strongly two categorical variables are associated (0 is none, 1 is complete), for example region by outcome in a cross-tabulation, on a scale that does not depend on the number of rows and columns.

Usage

cramerv(tbl, biascor = TRUE, minexp = 5, seed = 1L)

Arguments

Details

With \chi^2 the Pearson statistic (no continuity correction), n the total, r rows and k columns, \phi^2 = \chi^2/n and the uncorrected V = \sqrt{\phi^2 / \min(r-1, k-1)}. Bergsma's correction replaces \phi^2 with \max(0, \phi^2 - (r-1)(k-1)/(n-1)) and r, k with r - (r-1)^2/(n-1) and k - (k-1)^2/(n-1).

Empty rows and columns (an unused factor level) are dropped first, and method says so. When fewer than two rows or columns remain, every statistic is NA. When any expected count is below minexp, the p-value comes from stats::chisq.test(simulate.p.value = TRUE, B = 10000) under seed, and the caller's random number stream is restored afterwards.

Errors: tbl NULL, not numeric, or with a negative or non-finite count; biascor not TRUE or FALSE; minexp not a single non-negative number.

Value

A list with v (numeric, Cramer's V), chisq (the Pearson statistic), df (integer), p_value, n (the total count), minexp (the smallest expected count in the table), cellsbelow (integer), the number of cells whose expected count falls under the threshold. When any does, p_value comes from a Monte Carlo permutation instead of the chi-square approximation, and method says which was used.

References

Bergsma, W. (2013). A bias-correction for Cramer's V and Tschuprow's T. Journal of the Korean Statistical Society 42(3), 323-328. (Not in the local corpus; cited from the published paper.)

Examples

tbl <- rbind(c(120, 80), c(40, 160))
cramerv(tbl)

# No association gives a V near zero.
cramerv(rbind(c(100, 100), c(100, 100)))$v

# A sparse table's uncorrected V overstates the association.
sparse <- rbind(c(3, 1), c(1, 3))
c(raw = cramerv(sparse, biascor = FALSE)$v,
  corrected = cramerv(sparse)$v)

# small expected counts: a Monte Carlo p-value, and the method says so
cramerv(sparse)$method

# an unused level is dropped
cramerv(rbind(c(30, 10), c(0, 0), c(12, 28)))$method

# negative counts are refused
try(cramerv(rbind(c(-1, 2), c(1, 2))))

Rank correlation and midranks (C backend)

Description

cspr() is Spearman's rho: the Pearson correlation of the ranks, so it measures monotone association rather than linear association and is unaffected by any order-preserving transformation of either variable. cmidrk() exposes the ranks themselves; tied values share the average of the ranks they span, which is what makes the result agree with stats::cor() on tied data. Use Spearman's rho when two variables rise and fall together but not along a straight line (counts against rates, incomes against ranks), or when a few extreme values would dominate a Pearson correlation.

Usage

cspr(x, y)

cmidrk(x)

Arguments

Details

With r_i and s_i the midranks of x_i and y_i,

\rho = \frac{\sum_i (r_i - \bar r)(s_i - \bar s)} {\sqrt{\sum_i (r_i - \bar r)^2 \sum_i (s_i - \bar s)^2}},

the Pearson correlation of the ranks (Spearman, 1904). Because ties get midranks, the shortcut 1 - 6 \sum d_i^2 / (n(n^2 - 1)) is NOT what is computed; the full formula above is, and it stays correct with ties.

cspr() does not remove missing values: an NA in either vector gives NaN, and so does a constant vector (its ranks have no variance). cmidrk() leaves an NA in place and ranks the rest, as rank(x, na.last = "keep") does. Refusals: vectors of different lengths ("must have the same length") and a character, factor or list argument ("must be numeric"). Ranking sorts, so the cost is O(n \log n).

Value

cspr(): a length-1 numeric in [-1, 1] (NaN as described above). cmidrk(): a numeric vector the length of x, the rank of each value (1 for the smallest; tied values share the mean of their ranks; NA kept as NA).

References

Spearman C (1904). The proof and measurement of association between two things. The American Journal of Psychology 15(1), 72–101. DOI 10.2307/1412159

See Also

ccorr() for the Pearson correlation.

Examples

x <- c(1, 2, 3, 4, 5)
y <- c(2, 4, 9, 16, 25)

# Perfectly monotone but not linear: rho is 1 where Pearson is not.
cspr(x, y)
ccorr(x, y)

# Agrees with stats::cor(), ties included.
xt <- c(1, 2, 2, 2, 5, 5, 7)
yt <- c(3, 1, 1, 4, 4, 9, 2)
all.equal(cspr(xt, yt), stats::cor(xt, yt, method = "spearman"))

# Tied values share the average of the ranks they cover.
cmidrk(xt)
all.equal(cmidrk(xt), rank(xt))

# Invariant to any monotone rescaling.
all.equal(cspr(x, y), cspr(exp(x), log(y)))

# Missing values are kept in place by cmidrk().
cmidrk(c(3, NA, 1, 3))
try(cspr(1:3, 1:4))

Concentration-response functions for PM2.5 and NO2

Description

A concentration-response function (CRF) turns an air-pollution concentration into the relative risk (RR) of a health outcome compared with a counterfactual “clean air” level: an RR of 1.05 means 5% more deaths (or cases) than at the reference concentration. These are the first step of any burden estimate; polbrdn and vrfpol call them for you.

Usage

crfpm25(expo, y = "allcause", refconc = 5.8)

crfno2(expo, y = "allcause", refconc = 10, betaper10 = NULL)

Arguments

Details

crfpm25() is log-linear for all-cause mortality with the pooled cohort estimate of the WHO 2021 guideline review (Chen and Hoek 2020: RR 1.08, 95% CI 1.06-1.09, per 10 micrograms per cubic metre),

RR(z) = \exp(\ln(1.08) (z - z_{cf}) / 10)

for z > z_{cf} and 1 otherwise; for the cause-specific outcomes (ischaemic heart disease, stroke) it is the Integrated Exposure-Response curve of Burnett et al. (2014, equation 1),

RR(z) = 1 + \alpha (1 - e^{-\gamma (z - z_{cf})^{\delta}}),

with the GBD 2013 parameter triples; the IER was fitted per cause and has no all-cause form. crfno2() is the log-linear function of the WHO 2021 review (Huangfu and Atkinson 2020: RR 1.02, 95% CI 1.01-1.04, per 10 micrograms per cubic metre for all-cause mortality),

RR(z) = \exp(\beta (z - z_{cf}) / 10).

A vector of exposures gives the mean RR and log-RR. The per-unit relative risks are then in extra$rrperunit.

Every function has zero excess risk at and below the reference concentration (z \le z_{cf} gives RR 1). The published 95% limits of the log-linear coefficients are carried through the same curve into rrlo and rrup; the IER parameters come without an interval, so those are NA for "ihd" and "stroke", and also when betaper10 is supplied.

Refusals: an unknown y (the message lists the available ones); an exposure that is missing, infinite, negative, logical or empty; and a refconc (or betaper10) that is not a single finite number, or a negative refconc. A supplied betaper10 may be negative (a protective association) and replaces the table, in which case y is only a label. The functions are pure computations: no file, network or option is used.

Value

A list of class rmbl_crf (plot() draws the curve) with: rr and logrr, the relative risk and its logarithm at the exposure (their means over a vector of exposures); rrlo and rrup, the RR at the published 95% limits (NA where there are none); refconc; expconc (the exposure, or the mean of a vector); pollutant ("PM2.5" or "NO2"); citation (the source as a string); and extra, a list with the coefficient (betaper10 and betaper10ci for the log-linear forms, alpha, gamma, delta for the IER), form (PM2.5 only), outcome, and for a vector of exposures rrperunit, the RR at each one.

References

Chen, J. and Hoek, G. (2020). Long-term exposure to PM and all-cause and cause-specific mortality: a systematic review and meta-analysis. Environment International 143, 105974. Huangfu, P. and Atkinson, R. (2020). Long-term exposure to NO2 and O3 and all-cause and respiratory mortality: a systematic review and meta-analysis. Environment International 144, 105998. Khreis, H. et al. (2017). Exposure to traffic-related air pollution and risk of development of childhood asthma: a systematic review and meta-analysis. Environment International 100, 1-31. GBD 2019 Risk Factors Collaborators (2020). The Lancet 396(10258), 1223-1249. Burnett, R. T. et al. (2014). An integrated risk function for estimating the global burden of disease attributable to ambient fine particulate matter exposure. Environmental Health Perspectives 122(4), 397-403. WHO (2021). Global Air Quality Guidelines.

Examples

# All-cause mortality at an annual PM2.5 of 12 micrograms per cubic metre.
r <- crfpm25(12)
c(rr = r$rr, lower = r$rrlo, upper = r$rrup)

# The cause-specific IER curve for ischaemic heart disease.
crfpm25(12, y = "ihd")$rr

# NO2 and respiratory deaths; at or below the reference there is no
# excess risk.
crfno2(25, y = "respiratory")$rr
crfno2(8)$rr

# A vector of exposures (one per area): the mean RR and each area's RR.
crfpm25(c(5, 12, 20))$extra$rrperunit

# Your own coefficient: RR 1.04 per 10 micrograms per cubic metre.
crfno2(30, betaper10 = log(1.04))$rr

# An outcome without a published function is refused.
try(crfpm25(12, y = "asthma"))

Where incidents cluster: a kernel density surface

Description

Estimates the density of incident locations – police-reported occurrences, calls for service, addresses of people in a program – as a smooth surface in incidents per square kilometre, and the smallest areas holding a given share of them (the highest-density regions). The plot draws the surface as a heat map with those regions outlined.

Usage

crmdens(
  lon,
  lat,
  w = NULL,
  bw = NULL,
  ngrid = 200L,
  cover = c(0.5, 0.9),
  xy = c("lonlat", "planar")
)

Arguments

Details

The estimate is the Gaussian product-kernel density, computed on the grid by linear binning followed by an exact discrete convolution, which agrees with the direct sum to within the binning error (of order the squared grid spacing over the squared bandwidth). The grid extends three bandwidths past the outermost incident so the surface integrates to the total weight.

A highest-density region for share p is the set of grid cells with the largest densities whose estimated incidents add up to p of the total. Its area is the spatial counterpart of a concentration index: when half the incidents fall in a few per cent of the mapped area, they are concentrated. The bandwidth decides how sharp the surface is, so it should be reported with any map; the default suits a first look, not a definitive one.

Errors: lon and lat of different lengths, or (with coords = "lonlat") a coordinate that is not a longitude and latitude in degrees; fewer than two located incidents with positive weight; ngrid not a whole number from 10 to 2000; cover not shares strictly between 0 and 1; incidents all on one line with no bw given; a bw that is not one or two positive numbers, or longer than half the earth's circumference. Incidents with a missing coordinate are dropped.

Value

A list of class rmbl_density: field (incidents per square kilometre on the grid), x and y (grid coordinates in kilometres from the centre), lon and lat (the same grid in degrees, for "lonlat"), bandwidth, n, dropped, points (the projected incidents), and regions: one row per coverage share with the density level that bounds it, its areakm2, and its share of the mapped area.

References

Scott, D. W. (2015). Multivariate Density Estimation: Theory, Practice, and Visualization, 2nd ed. Wiley, Hoboken, NJ.

Wand, M. P. (1994). Fast computation of multivariate kernel estimators. Journal of Computational and Graphical Statistics, 3(4), 433-445.

Hyndman, R. J. (1996). Computing and graphing highest density regions. The American Statistician, 50(2), 120-126.

Chainey, S., Tompson, L. and Uhlig, S. (2008). The utility of hotspot mapping for predicting spatial patterns of crime. Security Journal, 21(1-2), 4-28.

See Also

knoxtst, plmfld (the same kind of surface, for an air pollutant).

Examples

b <- read.csv(system.file("extdata", "tps_break_enter_d14_2024.csv",
                          package = "rmoriebricklayer"))
kd <- crmdens(b$lon, b$lat, ngrid = 120)
kd
plot(kd)

# a fixed 300 m bandwidth on a coarser grid: a sharper surface
kd3 <- crmdens(b$lon, b$lat, bw = 0.3, ngrid = 60)
kd3$regions

# planar coordinates in metres
set.seed(2)
xy <- cbind(rnorm(200, sd = 500), rnorm(200, sd = 500))
crmdens(xy[, 1], xy[, 2], xy = "planar", ngrid = 40)$bandwidth

# longitudes and latitudes are checked
try(crmdens(c(-79.4, 200), c(43.7, 43.8)))

Quantiles, median and robust spread (C backend)

Description

cqnt() is the type-7 quantile, which is R's default, so it agrees with stats::quantile(x, probs, type = 7). cmed() is the 50% point. cmad() is the median absolute deviation, scaled by constant so that it estimates the standard deviation of a normal sample. ciqr() is the interquartile range and ctkyfnc() the outlier fences drawn at k IQRs beyond the quartiles.

Usage

cqnt(x, probs = c(0, 0.25, 0.5, 0.75, 1))

cmed(x)

cmad(x, const = 1.4826)

ciqr(x)

ctkyfnc(x, k = 1.5)

Arguments

Details

These are the robust counterparts of cmom(): a single corrupted row can move a mean or a variance arbitrarily far, but moves a median or a MAD hardly at all – which is what you want when deciding whether a freshly fetched column is still the column a capsule was pinned against, or when screening a column for values to inspect by hand.

The type-7 quantile at probability p interpolates linearly between the order statistics x_{(j)} and x_{(j+1)} with h = (n - 1)p + 1, j = \lfloor h \rfloor: Q(p) = x_{(j)} + (h - j)(x_{(j+1)} - x_{(j)}) (Hyndman and Fan, 1996). The MAD is c \cdot \mathrm{median}_i |x_i - \mathrm{median}(x)| with c = constant; the IQR is Q(0.75) - Q(0.25); and the fences are Q(0.25) - k \cdot \mathrm{IQR} and Q(0.75) + k \cdot \mathrm{IQR} (Tukey, 1977).

Missing values are not removed: an NA or NaN anywhere in x makes every result NaN (unlike stats::quantile(), which stops). An empty x gives NaN. cqnt() refuses an empty probs ("probs is empty") and any probability outside [0, 1] or missing ("must lie in [0, 1]"); every function refuses a character, factor or list x ("must be numeric"). The quantiles sort a copy of x, so the cost is O(n \log n).

Value

cqnt(): a numeric vector the length of probs, in the units of x, named by percentage ("25%", ...) as stats::quantile() names them. cmed(), cmad() and ciqr(): a length-1 numeric in the units of x. ctkyfnc(): a named length-2 numeric, lower and upper.

References

Hyndman RJ, Fan Y (1996). Sample quantiles in statistical packages. The American Statistician 50(4), 361–365. DOI 10.1080/00031305.1996.10473566

Tukey JW (1977). Exploratory Data Analysis. Addison-Wesley.

See Also

cmom() for the non-robust moments, ctrmmean() for a robust centre that still averages.

Examples

x <- c(2, 4, 4, 4, 5, 5, 7, 9)
cqnt(x, c(0.25, 0.5, 0.75))
all.equal(cqnt(x, c(0.1, 0.9)),
          as.numeric(stats::quantile(x, c(0.1, 0.9))))

cmed(x)
cmad(x)
all.equal(cmad(x), stats::mad(x))
cmad(x, const = 1)                  # unscaled median deviation
cmad(x, const = 1 / stats::qnorm(3/4))  # unrounded constant

ciqr(x)
ctkyfnc(x)

# Robustness: one wild value barely moves the median, but moves the
# mean a long way.
wild <- c(x, 1000)
c(mean = mean(wild), median = cmed(wild))

# Values outside the fences are the candidates to inspect.
f <- ctkyfnc(wild)
wild[wild < f[["lower"]] | wild > f[["upper"]]]

# A missing value is not dropped: the answer is NaN.
cmed(c(1, NA, 3))

# A probability outside [0, 1] is refused.
try(cqnt(x, 1.5))

SHA-256 hex digest (C backend)

Description

Computes the SHA-256 hash (FIPS 180-4) of text or raw bytes with the self-contained implementation in the package's compiled core, and returns it as a 64-character lowercase hexadecimal string. A hash is a short fingerprint: any change to the input changes it, so it is used throughout the package to record exactly which data, mapping or file a result came from. Sibling packages use the same routine through LinkingTo: rmoriebricklayer.

Usage

csha256(x)

Arguments

Details

For a raw vector the bytes are hashed as one message and one digest is returned. Any other input is converted with as.character() and each element is hashed separately, as the bytes of the string as stored in memory (no re-encoding). In a UTF-8 session those are the UTF-8 bytes; for a string marked latin1 they are the latin1 bytes, so use enc2utf8() first when the digest must be the same on every platform. NA gives NA. Numbers are hashed through their printed form (1 hashes as "1"), and a factor through its labels.

For a digest of a whole R object, hash its serialization (see the examples) or use dgstobj; for BLAKE2b use cblake2b.

Value

A character vector of 64-character lowercase hex digests: one per element for character input, length 1 for raw input (character(0) for an empty character vector).

References

National Institute of Standards and Technology (2015). Secure Hash Standard (SHS). FIPS PUB 180-4. DOI 10.6028/NIST.FIPS.180-4

Examples

# Hash a character scalar (NIST test vector for "abc").
csha256("abc")
# ba7816bf8f01cfea414140de5dae2223b00361a396177a9cb410ff61f20015ad

# Vectorised over character input: one digest per element.
csha256(c("abc", "def", NA))

# Raw input hashes the bytes directly; identical to the character form.
identical(csha256("abc"), csha256(charToRaw("abc")))

# The empty message.
csha256(raw(0))

# Fingerprint an arbitrary object via its serialization.
csha256(serialize(list(a = 1L, b = "x"), NULL))

SHA-512 hex digest (C backend)

Description

Hashes character or raw input with the self-contained SHA-512 (FIPS 180-4) in the compiled core and returns the digest as hex. A digest is a fingerprint: any change to the input gives a different value, so it pins a data file, a code string or a manifest. Use it over csha256() when a pin has to outlive the capsule by decades: the wider digest leaves more margin, including against a quantum adversary, for whom Grover's algorithm halves the effective preimage exponent.

Usage

csha512(x)

Arguments

Details

Character input is hashed element by element, each string as its bytes (UTF-8 text gives the UTF-8 bytes; no newline is added). NA elements give NA. Other atomic input (numbers, logicals) is first converted with as.character(), so csha512(1) is the digest of the text "1". A raw vector is hashed as one byte string. For files use sha512file(). Pure computation; no side effects.

Value

A character vector of 128-character lowercase hex digests: one per element for character input (NA for NA), length 1 for raw input.

See Also

csha256(), ccrc32(), sha512file()

Examples

# FIPS 180-4 test vector for "abc".
csha512("abc")

# Vectorised over character input; NA stays NA.
csha512(c("abc", "def", NA))

# Raw input hashes the bytes directly, and agrees with the character
# form for the same bytes.
identical(csha512("abc"), csha512(charToRaw("abc")))

# Twice the digest width of SHA-256.
c(sha256 = nchar(csha256("abc")), sha512 = nchar(csha512("abc")))

# numbers are hashed as their text
identical(csha512(1), csha512("1"))

Standard deviation and Euclidean distance (C backend)

Description

csdev() returns the standard deviation of a numeric vector and cdistm() the straight-line (Euclidean) distance between two vectors of the same length. Both run in the package's compiled core, the same code that rmorie and rmoriedata reach through LinkingTo: rmoriebricklayer, so a figure computed here is the figure those packages compute. Use them when you want the shared core's answer rather than base R's, for example when checking a capsule's recorded summary statistics against a fresh computation.

Usage

csdev(x, ddof = 1L)

cdistm(a, b)

Arguments

Details

csdev() computes

s = \sqrt{\frac{1}{n - \mathrm{ddof}} \sum_{i=1}^n (x_i - \bar x)^2},

so ddof = 1 is the sample standard deviation of stats::sd() and ddof = 0 the population figure. cdistm() computes \sqrt{\sum_i (a_i - b_i)^2}.

Missing values propagate: there is no na.rm, and one NA makes the result NA. Remove them first (for example with stats::na.omit()) when that is what you want. A vector with no more than ddof values has no defined spread and gives NaN (one value with the default ddof = 1, or an empty vector). Two empty vectors are at distance 0.

Refusals: a character, factor, list or data-frame x, a or b ("must be numeric"), and for cdistm() vectors of different lengths ("must have the same length"). Logical input is accepted and counted as 0/1. Both run in one pass over the data and allocate nothing beyond the coerced input.

Value

csdev(): a length-1 numeric, the standard deviation in the units of x (NA when x holds a missing value, NaN when length(x) <= ddof). cdistm(): a length-1 non-negative numeric, the distance in the units of a and b.

See Also

cmom() for the mean, variance, skewness and kurtosis in a single pass; cmad() for a spread that one wild value cannot move.

Examples

# Sample standard deviation, matching stats::sd().
csdev(c(2, 4, 4, 4, 5, 5, 7, 9))
all.equal(csdev(1:10), stats::sd(1:10))

# ddof = 0 divides by n instead of n - 1.
csdev(1:10, ddof = 0)
all.equal(csdev(1:10, ddof = 0), sqrt(mean((1:10 - mean(1:10))^2)))

# Missing values propagate; remove them first if they should be ignored.
csdev(c(3, NA, 5))
csdev(stats::na.omit(c(3, NA, 5)))

# Euclidean distance between two points.
cdistm(c(0, 0), c(3, 4))     # 5
cdistm(1:5, 1:5)             # 0 -- a point is zero from itself

# Vectors of different lengths are refused rather than recycled.
try(cdistm(1:3, 1:2))

Fast summary statistics (C backend)

Description

cmn(), cvari() and ccorr() compute the mean, the sample variance and the Pearson correlation in compiled code. They are thin R wrappers over the kernels in the package's C/C++ core, the same ones that other packages call through LinkingTo: rmoriebricklayer, so an R result and a compiled result agree exactly. Use them where you want those shared kernels; for everyday work mean(), var() and cor() give the same numbers.

Usage

cmn(x)

cvari(x)

ccorr(x, y)

Arguments

Details

With n values,

\bar x = \frac{1}{n}\sum_i x_i, \qquad s^2 = \frac{1}{n - 1}\sum_i (x_i - \bar x)^2, \qquad r = \frac{\sum_i (x_i - \bar x)(y_i - \bar y)} {\sqrt{\sum_i (x_i - \bar x)^2 \sum_i (y_i - \bar y)^2}}.

The variance uses the n - 1 (sample) denominator, as var does. The mean is accumulated so that a sum of finite values that would overflow a double (for example two values of 1e308) still gives the finite mean.

There is no na.rm: any NA or NaN gives NA or NaN; remove missing values first, for example with na.omit. Edge cases return NaN rather than an error: an empty vector for the mean, fewer than two values for the variance, and a constant vector for the correlation.

Errors: input that is a character vector, factor, list or data frame (logical input is accepted and read as 0/1); x and y of different lengths in ccorr().

Value

cmn(), cvari() and ccorr() each return a length-1 numeric.

See Also

cnormpdf, csha256.

Examples

## cmn(): sample mean (NA/NaN propagate; no na.rm)
cmn(1:10)                 # 5.5
cmn(c(1, 2, NA))          # NA -- call stats::na.omit() first if needed
cmn(stats::na.omit(c(1, 2, NA)))
cmn(c(1e308, 1e308))      # no overflow

## cvari(): n-1 (sample) variance, matching stats::var()
cvari(c(2, 4, 4, 4, 5, 5, 7, 9))
all.equal(cvari(1:10), stats::var(1:10))   # agrees with base R

## ccorr(): Pearson correlation of two equal-length vectors
ccorr(1:10, (1:10)^2)        # strong positive, near 0.97
ccorr(1:10, 10:1)            # perfect negative: -1

## refusals
try(cmn(c("1", "2")))     # character input is not coerced
try(ccorr(1:3, 1:4))         # lengths differ

Trimmed and winsorized means (C backend)

Description

Two ways to stop a handful of extreme rows dominating a column's centre. ctrmmean() DISCARDS the floor(n * trim) largest and smallest values, as mean(x, trim = ) does. cwinsmean() instead PULLS THEM IN to the most extreme surviving values, so every observation still contributes weight – usually the better choice when the extremes are real measurements rather than errors.

Usage

ctrmmean(x, trim = 0.1)

cwinsmean(x, trim = 0.1)

Arguments

Details

With g = \lfloor n \cdot \mathrm{trim} \rfloor and the sorted values x_{(1)} \le \dots \le x_{(n)}, the trimmed mean is \frac{1}{n - 2g} \sum_{i=g+1}^{n-g} x_{(i)} and the winsorized mean replaces x_{(1)}, \dots, x_{(g)} by x_{(g+1)} and x_{(n-g+1)}, \dots, x_{(n)} by x_{(n-g)} before averaging all n values. At trim = 0.5 both return the median, as mean() does.

A trim that is not a single number in [0, 0.5] is an error. An empty x, or any NA/NaN in x, gives NaN (missing values are not removed). A character, factor or list x is refused ("must be numeric"). Both sort a copy of x, so the cost is O(n \log n).

Value

A length-1 numeric in the units of x.

See Also

cmed() and cmad() for the other robust summaries; cmom() for the ordinary mean.

Examples

x <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 100)

mean(x)                               # dragged up by the 100
ctrmmean(x, 0.1)             # the 100 and the 1 dropped
cwinsmean(x, 0.1)          # the 100 pulled back to 9

# Agrees with base R's own trimming.
all.equal(ctrmmean(x, 0.2), mean(x, trim = 0.2))

# trim = 0 is the plain mean; trim = 0.5 is the median.
all.equal(ctrmmean(x, 0), mean(x))
all.equal(ctrmmean(x, 0.5), cmed(x))

# A trim past one half is not a trim.
try(cwinsmean(x, 0.7))
try(ctrmmean(letters))

Reproducible uniforms shared with morie's Python arm (splitmix64)

Description

Returns n pseudo-random numbers on [0, 1) from a 64-bit seed with the splitmix64 generator, the same numbers morie's Python computes for the same seed. Use it where the R and Python arms must draw identical "random" values so that their results agree: the within-day jitter of tied event dates in the TPS Hawkes fits, or a subsample.

Usage

cunif(n, seed)

Arguments

Details

splitmix64 (Steele, Lea & Flood 2014) advances a 64-bit state by the constant 0x9E3779B97F4A7C15 and mixes it with two multiply-xorshift rounds; the top 53 bits of each output divided by 2^{53} give a double in [0, 1). The stream depends only on seed: R's own generator and set.seed() are neither used nor changed. It is not a statistical replacement for R's generators.

Errors: n not one non-negative whole number; seed not one whole number in [0, 2^{53}).

Value

A numeric vector of length n, values in [0, 1).

References

Steele GL, Lea D, Flood CH (2014). Fast splittable pseudorandom number generators. OOPSLA 2014, 453–472. DOI 10.1145/2660193.2660195

Examples

cunif(3, 42)
identical(cunif(5, 1), cunif(5, 1))

# a stream is a prefix of a longer one with the same seed
all.equal(cunif(3, 7), cunif(10, 7)[1:3])

# R's own generator is left alone
set.seed(1); a <- runif(1)
set.seed(1); invisible(cunif(100, 9)); b <- runif(1)
identical(a, b)

# negative counts and seeds past 2^53 are refused
try(cunif(-1, 1))
try(cunif(2, 2^53))

The Community Well-Being index of a community

Description

Indigenous Services Canada's Community Well-Being (CWB) index: income, education, housing and labour force activity, each scored from 0 to 100, and their mean, from census profile figures.

Usage

cwbscr(
  incpercap,
  highsch20p,
  uni25p,
  notcro,
  nomajrep,
  partic,
  empl,
  bounds = c(2650, 75000)
)

Arguments

Details

The components, as the index's methodology sets them out (Penney, O'Sullivan and Senecal 2012; O'Sullivan and McHardy 2013; Indigenous Services Canada's reports on trends, 1981 to 2016 and 2021):

income

100 (\log y - \log a) / (\log b - \log a) for income per capita y (total income over total population) between the bounds a and b (2,650 and 75,000 dollars in the 2016 and 2021 series; 2,000 and 40,000 in the earlier ones), kept between 0 and 100 here.

education

two thirds of the percentage aged 20 and over with at least a high school certificate, plus one third of the percentage aged 25 and over with a university degree at the bachelor's level or above.

housing

half the percentage of people in dwellings with no more than one person per room, plus half the percentage in dwellings that need no major repairs.

labour

half the participation rate, plus half the employment rate of those participating (ages 20 to 64).

The index is the mean of the four. Indigenous Services Canada computes it from unrounded census microdata, scores only communities of 65 people or more and suppresses the components of communities under 250 people or 40 households: scores from published, rounded profile figures come close to the published index but need not equal it, and a small community's should be treated as the Department would. The full current methodology is available from the Department on request. The index measures socio-economic conditions as the census counts them; it says nothing of culture, language or relationship to land, and is not a measure of wellness as First Nations and Inuit define it. The CWB is Indigenous Services Canada's index, with bounds in the Canadian dollars of its series: for another country, set bounds in its currency and read the result as a CWB-style score, not comparable with Canada's.

Value

A data frame, one row per community: income, education, housing, labour, cwb.

References

Penney C, O'Sullivan E, Senecal S (2012). The Community Well-Being Index (CWB): Examining Well-Being in Inuit Communities, 1981-2006. Aboriginal Affairs and Northern Development Canada.

O'Sullivan E, McHardy M (2013). The Community Well-Being (CWB) Index: Disparity in Well-Being Between First Nations and Other Canadian Communities Over Time. In White JP, Peters J, Beavon D, Spence N (eds), Aboriginal Well-Being: Canada's Continuing Challenge. Thompson Educational.

Indigenous Services Canada (2019). Report on trends in First Nations communities, 1981 to 2016 (Community Well-Being index).
sac-isc.gc.ca/eng/1345816651029/1557323327644

Examples

# made-up profile figures for two communities, to show the computation
cwbscr(incpercap = c(18000, 31000), highsch20p = c(58, 81),
          uni25p = c(6, 19), notcro = c(71, 95),
          nomajrep = c(62, 90),
          partic = c(55, 68), empl = c(78, 92))

Weighted mean and variance (C backend)

Description

The weighted mean and weighted variance of a numeric vector, for data where some observations count for more than others: an average over areas weighted by population, or measurements weighted by how precise each one is. Both come from the compiled core in one pass.

Usage

cwtd(x, w)

Arguments

Details

The mean is \bar x_w = \sum_i w_i x_i / \sum_i w_i. The weights are treated as RELIABILITY weights (how precisely each observation is known), so the variance is

s_w^2 = \frac{\sum_i w_i (x_i - \bar x_w)^2} {\sum_i w_i - \sum_i w_i^2 / \sum_i w_i},

the unbiased estimator under that reading. With every weight equal to 1 it reduces exactly to stats::var(). Frequency weights (a weight of 3 meaning three identical rows) call for the denominator \sum w - 1 instead, and survey design weights call for a design-based variance, so do not use this variance for those.

Refusals: x and w of different lengths ("must have the same length"), any negative weight ("must be non-negative"), and a character, factor or list argument ("must be numeric"). All-zero weights, or an NA in either vector, give NaN for both entries.

Value

A named length-2 numeric vector: mean (units of x) and variance (squared units of x).

See Also

cmom() for the unweighted moments, cipw() for inverse-probability weights to feed in here.

Examples

x <- c(10, 20, 30, 40)
w <- c(1, 1, 2, 4)
cwtd(x, w)

# The mean agrees with base R.
all.equal(cwtd(x, w)[["mean"]], stats::weighted.mean(x, w))

# Equal weights recover the unweighted variance.
all.equal(cwtd(x, rep(1, 4))[["variance"]], stats::var(x))

# A single dominant weight pulls the mean onto that observation.
cwtd(x, c(1, 1, 1, 1000))[["mean"]]

# A negative weight is refused.
try(cwtd(x, c(1, -1, 1, 1)))

Decode numeric category codes with an explicit code-to-label dictionary

Description

Turns a vector of category codes (for example 1, 2, 3 from an SPSS, Stata or SAS file) into the category labels they stand for, using a code-to-label dictionary that you write out. Every code that occurs in the data must be in the dictionary; a code without a label is an error rather than a silent NA or a pass-through. Use it whenever a column arrives as numbers that are really categories, before any table or model, so that a group can never be renumbered by treating codes (or as.numeric() on a factor of them) as the categories.

Usage

deccodes(x, vlabs, keepna = TRUE)

Arguments

Details

x is converted to character (a factor by its labels, a number by as.character(), so 1 and 1.0 both become "1") and every distinct non-missing code is looked up by name in vlabs. The lookup itself is done by grdrcd, so the result carries the same recode_audit attribute, including a SHA-256 of the dictionary.

Errors:

• vlabs not a character vector with non-empty names;

• a code in x that is not a name of vlabs (the message lists the unmapped codes);

• any NA in x when keepna = FALSE.

For a haven-style labelled vector, which carries its dictionary as a labels attribute, use declbld, which returns a factor with the levels in code order.

Value

A character vector of labels, the same length as x, with NA where x is NA. It carries a recode_audit attribute: a list with mapping (the dictionary), kept (character(0)) and checksum (64-character hexadecimal SHA-256 of the code=label pairs).

See Also

declbld, grdrcd, audcats.

Examples

# codes from a survey file with their code book
race <- c(1, 2, 2, 3, NA)
book <- c("1" = "White", "2" = "Black", "3" = "Indigenous")
deccodes(race, book)

# the dictionary that was applied is kept on the result
attr(deccodes(race, book), "recode_audit")$checksum

# a code with no label is refused, never passed through
try(deccodes(c(1, 2, 9), book))

# missing codes can be refused too
try(deccodes(race, book, keepna = FALSE))

Decode a labelled import by its value labels, by code

Description

Converts a vector read from SPSS, Stata or SAS (for example by haven::read_sav()), which stores categories as numeric codes with a labels attribute, into a factor of the labels. Each code is looked up by its value, so neither the order of the codes nor the alphabetical order of the labels can affect which label a row gets. Codes with no label are refused. The factor's levels are in code order, which keeps the source program's reference category.

Usage

declbld(x, vlabs = NULL)

Arguments

Details

When x has a labels attribute (a named vector, label = code) it is used and vlabs is ignored; otherwise vlabs (code = label) is required. Decoding is done by deccodes. The levels are the labels ordered by their numeric code (codes that are not numbers are ordered after the numeric ones, by string), and include labels of codes that do not occur.

Errors: a labels attribute without names; no attribute and no vlabs; a code with no label.

Value

A factor whose levels are the labels in code order (not alphabetical), with a recode_audit attribute as described in grdrcd.

See Also

deccodes, relbl, xfrvrf

Examples

x <- structure(c(1, 2, 2, 4),
               labels = c(White = 1, Black = 2, Other = 3, Unknown = 4))
declbld(x)

# plain codes with a code book; levels follow the codes 1, 2, 10
declbld(c(2, 1, 10), vlabs = c("10" = "Ten", "2" = "Two", "1" = "One"))

# codes alone are not categories
try(declbld(c(1, 2, 2)))

Digest an arbitrary R object

Description

Fingerprints any R object by hashing its serialization, so a list, a data frame, a function or a fitted model all get a stable digest. The counterpart of digest::digest(), computed with this package's own hashes.

Serialization is pinned to version 2 with XDR byte order, so the digest is the same on a big-endian machine as on a little-endian one. Two objects that are identical() produce the same digest; two that merely print the same need not, because attributes are part of the serialization.

Usage

dgstobj(x, algo = c("sha256", "sha512", "blake2b", "crc32"), key = NULL)

Arguments

Details

The object is serialized with serialize(x, NULL, version = 2, xdr = TRUE) and the bytes are hashed with csha256(), csha512(), cblake2b() (32 bytes) or ccrc32(). CRC-32 is a checksum for accidental corruption, not a cryptographic hash. With a key the digest is HMAC-SHA-256 for "sha256" and keyed BLAKE2b for "blake2b"; a key with "sha512" or "crc32" is an error rather than a silently unkeyed digest. Objects that contain environments (functions, many fitted models) serialize their enclosing environments too, so their digest can change between sessions even when the code does not.

Value

A length-1 character vector of lowercase hex (64 characters for "sha256" and "blake2b", 128 for "sha512"), or a length-1 numeric for "crc32".

See Also

csha256() for hashing text or bytes directly, jsser() for a readable form that keeps types and attributes.

Examples

dgstobj(list(a = 1L, b = "x"))

# Stable across calls, and sensitive to any change.
identical(dgstobj(1:10), dgstobj(1:10))
dgstobj(1:10) == dgstobj(1:11)

# Attributes are part of the object, so they are part of the digest.
dgstobj(matrix(1:6, nrow = 2)) == dgstobj(1:6)

# Any of the hashes, and a keyed fingerprint.
dgstobj(mtcars, algo = "sha512")
dgstobj(mtcars, algo = "crc32")
dgstobj(mtcars, algo = "blake2b", key = "secret")

# A data frame's digest pins the data, so it can go in a manifest.
dgstobj(data.frame(x = 1:3))

# a key with an algorithm that has no keyed form is refused
try(dgstobj(1:3, algo = "crc32", key = "secret"))

The mix of outcomes over time, with simultaneous intervals and trend tests

Description

For a table of counts – periods in rows, outcomes in columns, such as cases withdrawn, pleaded guilty, or tried – gives each outcome's share of its period with Goodman's simultaneous confidence intervals, tests whether the mix is the same in every period, and tests each outcome for a linear trend in its share.

Usage

dispmix(cnts, prd = NULL, scores = NULL, conflvl = 0.95)

Arguments