dcon

Stats::LikeR: R’s statistics toolkit, native in Perl

I got tired of writing CSV files just so I could run a t-test in R.

My data was nearly always in Perl already. I’d parse a log or a lab export, clean it up, dump it to disk, open R, run t.test, and paste the p-value back into the Perl script that built the report. Every one of those hops is a chance for a column to get reordered or an NA to turn into a zero without anyone noticing.

Before: Perl, CSV, R, copy by hand, Perl. After: one Perl script.

Stats::LikeR is what I wrote so I could stop doing that. Version 0.34 is on CPAN now.

In short:

• R’s functions, Perl’s data. mean, sd, quantile, t_test, lm, aov, lmer, coxph, chisq_test and over a hundred more, plus pandas-style data-frame tools: group_by, agg, merge, pivot_table, read_table, write_table. Everything takes and returns plain arrays and hashes.
• No second language. Unlike Statistics::R, there’s no background R session. Unlike PDL, there’s no object system or matrix class to convert into. The heavy lifting is in XS.
• Checked against R. The test suite compares results with R itself, or with pandas 2.2.3 and SciPy 1.15.2 where those are the reference.
• Fast at everyday sizes. Usually the quickest of Perl, R and Python up to a few hundred thousand rows, and the quickest at writing files at every size.

cpanm Stats::LikeR
perl -MStats::LikeR -e 'h()'      # list every documented function

One small analysis, three languages

A typical job: read a CSV, summarize weight and age by sex, compare the two groups with a t-test, then fit a linear model. The file is made up; the calls are real.

Perl

use Stats::LikeR;

my $df = read_table('patients.csv');    # array of hashes, one per row

my $by_sex = agg($df,
    by  => 'sex',
    agg => { wt => [ 'mean', 'sd' ], age => [ 'mean', 'count' ] },
);

my @f = map { $_->{wt} } grep { $_->{sex} eq 'f' } @$df;
my @m = map { $_->{wt} } grep { $_->{sex} eq 'm' } @$df;
my $t = t_test(\@f, \@m);
printf "t = %.3f, p = %.4g\n", $t->{statistic}, $t->{p_value};

my $fit = lm(formula => 'wt ~ age * sex', data => $df);

R. Base aggregate takes one function, so several summaries mean writing an anonymous function that returns a vector:

df <- read.csv("patients.csv")
aggregate(cbind(wt, age) ~ sex, data = df,
          FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)))
t.test(df$wt[df$sex == "f"], df$wt[df$sex == "m"])
fit <- lm(wt ~ age * sex, data = df)

Python needs three libraries (pandas, SciPy and statsmodels) for the same four steps:

import pandas as pd
from scipy import stats
import statsmodels.formula.api as smf

df = pd.read_csv("patients.csv")
df.groupby("sex").agg({"wt": ["mean", "std"], "age": ["mean", "count"]})
stats.ttest_ind(df.wt[df.sex == "f"], df.wt[df.sex == "m"], equal_var=False)
fit = smf.ols("wt ~ age * sex", data=df).fit()

If you know R, the Perl will look familiar:

• Formulas use R’s syntax: * adds the interaction, - drops a term, . means every other column.
• t_test runs Welch’s test unless you pass var_equal => 1, the same default as R (Python makes you ask for it).
• The result hash has the fields R prints: statistic, df, p_value, conf_int, estimate.

And what comes back is ordinary Perl data. $by_sex is an array of hashes because that’s what went in; pass output_type => 'hoh' for a hash of hashes. Loop over it, dump it to JSON, or hand it to the rest of your program.

Less typing for everyday jobs

Most of what I wanted was for the obvious thing to happen without me remembering an argument.

Reading files: the file name says what to do

my $csv  = read_table('samples.csv');
my $tsv  = read_table('samples.tsv');                    # tab-separated; no sep needed
my $gz   = read_table('cohort.tsv.gz');                  # decompressed for you
my $xlsx = read_table('results.xlsx', sheet => 'Sheet2');

