Matrix model functions

OLS input and output

$result = Phit\fit_ols(
    x: [[1.0], [2.0], [3.0], [4.0], [5.0]],
    y: [2.1, 4.2, 5.8, 8.1, 10.2],
);
print_r($result);

Expected coefficients: [0.05, 2.01]; standard errors: [0.198074, 0.059722]; R-squared 0.997359; log-likelihood 2.516218. See examples/ols.php in the repository for formatted inference output, or the rendered comparison walkthrough for result tables.

Numeric inputs accept PHP integers or floats. null and NAN are missing values; rows containing them are dropped by default with a warning, or rejected with missing: ‘raise’. Strings, booleans and infinity are rejected. Arrays are consumed in iteration order; keys are ignored. Supply a rectangular n × k predictor matrix and n outcomes. The default prepends ones, giving p=k+1 coefficients; intercept: false gives p=k and uses x exactly. The retained sample must have n > p. An exact ones column enables centered R-squared and its adjusted version; other constant columns are not interpreted as intercepts. Without ones, R-squared is uncentered, matching R’s no-intercept models. Empty, ragged, infinite, numerically rank-deficient or overflowing inputs throw PHP exceptions. Rust returns Result<OlsResult, StatsError>.

Returned keys (all numeric arrays are zero-indexed):

Inference assumes independent homoskedastic Gaussian errors. Rank is checked by SVD after column scaling, using epsilon * max(n,p) * largest singular value. Very ill-conditioned full-rank designs may therefore be rejected.

Exact zero RSS yields zero standard errors and point confidence intervals; nonzero coefficients yield infinite t and zero p, while 0/0 gives NaN t and p. Intercept-only OLS reports R-squared and adjusted R-squared as zero, matching R. For other designs, zero total outcome variation gives NaN R-squared when RSS is also zero (otherwise negative infinity). Zero RSS gives positive infinite log-likelihood: the Gaussian likelihood is unbounded as variance approaches zero. PHP exposes NAN/INF, not null; these values cannot be encoded directly by ordinary JSON encoding. Mathematically perfect fits often retain tiny decomposition-dependent residuals; we do not round them to zero. Their t-statistics and likelihood can differ greatly from R despite coefficients and predictions agreeing to machine precision.

Binary logistic regression

$result = Phit\fit_logistic(x: $design, y: $binaryOutcomes, missing: 'drop', intercept: false);

Use the same plain rectangular arrays, automatic fixed intercept, ordering, numeric validation and missing-data policy as OLS. With default unit weights, outcomes must encode 0 or 1 under the documented success-count rounding tolerance. Integer trial weights also permit grouped proportions; successes and failures must both remain after deletion, and retained n must exceed p. Invalid numeric outcomes are rejected even if another value in their row is missing. The model uses independent binomial observations with linear log-odds offset + X beta. Unit trial weights give Bernoulli outcomes; positive integer trial weights allow grouped successes/trials supplied as proportions. General fractional likelihood weights and penalties are unsupported.

Bernoulli dispersion is fixed at one. Inference uses asymptotic normal z tests, not Student t tests or profile-likelihood confidence intervals. No sigma, R-squared or Gaussian variance estimate is returned. Other constant columns are not treated as intercepts for the null baseline. This baseline convention agrees with R glm’s no-intercept null deviance; statsmodels may choose a different null model, so our comparison constructs that baseline explicitly.

The dedicated fitter normally starts coefficients at zero. For an absolute logit offset >20, an explicit ones-column intercept starts from the marginal response link minus the mean offset. It then solves Newton/IRLS updates with column-scaled weighted SVD. Rank tolerance is epsilon × max(n,p) × largest singular value. Backtracking tries the full step and at most 30 halvings. Convergence requires each full scaled coefficient step to be at most 1e-8 × (1 + absolute next scaled coefficient), together with likelihood change at most 1e-10 × (1 + absolute likelihood). These fixed defaults are currently implementation choices, without a public tuning API.

Complete or quasi separation means the ordinary finite maximum-likelihood estimate may not exist. A separating direction encountered during fitting is rejected. This is a certificate check, not an exhaustive linear-programming separation detector: other difficult cases can fail the information-rank check, line search or iteration limit. Those failures throw exceptions rather than returning an unconverged model. Near separation can still produce very uncertain finite estimates; convergence is not assurance of reliable statistical inference. Overflowing inference also fails explicitly.

