# rStats.js

> R-compatible statistics for JavaScript and TypeScript: distributions, random numbers that match
> R's `set.seed()`, descriptive statistics, hypothesis tests, `lm`, `glm` and lme4's `lmer` and `glmer`, special functions,
> numerical methods and linear algebra, with R's function names, defaults and numerical results.
> Every function is tested against R 4.6.2. ESM only, no dependencies, runs in browsers, Node.js,
> Deno and Bun.

Install with `npm install @numero-uno-jpn/rstats`. All functions are named exports:
`import { pnorm, tTest } from "@numero-uno-jpn/rstats"`, or `import * as rs from
"@numero-uno-jpn/rstats"` and call `rs.tTest(x, y)`. Each area is also a subpath export
(`/distributions`, `/tests`, `/linalg`, `/descriptive`, `/numerics`, `/special`). Mixed models are only in the
subpath `/mixed`: `import { lmer, glmer } from "@numero-uno-jpn/rstats/mixed"` (not in the main entry point). For use from
a CDN without a bundler, see the README section "Using it in the browser".

## Translating R code

- Function names are R's names. A dot in R's name becomes camelCase: `t.test` → `tTest`,
  `chisq.test` → `chisqTest`, `p.adjust` → `pAdjust`, `set.seed` → `setSeed`,
  `sample.int` → `sampleInt`, `weighted.mean` → `weightedMean`, `pairwise.t.test` →
  `pairwiseTTest`. Methods of `glm` carry a prefix: `summary(fit)` → `glmSummary(fit)`,
  `predict` → `glmPredict`, `confint` → `glmConfint`, `anova` → `glmAnova`, `residuals` →
  `glmResiduals`. Other renamed functions: `var` → `variance`, `%*%` → `matMul`,
  `.Random.seed` → `getRandomSeed()` / `setRandomSeed(seed)`, `qr.Q` → `qrQ` (likewise `qrR`,
  `qrCoef`, `qrResid`, ...). Functions from other packages keep their own names (`zTest` from
  BSDA, `brunnerMunzelTest` from lawstat, `ginv` from MASS, `grad` and `hessian` from numDeriv).
- Required arguments are positional; optional ones go in one options object with camelCase keys:
  `pt(q, df, lower.tail = FALSE, log.p = TRUE)` → `pt(q, df, { lowerTail: false, logP: true })`,
  `t.test(x, y, var.equal = TRUE, conf.level = 0.99)` →
  `tTest(x, y, { varEqual: true, confLevel: 0.99 })`. When an optional second sample is
  omitted, the options take its place: `tTest(x, { mu: 1 })`.
- Option values are R's strings: `alternative: "less"` (unambiguous abbreviations such as `"l"`
  are accepted, as with R's `match.arg`), `method: "kendall"`, `family: "binomial"`.