One function for all four. pandas needs sep="\t" for the TSV and a separate pd.read_excel (plus openpyxl) for the spreadsheet. Base R needs read.csv, read.delim and the readxl package. (To be fair, data.table::fread guesses separators.) You can still pass sep when a file is named something odd.

It also reads a commented-out header, which many bioinformatics tools write:

# PDB   score
1a2b    10
3c4d    20

read_table uses PDB and score as the column names. pandas and R skip that line as a comment, and you name the columns by hand.

Writing files works the same way

write_table($df, 'out.csv');
write_table($df, 'out.xlsx');    # a real Excel workbook
write_table($df, 'out.tex');     # a LaTeX tabular

pandas spreads this over to_csv, to_excel (with openpyxl or XlsxWriter installed) and to_latex. R uses write.csv, writexl and xtable.

Every column against every other column

my %data = (height => \@height, weight => \@weight, age => \@age);   # hash of arrays
my $r = col2col(\%data, 'cor');
print $r->{height}{weight};

col2col runs any two-column function over every pair of columns. A correlation matrix is easy anywhere, but swap 'cor' for 't_test' and you get a t-test for every pair. In R or pandas that’s a loop you write yourself.

Table 1 without listing your factors

my $t1 = table_one(\@cohort, by => 'arm');
write_table($t1, 'table1.xlsx');

table_one works out which columns are continuous (all numbers) and which are categorical, and picks a t-test, ANOVA or chi-squared test to match. R’s tableone wants the categorical columns listed in factorVars.

ANOVA straight from groups

my $fit = aov({
    ctrl => [4.17, 5.58, 5.18, 6.11, 4.50],
    trt1 => [4.81, 4.17, 4.41, 3.59, 5.87],
    trt2 => [6.31, 5.12, 5.54, 5.50, 5.37],
});

With no formula, aov stacks the groups into long form itself. In R you’d write aov(values ~ ind, data = stack(groups)).

Help without leaving the terminal

perl -MStats::LikeR -e 'h(*agg)'

That prints the documentation for agg. Passing the glob *agg saves you escaping quotes in a one-liner; inside a script, h('agg') works too.

Resource use

I timed each function against R 4.6.1 and Python 3.14.2, from 1,000 elements up to 300,000 rows for data frames and files, and up to 1,000,000 elements for vectors, seven runs each.

The short version: at the sizes most scripts handle, Stats::LikeR is usually the fastest of the three. At the top end, NumPy and pandas catch up on plain arithmetic and pull well ahead on joins and filtering. Writing files is where it wins at every size.

Small vectors: overhead is everything

sum over 1,000 numbers: Stats::LikeR 1.2 µs, R 2.9 µs, List::Util 4.5 µs, Python 24 µs

On a 1,000-element vector, the cost of the call matters more than the arithmetic, and an XS function given an array ref has very little overhead. That’s why sum, min and max are 3.6 to 5.7 times faster than List::Util 1.70 at this size.

By a million elements the picture changes: R’s sum is faster (0.8 ms against 1.4 ms) and NumPy’s mean is about seven times faster than mine. skew and kurtosis are the exception and stay ahead of SciPy all the way up.

Pass a reference, save the memory

Peak memory: writing CSV under 10 kB vs 6 MB for pandas CSV and 562 MB for pandas XLSX; sum over an array ref 4 kB vs 8 MB for a flattened list

The reductions also accept a flat list, List::Util-style: sum(@x). That makes Perl copy the whole array onto the stack. Pass \@x instead and sum, sd, var, skew and friends stay at about 4 kB however long the vector is. A flattened million-element list costs 8 MB, whether it’s List::Util, Python or my own list form.

write_table streams, so it used under 10 kB for every CSV I tried. pandas reaches 6 MB for CSV and 562 MB for XLSX.

Reading and writing files

300,000 rows. Write CSV: 55 ms vs 346 ms (R) and 483 ms (pandas). Write XLSX: 0.41 s vs 17.3 s. Read mixed CSV: 174 ms vs 376 ms and 409 ms.

This is the clearest win. A 300,000-row CSV takes 55 ms to write, against 346 ms in R and 483 ms in pandas. The same table as XLSX takes 0.41 seconds, against 17.3 seconds in Python.