The PHP-free Rust entry point is fit_logistic(x: Vec<Vec<f64>>, y: Vec<f64>) -> Result<LogisticResult, StatsError>. It rejects nonfinite inputs; missing deletion belongs to the binding. Numeric field types match OlsResult, with usize for iterations and degrees of freedom. converged and sample deletion metadata are binding-only. See the executed derivation and verification record.

Probit regression

Phit\fit_probit(array $x, array $y, string $missing = 'drop', ?array $weights = null, ?array $offset = null, bool $intercept = true): array;

The binomial input contract and output keys match logistic regression: binary outcomes or grouped proportions with integer trial weights, both successes and failures after deletion, full rank, n > p after fixed-intercept insertion, probabilities in fitted, response residuals, coefficients/covariance/SE/z/p/Wald intervals, likelihood, deviance, null_deviance, AIC, df_residual, iterations, converged and sample metadata. Its link is Phi(offset + X beta), not sigmoid(offset + X beta). No latent variance parameter is estimated; it is fixed to one for identification.

The dedicated Fisher-scoring fitter uses expected information X’ diag(trials × phi(eta)^2/[Phi(eta)Phi(-eta)]) X, recalculated at the final fit. This convention matches R glm(probit) and statsmodels GLM(Binomial(Probit)). An observed-Hessian covariance from a discrete Probit implementation may differ. Confidence intervals and tests remain normal Wald, not profile likelihood. Coefficients are latent-normal effects; exponentiating them does not give a constant odds ratio. Null deviance uses an intercept-only fit retaining offsets with an exact ones column, otherwise the offset-only baseline (probability 0.5 with zero offsets).

The limit is 100 updates with up to 30 step halvings. Probit uses tighter full scaled-coordinate step tolerance 1e-11 × (1 + absolute next coordinate), and likelihood-change tolerance 1e-13 × (1 + absolute likelihood), because Fisher scoring can converge more slowly than canonical Newton updates. Rank tolerance and separation-certificate behavior match logistic regression. Failed convergence, rank failure or overflowing calculations throw rather than returning a partial fit. Tail likelihood uses erfc and a Mills-ratio continued fraction beyond -8.

The PHP-free entry point is fit_probit(x: Vec<Vec<f64>>, y: Vec<f64>) -> Result<ProbitResult, StatsError>. Fields have the same Rust types as LogisticResult. Missing filtering and converged/sample metadata remain in the PHP binding. See the walkthrough.

Poisson regression

Phit\fit_poisson(array $x, array $y, string $missing = 'drop', ?array $weights = null, ?array $offset = null, bool $intercept = true): array;

The raw design and missing-data contracts are the same. Outcomes must be finite, nonnegative integer counts (PHP integers or floats with an integer value). Fractions, negatives, strings and booleans are rejected, even in rows with a missing predictor. The retained sample must have at least one positive count, full numerical rank and n > p. All-zero samples are conservatively rejected, including no-intercept designs where a finite optimum might otherwise exist.

The model is independent Poisson observations with log mean offset + X beta. Returned keys and types match the logistic table, but fitted is a mean count, not a probability, and residuals are y - fitted. Confidence intervals/tests are normal Wald with dispersion fixed at one. There is no estimated sigma, RSS-based variance or robust covariance. Positive prior likelihood weights and additive offsets are supported. There is no separate exposure argument: supply offset = log(exposure) for positive known exposure. A coefficient’s exponential is a conditional mean multiplier per unit.

Likelihood includes -log Gamma(y+1). Deviance is twice saturated minus fitted log-likelihood, including the limiting zero-count convention; it is not -2 × log_likelihood. Null deviance uses mean count y-bar with a ones column, otherwise fixed mean one. AIC is 2p - 2 × log_likelihood, and df_residual is n-p. Zero individual counts are valid. Conditional variance equals the conditional mean; overdispersion can invalidate conventional standard errors.

The dedicated weighted-SVD Newton/IRLS fitter uses information weights equal to prior likelihood weights times the fitted means and final inverse information for covariance. It uses the logistic step and likelihood tolerances (1e-8 and 1e-10), at most 100 updates and 30 halvings. Certain predictor patterns have no finite MLE; numerical rank failure or iteration exhaustion throws. There is no exhaustive boundary-MLE existence check. Converged is always true in returned results; no partial fits are returned.

