R compatibility: target, coverage and remaining work

Distribution and statistical utility compatibility

The v0.2 probability APIs cover the planned distribution families with scalar/batch density, mass, tail and quantile operations. Descriptive statistics follow sample variance/covariance and R type-7 quantiles, with explicit frequency/reliability/population weight definitions. R and independent 60-digit mpmath fixtures verify ordinary and extreme values; the comparison record documents intentional point-mass/endpoint differences, R inaccuracies and SciPy peer discrepancies. This numerical coverage does not implement R recycling, coercion, vector names or formulas. Seeded sampling uses local ChaCha20 state; its sequences are independent of R’s RNG algorithms/state. Existing model paths and their behavior remain unchanged; the estimator reuse experiment is recorded in the v0.2 assessment.

Implemented and verified

The new mixed call is fit_mixed(). The existing fit_mixed_intercept() remains the optimized one-group random-intercept ML call used by the million-row benchmarks. The general implementation caches a dense QR square-root summary and then solves penalized least squares; it does not allocate an observation covariance matrix. It has explicit dense-design limits: at most 512 random coefficients and 50,000,000 entries in the initial augmented matrix. This is not yet lme4’s sparse implementation and has not been benchmarked on large general designs.

An intercept-only OLS fit now reports R-squared and adjusted R-squared as zero, matching summary.lm(). Other undefined statistics retain their documented NAN/INF representation. This is a documented change to the previous intercept-only zero-response R-squared convention.

Optimizers: runtime versus reference

BOBYQA implementations can share an algorithm without producing identical optimization trajectories. Our Nelder–Mead bounds and stopping rules also differ from lme4’s implementation. We compare final estimates and likelihoods and name the runtime implementation in each result. Unsupported R wrapper names are rejected; they do not silently select another optimizer. The existing fast random-intercept call keeps its specialized scalar profile search.

The lme4 control documentation and wrapper documentation describe their controls. Supporting all those controls and backends is still open.

Evidence and tolerance policy

The executed comparison notebook fits 12 fixed-model cases and 15 mixed-model cases, including the published sleepstudy random-slope model under both ML and REML, plus three- and four-grouping-level nested designs. It checks coefficients, covariance, inference, every observation’s fitted/residual values, random effects/covariances, scale, likelihood and information criteria. Both native optimizers are compared against all four R wrappers. sleepstudy ML/REML also have independent statsmodels fits. Older fixtures and the four published random-intercept datasets remain tested.

Fixed-model comparisons use abs_error <= 1e-7 + 1e-6 * abs(reference). Native BOBYQA versus R BOBYQA uses 1e-5 + 1e-5 * abs(reference) for mixed estimates; cross-optimizer comparisons use 1e-3 + 3e-4 * abs(reference). All mixed likelihood/AIC/BIC comparisons use the stricter 5e-6 + 1e-9 * abs(reference). Shapes, sample metadata, grouping order, parameter counts and singular flags are exact where compared. The notebook reports actual errors, not only pass/fail.

GLM covariances are the inverse Fisher information at the final coefficients. R’s vcov() and summary() for a glm use the working weights of the last IRLS iteration, which were computed from the coefficients one step before the final update. Newton steps converge quadratically, so R’s coefficients are already at the maximum to many digits, but its weights lag one step behind. Its covariance therefore differs from the information at its own coefficients by about the size of that last step. This is not a difference in the convergence tolerance, in Fisher versus observed information (they are equal for the canonical logit and log links), or in dispersion (fixed at one for the binomial and Poisson families).

For am ~ wt + hp on mtcars, the coefficients agree to 12 digits. R’s last step changed them by 3.4e-8 relative. The variance of the wt coefficient is 9.41676694948 from vcov() and 9.41676794579 from phit, 1.1e-7 relative. Recomputing the information in R at coef(fit), or refitting with start = coef(fit), gives 9.41676794579. On designs near separation the last step is larger and the two can differ by about 1e-4 relative. phit’s R fixtures restart glm() from its own solution until the coefficients stop changing, so they compare against the information at the converged coefficients.