- Results are plain objects with camelCase fields. Hypothesis tests return R's `htest` fields:
  `statistic`, `parameter`, `pValue` (R's `p.value`), `confInt` (`conf.int`) with `confLevel`,
  `estimate`, `nullValue`, `method`, `alternative`. R's names of these values are separate
  fields: `statisticName`, `parameterName`, `estimateNames`, `nullValueName`.
- Data shapes: a sample is an array of numbers or a `Float64Array`. A matrix or contingency table
  is an array of rows. Grouped data is either an array of samples (`kruskalTest([a, b, c])`) or
  values with a parallel array of group labels (`kruskalTest(x, g)`). The first argument of the
  d/p/q functions may be a number or an array (the result has the same shape); distribution
  parameters are scalars and are not recycled as vectors.
- Models: `lm` and `aov` take an R formula string and a data frame given as an object of columns:
  `lm("y ~ A * B + x", { y, A, B, x })`, `aov("y ~ W + Error(id/W)", data)`; string columns are
  factors (levels sorted by code point unless `levels: { A: [...] }`), `contrasts: { A: "contr.sum" }`
  sets the coding. Formula functions such as `log(x)` or `I(x^2)` are not supported: compute the
  column first. Tables (`anova(fit)`, `Anova(fit, { type: 3 })` for car's `Anova`, `aovSummary(fit)`)
  are `{ colNames, rowNames, table }` with R's names. afex's `aov_ez(id, dv, data, between, within)`
  is `aovEz(data, { id, dv, between, within })`. `lm(x, y)` also takes a matrix, and `glm(x, y, {
  family })` takes a matrix (array of rows) or a vector of predictors (no formula); the intercept
  is added unless `intercept: false`. A binomial response can be rows of `[successes, failures]`,
  like `cbind(s, f)` in R.
- Mixed models: lme4's `lmer(y ~ x + (x | subject), data = d, REML = FALSE)` →
  `lmer("y ~ x + (x | subject)", d, { REML: false })` from `@numero-uno-jpn/rstats/mixed`, with `(1 | a/b)`,
  `(1 | a:b)`, `(x || g)` and crossed terms. The result has the fields `fixef`, `coefTable` (`{ estimate, stdError,
  tValue }`), `vcov`, `sigma`, `varCorr` (one `{ group, columns, sd, corr }` per term, as
  `VarCorr`), `ranef` (one `{ group, levels, columns, values }` per grouping factor), `fitted`, `residuals`,
  `logLik`, `AIC`, `BIC`, `deviance`, `REMLcrit`, `theta` and `singular`. The estimates are the exact optimum, so
  the variance components differ from lme4's in the 5th to 9th significant digit.
- lmerTest: `summary(fit)` of `lmerTest::lmer` → `lmerSummary(fit)` (`{ names, estimate, stdError, df, tValue,
  pValue }`, Satterthwaite's df); `contest(fit, L)` → `lmerContest(fit, L, { joint, rhs, level })`;
  `anova(fit, type = "II")` → `lmerAnova(fit, { type: "II" })` (`{ terms, sumSq, meanSq, numDF, denDF, F, pValue }`,
  Type III by default); lme4's `anova(m1, m2)` → `lmerAnova([m1, m2], { names })` (likelihood-ratio tests; REML fits
  are refitted by ML). Pass the object `lmer` returned (not a copy). Kenward-Roger is not available.
- Generalized linear mixed models: lme4's `glmer(y ~ x + (1 | g), data = d, family = poisson, nAGQ = 1)` →
  `glmer("y ~ x + (1 | g)", d, { family: "poisson", nAGQ: 1 })` (families binomial, poisson, Gamma, inverse.gaussian,
  gaussian with a link other than identity; `nAGQ` 2 to 100 for one scalar term, binomial and poisson only; a binomial
  response may be `cbind(successes, failures)` on the left-hand side). `summary(fit)` → `glmerSummary(fit)`;
  `anova(fit)` → `glmerAnova(fit)`, `anova(m1, m2)` → `glmerAnova([m1, m2], { names })` (glmer, lmer and glm fits);
  `predict(fit, newdata, re.form, type, se.fit)` → `glmerPredict(fit, { newdata, reForm, type, seFit })` and
  `lmerPredict(fit, ...)`. Pass the object `glmer` returned (not a copy). For `nAGQ > 1` the log-likelihood is the
  true one (R leaves out a constant).
- Indices are 0-based where R returns indices: `sampleInt(n)` returns `0..n-1` (R's
  `sample.int(n) - 1` from the same random stream), and the pivots of `qr` and `chol`.
- Contingency tables and proportions: DescTools' `BinomCI(x, n, method = "wilson")` → `binomCI(x, n, { method:
  "wilson" })` (also `BinomDiffCI` → `binomDiffCI(x1, n1, x2, n2, ...)`, `BinomRatioCI` → `binomRatioCI`, `GTest` →
  `gTest`, `StuartMaxwellTest` → `stuartMaxwellTest`, `GoodmanKruskalGamma` → `goodmanKruskalGamma`), epitools'
  `oddsratio(m)` → `oddsRatio(m)`, exact2x2's `fisher.exact(m, midp = TRUE)` → `fisherTest(m, { midp: true })`,
  Exact's `exact.test(m, method = "boschloo")` → `boschlooTest(m)`, coin's `lbl_test(table)` → `lblTest(table)`,
  equivalence's `tost(x, y, epsilon = 1)` → `tost(x, y, { epsilon: 1 })`; the 'N − 1' chi-squared test is
  `chisqTest(m, { nMinus1: true })` and `binomDiffEquivalenceTest(x1, n1, x2, n2, { margin })` is the score TOST of
  two proportions. Where these packages' root searches are inexact, rStats.js returns the true roots.
- R's warnings: every test, model and fit returns them in its result (`warnings: string[]`, R's text, each distinct
  warning once, `[]` when there are none). Nothing is printed by default. `setWarningHandler((message, source) =>
  console.warn(source + ": " + message))` prints each warning of a call once, as R does (`source` is the function
  called; a fit's warnings are passed when it is fitted, not again by `anova(fit)`); a handler that throws plays the
  part of `options(warn = 2)`; `setWarningHandler(null)` stops it and the call returns the previous handler.
- Errors: input for which R returns NaN, NA or a meaningless value because the result is not defined (`mcnemarTest`
  on a table without discordant pairs, a non-finite value where a finite one is needed, too few observations or
  groups) throws an `Error` or `RangeError`; wrap such calls in `try`/`catch` where R code checks `is.na()` on the result. The
  distribution functions (d/p/q/r) return NaN for invalid parameters, as R. A computation that would allocate more
  than 2^24 values in one array or run for more than about 10 seconds (large exact tests, huge simulations) throws a
  RangeError before it starts; for exact tests the message suggests `exact: false`.
- Power calculations take one options object and solve for the quantity given as `null`:
  `power.t.test(delta = 0.5, power = 0.8)` → `powerTTest({ delta: 0.5, power: 0.8 })`, pwr's
  `pwr.anova.test(k = 3, f = 0.25, power = 0.8)` → `pwrAnovaTest({ k: 3, f: 0.25, power: 0.8 })`
  (`sig.level` → `sigLevel`). effectsize's `cohens_d(x, y)` → `cohensD(x, y)` (`{ name, value, ciLow, ciHigh }`),
  `eta_squared(model)` → `etaSquared(anova(fit))` or `etaSquared(Anova(fit, { type: 3 }))`.
- Cox regression takes the covariates as a matrix (array of rows): `coxph(Surv(time, status) ~ x1 + x2, ties =
  "breslow", weights = w)` → `coxph(time, status, x, { ties: "breslow", weights: w })`, with `strata` as labels.
  psych's `omega(m, nfactors = 3)` → `omega(x, { nfactors: 3 })` (`omega(null, { covmat: R })` for a correlation
  matrix).
- Robust standard errors of `lm` and `glm` fits: sandwich's `vcovHC(fit, type = "HC0")` → `vcovHC(fit, { type:
  "HC0" })` (HC3 by default), `vcovCL(fit, cluster = ~ id)` → `vcovCL(fit, { cluster: d.id })` (one label per input
  row; several clustering dimensions as an array of label arrays); both return `{ vcov, names, warnings }`. lmtest's
  `coeftest(fit, vcov. = vcovHC(fit, type = "HC0"))` → `coeftest(fit, { vcov: vcovHC(fit, { type: "HC0" }) })` (`vcov`
  is a matrix, such a result or a function of the fit) and `coefci(fit, parm, level, vcov.)` → `coefci(fit, { parm,
  level, vcov })` (`parm` 0-based); a glm gets z tests, as lmtest. A modified Poisson regression is `glm(x, y, {
  family: "poisson" })` of a 0/1 outcome with `coefci(fit, { vcov: vcovHC(fit, { type: "HC0" }) })`, exponentiated.
  survival's `coxph(..., cluster = id)` → `coxph(time, status, x, { cluster: id })` (also `robust: true`): `var` is
  then the robust variance and `se` stays R's `se(coef)`.
- Diagnostic accuracy and agreement: pROC's `roc(response, predictor)` → `roc(response, predictor)` (also `roc("y ~
  x", d)` and `roc({ controls, cases })`), `auc(r, partial.auc = c(1, 0.8))` → `auc(r, { partialAuc: [1, 0.8] })`,
  `ci.auc(r)` → `ciAuc(r)` (DeLong; `{ method: "bootstrap" }`), `roc.test(r1, r2)` → `rocTest(r1, r2)`, `coords(r,
  "best")` → `coords(r, "best")` (column arrays); epiR's `epi.tests(tab)` → `diagnosticTests(tab)`; BlandAltmanLeh's
  `bland.altman.stats(x, y)` → `blandAltman(x, y)`; mcr's `mcreg(x, y, method.reg = "PaBa")` → `mcreg(x, y, {
  methodReg: "PaBa" })`; psych's `ICC(m)` → `ICC(m)` (irr's `icc` is a separate function); TOSTER's `t_TOST(x, y, eqb =
  0.5)` → `tTOST(x, y, { eqb: 0.5 })` and `tsum_TOST(m1, sd1, n1, m2, sd2, n2, eqb = 0.5)` → `tsumTOST({ m1, sd1, n1,
  m2, sd2, n2 }, { eqb: 0.5 })`.
- Data frames are objects of columns (`null` or `NaN` for NA); `fromRows(rows)` converts an array of row objects
  (`[{ age: 30, score: 5 }, ...]`) and `toRows(d)` converts back. Tests take formulas like R's:
  `t.test(y ~ g, data = d)` → `tTest("y ~ g", d)`, `friedman.test(y ~ w | b, d)` → `friedmanTest("y ~ w | b", d)`,
  `cor.test(~ x + y, d)` → `corTest("~ x + y", d)`. `na.omit(d)` → `naOmit(d)`, `complete.cases(d)` →
  `completeCases(d)`, `subset(d, age >= 20, select = c(y, g))` → `subset(d, (r) => r.age >= 20, { select: ["y",
  "g"] })`, `tapply(x, g, mean)` → `tapply(x, g, mean)` (`{ dim, dimnames, values }`), `split(x, g)` → a `Map`,
  `aggregate(y ~ g + sex, d, mean)` → `aggregate("y ~ g + sex", d, mean)`, `merge(a, b, by = "id", all.x = TRUE)`
  → `merge(a, b, { by: "id", allX: true })`, `rbind(d1, d2)` → `rbind(d1, d2)`. tidyr's `pivot_longer(d, cols =
  c(t1, t2), names_to = "time", values_to = "score")` → `pivotLonger(d, { cols: ["t1", "t2"], namesTo: "time",
  valuesTo: "score" })`, `pivot_wider(d, names_from = time, values_from = score)` → `pivotWider(d, { namesFrom:
  "time", valuesFrom: "score" })`.
- Random numbers: after `setSeed(42)`, `rnorm`, `runif`, `sample` and the other `r*` functions
  give the same numbers as R after `set.seed(42)`. `RNGkind` supports Mersenne-Twister (R's
  default) and L'Ecuyer-CMRG.

## Examples

```ts
import { pnorm, pt, setSeed, rnorm, tTest, glm, glmSummary } from "@numero-uno-jpn/rstats";