The PHP-free entry point is fit_poisson(x: Vec<Vec<f64>>, y: Vec<f64>) -> Result<PoissonResult, StatsError>. Its numeric Rust field types match LogisticResult. See the walkthrough for the likelihood, score, matrices and deviance.

Gaussian random-intercept mixed model

Phit\fit_mixed_intercept(array $x, array $y, array $groups, string $missing = 'drop', bool $intercept = true): array;

This first mixed model has continuous Gaussian outcomes and one random intercept per group. The fixed intercept is prepended by default; use intercept: false for a literal fixed design. This specialized call has no random slopes, crossed/nested grouping arguments, random-effects covariance matrices, offsets or observation weights. The estimation method is marginal maximum likelihood (ML), not REML.

Groups is a list of integer or UTF-8 string IDs aligned with x/y in iteration order. Labels are identities, not numbers used as predictors: integer 1 and string ‘1’ identify different groups. null is missing; floats (including NAN), booleans, arrays and objects are invalid IDs. Validation occurs before deletion. Listwise deletion includes all three inputs, warns once on a successful fit, and preserves original row accounting. Groups with no retained observations vanish. Labels and group-level outputs follow first retained appearance.

The retained design must be finite, full rank and have n > p. At least two groups and at least one repeated group are required; singleton groups are allowed when other groups repeat. One group or all-singleton data cannot identify the two variance components under this first contract. Numerically perfect fits, rank failures, unsupported variance-ratio ranges and overflowing calculations throw exceptions. No partial failed fit is returned.

Fixed-effect covariance is the plug-in GLS covariance (X’ V^-1 X)^-1 with estimated variance components. This matches lme4’s conventional covariance. It does not account fully for variance-component uncertainty. Normal Wald inference can be optimistic with few groups; there are no Satterthwaite or Kenward–Roger degrees of freedom, cluster-robust errors, variance-component p-values, or uncertainty intervals for predicted random effects. Our Python comparison recomputes this declared GLS convention from statsmodels estimates; statsmodels’ joint observed-Hessian standard errors may differ. df_residual n-p is a project convention; R merMod methods may count variance parameters too.

The dedicated fitter profiles sigma² and beta for rho=tau²/sigma², using group whitening and column-scaled SVD. A one-time QR square-root summary of within-group centered [X, y] rows, plus cached group means, reduces each trial to p+1+G rows. This is the random-intercept specialization of penalized least-squares elimination; no full observation-level SVD or covariance calculation is repeated in the search. The final covariance uses the compact SVD. It scans t=log(1+rho) from 0 to 32 in steps of 0.5 and refines every sampled negative-to-positive objective-score crossing with 40 bisections. It compares those candidate minima with the exact tau²=0 boundary. A decreasing objective at the upper limit is rejected. The finite grid is not a proof of a global optimum for every possible sample; very narrow unbracketed stationary points or extreme ratios are outside this initial guarantee. This is one dedicated profile search, without a generic optimizer interface.

A zero random-intercept variance is a valid ML boundary fit, indicated by boundary=true. Positive residual variance remains required. At that boundary fixed coefficients match OLS, but covariance uses the ML scale RSS/n, rather than OLS’s unbiased RSS/(n-p). Missing deletion is the only automatic warning policy; callers should also inspect boundary.

The PHP-free Rust entry point is fit_random_intercept(x: Vec<Vec<f64>>, y: Vec<f64>, groups: Vec<usize>) -> Result<MixedResult, StatsError>. Rust group IDs are arbitrary usize values and group_ids records first appearance. The binding preserves PHP labels separately and adds method/convergence/deletion metadata. See the walkthrough.

Weights, offsets and grouped binomial responses

The four fixed-model signatures accept optional weights/offsets and append the fixed-intercept option after existing positional arguments:

Phit\fit_ols(array $x, array $y, string $missing = 'drop', ?array $weights = null, ?array $offset = null, bool $intercept = true): array;
Phit\fit_logistic(array $x, array $y, string $missing = 'drop', ?array $weights = null, ?array $offset = null, bool $intercept = true): array;
Phit\fit_probit(array $x, array $y, string $missing = 'drop', ?array $weights = null, ?array $offset = null, bool $intercept = true): array;
Phit\fit_poisson(array $x, array $y, string $missing = 'drop', ?array $weights = null, ?array $offset = null, bool $intercept = true): array;