Different optimizers can stop at slightly different variance estimates on a flat likelihood. The cross-optimizer tolerance accommodates that observed spread, particularly covariance entries and small residuals on sleepstudy. It does not replace the stricter BOBYQA check or likelihood check.

The R fixtures record each wrapper’s controls and warnings. They use tightened optimization controls and boundary.tol=1e-3, with lme4 checking whether snapping near-boundary Cholesky parameters to zero improves the objective. At the default boundary tolerance, one singular REML case was left just above zero by nloptwrap and classified nonsingular. Tightening nlminb excessively also produced singular convergence warnings. These fixtures validate explicit controls, not identity under every package default. The saved reference runs have convergence code zero and no warnings. The multi-level BOBYQA references also use a local R-only restart (rhobeg=0.001, rhoend=1e-9) to improve variance-component precision on flat profiles; no native estimates are used to seed R. Published historical benchmark reports retain their original model/optimizer definitions; their timings do not describe the new general fitter.

Remaining compatibility work

The following items are not implemented and are not covered by a blanket parity claim:

• Native NLopt and PORT/nlminb backends, full optimizer controls, restart policies, gradient/Hessian convergence diagnostics, and an allFit-style runtime audit.
• Sparse general mixed-model computation and large-design verification.
• R-style handling of zero weights, aliased/rank-deficient fixed columns, zero-column designs, and fit objects retained after nonconvergence. We currently reject these or fail with an explicit error rather than returning unreliable inference.
• Configurable confidence levels; profile and bootstrap intervals; likelihood-ratio comparisons and fixed-effect degrees-of-freedom corrections. Mixed Wald p-values are an explicit extension; lme4’s default Gaussian summary does not supply them.
• Simulation on new observations, prediction standard errors and intervals for mixed models, residual diagnostics, leverage/influence, and model-object methods such as update/refit/anova. (predict() covers point predictions, standard errors and lm intervals for new data, as described with the formulas; random-effect conditional variances are reported by the mixed fits.)
• Formula parsing, contrasts and factor encoding, subset/na.exclude semantics, multivariate lm responses, fractional binomial weights and additional GLM families/links.
• Broader generalized/nonlinear mixed models: binomial random slopes/multiple terms, mixed Poisson/probit/negative-binomial likelihoods, nlmer, and related methods. The new binomial random-intercept glmer counterpart is implemented separately, with Laplace/adaptive quadrature; it does not yet support the general Gaussian random-term designs.

Each addition needs an input/output contract, a mathematical walkthrough, saved R reference cases and PHP integration checks. This document is the compatibility ledger: move an item to implemented only when that evidence exists.

Reproduce

cargo test --workspace --locked
cargo build -p phit --release --locked
php -n -d "extension=$PWD/target/release/libphit.dylib" tests/r-compatibility.php
php -n -d "extension=$PWD/target/release/libphit.dylib" tests/lmm-expanded.php
PHP=/opt/homebrew/bin/php .context/notebooks-venv/bin/python docs/verifications/compare_r_compatibility.py
# Refresh expected R values explicitly (R and lme4 are development dependencies):
Rscript tests/generate-r-compatibility-fixtures.R
Rscript tests/generate-lmm-expanded-fixtures.R

Use .so on Linux and select matching PHP/php-config as described in development. Normal Rust/PHP fixture tests need no installed R. The live Python comparison requires the pinned notebook environment.

Binomial mixed-model evidence

The new executed verification adds 28 R/PHP binomial random-intercept comparisons, independent normalized-likelihood checks, QUADPACK approximation-error measurements, an independently optimized cbpp marginal likelihood and statsmodels GLM boundary reduction. R GLMM references currently use BOBYQA; this does not extend the four-optimizer Gaussian audit into a claim about all GLMM optimizer controls.

The v0.1 pass stabilizes this tested subset; full R compatibility is not a v0.1 objective or promise. This ledger describes coverage, not a release guarantee.