Package 'crimlab'

Description

Spreads each event known only to have happened in an interval ⁠[from, to)⁠ with uniform density over that interval and counts each instant in the bin of its local clock time (Ratcliffe 2002): with the defaults, the 168 hours of the week from Sunday 00:00. A change of clock (daylight saving) moves the remaining time to the right hours: skipped hours get nothing and a repeated hour gets both. Missing values follow stated rules: an event with no from is left out; with no to, a reversed interval or one shorter than lo seconds, its whole weight goes to the start's bin; one of at least hi seconds is spread evenly over all bins.

Usage

aorwt(
  from,
  to,
  bw = 3600,
  nb = 168,
  origin = 259200,
  off = 0,
  tt = NULL,
  toff = NULL,
  lo = 120,
  hi = 604800,
  each = FALSE
)

Arguments

Value

A list: tot (the summed weights by bin), w (NULL, or a matrix with one row per event and one column per bin, rows summing to 1, 0 for events left out) and cnt (events left out, without end, reversed, short and long).

References

Ratcliffe, J. H. (2002). Aoristic signatures and the spatio-temporal analysis of high volume crime patterns. Journal of Quantitative Criminology, 18(1), 23-43. doi:10.1023/A:1013240828824

Examples

f <- as.POSIXct(c("2024-03-04 22:00", "2024-03-05 09:30"), tz = "UTC")
t <- as.POSIXct(c("2024-03-05 02:00", "2024-03-05 10:30"), tz = "UTC")
a <- aorwt(f, t, bw = 3600, nb = 24, origin = 0)
a$tot

Description

bdig() takes the first nd significant digits of each value. With dp = -1 they are the digits of the shortest decimal that gives back the double, decided exactly, so 0.29 gives 29 where truncating x times a power of ten gives 28; with dp from 0 to 22 the data are taken as recorded to dp decimals. bnft() tests the digits against Benford's law, P(d) = log10(1 + 1 / d), by the chi-square statistic with its exact upper tail, a Monte Carlo p-value from samples of Benford digits of the same size, and Nigrini's mean absolute deviation. A departure from Benford's law is a reason to look, not a finding of fabrication: many honest data sets (bounded amounts, assigned numbers) do not follow it.

Usage

bdig(x, nd = 1, dp = -1, sgn = "pos")

bnft(dig, nd = 1, nsim = 0, seed = NULL)

Arguments

Value

bdig(): a list of dig (an integer per value, 0 when set aside) and cnt (zeros, non-finite, excluded sign, out of range). bnft(): a list of obs (a data frame of digit, observed count and Benford expected count) and stat (m digits used, chisq, df, p, pmc (Monte Carlo, NaN without samples) and mad).

References

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

Nigrini, M. J. (2012). Benford's Law: Applications for Forensic Accounting, Auditing, and Fraud Detection. Wiley. doi:10.1002/9781119203094

Examples

set.seed(1)
x <- exp(runif(500, 0, 10))
d <- bdig(x)$dig
bnft(d, nsim = 199, seed = 1)$stat

Description

Combines a journey-to-crime surface over origin zones with an origin-destination matrix of solved offences (Levine 2009; Levine and Block 2011): P(Jtc), the surface; P(O), the origins of all solved offences; P(O|Jtc), the origins of offences in the zones where the series happened; their product, Bayes' numerator with the surface as prior; and the risk, the product over P(O), renormalised. Each block sums to 1. The matrix describes past solved offences and inherits their selection.

Usage

bjtc(pj, od, dest)

Arguments

Value

A data frame with one row per origin zone: pj, po, pod, prod and risk (0 where po is 0).

References

Levine, N. and Block, R. (2011). Bayesian journey-to-crime estimation: an improvement in geographic profiling methodology. The Professional Geographer, 63(2), 213-229. doi:10.1080/00330124.2010.547152

Examples

od <- matrix(c(5, 1, 0, 2, 4, 1, 0, 1, 6), 3)
bjtc(pj = c(0.5, 0.3, 0.2), od = od, dest = c(1, 2, 2))

Description

Does the order in which an algorithm ranks areas by predicted risk match the order of what was later recorded there, and does the disagreement follow the areas' demographic make-up? calaud() ranks the areas both ways, measures the agreement by Spearman's correlation, and gives each group the mean gap between its areas' two ranks, with a permutation test. It generalises the district-level analysis of the SciencesPo predictive-policing study of Chicago (Lacherade et al. 2021) to any place: the five columns of the audit are area, risk, outcome, group and population.

Usage

calaud(
  area,
  risk = NULL,
  outcome = NULL,
  group = NULL,
  population = NULL,
  nperm = 999,
  seed = NULL
)

Arguments

Details

Records are rolled up to areas first: an area's risk is the mean over its records, its rate the mean outcome (or, with population, the summed outcome per head), and its group the most frequent label among its records (the first in order on ties). One record per area is the per-area audit itself. Rank 1 is the highest; tied values share the mean of their ranks. An area's gap is its outcome rank minus its risk rank: positive when the algorithm ranks it as riskier than its recorded outcome warrants, negative when less risky.

Spearman's rho is Pearson's correlation of the two rankings (Spearman 1904). pex is its two-sided p-value by algorithm AS 89 (Best and Roberts 1975): exact for nine areas or fewer, from the Edgeworth series up to 1289, and NA when either ranking has ties; pt is the large-sample one from the t distribution, as stats::cor.test() gives them. A group's mean gap is tested against shuffling the groups across the areas, nperm times; between tests all groups at once by the between-group sum of squares of the gaps (the gaps sum to 0).

The audit compares the algorithm with the records, not with crime: a recorded outcome follows where the police looked as well as what happened, so a group's areas can look over-predicted because their incidents were recorded less, or under-predicted because they were recorded more (see fairsim()). A gap says the two orders disagree for that group's areas; it does not say which order is right or why.

Value