Optional vectors must have one value per original observation. null means unit weights or zero offsets. An explicitly supplied empty array is an invalid length, not the default. Weights are positive finite numbers; zero/negative weights are rejected. Offsets are known additive terms: eta = offset + X beta. For OLS they add directly to the outcome mean; Poisson log exposure is supplied as offset = log(exposure) with exposure positive. Numeric strings are not coerced.

OLS weights specify residual variance sigma²/weight. OLS RSS is weighted RSS; returned fitted values/residuals are on the original response scale. R-squared removes the known offset before calculating fitted variation, following R summary.lm. Gaussian likelihood includes sum(log(weights))/2.

Logistic/probit weights are positive integer trial counts, with y a proportion in [0,1] and weights[i] * y[i] an integer success count (absolute rounding tolerance 1e-8). Unit weights preserve binary-only behavior. For ten trials with three successes, supply y: [0.3, ...], weights: [10, ...]. Do not repeat a grouped response as ten fractional Bernoulli outcomes. The likelihood includes the binomial coefficient; deviance compares to the grouped saturated model. Fitted values and response residuals are proportions. Residual degrees of freedom count rows, not the sum of trials. Poisson weights are positive prior likelihood weights, not an exposure substitute.

Missing values in any supplied optional vector participate in listwise deletion with x and y. One warning and the same retained/dropped row metadata apply. Invalid values and original shapes are checked before deleting rows.

$exposure = [1, 2, 4, 8, 16, 3, 6, 12];
$counts = [1, 4, 3, 12, 23, 2, 11, 14];
$r = Phit\fit_poisson(
    x: array_fill(0, count($counts), []), // Intercept-only rate model.
    y: $counts,
    offset: array_map(fn($e) => log($e), $exposure),
);
$rate = exp($r['coefficients'][0]);
// fitted[i] = exposure[i] * rate; the result contains expected counts.

General Gaussian mixed-model call

Phit\fit_mixed(
    array $x, array $y, array $random,
    string $missing = 'drop', bool $reml = false,
    ?array $weights = null, ?array $offset = null, ?array $start = null,
    int $maxEvaluations = 20000, float $tolerance = 1e-7,
    string $optimizer = 'bobyqa', float $objectiveTolerance = 1e-12,
    bool $intercept = true,
): array;

Each random term is an array with z (n×k explicit random-design matrix), groups (n integer/string labels), and optional boolean correlated (default true). For (1 + slope | group), supply z columns [1, slope]; for independent intercept and slope use the same z with correlated: false. Multiple terms are independent of one another. Repeated labels across terms do not connect their random effects. For nested groups, supply explicit identifiers for each parent and parent/child combination. There is no two-level limit: supply three, four or more random terms within the documented dense-design limits. Groups and output ordering follow first appearance after deletion; integer 1 and string “1” remain distinct.

$r = Phit\fit_mixed(
    x: $fixedPredictors, // Fixed intercept is added automatically.
    y: $outcomes,
    random: [
        ['z' => $interceptAndSlope, 'groups' => $subjects, 'correlated' => true],
        ['z' => $interceptOnes, 'groups' => $sites],
    ],
    reml: true,
    optimizer: 'bobyqa', // or 'nelder_mead'
);

All original random-design/group/weight/offset rows participate in one deletion mask. Each term needs at least two groups and a repeated observation. For a correlated term, n must exceed its total random-coefficient count, matching lme4’s default check. Fixed columns must be numerically full rank, n>p, and residual scale must be positive. Singular random covariance is permitted and flagged with singular rather than an error. The engine rejects nonfinite results, exhausted budgets and invalid starts/controls.

BOBYQA is the pinned Rust PRIMA port; tolerance is its final parameter trust radius in (0, 0.1]. Nelder–Mead is the pinned argmin implementation; objectiveTolerance is the positive finite standard-deviation threshold of simplex objective values. The other control is unused by that optimizer. Both use a 0.1 initial step/radius, nonnegative Cholesky diagonals, off-diagonal bounds ±10,000 and diagonal upper bound 10,000. Starting values are in the supplied random-column units; omission uses diagonal factors one and off-diagonal factors zero. maxEvaluations limits optimizer cost calls, including rejected Nelder–Mead bound trials. evaluations also counts two pre/post likelihood evaluations. These are local optimizers; convergence does not establish global optimality or model identifiability.