Transforms

rank and auc beat both rivals at every size. auc is 30 times faster than R at 1,000 elements and still 3 times faster than Python at a million. uniq and seq lose everywhere; R’s seq is lazy and doesn’t build anything, so that one was never going to go my way.

Data frames

At 1,000 rows, Stats::LikeR wins most comparisons: group_by plus mean is 5 times faster than pandas, and select_cols takes about 2 microseconds whatever the frame size. aoh2hoa, which turns a list of rows into a hash of columns, is 20 times faster than the pandas equivalent even at 300,000 rows.

Against other Perl modules

The set functions (union, intersection, Lonly, Ronly, is_equivalent) are 1.5 to 3.3 times faster than List::Compare 0.55 at a million elements, and use 2 to 40% less memory.

Where to use something else

merge of 300,000 rows: 214 ms vs 2 ms in pandas. All-numeric CSV read: 148 ms vs 61 ms. sum over a million: 1.4 ms vs 0.8 ms in R.

I’d rather you hear this from me. At 300,000 rows pandas does a merge in 2 ms and I take 214 ms; filter, group_by and transpose fall behind too. pandas also reads all-numeric CSV faster. If you’re joining big tables, use pandas.

What’s in it

The README documents over 115 functions. Roughly:

Area	Functions	What you’d use otherwise
Descriptive statistics	mean, median, mode, sd, var, quantile, skew, kurtosis, summary, table_one, density, hist	R base/stats; NumPy, SciPy
Hypothesis tests	t_test, wilcox_test, chisq_test, fisher_test, ks_test, shapiro_test, kruskal_test, friedman_test, mcnemar_test, cmh_test, binom_test, prop_test, var_test, dunn_test, p_adjust, power_t_test	R stats; scipy.stats
Models	lm, glm, aov, anova, lmer, ivreg, svyglm, coxph, hurdle, zerotrunc, prcomp, predict, vif	lm, lme4, survival, survey; statsmodels
Effect sizes & diagnostics	cohen_d, cramers_v, eta_squared, smd, roc, auc, bedroc, hosmer_lemeshow, survfit, logrank_test	effectsize, pROC; scikit-learn
Data frames	read_table, write_table, filter, select_cols, group_by, agg, merge, pivot_table, melt, concat, drop_duplicates, dropna, fillna, interpolate, qcut, value_counts, view	dplyr, tidyr; pandas
Distributions & sampling	dnorm, pnorm and the rest of the d/p/q/r family, rnorm, runif, rbinom, sample, seq	R stats; numpy.random, scipy.stats

Correct first, then fast

A stats library that’s fast but wrong is worse than useless, so the tests compare against R’s output, not against the module’s own earlier results. Some specifics:

• lmer gives the Satterthwaite degrees of freedom that lmerTest adds.
• coxph handles ties with Efron’s method, like survival::coxph.
• interpolate matches pandas and SciPy to machine precision.
• min, max and sum pass List::Util’s own test suite, apart from a few deliberate R-like choices: a NaN anywhere in the list gives NaN, and an empty list dies instead of returning undef.

Where I chose not to copy R, the docs say so and explain why. For example, lm sorts factor levels by byte value, as patsy and pandas do, rather than by locale. The fit is identical, but the reference level can differ.

Should you use it?

If your data is already in Perl and fits in memory, I think so. That covers a lot: log analysis, lab pipelines, survey files, nightly reports, most things under a few hundred thousand rows. You drop the export step and the second language, and at those sizes it’s usually quicker too.

It runs on Perl 5.10.1 and later and needs a C compiler to build:

cpanm Stats::LikeR

use Stats::LikeR; exports the whole toolkit, which is handy for quick scripts. In a larger program, import just what you need:

use Stats::LikeR qw(read_table t_test lm);

Bug reports are welcome on GitHub, and so are benchmarks where it does badly. Those are the most useful ones.

Links

• Stats::LikeR on MetaCPAN
• Source and issue tracker on GitHub

Also posted on my WordPress blog.