A list: areas, a data frame with one row per area (sorted by label): area, n (records), risk, rate, group, riskrk, raterk and gap; groups, one row per group: group, n (areas), gap (mean) and p (NaN without areas or permutations); spear, a named vector: rho, s (Spearman's statistic, the sum of squared rank differences), pex, t and pt; ties (TRUE when either ranking has ties); between, ss and p; and nperm.

References

Spearman, C. (1904). The proof and measurement of association between two things. American Journal of Psychology, 15(1), 72-101. doi:10.2307/1412159

Best, D. J. and Roberts, D. E. (1975). Algorithm AS 89: the upper tail probabilities of Spearman's rho. Journal of the Royal Statistical Society, Series C, 24(3), 377-379. doi:10.2307/2347111

Lacherade, P., Szabo, B., Krikava, J. and Aeby, M. (2021). Predictive policing in the City of Chicago. Group project, Decoding Biases in AI, Sciences Po. https://github.com/marquito95/Predictive-policing-Chicago

Examples

set.seed(3)
d <- fairsim(4000, nareas = 24, bias = 0.5, seed = 1)
# the recorded outcome is `det`: the incidents that were detected
r <- calaud(d$area, d$risk, d$det, d$group, nperm = 499, seed = 1)
r$groups
r$spear

Description

Reports, for each variable, what the functions of this package would make of it: its storage type, its levels and their counts, missing values, unused levels, levels counted fewer than nmin times, labels that differ only by case or white space, a 0/1 or logical coding, integer codes that look like categories, value labels (as haven makes) that do not match the values, and which level a default would take as the reference, and why. Nothing is changed.

Usage

catchk(x, ref = NULL, nmin = 5)

Arguments

Details

Levels are those of a factor, in level order; FALSE and TRUE for a logical, and 0 and 1 for numbers that are all 0 or 1, both always listed; otherwise the distinct values as text, in order of first appearance, and for a labelled vector then the labelled values that do not occur. Doubles that are not all whole numbers, dates, times, lists and matrices are not categorical: they are reported with their type and missing values only.

The default reference: FALSE for a logical and 0 for 0/1 numbers (sntbnd() takes TRUE and 1 as the treatment), the first level of a factor (ref = 1 on a table of it), and none otherwise: character values and codes have no order the package uses, so the reference must be stated, and a labelled vector is refused until it is decoded. Labels are compared after trimming white space, making each run of it one space and ASCII letters lower case; other letters are compared as they are.

Value

An object of class catchk, a list: ok (TRUE when no variable has a problem), nmin, vars, a data frame with one row per variable, and levels, one row per level. vars: var, type, n (values), nna (missing), nlev (levels), unused (levels never seen), low (levels seen fewer than nmin times), twin (levels with a label that differs from another only by case or white space), bin (logical or 0/1), codes (whole numbers, not 0/1, with at most 10 distinct values: likely category codes), nolab (values seen with no value label), dlab (value labels shared by two values), ref (the default reference, NA for none), why, want (the stated reference), state ("default": the stated level is the default; "other": present, but must be passed explicitly; "low", "empty" (never seen), "case" (present only up to case or white space), "absent"; NA when none is stated) and ok. nlev, unused, low and twin are NA for a variable that is not categorical. levels: var, level, label (the value label, NA where none), n, low, twin (the level it differs from only by case or white space, NA where none) and ref (TRUE for the default reference). A variable has a problem when it has missing values, unused or low levels, twins, value labels, likely codes, fewer than two levels seen, or a stated reference in state "low", "empty", "case" or "absent".

Examples

d <- data.frame(sex = c("Male", "male", "Female", "Female", NA),
                custody = c(TRUE, FALSE, TRUE, TRUE, FALSE),
                court = c(1L, 2L, 2L, 3L, 1L))
catchk(d, nmin = 2)
catchk(factor(c("b", "a", "b"), levels = c("a", "b", "c")), ref = "b")

Description

The median of the x coordinates and the median of the y coordinates (the mean of the two middle values for an even number of points). It depends on the orientation of the axes: rotating the points does not rotate it. webpt() gives the rotation-invariant centre of minimum distance.

Usage

cgmedn(x, y = NULL, lonlat = FALSE)

Arguments

Value

A named numeric vector c(x, y), in degrees when lonlat = TRUE (the medians are taken of the projected coordinates).

References

Neft, D. S. (1966). Statistical Analysis for Areal Distributions. Regional Science Research Institute, Philadelphia.

Examples

cgmedn(c(1, 2, 3, 10), c(5, 1, 2, 2))

Description

The (weighted) mean centre of a point pattern, the standard deviations of the two coordinates, the standard distance deviation (the root mean squared distance from the mean centre), and the geometric and harmonic means of the coordinates (Bachi 1957; Neft 1962). The mean centre is the point that minimises the sum of squared distances to the points. The dispersions divide by the total weight less k: k = 1 makes the squared values unbiased for the variances; k = 2 is the convention of some crime-mapping texts; k = 0 divides by the total weight. That the mean centre minimises the weighted sum of squared distances is proved in Lean: see proof("cgstat").

Usage

cgstat(x, y = NULL, w = NULL, k = 1, lonlat = FALSE)

Arguments

Value

A named numeric vector: w (the total weight), mx, my (the mean centre, in degrees when lonlat = TRUE), sx, sy (standard deviations of the coordinates), sdd (the standard distance deviation), gx, gy (geometric means) and hx, hy (harmonic means). The last four are NaN for a coordinate with a non-positive value under positive weight, and for lonlat = TRUE (they depend on the origin of the coordinates). Lengths are in the unit of the coordinates, or metres.

References

Bachi, R. (1957). Statistical analysis of geographical series. Bulletin de l'Institut International de Statistique, 36(2), 229-240.

Neft, D. S. (1966). Statistical Analysis for Areal Distributions. Regional Science Research Institute, Philadelphia.

Examples

set.seed(1)
p <- cbind(x = rnorm(50, 100, 10), y = rnorm(50, 200, 5))
cgstat(p)
cgstat(p[, 1], p[, 2], w = rep(1:2, 25))
cgstat(p, k = 0)  # dispersions over n, not n - 1

Description

Compares two striated tool-mark profiles by the deterministic Chumbley score (Chumbley et al. 2010; Hadler and Morris 2018): each profile is trimmed and detrended by its lowess smooth, the pair of windows of wopt points with the largest correlation fixes the shift, and windows of wval points at that shift are compared with windows at other shifts by a standardised Wilcoxon rank-sum statistic.

Usage

chmb(a, b, wopt, wval, coarse = 0.25)

Arguments

Value

A named numeric vector: nsame, ndiff (window pairs used), u (the statistic) and p (its upper normal tail; NaN when either set is empty).

References

Chumbley, L. S., Morris, M. D., Kreiser, M. J. et al. (2010). Validation of tool mark comparisons obtained using a quantitative, comparative, statistical algorithm. Journal of Forensic Sciences, 55(4), 953-961. doi:10.1111/j.1556-4029.2010.01424.x

Hadler, J. R. and Morris, M. D. (2018). An improved version of a tool mark comparison algorithm. Journal of Forensic Sciences, 63(3), 849-855. doi:10.1111/1556-4029.13640

Examples

set.seed(1)
m <- cumsum(rnorm(700))
a <- m[1:600] + rnorm(600, sd = 0.2)
b <- m[51:650] + rnorm(600, sd = 0.2)
chmb(a, b, wopt = 120, wval = 30)

Description

Scores a set of likelihood ratios against the truth. Ratios from same-source comparisons should be large and ratios from different-source comparisons small. The cost is 1 for a system that always answers 1. It is above 0 for any set of ratios, and tends to 0 as the ratios of a system that separates the two kinds perfectly grow without bound. See proof("cllr") and cllrmin().

Usage

cllr(same, diff)

Arguments

Value

One number, the cost in bits.

References

Brummer, N. and du Preez, J. (2006). Application-independent evaluation of speaker detection. Computer Speech and Language, 20(2-3), 230-275. doi:10.1016/j.csl.2005.08.001

Examples

cllr(same = c(50, 200, 8), diff = c(0.01, 0.2, 1.5))

Description

Splits the cost of a set of likelihood ratios into min, their cost after recalibration with the pool-adjacent-violators algorithm, which keeps their order and gives each block of pooled ratios the value that costs that block least, and cal, the rest (cllr - min). That min is the least cost of every order-keeping recalibration, so that cal is never negative, is not proved. A system can separate the two kinds of comparison well and still report ratios that are too large or too small. See proof("cllrmin") and proof("lrcal").

Usage

cllrmin(same, diff)

Arguments

Value

A named numeric vector: cllr, min and cal, with cllr = min + cal.

References

Brummer, N. and du Preez, J. (2006). Application-independent evaluation of speaker detection. Computer Speech and Language, 20(2-3), 230-275. doi:10.1016/j.csl.2005.08.001

Examples

cllrmin(same = c(50, 200, 8), diff = c(0.01, 0.2, 1.5))

Description

A preset says how the columns of one published dataset (an Ontario segregation table, the Special Investigations Unit's cases, a use-of-force file) feed the arguments of one of crimlab's generic functions, so that generic(clprep(data, preset)) gives the published analysis while the same generic runs on anyone's data given as plain vectors. A preset describes a layout and ships no data: the user supplies the table.

Usage

clprep(data = NULL, preset = NULL)

Arguments

Details

The presets are rows of presets.tsv (in the installed package's presets folder), one per preset and argument: preset, generic, arg, column (the dataset's column), recode, filter, units, note, source. Python reads the same file. A recode is one of: empty (the column as it is); fac (a factor, one effect per value); num (numbers); one (a 1 for every row); ⁠sub:a,b⁠ (characters a to b, as a year from an ISO date); const:v (the argument is v, no column); or a mapping ⁠x|y=1;has:abc=1;*=0⁠ from labels, compared without case or surrounding space (⁠has:⁠ a label containing the text), to values, * covering every other label and ⁠*=NA⁠ dropping those rows. A filter col=value keeps the rows whose column has that value. An arg of recipe derives columns first: long stacks the two count columns named in recode as count with with 1 for the first; acvm and pyr reduce one row per placement to one per person and year (see Details of the returned attributes).

Rows with a missing value in any column the preset reads are dropped, as the source analyses drop them; their number is the attribute drop.

Value

With a preset, an object of class clpre: a list of the generic's arguments, which the generic takes as its first argument, with attributes preset, generic, source, units, note and drop. Without one, a data frame of the presets: preset, generic, args and source.

Examples

head(clprep())
# a table laid out as the published OTIS c01 table (numbers made up)
d <- data.frame(EndFiscalYear = rep(2020:2022, each = 2),
                Gender = rep(c("Female", "Male"), 3),
                NumberIndividuals_RestrictiveConfinement =
                  c(40, 310, 52, 335, 47, 360))
p <- clprep(d, "otis-c01-cnt")
cntreg(p)
# the same model on plain vectors
cntreg(d$NumberIndividuals_RestrictiveConfinement, d$Gender == "Female",
       adj = factor(d$EndFiscalYear))

Description

Splits the change in an aggregate clearance rate between two periods, areas or services (A and B) into a rate effect, from each offence type clearing more or less, and a mix effect, from the workload moving between types that clear at different rates (Kitagawa's decomposition). The two always add up to the change. When every type clears better in B and the aggregate still falls, the fall is all mix: Simpson's paradox, and the reason an aggregate clearance rate cannot be read as performance without the split by type. See proof("clrsmp").

Usage

clrsmp(clra, reca, clrb, recb)

Arguments

Details

Each part's mix term is measured from the mean of the two aggregate rates, so a part's term is negative when work moved toward a part that clears below average; the centre changes how the mix effect is shared among the parts, never its total.

Value

A list. parts, a data frame with one row per type, in the order of clra: part, ratea, rateb (clearance rates), sharea, shareb (shares of the recorded workload), rateeff and mixeff. total, a named numeric vector: ratea, rateb (aggregate rates), change, rateeff, mixeff (with change = rateeff + mixeff), up and down (types whose rate rose or fell) and simpson: 2 when every type moved one way and the aggregate the other, 1 when the aggregate change and the rate effect differ in sign, 0 otherwise.

References

Kitagawa, E. M. (1955). Components of a difference between two rates. Journal of the American Statistical Association, 50(272), 1168-1194. doi:10.1080/01621459.1955.10501299

Examples

# both types clear better in B, yet the aggregate falls
clrsmp(clra = c(minor = 90, serious = 20), reca = c(100, 50),
       clrb = c(19, 60), recb = c(20, 140))

Description

Fits a log-linear model of a count on a group, with adjustment columns, an exposure and clusters, by Poisson and negative binomial maximum likelihood, and reports the incidence rate ratio (IRR) of each group against the reference with its interval. The same model serves counts of people, placements, days or deaths in the cells of a published table and counts per person or per area.

Usage

cntreg(
  y,
  grp = NULL,
  adj = NULL,
  expo = NULL,
  clus = NULL,
  ref = NULL,
  gee = FALSE,
  cadj = TRUE,
  level = 0.95
)

Arguments

Details

The Poisson fit takes the steps of stats::glm() (iteratively reweighted least squares from y + 0.1, stopping when the relative change of the deviance is below 1e-8), and the negative binomial fit (variance mu + mu^2 / theta) those of MASS::glm.nb(): theta by maximum likelihood and the coefficients in turn. A column of the design that the columns before it explain (an adjustment level the group determines, say) is aliased and dropped, as R does: its coefficients are NA.

With clus, each fit also gets cluster-robust (sandwich) standard errors, rows "poiscl" and "nbcl": ⁠(X'WX)^-1 (sum_g U_g U_g') (X'WX)^-1⁠ over the clusters' scores U_g, times G / (G - 1) for G clusters when cadj (the default of sandwich::vcovCL(); cadj = FALSE gives the independence GEE variance of 'geepack'). With gee = TRUE the Poisson model is also fitted by generalised estimating equations with an exchangeable working correlation (row "gee"), as geepack::geeglm(corstr = "exchangeable") estimates it.

Intervals are Wald intervals on the log scale, ⁠exp(est +- z se)⁠, with z the normal quantile for level. sep flags a ratio whose group or reference has no counts at all: the estimate then runs off towards 0 or infinity and its interval means nothing.

Counts need not be whole numbers (a median duration, say): the estimating equations are those of the Poisson pseudo-likelihood, and the log-likelihood uses lgamma(y + 1), where stats::glm() reports an infinite AIC.

Value

An object of class cntreg, a list: irr, a data frame with one row per fit and group (fit, term, irr, lo, hi, p, sep); coef, every coefficient of every fit on the log scale (fit, term, est, se, z, p); fit, one row per model (fit, loglik, aic, theta, setheta, alpha, iter, conv, rank); ref, the reference group; n, the number of counts; ncl, the number of clusters (0 without clus); and level.

References

Venables, W. N. and Ripley, B. D. (2002). Modern Applied Statistics with S, 4th ed. Springer. doi:10.1007/978-0-387-21706-2

Liang, K.-Y. and Zeger, S. L. (1986). Longitudinal data analysis using generalized linear models. Biometrika, 73(1), 13-22. doi:10.1093/biomet/73.1.13

Cameron, A. C. and Miller, D. L. (2015). A practitioner's guide to cluster-robust inference. Journal of Human Resources, 50(2), 317-372. doi:10.3368/jhr.50.2.317

Examples

# placements by group in eight institutions over three years
d <- data.frame(
  n = c(12, 30, 8, 25, 14, 41, 9, 22, 11, 35, 7, 28),
  female = rep(c("no", "yes"), 6),
  year = factor(rep(c(2021, 2022, 2023), each = 4)),
  inst = rep(c("A", "B", "C", "D"), each = 2, length.out = 12)
)
cntreg(d$n, d$female, adj = d["year"], clus = d$inst)
# rates: a count per person-days of exposure
cntreg(c(5, 12), c("a", "b"), expo = c(1000, 1500))

Description

The one-mode projection of offenders onto incidents: two offenders are joined when they are recorded in the same incident, and the weight of the tie is the number of incidents they share. An offender recorded twice in one incident counts once there. Offenders who never offend with anyone appear in nodes but in no edge; give nodes to the network functions to keep them as isolated nodes.

Usage

coofg(inc, off)

Arguments

Details

Ties are listed in the order each pair is first met: incidents in the order they first appear, and within an incident the pairs of its offenders in the order of the rows; from is the offender who appears first in off. Labels are compared as text (as.character(); a factor by its labels).

Value

A list: edges, a data frame with from, to (offender labels) and w (incidents shared), and nodes, every offender in order of first appearance.

References

Reiss, A. J. (1988). Co-offending and criminal careers. Crime and Justice, 10, 117-170. doi:10.1086/449145

Examples

g <- coofg(inc = c(1, 1, 1, 2, 2, 3), off = c("a", "b", "c", "a", "b", "d"))
g$edges
ntdeg(g$edges$from, g$edges$to, w = g$edges$w, nodes = g$nodes)

Description

Treats the successive locations of a series (incidents attributed to one offender, in time order) as a correlated random walk (Kareiva and Shigesada 1983): steps of independent lengths and changes of heading (turning angles) independent of each other and of the steps. The fit gives the mean and mean squared step and the mean cosine and sine of the turns; the expected squared displacement after k moves is Kareiva and Shigesada's for symmetric turning (mean sine 0), ⁠k E l^2 + 2 (E l)^2 c / (1 - c) (k - (1 - c^k) / (1 - c))⁠, next to the observed mean over every start. The walk is then continued nsim times for h steps from the last location and heading, each turn and step drawn independently and uniformly from the observed ones (an empirical correlated random walk), and the simulated positions summarise where the next events may fall. It describes movement; it does not identify anyone, and a series with few events gives a rough walk.

Usage

crwalk(x, y = NULL, h = 1, nsim = 999, seed = NULL, lonlat = FALSE)

Arguments

Value

A list. fit, a named vector: steps, ml (mean step), ml2 (mean squared step), c and s (mean cosine and sine of the turning angles) and head (the last heading, radians anticlockwise from east). msd, a data frame with one row per step ahead k: obs (the observed mean squared displacement over k moves, NA beyond the series), exp (its expectation; it assumes s near 0), and x, y, sdx, sdy (the mean and standard deviations of the simulated positions; the means in degrees when lonlat = TRUE, the rest in metres). sim, the simulated positions: matrices x and y, one row per walk and one column per step (planar, metres when lonlat = TRUE).

References

Kareiva, P. M. and Shigesada, N. (1983). Analyzing insect movement as a correlated random walk. Oecologia, 56(2-3), 234-238. doi:10.1007/BF00379695

Examples

set.seed(1)
a <- cumsum(rnorm(20, 0, 0.5))
p <- cbind(x = cumsum(cos(a)), y = cumsum(sin(a)))
r <- crwalk(p, h = 3, nsim = 200, seed = 1)
r$fit
r$msd

Description

Builds a two-way table from labels (counting rows, summing a count column of an aggregate table, or counting distinct ids per cell) and gives Pearson's chi-square test, the likelihood-ratio test, Cramer's V with an interval and its bias-corrected form, the contingency coefficient, phi for a 2 x 2 table, the share on the diagonal (the rows whose label is also the column's, as a region at placement against the region now), and on request Fisher's exact test with, for a 2 x 2 table, the conditional odds ratio and its exact interval.

Usage

ctabx(
  row,
  col = NULL,
  w = NULL,
  agg = NULL,
  correct = FALSE,
  exact = FALSE,
  level = 0.95
)

Arguments

Details

Rows or columns with no counts are dropped before the statistics (their number is drop), as an unused factor level carries no information; the table tab keeps them. Cramer's V is sqrt(chisq / (n (k - 1))), k the smaller dimension, from the uncorrected statistic even when correct is TRUE; its interval inverts the noncentral chi-square distribution, sqrt(ncp / (n (k - 1))), as the default of DescTools::CramerV(); the bias-corrected V is Bergsma's. The exact test enumerates every table with the observed margins and refuses one that would take more than 1e7 tables. With agg other than "n", "sum" or "dist" the table holds values, not counts, and the statistics are NA.

Value

An object of class ctabx, a list: tab, the table (all labels); expd and res, the expected counts and Pearson residuals of the table without empty rows and columns; stat, a one-row data frame: n, chisq, df, p, g2, g2p (the likelihood-ratio statistic and its p-value), v, vlo, vhi, vbc (bias-corrected V), cc (the contingency coefficient), phi, diag, fp (Fisher's p-value), or, orlo, orhi; and drop, agg, correct and level.

References

Agresti, A. (2013). Categorical Data Analysis, 3rd ed. Wiley. ISBN 978-0-470-46363-5.

Bergsma, W. (2013). A bias-correction for Cramer's V and Tschuprow's T. Journal of the Korean Statistical Society, 42(3), 323-328. doi:10.1016/j.jkss.2012.10.002

Examples

# force type by weapon, one row per incident
wp <- c("none", "none", "knife", "gun", "none", "knife", "gun", "none")
fc <- c("hands", "hands", "CEW", "firearm", "CEW", "hands", "firearm",
        "hands")
ctabx(wp, fc)
# region at placement against region now: the diagonal share
ctabx(c("N", "N", "E", "W", "E"), c("N", "E", "E", "W", "E"))
# a published aggregate table: sum its count column
d <- data.frame(race = c("A", "A", "B", "B"), sex = c("F", "M", "F", "M"),
                n = c(12, 40, 9, 51))
ctabx(d$race, d$sex, w = d$n, exact = TRUE)

Description

The smallest convex polygon holding every point, by Andrew's monotone chain with an exact orientation predicate (Shewchuk 1997), so near-collinear points never make a wrong turn. Points on an edge are not vertices; of coincident points the first is reported.

Usage

cvhull(x, y = NULL, lonlat = FALSE)

Arguments

Value

A list: idx, the row numbers of the vertices, anticlockwise from the point with the lowest x (then lowest y); area and perim, the hull's area and perimeter (in the square and the unit of the coordinates, or square metres and metres when lonlat = TRUE). One distinct point gives one vertex and zeros; collinear points give their two ends, area 0 and twice the segment's length. Each coordinate must be 0 or of magnitude from 1e-140 to 1e150, where the predicate is exact.

References

Andrew, A. M. (1979). Another efficient algorithm for convex hulls in two dimensions. Information Processing Letters, 9(5), 216-219. doi:10.1016/0020-0190(79)90072-3

Shewchuk, J. R. (1997). Adaptive precision floating-point arithmetic and fast robust geometric predicates. Discrete and Computational Geometry, 18(3), 305-363. doi:10.1007/PL00009321

Examples

cvhull(c(0, 1, 1, 0, 0.5), c(0, 0, 1, 1, 0.5))

Description

The rate each group would have with the stratum make-up (age bands, say) of a standard population: the sum over strata of the standard's weights times the group's stratum rates. The interval is Fay and Feuer's gamma interval; the ratio of two standardised rates has a log-normal interval and their difference a normal one.

Usage

dastd(x, n, strat, grp = NULL, std, ref = NULL, per = 1, level = 0.95)

Arguments

Value

An object of class dastd, a list: rates (group, x the events, rate, lo, hi, var the variance of the rate), cmp (one row per group but the reference: group, ratio, lo, hi, diff, dlo, dhi), ref, per and level.

References

Fay, M. P. and Feuer, E. J. (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-#

Examples

d <- data.frame(age = rep(c("18-29", "30-49", "50+"), 2),
                grp = rep(c("A", "B"), each = 3),
                x = c(30, 22, 5, 41, 30, 12),
                n = c(5000, 9000, 6000, 3000, 8000, 7000))
dastd(d$x, d$n, d$age, d$grp, std = c(30, 45, 25), per = 1e4)

Description

What a match found by trawling a database is worth. np is the bound of the second National Research Council report: the match probability times the number of profiles searched. any is the exact chance of at least one chance match among that many unrelated profiles, never above np. post takes the other view (Balding and Donnelly): everyone else in the database was excluded, so among pop possible sources with equal prior the matching person's probability of being the source goes up, not down. See proof("dbsrch").

Usage

dbsrch(mp, n, pop = NULL)

Arguments

Value

A named numeric vector: np, any and, when pop is given, post.

References

Balding, D. J. and Donnelly, P. (1996). Evaluating DNA profile evidence when the suspect is identified through a database search. Journal of Forensic Sciences, 41(4), 603-607.

Examples

dbsrch(mp = 1e-9, n = 5e6, pop = 3e7)

Description

The mean direction of a set of directions, the mean resultant length R (from 0, directions spread evenly, to 1, all the same), the circular variance 1 - R and the circular standard deviation ⁠sqrt(-2 log R)⁠ (Mardia and Jupp 2000). Directions are given as angles, or as points seen from an origin (the directions of incidents from an offender's base, say). The mean direction of points depends on the origin chosen. That R lies in ⁠[0, 1]⁠ is proved in Lean: see proof("dirmn").

Usage

dirmn(
  x = NULL,
  y = NULL,
  origin = c(0, 0),
  w = NULL,
  angle = NULL,
  lonlat = FALSE
)

Arguments

Value

A named numeric vector: theta (the mean direction in radians anticlockwise from east, in ⁠[0, 2 pi)⁠; NaN when R is 0), bearing (the same as a compass bearing, degrees clockwise from north), r, var, sd, w (the total weight) and n.

References

Mardia, K. V. and Jupp, P. E. (2000). Directional Statistics. Wiley. doi:10.1002/9780470316979

Examples

dirmn(c(1, 2, 0), c(1, 1, 3), origin = c(0, 0))
dirmn(angle = c(0.1, 0.2, 6.2))

Description

Estimates a causal parameter from orthogonal scores whose nuisance functions (outcome regressions, propensities) are learned out of fold: each fold's predictions come from learners fitted on the other folds, so that the estimate keeps root-n inference while the learners are flexible (Chernozhukov et al. 2018). The learning, the scores, the estimate and its standard error are all compiled; the folds are the caller's, so results can be matched with other software fold for fold.

Usage

dmlx(
  y,
  d,
  x,
  model = "plr",
  fold,
  z = NULL,
  s = NULL,
  t = NULL,
  score = NULL,
  lrn = "lasso",
  lrnd = lrn,
  alpha = 0.5,
  nlam = 100,
  nfin = 5,
  tol = 1e-07,
  scr = TRUE,
  hon = FALSE,
  ntree = 500,
  mtry = NULL,
  trim = NULL,
  clu = NULL,
  pred = NULL,
  seed = 1
)

Arguments

Details

The models are those of 'DoubleML', with its scores and its defaults:

• "plr", the partially linear regression ⁠y = theta d + g(x) + e⁠, score "po" (partialling out: nuisances l = E[y | x], m = E[d | x]) or "iv" (IV-type: also g, learned on ⁠y - theta0 d⁠ with theta0 from the partialling-out score);

• "pliv", the partially linear IV model with one instrument z (nuisances l, m = E[z | x], r = E[d | x], and g for "iv");

• "irm", the interactive regression model for a 0/1 treatment, score "ate" or "atte" (nuisances g0, g1, the outcome regressions among the untreated and the treated, and the propensity m);

• "iivm", the interactive IV model (local average effect) for a 0/1 treatment and a 0/1 instrument (g0, g1, m = P(z = 1 | x), r0, r1, the treatment probabilities given z);

• "did", difference in differences on a panel: y is the change in the outcome, d the 0/1 group (nuisances g0, the change among the untreated, and m), with the in-sample normalised observational score of 'DoubleML' for 'Python' (Sant'Anna and Zhao 2020; Chang 2020);

• "didcs", the same on repeated cross-sections, t the 0/1 period (nuisances g00, g01, g10, g11, the outcome regressions in each group and period, and m);

• "ssm", the average effect when the outcome is observed only where s is 1, missing at random given d and x (Bia, Huber and Laffers 2023; nuisances pi = P(s = 1 | d, x), m, g0, g1);

• "aipw" and "tmle", the average effect of a 0/1 treatment from the nuisances of "irm": augmented inverse probability weighting, with the standard error of the 'AIPW' package (the standard deviation of the influence function), and targeted maximum likelihood with the fluctuation of the 'tmle' package (Gruber and van der Laan 2012).

The partialling-out score is Neyman-orthogonal, the AIPW term is unbiased when either nuisance is right, and the estimate is the root of the empirical score: see proof("dmlx").

The estimate solves the mean of the score over all folds (DoubleML's "dml2"); its variance is mean(psi^2) / mean(psi_a)^2 / n, or, with clusters (one-way, models "plr" and "pliv"), the cluster-robust form of Chiang et al. (2022). With several repetitions of the folds (columns of fold), the estimate is the median over repetitions and the variance the median of var + (theta_r - theta)^2.

Value

A list: coef and se, the estimate and its standard error; reps, a data frame with each repetition's coef and se; psia and psib, the score elements (one column per repetition; the estimate solves mean(psia * coef + psib) = 0); pred, the out-of-fold predictions, an array named by nuisance.

References

Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W. and Robins, J. (2018). Double/debiased machine learning for treatment and structural parameters. The Econometrics Journal, 21(1), C1-C68. doi:10.1111/ectj.12097

Chiang, H. D., Kato, K., Ma, Y. and Sasaki, Y. (2022). Multiway cluster robust double/debiased machine learning. Journal of Business and Economic Statistics, 40(3), 1046-1056. doi:10.1080/07350015.2021.1895815

Sant'Anna, P. H. C. and Zhao, J. (2020). Doubly robust difference-in-differences estimators. Journal of Econometrics, 219(1), 101-122. doi:10.1016/j.jeconom.2020.06.003

Chang, N.-C. (2020). Double/debiased machine learning for difference-in-differences models. The Econometrics Journal, 23(2), 177-191. doi:10.1093/ectj/utaa001

Bia, M., Huber, M. and Laffers, L. (2023). Double machine learning for sample selection models. Journal of Business and Economic Statistics, 42(3), 958-969. doi:10.1080/07350015.2023.2271071

Gruber, S. and van der Laan, M. J. (2012). tmle: an R package for targeted maximum likelihood estimation. Journal of Statistical Software, 51(13), 1-35. doi:10.18637/jss.v051.i13

Examples

set.seed(1)
n <- 500
x <- matrix(rnorm(n * 5), n)
d <- x[, 1] + rnorm(n)
y <- 0.5 * d + x[, 1] - x[, 2] + rnorm(n)
fold <- sample(rep(1:5, length.out = n))
r <- dmlx(y, d, x, "plr", fold, lrn = "glm", lrnd = "glm")
c(r$coef, r$se)

Description

Splits each group's disparity in force per resident, against a reference group, into a disparity in contact per resident times a disparity in force per contact. The split is a product, never a sum. See proof("dspbm").

Usage

dspbm(pop, contact, force, ref = 1)

Arguments

Value

A data frame with one row per group, in the order of pop: group, ref (TRUE for the reference group), cd (contact per resident), fd (force per contact) and rd (force per resident), each against the reference group, with rd = cd * fd.

Examples

dspbm(pop = c(a = 9000, b = 1000), contact = c(a = 900, b = 300),
      force = c(a = 9, b = 6))

Description

Offences by one group against another, as a rate per member of the offending group and as a rate per offender-victim pair exposed. The rate per offender rises with the size of the victim group even when offenders choose victims at random; the rate per pair does not. See proof("dyadrt").

Usage

dyadrt(count, pop)

Arguments

Value

A data frame with one row per dyad: off, vic, count, rtoff, rtpair, nulloff (the rate per offender under random mixing) and ratio (rtpair over the pooled pair rate; 1 under random mixing).

Examples

pop <- c(a = 800, b = 200)
dyadrt(count = outer(pop, pop) / 1000, pop = pop)

Description

The information a set of likelihood ratios still leaves missing, across prior odds: for the ratios as reported, for the same ratios after calibration, and for a system that always reports 1. Reported ratios above the last curve do harm at that prior. At even prior odds the first two curves pass through the values of cllrmin(). proof() records no theorem for this plot.

Usage

eceplt(same, diff, lo = seq(-2.5, 2.5, by = 0.05), ...)

Arguments

Value

Invisibly, a data frame: lo, ece, min and ref.

References

Ramos, D. and Gonzalez-Rodriguez, J. (2013). Reliable support: measuring calibration of likelihood ratios. Forensic Science International, 230(1-3), 156-169. doi:10.1016/j.forsciint.2013.04.014

Examples

eceplt(same = c(50, 200, 8, 0.5), diff = c(0.01, 0.2, 1.5, 0.001))

Description

Fits the lasso, ridge and elastic net path by cyclical coordinate descent, for a numeric (gaussian) or a 0/1 (binomial) outcome, and cross-validates it over folds the caller supplies. The objective, the standardisation of the columns and the scale of lambda are those of 'glmnet', so that with the same lambda and the same folds the two give the same coefficients, up to their convergence thresholds: the gaussian family minimises ⁠(1 / 2n) |y - a0 - x b|^2 + lambda pen(b)⁠, the binomial family the mean negative log-likelihood plus ⁠lambda pen(b)⁠, with ⁠pen(b) = sum((1 - alpha) / 2 (s b)^2 + alpha |s b|)⁠, s the population standard deviations of the columns. For the gaussian family the outcome is scaled to unit variance while fitting, as 'glmnet' does. The intercept is never penalised; a constant column keeps a zero coefficient.

Usage

enpath(
  x,
  y,
  fam = "gaussian",
  alpha = 1,
  lambda = NULL,
  nlam = 100,
  lmr = NULL,
  fold = NULL,
  tol = 1e-07,
  maxit = 1e+05,
  scr = TRUE
)

Arguments

Details

Speed against accuracy: tol and maxit set how far the descent goes. The default tol = 1e-7 is glmnet's; ⁠tol = 1e-20, maxit = 1e7⁠ is the exact setting, which reproduces glmnet run at the same threshold to 1e-8 in every coefficient (tests). Between the two, the penalised objective at each penalty differs by less than 1e-6 of its value (tested on the gaussian and binomial paths). nlam and lmr set the length and the end of a path chosen here, and scr the screening of a binomial path.

Each coordinate step is the minimiser of its one-dimensional objective, by soft-thresholding, the proximal map of the absolute value: see proof("enpath").

The binomial fit is a proximal Newton method: iteratively reweighted least squares around coordinate descent over the active columns, with the optimality conditions checked over every column before a lambda is left. The gaussian fit updates the gradient through the column inner products, formed only for the columns that enter the model.

Value

A list: lambda, the penalties fitted; a0, the intercepts; beta, a matrix of coefficients with one column per penalty and one row per column of x; and, when fold is given, cvm and cvsd, the mean held-out loss (squared error, or the binomial deviance with probabilities kept within 1e-5 of 0 and 1) and its standard error over folds, computed as 'glmnet' computes them, lmin, the penalty at the minimum, and l1se, the largest penalty within one standard error of it; these four are NULL without folds.

References

Friedman, J., Hastie, T. and Tibshirani, R. (2010). Regularization paths for generalized linear models via coordinate descent. Journal of Statistical Software, 33(1), 1-22. doi:10.18637/jss.v033.i01

Tibshirani, R., Bien, J., Friedman, J., Hastie, T., Simon, N., Taylor, J. and Tibshirani, R. J. (2012). Strong rules for discarding predictors in lasso-type problems. Journal of the Royal Statistical Society: Series B, 74(2), 245-266. doi:10.1111/j.1467-9868.2011.01004.x

Examples

set.seed(1)
x <- matrix(rnorm(200 * 5), 200)
y <- x[, 1] - 0.5 * x[, 2] + rnorm(200)
f <- enpath(x, y, nlam = 20, fold = rep(1:5, 40))
f$lmin
f$beta[, f$lambda == f$lmin]

Description

eqproj() maps longitude and latitude to planar metres by the Lambert azimuthal equal-area projection of the sphere about a centre; eqinv() maps them back. Areas are exact, so intensities, nearest-neighbour expectations and K functions computed on the result are right; lengths are scaled by cos(c / 2) along the radius from the centre and its inverse across it (c the central angle), so within 100 km of the centre every length is right to 3.1e-5 relative. Every planar function of the package takes the projected coordinates; those with a lonlat argument call this projection themselves.

Usage

eqproj(x, y = NULL, lon0 = NULL, lat0 = NULL, r = 6371008.8)

eqinv(x, y = NULL, lon0, lat0, r = 6371008.8)

Arguments

Value

A data frame with one row per point, in input order: x and y (east and north, in the unit of r) from eqproj(), with the centre used as the attribute centre; lon and lat from eqinv().

References

Snyder, J. P. (1987). Map Projections: A Working Manual. U.S. Geological Survey Professional Paper 1395, chapter 24. doi:10.3133/pp1395

Examples

p <- eqproj(c(-79.40, -79.38, -79.35), c(43.64, 43.66, 43.65))
p
a <- attr(p, "centre")
eqinv(p, lon0 = a[1], lat0 = a[2])

Description

Recorded counts (stops, searches, uses of force) are rates per unit of exposure to police, and exposure is seldom observed. If each group's true exposure lies within a factor gamma of an observed proxy (residents, drivers on the road, calls for service), each group's true rate lies within a factor gamma of its proxy rate, and the ratio of two groups' true rates within a factor gamma^2 of the proxy ratio. The direction of a disparity is settled only when the proxy ratio lies outside that band. See proof("expbnd").

Usage

expbnd(y, m, gamma, ref = 1)

Arguments

Value

A data frame with one row per group, in the order of y: group, ref (TRUE for the reference group), rate (y / m), lo and hi (the bounds on the true rate), ratio (the proxy rate against the reference group), rlo and rhi (the bounds on the true ratio; every end is reached by some admissible exposure), and sign: 1 when the true ratio exceeds 1 for every admissible exposure, -1 when it is below 1 for every one, 0 when the data leave it open.

References

Knox, D., Lowe, W. and Mummolo, J. (2020). Administrative records mask racially biased policing. American Political Science Review, 114(3), 619-637. doi:10.1017/S0003055420000039

Examples

expbnd(y = c(white = 100, black = 300), m = c(1000, 1000), gamma = 1.5)

Description

The confusion matrix of a risk score cut at cut within each group, and the usual rates computed from it: accuracy, sensitivity (true positive rate), specificity, false positive and false negative rates, positive and negative predictive values, the share and number predicted positive and the Matthews correlation, each with its ratio to the base group and, for the first eight, a Wilson score interval (Wilson 1927). A case is predicted positive when its score exceeds cut. No single metric defines fairness, and when base rates differ between groups the error rates and the predictive values cannot all be equal at once (Chouldechova 2017); the table shows which differ, not which should not.

Usage

fair(score, y, grp, cut = 0.5, base = 1, conf = 0.95)

Arguments

Value

A list of data frames, one row per group with its label in group: cnt (tp, fp, tn, fn, n), met (acc, tpr, tnr, fpr, fnr, ppv, npv, ppos, npos, mcc; NaN where a rate has an empty denominator), rat (each metric over the base group's), lo and hi (the interval ends of the first eight metrics; the last two columns NaN).

References

Wilson, E. B. (1927). Probable inference, the law of succession, and statistical inference. Journal of the American Statistical Association, 22(158), 209-212. doi:10.1080/01621459.1927.10502953

Chouldechova, A. (2017). Fair prediction with disparate impact: a study of bias in recidivism prediction instruments. Big Data, 5(2), 153-163. doi:10.1089/big.2016.0047

Examples

set.seed(1)
g <- rep(c("a", "b"), each = 100)
y <- rbinom(200, 1, ifelse(g == "a", 0.3, 0.5))
s <- plogis(rnorm(200, ifelse(y == 1, 0.8, -0.8)))
fair(s, y, g)$met

Description

How the share of cases a tool flags (predicts positive) differs between groups, as a difference and as a ratio against a base group, each with an interval; with the outcomes, how its true and false positive rates differ; and the bias amplification score, the largest parity gap times the Gini coefficient of the groups' shares.

Usage

fairdf(
  pred,
  grp = NULL,
  y = NULL,
  base = NULL,
  conf = 0.95,
  nboot = 0,
  seed = NULL
)

Arguments

Details

fair() reports each group's rates and their ratios; fairdf() reports the comparisons the fairness literature names. The parity gap dp is a group's share flagged minus the base group's (demographic parity holds when every gap is 0), with Newcombe's (1998) hybrid score interval. The ratio di is the disparate impact ratio, with the log interval of Katz et al. (1978); under the four-fifths rule a ratio below 0.8 against the group with the highest selection rate is the usual sign of adverse impact, so by default that group is the base. With the outcomes y, dtpr and dfpr are the gaps in true and false positive rates (both 0 is equalized odds), aod their mean (the average odds difference) with the interval of Zou and Donner (2008) built from the two Newcombe intervals, and eo the larger of the two gaps in size. The bias amplification score of Barman and Barman (2026) is large only when a sizeable gap comes with unequal shares across all the groups; its bootstrap interval, and the Gini coefficient's, are percentile intervals from redrawing each group's decisions with replacement (Efron and Tibshirani 1993), the base group held.

These numbers describe what the tool did to each group. They do not say why: a gap in shares flagged can come from a gap in the outcomes the tool was built to predict, from how those outcomes were recorded, or from the tool; and when the groups' base rates differ, equal error rates and equal predictive values cannot all hold at once (Chouldechova 2017, see impfpr()). A ratio's interval is NaN when either share is 0, and a ratio against a base group with no case flagged is Inf (or NaN for 0 over 0).

Value

A list: tab, a data frame with one row per group (in the order of grp): group, n, npos (cases flagged), ppos (their share) with lo and hi, dp with dplo and dphi, di with dilo and dihi, and with y also tpr, fpr, dtpr (dtprlo, dtprhi), dfpr (dfprlo, dfprhi), aod (aodlo, aodhi) and eo; the base group's own intervals are NaN. sum, a named vector: di (the smallest finite ratio of another group), dp (the gap of largest size, signed), gini, bas, ginilo, ginihi, baslo, bashi (NaN when nboot is 0) and with y also aod and eo (those of largest size, signed). base, the base group's label.

References

Newcombe, R. G. (1998). Interval estimation for the difference between independent proportions: comparison of eleven methods. Statistics in Medicine, 17(8), 873-890. doi:10.1002/(SICI)1097-0258(19980430)17:8<873::AID-SIM779>3.0.CO;2-I

Katz, D., Baptista, J., Azen, S. P. and Pike, M. C. (1978). Obtaining confidence intervals for the risk ratio in cohort studies. Biometrics, 34(3), 469-474. doi:10.2307/2530610

Zou, G. Y. and Donner, A. (2008). Construction of confidence limits about effect measures: a general approach. Statistics in Medicine, 27(10), 1693-1702. doi:10.1002/sim.3095

Barman, P. K. and Barman, P. K. (2026). Spatial disparities in simulated patrol detection: a reproducible generative-model audit across Baltimore and Chicago. doi:10.48550/arXiv.2603.18987

Efron, B. and Tibshirani, R. J. (1993). An Introduction to the Bootstrap. Chapman and Hall. doi:10.1007/978-1-4899-4541-9

Examples

set.seed(1)
g <- rep(c("a", "b"), each = 200)
y <- rbinom(400, 1, 0.3)
p <- rbinom(400, 1, ifelse(g == "a", 0.4, 0.25))
r <- fairdf(p, g, y, nboot = 199, seed = 1)
r$tab[, c("group", "ppos", "dp", "dplo", "dphi", "di", "dilo", "dihi")]
r$sum

Description

Incidents in segregated areas, whose true occurrence is the same for every group but whose recording, and the risk an algorithm predicts for them, differ by group by a set amount: data on which an audit's answer is known (after Barman and Barman 2026). The first group is the reference.

Usage

fairsim(
  n = 2000,
  groups = c("A", "B"),
  props = NULL,
  nareas = 20,
  base = 0.3,
  bias = 0.5,
  seed = NULL
)

Arguments

Details

Each record draws its group by props, then an area among those its group lives in (with g groups, areas 1, g + 1, ⁠2 g + 1⁠, ... hold the first group only, areas 2, g + 2, ... the second, and so on), then truth, 1 with probability base whatever the group; det, the recorded outcome, 1 with probability base in the reference group and base (1 - bias) in the others (clamped to 0 and 1), so that the other groups' incidents are recorded less when bias is positive; and risk, normal with mean 250 in the reference group and ⁠250 + 100 bias⁠ in the others and standard deviation 40, clamped to 0 and 500. So fairdf() on det should find a ratio near 1 - bias, scrdsp() on risk a gap near ⁠100 bias⁠, and calaud() the other groups' areas ranked as riskier than their recorded outcomes.

Value

A data frame with one row per record: area (a number from 1), group, truth, det and risk.

References

Barman, P. K. and Barman, P. K. (2026). Spatial disparities in simulated patrol detection: a reproducible generative-model audit across Baltimore and Chicago. doi:10.48550/arXiv.2603.18987

Examples

d <- fairsim(2000, seed = 1)
head(d)
# the injected ratio is 1 - 0.5
fairdf(d$det, d$group, base = "A")$tab[, c("group", "di", "dilo", "dihi")]

Description

For each point, the number of points (itself and coincident ones included) within distance r of it. With r = 0 it counts coincident points, the mode of exact locations; with r > 0 it is the fuzzy mode of the crime-mapping literature, a fixed-radius scan centred on each incident. Counts of nearby points overlap by design: each point counts its own neighbourhood.

Usage

fzmode(x, y = NULL, r, count = NULL, metric = "euclid", lonlat = FALSE)

Arguments

Value

An integer vector, one count per point, in input order (with count, a numeric vector of zone totals).

References

Eck, J. E., Chainey, S., Cameron, J. G., Leitner, M. and Wilson, R. E. (2005). Mapping Crime: Understanding Hot Spots. National Institute of Justice, NCJ 209393.

Examples

fzmode(c(0, 0, 1, 5), c(0, 0, 0, 5), r = 1)

Description

Geary's (1954) contiguity ratio of values on units with spatial weights w: ⁠C = (n - 1) sum w_ij (x_i - x_j)^2 / (2 S0 sum (x_i - xbar)^2)⁠, 1 under no autocorrelation, below 1 when neighbours are alike. Its variance is given under normality and under randomisation (Cliff and Ord 1981), with z statistics (1 - C) / sd and their upper normal tails, and, with nsim > 0, a permutation p-value against positive autocorrelation. That C has expectation 1 under randomisation is proved in Lean: see proof("gearc").

Usage

gearc(x, w, nsim = 0, seed = NULL)

Arguments

Value

A named numeric vector: c, vn and vr (its variance under normality and randomisation), zn, zr, pn, pr, pperm (⁠(1 + #{C* <= C}) / (nsim + 1)⁠, NaN without permutations) and s0 (the sum of the weights).

References

Geary, R. C. (1954). The contiguity ratio and statistical mapping. The Incorporated Statistician, 5(3), 115-145. doi:10.2307/2986645

Cliff, A. D. and Ord, J. K. (1981). Spatial Processes: Models and Applications. Pion, London.

Examples

set.seed(1)
p <- cbind(x = runif(50), y = runif(50))
gearc(p[, 1] + p[, 2] + rnorm(50, sd = 0.1), spwt(p, k = 4))

Description

The global concentration statistic of Getis and Ord (1992) for non-negative values: ⁠G = sum_ij w_ij x_i x_j / sum_{i != j} x_i x_j⁠, its expectation S0 / (n (n - 1)) and variance under randomisation, the z statistic, its upper normal tail (high values near high values) and, with nsim > 0, a permutation p-value. Unlike Moran's I and Geary's C it tells high-value clusters from low-value ones. Its expectation under randomisation is proved in Lean: see proof("gorg").

Usage

gorg(x, w, nsim = 0, seed = NULL)

Arguments

Value

A named numeric vector: g, e, v, z, p and pperm (⁠(1 + #{G* >= G}) / (nsim + 1)⁠).

References

Getis, A. and Ord, J. K. (1992). The analysis of spatial association by use of distance statistics. Geographical Analysis, 24(3), 189-206. doi:10.1111/j.1538-4632.1992.tb00261.x

Examples

set.seed(1)
p <- cbind(x = runif(50), y = runif(50))
gorg(exp(3 * p[, 1]), spwt(p, type = "band", d = 0.2, style = "raw"))

Description

The rate (events over exposure) or the proportion (events over trials) of each group, and its ratio to and difference from a reference group, with approximate and exact intervals. Rows that share a group are pooled: give one row per group, or one per case (0/1 outcomes with n = 1, the default for proportions), or anything between.

Usage

grprr(
  x,
  n = NULL,
  grp = NULL,
  ref = NULL,
  type = "rate",
  per = 1,
  level = 0.95
)

Arguments

Details

For a rate x / n the approximate interval is rate * exp(-+z / sqrt(x)) and the exact one the Poisson interval through the gamma quantiles, as stats::poisson.test(); the ratio's approximate interval is log-normal with variance 1 / x + 1 / x0, its exact interval the one conditional on the total count (the Clopper-Pearson interval of x / (x + x0), as the two-sample stats::poisson.test()), and the difference's is Wald. For a proportion the intervals are Wilson's score interval and the exact Clopper-Pearson one (stats::binom.test()); the ratio's is log-normal (Katz et al.), with no exact interval (NA), and the difference's is Newcombe's hybrid score interval. With a zero count the log-normal intervals are NA: no continuity correction is made.

rep is a group's share of the events over its share of the exposure (or the trials): over 1 when the group is over-represented among the events, as in a comparison of who is placed in segregation with who is in custody.

Value

An object of class grprr, a list: rates, a data frame with one row per group (group, x, n, rate, lo, hi the approximate interval, elo, ehi the exact one, rep); cmp, one row per group but the reference (group, ratio, lo, hi, elo, ehi, diff, dlo, dhi); ref, type, per and level.

References

Katz, D., Baptista, J., Azen, S. P. and Pike, M. C. (1978). Obtaining confidence intervals for the risk ratio in cohort studies. Biometrics, 34(3), 469-474. doi:10.2307/2530610

Newcombe, R. G. (1998). Interval estimation for the difference between independent proportions: comparison of eleven methods. Statistics in Medicine, 17(8), 873-890. doi:10.1002/(SICI)1097-0258(19980430)17:8<873::AID-SIM779>3.0.CO;2-I

Examples

# use-of-force incidents per 1,000 contacts by perceived race
grprr(c(white = 120, black = 95, other = 30),
      c(white = 60000, black = 21000, other = 15000), per = 1000)
# proportions from one row per case
grprr(c(1, 0, 0, 1, 1, 0, 1, 1), grp = c(1, 1, 1, 1, 2, 2, 2, 2),
      type = "prop")

Description

The local statistic Gi* of Getis and Ord (1992; Ord and Getis 1995) for each cell of a grid of counts, with binary weights on the cells within rad cell widths (rad = 1.5: the 3 by 3 block, cell included), its normal p-value and Benjamini-Hochberg q-value, a conditional permutation p-value, and an exposure-aware version: given the total, the count of a neighbourhood is binomial with the neighbourhood's share of the exposure under complete spatial randomness, so cells overhanging the study area neither inflate the mean nor get flagged. No threshold is applied.

Usage

gstr(cnt, ex = NULL, rad = 1.5, nsim = 0, seed = NULL)

Arguments

Value

A data frame with one row per cell (down the columns of cnt): row, col, gz (the classic Gi* z), pz, qz, pp (the permutation p-value, NaN without permutations), ze (the exposure-aware z), pe (its exact two-sided binomial p-value) and qe.

References

Getis, A. and Ord, J. K. (1992). The analysis of spatial association by use of distance statistics. Geographical Analysis, 24(3), 189-206. doi:10.1111/j.1538-4632.1992.tb00261.x

Ord, J. K. and Getis, A. (1995). Local spatial autocorrelation statistics: distributional issues and an application. Geographical Analysis, 27(4), 286-306. doi:10.1111/j.1538-4632.1995.tb00912.x

Examples

set.seed(1)
m <- matrix(rpois(100, 2), 10)
m[4:5, 4:5] <- 9
g <- gstr(m)
g[g$qe < 0.05, ]

Description

The chance that an unrelated person from the same population shares a DNA profile, locus by locus and overall. With theta = 0 this is the product rule (p^2 for a homozygote, ⁠2pq⁠ for a heterozygote). A positive theta applies the Balding-Nichols correction for shared ancestry. That never makes a homozygote match rarer than the product rule says; it can make a heterozygote match rarer when both alleles are common. See proof("gtmp").

Usage

gtmp(p, q = NULL, hom = FALSE, theta = 0)

Arguments

Value

A list: locus, a data frame with locus, hom and mp; log10mp, the base-10 logarithm of the profile match probability; and log10lr, the same with the sign changed, the weight of a single-source match against an unrelated alternative.

References

Balding, D. J. and Nichols, R. A. (1994). DNA profile match probability calculation. Forensic Science International, 64(2-3), 125-140. doi:10.1016/0379-0738(94)90222-4

Examples

gtmp(p = c(D3 = 0.1, vWA = 0.2), q = c(0.3, NA), hom = c(FALSE, TRUE),
     theta = 0.03)

Description

Smooths values at points or areal units (rates by area centroid, say) by Hansen's (1991) head-bang rule, which removes isolated peaks and troughs while keeping edges between regions of different level. For each unit i the rule takes its nn nearest neighbours and every pair ⁠(j, k)⁠ of them whose angle at i is at least ang degrees, so that i lies roughly between j and k (a triple); with ntrip > 0 only the ntrip triples whose ends are nearest (smallest ⁠|p_j - p_i| + |p_k - p_i|⁠) are kept. Each triple gives a low value min(v_j, v_k) and a high value max(v_j, v_k); the low screen is the median of the lows, the high screen the median of the highs, and the new value is the median of the low screen, v_i and the high screen, so v_i moves only if it lies outside both. With weights w (populations at risk, say) the screens are weighted medians, each end weighted by its own unit's weight, after Mungiole, Pickle and Simonson (1999): the smallest value whose cumulative weight reaches half the total. A unit with no triple (on the outer edge of the points) keeps its value; Hansen's extrapolated end triples are not formed. Each of iter passes starts from the previous pass's values.

Usage

hdbang(
  x,
  y = NULL,
  v,
  w = NULL,
  nn = 8,
  ang = 135,
  ntrip = 0,
  iter = 1,
  lonlat = FALSE
)

Arguments

Value

A numeric vector of smoothed values, one per point, in input order.

References

Hansen, K. M. (1991). Head-banging: robust smoothing in the plane. IEEE Transactions on Geoscience and Remote Sensing, 29(3), 369-378.

Mungiole, M., Pickle, L. W. and Simonson, K. H. (1999). Application of a weighted head-banging algorithm to mortality data maps. Statistics in Medicine, 18(23), 3201-3209.

Examples

g <- expand.grid(x = 1:6, y = 1:6)
v <- rep(1, 36)
v[15] <- 50
hdbang(g, v = v)[15]

Description

The distance along the sphere between pairs of points given in longitude and latitude. The central angle is computed as the atan2 of the cross and dot products of the two unit vectors, written so that nothing cancels: accurate to a few units in the last place from coincident points to antipodes (the spherical case of Vincenty 1975). The Earth is taken as a sphere of mean radius r; distances on the ellipsoid differ by up to about 0.5 per cent.

Usage

hvdist(from, to, r = 6371008.8)

Arguments

Value

A numeric vector of distances, one per pair, in the order of the input.

References

Vincenty, T. (1975). Direct and inverse solutions of geodesics on the ellipsoid with application of nested equations. Survey Review, 23(176), 88-93. doi:10.1179/sre.1975.23.176.88

See Also

eqproj() to project longitude and latitude to metres.

Examples

# Toronto to Montreal, metres
hvdist(c(-79.38, 43.65), c(-73.57, 45.50))
pts <- data.frame(lon = c(0, 90, 180), lat = c(0, 0, 0))
hvdist(pts, c(0, 0))

Description

A share s of a cohort is at high risk of reconviction, h a period, and the rest at low risk, l. If a programme changes who survives the first period without reconviction, the period-two hazard among the survivors is a weighted average of h and l with weight the share of high-risk survivors, and it rises with that share. An arm that keeps more of its high-risk people out of reconviction in period one therefore shows the higher period-two hazard even if it has no effect on anyone in period two: a period-by-period hazard ratio is not a causal contrast; a risk difference at a fixed follow-up is. See proof("hzdsel").

Usage

hzdsel(s, h, l, surh, surl)

Arguments

Value

A named numeric vector: share0, share1 (the high-risk share among the survivors of each arm), haz0, haz1 (the period-two hazard in each arm) and ratio (haz1 / haz0).

References

Hernan, M. A. (2010). The hazards of hazard ratios. Epidemiology, 21(1), 13-15. doi:10.1097/EDE.0b013e3181c1ea43

Examples

hzdsel(s = 0.5, h = 0.5, l = 0.1, surh = c(0.5, 0.8), surl = c(0.9, 0.9))

Description

Chouldechova's (2017) identity: in any confusion table, the false positive rate is p / (1 - p) (1 - ppv) / ppv (1 - fnr), from the base rate p, the positive predictive value ppv and the false negative rate fnr. So two groups whose scores are equally calibrated (the same ppv) and miss equally (the same fnr) but whose base rates differ must have different false positive rates: no score, however it is built, can equalise all three. The identity is arithmetic, true of any tool on any data; it says which equalities cannot hold together, not which one should give way. See proof("impfpr").

Usage

impfpr(p, ppv, fnr)

Arguments

Value

The false positive rates, one per triple.

References

Chouldechova, A. (2017). Fair prediction with disparate impact: a study of bias in recidivism prediction instruments. Big Data, 5(2), 153-163. doi:10.1089/big.2016.0047

Examples

# the same ppv and fnr, different base rates: different fpr
impfpr(p = c(a = 0.3, b = 0.5), ppv = 0.6, fnr = 0.25)

Description

An offender commits crimes at rate lam a year while free, and each crime leads to a custodial sentence of s years with probability q. Over a career the crime rate is ⁠lam / (1 + lam q s)⁠: the share of offending prevented, ⁠lam q s / (1 + lam q s)⁠, never reaches 1, longer or surer sentences lower the rate, and one more year of sentence prevents ⁠lam^2 q / ((1 + lam q s) (1 + lam q (s + 1)))⁠ crimes a year, less for each year added. A sentence uniform across offenders prevents a larger share of the offending of those with the higher rate. The model holds the rate fixed over the career and replaces nobody; see incapr() for both. See proof("incap").

Usage

incap(lam, q, s, share = NULL)

Arguments

Value

A data frame with one row per group: group, lam, q, s, rate (the long-run crime rate), prev (the share prevented) and marg (the crimes a year one more year of sentence prevents). With share, a last row all holds the weighted free rate in lam, the weighted rate in rate and the share of all offending prevented in prev.

References

Avi-Itzhak, B. and Shinnar, R. (1973). Quantitative models in crime control. Journal of Criminal Justice, 1(3), 185-217. doi:10.1016/0047-2352(73)90061-5

Examples

incap(lam = c(low = 2, high = 10), q = 0.1, s = 1, share = c(0.8, 0.2))

Description

Crimes prevented by a sentence of s periods served from career age t0, when the offending rate falls with career age and a share r of the crimes prevented is committed by someone else instead. The count is at most s times the rate at the start of the sentence and at least s times the rate at its end; the same sentence served one period later prevents no more; and, net of replacement, the estimate from the rate at the start is still an upper bound. The usual constant-rate estimate therefore overstates what a sentence buys under both assumptions. See proof("incapr").

Usage

incapr(lam, t0, s, r = 0)

Arguments

Value

A named numeric vector: prev (crimes prevented), net (net of replacement), upper (s * lam at age t0), lower (s * lam at age t0 + s) and later (crimes prevented by the same sentence from age t0 + 1).

References

Shinnar, S. and Shinnar, R. (1975). The effects of the criminal justice system on the control of crime: a quantitative approach. Law and Society Review, 9(4), 581-611. doi:10.2307/3053340

Examples

incapr(lam = c(12, 10, 8, 6, 5, 4, 3, 2, 2, 1), t0 = 2, s = 3, r = 0.25)

Description

Given a surface over the grid of jtsurf() and the offender's known base: the surface's value in the base's cell, the hit score (the share of the grid searched, in decreasing order of the surface, up to and including the base's cell), in the best and the worst order of ties (Rossmo 2000; Canter et al. 2000), and the error distance from the peak cell's centre to the base.

Usage

jtacc(prob, win, base)

Arguments

Value

A named numeric vector: p (the value in the base's cell), hitbest and hitworst (shares from 0 to 1), err (the error distance), peak and cell (the peak's and the base's cell, numbered down the columns).

References

Rossmo, D. K. (2000). Geographic Profiling. CRC Press.

Examples

s <- jtsurf(c(2, 3, 4), c(2, 4, 3), win = c(0, 6, 0, 6), nx = 30,
            ny = 30, type = "negexp", par = c(1, 1), mode = "lik")
jtacc(s$prob, win = c(0, 6, 0, 6), base = c(3.2, 2.9))

Description

A kernel density estimate of trip lengths from known offender trips, to calibrate a decay function (Levine 2009; Silverman 1986). With refl = TRUE the kernels are reflected at 0, so no mass falls below 0 and the estimate integrates to 1 over ⁠[0, Inf)⁠. Tabulated at evenly spaced distances it is a "table" decay for jtdecay() and jtsurf().

Usage

jtcal(d, t, h, w = NULL, kern = "normal", refl = TRUE)

Arguments

Value

A numeric vector of densities, one per value of t.

References

Silverman, B. W. (1986). Density Estimation for Statistics and Data Analysis. Chapman and Hall. doi:10.1201/9781315140919

Levine, N. (2009). Introduction to the special issue on Bayesian journey-to-crime modelling. Journal of Investigative Psychology and Offender Profiling, 6(3), 167-185. doi:10.1002/jip.107

Examples

set.seed(1)
trips <- rexp(200, 1 / 2)
jtcal(trips, t = seq(0, 10, 0.5), h = 0.5)

Description

The decay of the chance of offending with distance from an offender's base, in the forms of Brantingham and Brantingham (1981), Rossmo (2000) and Levine (2009): "linear", ⁠f = max(0, A + B d)⁠, par = c(A, B); "negexp", ⁠f = A exp(-b d)⁠, par = c(A, b) with b >= 0; "normal", ⁠f = A dnorm(d, mu, s)⁠, par = c(A, mu, s); "lognormal", ⁠f = A dlnorm(d, mu, s)⁠, par = c(A, mu, s) (the lognormal density of d; some texts use ⁠log d^2⁠ and 1 / d^2, which is not a density); "trunc", the truncated negative exponential, ⁠f = P d / dp⁠ up to the peak dp and ⁠P exp(-C (d - dp))⁠ beyond, par = c(P, dp, C); and "table", values tab at distances ⁠0, delta, 2 delta, ...⁠, linear in between and 0 beyond, par = delta (as jtcal() tabulates).

Usage

jtdecay(d, type, par, tab = NULL)

Arguments

Value

A numeric vector of f(d), one per distance.

References

Brantingham, P. J. and Brantingham, P. L. (1981). Environmental Criminology. Sage.

Levine, N. (2009). Introduction to the special issue on Bayesian journey-to-crime modelling. Journal of Investigative Psychology and Offender Profiling, 6(3), 167-185. doi:10.1002/jip.107

Examples

jtdecay(c(0, 0.5, 1, 2, 5), type = "trunc", par = c(1, 0.5, 0.8))

Description

Ranks the cells of a grid as candidate bases of an offender from the locations of a series of incidents. mode = "sum" gives the classic aggregate surface, the sum over incidents of the decay of their distance to the cell (Canter et al. 2000; Levine 2009); mode = "lik" gives the sum of the logarithms, the log-likelihood of the cell as base when the decay function is the planar density of an incident about the base and the incidents are independent (O'Leary 2009), so that prob is the posterior over cells under a uniform prior. The sum surface is not a likelihood.

Usage

jtsurf(
  x,
  y = NULL,
  win,
  nx,
  ny,
  type,
  par,
  tab = NULL,
  metric = "euclid",
  mode = "sum"
)

Arguments

Details

What this does not license: the surface ranks places to search under the assumptions that the series is one offender's, that the offender travelled from one fixed base inside the grid, and that the decay function, which describes offenders in general, fits this one. It is a search priority, not evidence that any person lives anywhere. Accuracy on solved series in the same place (jtacc()) is the only guide to its worth.

Value

A list. raw and prob, matrices of nx rows (across) and ny columns (up), cell ⁠[i, j]⁠ centred at ⁠(xmin + (i - 1/2) dx, ymin + (j - 1/2) dy)⁠: the surface and its normalised version (summing to 1). peak, a named vector: cell (the first cell of maximum value, numbered down the columns of the matrix), x and y (its centre).

References

Canter, D., Coffey, T., Huntley, M. and Missen, C. (2000). Predicting serial killers' home base using a decision support system. Journal of Quantitative Criminology, 16(4), 457-478. doi:10.1023/A:1007551316253

O'Leary, M. (2009). The mathematics of geographic profiling. Journal of Investigative Psychology and Offender Profiling, 6(3), 253-265. doi:10.1002/jip.111

Examples

s <- jtsurf(c(2, 3, 4), c(2, 4, 3), win = c(0, 6, 0, 6), nx = 30,
            ny = 30, type = "negexp", par = c(1, 1), mode = "lik")
s$peak

Description

The weight of two DNA profiles for a stated relationship against no relationship, locus by locus and overall. Loci are taken as independent; there is no allowance for mutation or coancestry. See proof("kinidx").

Usage

kinidx(g1, g2, freq, rel = "parent")

Arguments

Value

A list: locus, a data frame with locus and idx; log10idx, the base-10 logarithm of the combined index; and post, the probability of the relationship at even prior odds.

Examples

kinidx(g1 = c(12, 14), g2 = c(12, 14), freq = c(0.1, 0.2), rel = "sibling")

Description

Partitions the points into k groups by Lloyd's iterations: each point joins its nearest centre (the lowest-numbered on ties), then each centre moves to its points' mean (upd = "mean", K-means, minimising the within-group sum of squares) or to their centre of minimum distance (upd = "weber", a spatial K-medians minimising the sum of distances), until no point changes group or after maxit moves. An emptied group keeps its centre. The result depends on the start: k-means++ seeding (Arthur and Vassilvitskii 2007), the centres of the busiest grid cells at least a given distance apart, or centres given.

Usage

kmpart(
  x,
  y = NULL,
  k,
  init = "pp",
  cell = NULL,
  upd = "mean",
  maxit = 100,
  seed = NULL,
  count = NULL,
  lonlat = FALSE
)

Arguments

Value

A list. cl, the group of each point (1 to k, input order). cen, a data frame with one row per group: x, y (its centre, in degrees when lonlat = TRUE), n (points), ss (sum of squared distances to the centre) and sd (sum of distances). fit, a named vector: ss and sd (totals), iter (moves made) and conv (1 if no point changed group at the end).

References

Lloyd, S. (1982). Least squares quantization in PCM. IEEE Transactions on Information Theory, 28(2), 129-137. doi:10.1109/TIT.1982.1056489

Arthur, D. and Vassilvitskii, S. (2007). k-means++: the advantages of careful seeding. Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms, 1027-1035.

Examples

set.seed(1)
x <- c(rnorm(30, 0), rnorm(30, 5))
y <- c(rnorm(30, 0), rnorm(30, 5))
kmpart(x, y, k = 2, seed = 1)$cen

Description

Maximum likelihood for a 0/1 outcome and a logit link, by Newton's method on the design with centred columns, the method glm() uses with family = binomial() when the columns are not collinear. Standard errors come from the inverse Fisher information.

Usage

lgtreg(x, y, tol = 1e-10, maxit = 50)

Arguments

Value

A list: coef and se, named numeric vectors with the intercept first, and dev, the residual deviance.

References

McCullagh, P. and Nelder, J. A. (1989). Generalized Linear Models (2nd ed.). Chapman and Hall. doi:10.1007/978-1-4899-3242-6

Examples

set.seed(1)
x <- matrix(rnorm(300 * 2), 300)
y <- rbinom(300, 1, plogis(0.5 + x[, 1] - x[, 2]))
lgtreg(x, y)$coef

Description

Dwork et al. (2012) ask that similar individuals be treated similarly: that the gap between two people's scores be at most lip times the distance between them, for a distance that says who is alike for the task. lipaud() checks the given pairs against that rule and reports the smallest constant the scores satisfy on them.

Usage

lipaud(h, d = NULL, pairs = NULL, lip = 1)

Arguments

Details

A pair at distance 0 is two people the distance cannot tell apart: any gap between their scores breaks the rule, and no constant repairs it (its ratio is Inf). The audit is only as good as the distance: it says whether the scores treat alike those the distance calls alike, not whether the distance is the right one.

Value

A list: lreq (the smallest constant the scores satisfy on these pairs), nviol and rate (the pairs breaking the rule at lip, and their share), ok (TRUE when none does), maxgap (the largest gap), worst (the first pair with the largest ratio), and pairs, a data frame with one row per pair: i, j, d, gap, ratio and viol.

References

Dwork, C., Hardt, M., Pitassi, T., Reingold, O. and Zemel, R. (2012). Fairness through awareness. Proceedings of the 3rd Innovations in Theoretical Computer Science Conference, 214-226. doi:10.1145/2090236.2090255

Examples

x <- c(1, 2, 2.5, 6)
h <- c(0.10, 0.20, 0.60, 0.65)
# scores on a scale of 0 to 1 against distances in years
lipaud(h, dist(x), lip = 0.2)

Description

Anselin's (1995) local indicator of spatial association for each unit: with z = x - mean(x) and m2 = mean(z^2), ⁠I_i = z_i / m2 sum_j w_ij z_j⁠, its expectation and variance under total randomisation, the two-sided normal p-value, and a conditional permutation p-value: unit i keeps its value and its neighbours' values are drawn without replacement from the others, nsim times from seed. The quadrant tells high-high and low-low clusters from high-low and low-high outliers. The p-values are per unit and not adjusted for the many tests made (use stats::p.adjust()); with row-standardised weights the sum of the I_i is S0 times global Moran's I. That identity is proved in Lean: see proof("lmorn").

Usage

lmorn(x, w, nsim = 0, seed = NULL)

Arguments

Value

A data frame with one row per unit: ii (I_i), e, v, z, p, pperm (⁠(1 + min(#{I* >= I}, #{I* <= I})) / (nsim + 1)⁠, NaN without permutations) and quad (1 high-high, 2 low-low, 3 low-high, 4 high-low, 0 neither).

References

Anselin, L. (1995). Local indicators of spatial association: LISA. Geographical Analysis, 27(2), 93-115. doi:10.1111/j.1538-4632.1995.tb00338.x

Examples

set.seed(1)
p <- cbind(x = runif(60), y = runif(60))
v <- ifelse(p[, 1] > 0.7 & p[, 2] > 0.7, 5, 1) + rnorm(60, sd = 0.3)
l <- lmorn(v, spwt(p, k = 5), nsim = 199, seed = 1)
head(l)

Description

The local statistics of Getis and Ord (1992) as z values (Ord and Getis 1995): G_i compares the weighted sum of the neighbours' values with its expectation from the other n - 1 values; ⁠G*_i⁠ (star = TRUE) includes the unit itself with weight 1 and uses all n. Large positive values mark concentrations of high values, large negative ones of low values. With star = TRUE the weights are best binary (style = "raw").

Usage

locgi(x, w, star = TRUE)

Arguments

Value

A numeric vector of z values, one per unit; NaN for a unit without neighbours.

References

Ord, J. K. and Getis, A. (1995). Local spatial autocorrelation statistics: distributional issues and an application. Geographical Analysis, 27(4), 286-306. doi:10.1111/j.1538-4632.1995.tb00912.x

Examples

set.seed(1)
p <- cbind(x = runif(60), y = runif(60))
v <- rpois(60, ifelse(p[, 1] > 0.7, 6, 2))
locgi(v, spwt(p, type = "band", d = 0.2, style = "raw"))

Description

The mean length of stay of the people released (or of the cases disposed) leaves out everyone still inside. Their elapsed time is a lower bound on their stay, so the mean over everyone admitted is at least lower, the mean of the completed stays and the elapsed times together; the released-only mean falls short of it by m / (n + m) times the excess of the mean elapsed time over the released mean, and understates the truth for certain when that excess is positive. No upper bound exists without a cap on how long the current stays will last. See proof("losbnd").

Usage

losbnd(done, open)

Arguments

Value

A named numeric vector: mean (of the completed stays), age (the mean elapsed time of the open ones), lower (the lower bound on the mean over all) and bias (lower - mean, the least by which the completed mean understates when it is positive).

References

Manski, C. F. (2003). Partial Identification of Probability Distributions. Springer. doi:10.1007/b97478

Examples

losbnd(done = c(30, 45, 60, 90, 120), open = c(100, 150, 200))

Description

Recalibrates a set of likelihood ratios against the truth with the pool-adjacent-violators algorithm, keeping their order. The result can hold 0 and Inf: the data never contradicted those ratios. See proof("lrcal").

Usage

lrcal(same, diff)

Arguments

Value

A list: same and diff, the calibrated ratios in the order given.

Examples

lrcal(same = c(50, 200, 8), diff = c(0.01, 0.2, 1.5))

Description

Turns a likelihood ratio into the probability of the first hypothesis, which needs a prior. A large ratio with a small prior can still leave that probability low; reading the ratio as the probability is the prosecutor's fallacy. See proof("lrpost").

Usage

lrpost(lr, prior)

Arguments

Value

The posterior probabilities.

Examples

lrpost(lr = 1e6, prior = 1 / 5e6)

Description

Robust locally weighted regression (Cleveland 1979) with the conventions of the classic implementation, so that the fit equals that of stats::lowess(); the fitted values come back in the order of the input, not sorted.

Usage

lwss(x, y, f = 2/3, iter = 3, delta = NULL)

Arguments

Value

A numeric vector of fitted values, one per point in input order.

References

Cleveland, W. S. (1979). Robust locally weighted regression and smoothing scatterplots. Journal of the American Statistical Association, 74(368), 829-836. doi:10.1080/01621459.1979.10481038

Examples

lwss(c(3, 1, 2, 5, 4), c(2, 1, 3, 6, 4))

Description

The Pearson correlation, over all pairs of events, between their spatial dissimilarity (the distance, or 1 / (distance + cs) when cs > 0, Mantel's reciprocal transform) and their temporal one (⁠|t_i - t_j|⁠, or ⁠1 / (|t_i - t_j| + ct)⁠), with a permutation test that permutes the times over the fixed locations (Mantel 1967). Drawing distances and intervals independently instead, as some descriptions do, does not keep the dependence among the pairs of one distance matrix and gives no test of stated size. Memory grows with the number of events, not pairs; time with pairs times permutations.

Usage

mantst(
  x,
  y = NULL,
  t,
  cs = 0,
  ct = 0,
  nsim = 999,
  seed = NULL,
  metric = "euclid",
  lonlat = FALSE
)

Arguments

Value

A named numeric vector: r (the correlation), pairs, and the permutation p-values pup (⁠(1 + #{r_perm >= r}) / (nsim + 1)⁠, large r: events close in space are close in time when both are raw distances), plo (r_perm <= r) and p2 (⁠|r_perm| >= |r|⁠); NaN without permutations. 3 to 46340 events.

References

Mantel, N. (1967). The detection of disease clustering and a generalized regression approach. Cancer Research, 27(2), 209-220.

Examples

set.seed(1)
x <- runif(60)
y <- runif(60)
mantst(x, y, t = x + rnorm(60, sd = 0.1), nsim = 199, seed = 1)

Description

When the group a person is recorded in (an officer's perception of race, say) differs at random from the true one, with the same error rates whatever the outcome, the measured groups are mixtures of the true ones. The measured risk difference is then the true one times the factor ⁠(se + sp - 1) p (1 - p) / (d1 d0)⁠, where d1 and d0 are the shares recorded in each group. For any sensitivity and specificity it is never larger in size than the true difference, so a disparity in recorded data understates the true one; with se + sp >= 1 it never has the opposite sign, while with se + sp < 1 a difference other than 0 has it. See proof("misrd").

Usage

misrd(r1, r0, prev, se, sp)

Arguments

Value

A named numeric vector: true (r1 - r0), measured (the difference between the recorded groups), factor (measured / true whenever true is not 0), d1 and d0.

References

Bross, I. (1954). Misclassification in 2 x 2 tables. Biometrics, 10(4), 478-486. doi:10.2307/3001619

Examples

misrd(r1 = 0.3, r0 = 0.1, prev = 0.4, se = 0.8, sp = 0.9)

Description

The weight of a DNA mixture at one locus for one account of who contributed against another. Each account names the alleles of the mixture its known contributors do not explain and the number of unknown contributors who must explain them. Peak heights, drop-out and coancestry are not modelled. See proof("mixlr").

Usage

mixlr(p, np, xp, nd, xd)

Arguments

Details

With no alleles to explain and one unknown, the denominator is the chance a random person is not excluded from the mixture.

Value

A named numeric vector: num, den and lr.

References

Weir, B. S., Triggs, C. M., Starling, L., Stowell, L. I., Walsh, K. A. J. and Buckleton, J. (1997). Interpreting DNA mixtures. Journal of Forensic Sciences, 42(2), 213-222.

Examples

# victim explains a and b; suspect would explain c and d
p <- c(a = 0.1, b = 0.2, c = 0.05, d = 0.15)
mixlr(p, np = character(0), xp = 0, nd = c("c", "d"), xd = 1)

Description

The group-based trajectory model with several normal indicators per unit (Nagin 1999; Magrini 2022): given group k, the r indicators of a unit at period t are multivariate normal with mean ⁠B_k' x_t⁠ and covariance S_k, independently over periods. mvtll() evaluates the log-likelihood and posterior memberships; mvtem() fits by EM with closed-form steps from the values given. A period with any missing indicator is left out for that unit.

Usage

mvtll(y, x, beta, sig, prob)

mvtem(y, x, beta, sig, prob, maxit = 500, tol = 1e-08)

Arguments