pnorm(1.96);                                  // ≈ 0.9750021049
pt(2, 10, { lowerTail: false, logP: true });  // ≈ -3.305141551

setSeed(42);
rnorm(3);  // [1.3709584471466685, -0.5646981713960887, 0.3631284113373392], as in R

const r = tTest([5.1, 4.9, 6.2, 5.8, 6.0, 5.5], [4.1, 4.5, 5.0, 4.8, 4.3]);
r.pValue;  // ≈ 0.003536530888 (Welch, R's default)
r.confInt; // ≈ [0.4423253335, 1.644341333]

// R: glm(cbind(s, f) ~ dose, family = binomial)
const fit = glm([1, 2, 3, 4, 5, 6], [[1, 19], [4, 16], [9, 11], [13, 7], [18, 2], [19, 1]], { family: "binomial" });
fit.coefficients;        // ≈ [-3.77413, 1.14807]
glmSummary(fit).pValue;  // ≈ [5.357e-8, 3.579e-9]
```

## Docs

- [README](https://raw.githubusercontent.com/numero-uno-jpn/rStats.js/main/README.md): install, examples, the full list of functions by area
- [Differences from R](https://raw.githubusercontent.com/numero-uno-jpn/rStats.js/main/docs/differences-from-r.md): every known difference from R, function by function
- [Accuracy](https://raw.githubusercontent.com/numero-uno-jpn/rStats.js/main/docs/accuracy.md): how results are checked against R and high-precision references
- [API reference](https://numero-uno-jpn.github.io/rStats.js/api/): every function with its options, results and an example
- [R functions](https://raw.githubusercontent.com/numero-uno-jpn/rStats.js/main/docs/r-functions.md): the R function each rStats.js function corresponds to
- [Guide for JavaScript developers](https://raw.githubusercontent.com/numero-uno-jpn/rStats.js/main/docs/guide-for-js-developers.md): data frames, missing values, factors, formulas and a worked analysis
- [Migrating from jStat](https://raw.githubusercontent.com/numero-uno-jpn/rStats.js/main/docs/migrating-from-jstat.md): the replacement for every jStat function

## Optional

- [Performance](https://raw.githubusercontent.com/numero-uno-jpn/rStats.js/main/docs/performance.md): speed and bundle sizes
- [Roadmap](https://raw.githubusercontent.com/numero-uno-jpn/rStats.js/main/ROADMAP.md): what is planned
- [Changelog](https://raw.githubusercontent.com/numero-uno-jpn/rStats.js/main/CHANGELOG.md)
- [Demo](https://numero-uno-jpn.github.io/rStats.js/demo/): every function in the browser, compared with R