Theta order is term order, then lower-triangular column-major order; independent terms include diagonals only. Each term’s covariance is sigma² L Lᵀ. Random effects are conditional means (BLUPs), ordered by retained group then z column. Each term’s effects_cond_var gives, per group, the k×k conditional covariance of that group’s effect vector given y at the estimated theta and sigma. It is the group’s diagonal block of sigma² Λ (ΛᵀZᵀWZΛ + I)⁻¹ Λᵀ, the matrix lme4 reports through ranef(condVar = TRUE) (the postVar attribute). The variance is conditional on the fixed effects, as in lme4. It uses the REML or ML sigma according to method. Blocks between different terms are not returned. A zero relative factor gives a zero block. lme4 splits (x || g) into one-column terms and reports only the diagonal; here the term is one block, so its off-diagonal conditional covariances are also returned. Dense limits are q≤512 random coefficients and n(q+p+1)≤50,000,000 entries; exceeding them fails explicitly. The existing fit_mixed_intercept remains the faster random-intercept ML path.

The result contains:

The covariance is plug-in GLS at estimated variance components. Wald intervals cover fixed effects only; random effects and variance estimates do not yet have intervals. REML likelihoods should not be used to compare different fixed-effect designs. Compatibility coverage and remaining work describes the exact boundary of the R/lme4 parity claim.

Binomial random-intercept ML

Phit\fit_logistic_mixed(
    array $x, array $y, array $groups, string $missing = 'drop',
    ?array $weights = null, ?array $offset = null,
    int $quadraturePoints = 1, int $maxEvaluations = 20000,
    float $tolerance = 1e-7, bool $intercept = true,
): array;

The conditional model is successes_i ~ Binomial(trials_i, logistic(offset_i + x_i beta + b_group)), with independent Gaussian group intercepts of variance tau². Unit weights mean one trial per row. For grouped responses, supply y as successes/trials and positive integer weights as trials. Success counts must be integer within absolute tolerance 1e-8; accepted floating-point count encodings are rounded to integer successes before fitting. This prevents near-zero input rounding from creating fictitious events. Residuals use the canonical count proportions. The first implementation limits trials per row to 1,000,000 to protect direct binomial likelihood arithmetic; larger counts need another stable log-mass implementation. R has broader input support.

The fixed intercept is prepended by default; intercept: false uses x literally. Require finite full-rank fixed design, n>p, at least two groups, a repeated group, and both successes/failures after count canonicalization and deletion. Some groups may contain only successes or only failures. Fixed-design separation prevents a finite unpenalized MLE and fails explicitly. Group labels use the existing integer/string contract and first-retained-appearance ordering. All original x/y/groups/weights/offset values and lengths are checked before one shared listwise deletion mask.

quadraturePoints: 1 uses Laplace, matching lme4’s default integration method. Values 2..61 use adaptive Gauss-Hermite quadrature for each group. Begin with 25 or 61 nodes for accuracy checks on small groups or rare events, and inspect sensitivity; no node count guarantees accuracy for every dataset. This function has one grouping factor and a random intercept only. Random slopes, multiple random terms, Poisson/probit mixed likelihoods, and REML are not implemented here. The general Gaussian mixed call continues to support its broader designs.

The optimizer is the pinned native Rust BOBYQA port. It starts from the ordinary binomial GLM coefficients and group SD one, and compares the solution with an exact zero-variance GLM candidate. Fixed columns are RMS-scaled for optimization. Scaled fixed coefficients are bounded ±100; group SD is bounded [0,20]. Reaching a bound fails. Tolerance is the final trust radius in (0,0.1]; maxEvaluations limits optimizer calls. The group mode solver uses bracketed Newton steps and positive posterior curvature. Nonconvergence, nonfinite results, nonstationary solutions or nonpositive observed information fail with an exception. A near-zero SD below 1e-4 is snapped to exact zero only when the likelihood gain is within floating-point resolution. These controls establish a local fit, not global optimality or scientific model adequacy.

At an exact zero variance boundary, coefficient inference and likelihood reduce to the ordinary binomial GLM; AIC/BIC still count the additional variance parameter, as lme4 does. Otherwise covariance comes from the beta block of inverse joint information in beta/tau, computed with centered finite differences and Richardson extrapolation. Normal Wald inference is asymptotic; no random-effect or variance-component confidence intervals are returned. neg2_log_likelihood is deliberately explicit: it is not lme4’s conditional deviance extractor.

For grouped data, lme4 nAGQ>1’s reported likelihood uses a different additive constant. Our verification retains that raw value and adds the saturated binomial log likelihood before comparing LL/AIC/BIC. See the executed verification and math walkthrough.