Changes for version 0.319 - 2026-09-22
- glm(): offsets, prior weights, robust covariance, absorbed factors
- glm() took only formula, data, family, theta and conf.level, and died on anything else, which left count models unable to put a rate on person-time. It now takes offset() terms in the formula and an offset argument (a column, an expression such as 'log(t)', or an array ref), and the negative-binomial theta search and the null deviance both see the offset, as MASS::glm.nb()'s do.
- weights are R's prior weights. A binomial fit whose weights make a non-integer number of successes warns as R does.
- vcov => 'HC0' .. 'HC3' gives sandwich::vcovHC()'s covariance and cluster gives vcovCL()'s, one to four ways ('firm + year'), with every standard error, z, p-value and interval recomputed from it. A poisson fit on a 0/1 outcome with HC0 is the modified-Poisson risk ratio.
- A factor after '|' in the formula, or in absorb, is absorbed by weighted within-group demeaning instead of expanded into dummy columns, so a factor with thousands of levels costs a vector per level rather than a column. Groups whose outcome sits at a boundary are dropped, as fixest::feglm() drops them, and fe.removed counts their rows; only rows with a positive weight decide that, so a zero-weight count cannot keep an all-zero group in the fit.
- The IRLS loop was rewritten around one core shared with the new models below, with maxit and epsilon exposed and R's convergence rule and penultimate-iterate standard errors. New results: loglik, vcov, vcov.type, dispersion, nobs and, where they apply, n.clusters, absorb, fe.removed, twologlik and SE.theta. predict() re-evaluates an offset on new rows, and croaks on a model it cannot predict from (an offset given as an array, or absorbed factors) rather than leaving the term out.
- The formula reader now evaluates log(), exp(), sqrt(), log2(), log10(), log1p() and abs() of a column, which offsets and ivreg's examples need.
- Validated in t/glm_offset_weights.R.t, t/glm_vcov.R.t and t/glm_absorb.R.t against R's glm, MASS, sandwich and fixest test suites, statsmodels' and Stata's pinned results, and a full-dummy fit.
- glm(): Inf, a NaN in step-halving, and large negative-binomial counts
- An infinite response, covariate, offset or weight is not missing, so its row went into the fit and every coefficient came back NaN without a word. glm() now croaks "NA/NaN/Inf in 'y'" (or 'x', 'offset', 'weights') with the row name, as R's Cdqrls() stops.
- A coefficient aliased on one IRLS iteration was carried into the next as NaN, so if the column stopped being aliased and the step was halved, the NaN reached the linear predictor. It is carried as 0, as glm.fit() does.
- The negative-binomial log-likelihood summed y logs per row for every count below 1e6, so a fit cost O(sum y): 20000 counts near 34000 took 11.7 s against 0.01 s for the Poisson fit. Counts from 64 up now go through Loader's stirlerr(), and the fit takes 0.03 s. It is also more accurate: 3e-13 relative against mpmath at 60 digits, where the long sum was 1.9e-9 out. t/glm.t pins it.
- X'WX is accumulated over its upper triangle only and skips zero-weight rows, halving the IRLS loop's main cost and leaving the matrix exactly symmetric; coefficients move by about 1e-12. Absorbed factors are demeaned all columns in one pass over the rows, bit for bit as before, and vcovHC() no longer holds an n x p copy of the scores.
- coxph(): counting-process data, strata, robust variance, a formula
- coxph() took only (\@time, \@status, covariates), which cannot express a time-varying covariate or late entry. It now takes (start, stop] data, strata, a cluster with the grouped-dfbeta robust variance, case weights and an offset, either as options to the positional form or through formula => 'Surv(start, stop, event) ~ x + strata(g) + cluster(id)' over a data set. Intervals that span no event are skipped, as survival's agreg.fit skips them. New results: var and vcov, the score and Wald tests, and naive.se, naive.var, robust.score.test and n.clusters under a robust variance.
- Validated in t/coxph_extended.R.t against survival's own tests (bladder, the phreg corpora) and statsmodels'.
- New models: zerotrunc, hurdle, svyglm, ivreg, lmer
- zerotrunc() is countreg::zerotrunc(), a Poisson, negative binomial or geometric regression truncated at zero, fitted by damped Newton on the exact likelihood and its analytic Hessian. hurdle() is pscl::hurdle(), with a logit or a censored count zero part and the regressors after '|' for it. t/zerotrunc_hurdle.R.t pins both to countreg and pscl, with an mpmath third opinion at 60 digits (t/zerotrunc_hurdle.mpmath.py).
- svyglm() is survey::svyglm() on a one-stage design: sampling weights, strata, PSUs and a finite-population correction, with the Taylor linearisation variance and design degrees of freedom. t/svyglm.R.t is taken from survey's tests on the api data.
- ivreg() is ivreg::ivreg(): two-stage least squares by Householder QR, two- or three-part formulas, weights, HC0/HC1 and clustering, and the weak-instrument, Wu-Hausman and Sargan diagnostics. t/ivreg.R.t follows ivreg's tests and examples and statsmodels' Stata ivreg2 results.
- lmer() is lme4::lmer() by REML or ML, with random intercepts and slopes, correlated or not, crossed or nested grouping factors, and lmerTest's Satterthwaite degrees of freedom. The fit is held to a tightly converged lme4 fit rather than to lme4's default one, which stops about 1e-6 short in theta. t/lmer.R.t covers lme4's and lmerTest's examples and statsmodels' mixed-model corpora.
- anova() compares fitted models
- anova($m0, $m1, ...) with lm or glm fits is R's anova.lmlist / anova.glmlist, with test => 'F', 'Chisq' or 'LRT' and dispersion, and with negbin fits MASS's anova.negbin likelihood-ratio table; data and formulas still go to the XS anova() as before. The dispatching wrapper has no prototype, because the XS one's ($@) would put anova(@fits) in scalar context. Validated in t/anova_fits.R.t against R and MASS.
- read_table: five bugs
- The CSV parser split lines on $/, not on newlines, because sv_gets() reads PL_rs. Under a `local $/;` anywhere up the call stack the whole file was one "line" and read_table returned [] without a word; a record length ($/ = \N) cut rows at N bytes. The parser, and the perl-side peek that recovers a commented-out header, now split on "\n" whatever $/ is, and a filter still sees the caller's $/. t/read_table.input_record_separator.t.
- A UTF-8 byte-order mark, which Excel's "CSV UTF-8" export writes, was read as part of the first header name ("\xEF\xBB\xBFid"), and hid a comment or commented-out header behind it. It is dropped from the start of the file, as pandas drops it; t/read_table.bom.pandas.t takes its cases from pandas' test_utf8_bom and test_first_row_bom.
- In an .xlsx, a cell with formatting and no value (<c r="D1" s="1"/>) past a row's last value widened every row to reach it, so a shaded column came back as unnamed columns of undef and a duplicate-name warning. Such cells no longer count toward the width, which is what readxl 1.5.0 and pandas 2.2.3 give for the same workbook.
- A hoh read of a file whose only column is the row name came back as {}: the per-row hash was only made when a value was stored in it. Each row is now an empty hash, as R's read.table gives n rows of 0 columns. A read error part-way through a CSV is now an error rather than a truncated table.
- read_table: faster
- 'hoh' now goes through the XS fast path that aoh and hoa already took, with the perl closure's duplicate-row-name warning and undefined-row-name error spelled the same; a 300,000 x 5 CSV reads in 0.20 s rather than 0.91 s. Only a 'filter' still needs the closure.
- An unquoted field is copied straight from the line buffer instead of through the field accumulator, and the row hashes are filled through shared-hash-key copies of the column names, so no key is hashed again per cell. The same file reads as an aoh in 0.096 s rather than 0.120 s.
- write_table: wide characters in row.names and tex.longtable.head
- write_table() croaked "Wide character in subroutine entry" when row.names named a column outside Latin-1, or when tex.longtable.head was such a caption: the XS check for a non-numeric value read it with SvPVbyte. It now reads the bytes as stored, and no longer downgrades the caller's string in place or hands isdigit() a negative char. Found by XS::Check; t/write_table.t and t/write_table.longtable.t pin both.
- dunn_test, p_adjust: method names fold by ASCII rules
- The method name was lowercased with tolower(), which is undefined for a byte >= 0x80 in a signed char and follows LC_CTYPE, so under a Turkish single-byte locale 'BONFERRONI' folded its I outside ASCII and was rejected. Both now fold A-Z only. t/dunn_test.t and t/p_adjust.R.t.
Modules
Get basic statistical functions, like in R, but with Perl using XS for performance