Model formulas

Example

$data = [
    'rent' => [1879.25, 1805.03, 2022.29, 1267.85, 1509.24, 2240.10, 959.51, 2386.54],
    'area' => [74.4, 57.4, 66.4, 49.1, 45.7, 92.4, 37.0, 93.9],
    'district' => ['north', 'south', 'north', 'east', 'south', 'north', 'east', 'south'],
];

$fit = Phit\fit_lm('log(rent) ~ log(area) + district', $data);
$fit['columns'];  // ['(Intercept)', 'log(area)', 'districtnorth', 'districtsouth']
$fit['named'];    // ['(Intercept)' => ..., 'log(area)' => ..., ...]
$fit['p_values']; // Aligned with columns, as in Phit\fit_ols.

$model = Phit\fit_glm('passed ~ study_hours * practice_score', $learning, family: 'binomial', link: 'probit');
$mixed = Phit\fit_lmer('score ~ study_hours + prior_score + (1 | school)', $schools);

Data

$data is an associative array of columns, name => list of values. An object is also accepted. The extension calls its public toModelData() method, or toColumns() when the first method is absent, and uses the returned array. Exceptions thrown by those methods propagate unchanged. phibble’s DataFrame::toModelData() returns this shape.

Each column is one of the following.

• Numbers: integers and floats. null and NAN are missing. Infinity raises 1003.
• Booleans: a logical column. R codes it as a two-level factor named {col}TRUE.
• Strings: a character column. It becomes a factor whose levels are the distinct retained values in byte (C-locale) order, which is R’s order under LC_COLLATE=C. Other R locales can sort mixed-case or accented strings differently.
• A factor descriptor: an array with exactly the keys levels (list of unique strings), codes (list of 0-based integers, null for missing), and ordered (bool). Levels keep the caller’s order and the first level is the reference.
• Dates arrive as numbers, for example days or seconds since the epoch.

Every column must have the same length (code 1002). Only the columns a formula uses are type-checked. An unknown column, a malformed factor descriptor, or a referenced list-column or mixed-type column raises 1001. A list of row arrays raises 1001.

Formula syntax

The formula grammar and operator precedence follow R.

Terms are ordered as R orders them: main effects first, then two-way interactions, and so on, with ties in formula order. Within an interaction, variables appear in order of first appearance in the formula, and the first variable varies fastest in the column order. y ~ b:a + a therefore gives bv:ap, bv:aq, and so on.

Missing values and unused levels

Listwise deletion uses only the variables that the formula references. These include the response, transforms, offsets, weights, grouping factors, and random-effect columns. Other columns may contain missing values. As in the matrix functions, missing: 'drop' removes incomplete rows with one E_WARNING, and missing: 'raise' raises 1001. retained_rows and dropped_rows are original zero-based positions.

model_matrix keeps unused factor levels, as R’s model.matrix() does, and their indicator columns are all zero. The fitting functions drop levels that are absent after deletion, as lm(), glm(), lmer(), and glmer() do through drop.unused.levels = TRUE. Character levels always come from the retained rows.

model_matrix

The result has these keys.

x includes the intercept column. Pass it to a matrix function with intercept: false. Each random entry has the shape that fit_mixed expects:

$mm = Phit\model_matrix('y ~ x + (1 + x | g)', $data);
$fit = Phit\fit_mixed(x: $mm['x'], y: $mm['y'], random: $mm['random'], intercept: false);

fit_lm and fit_glm

fit_lm calls the fit_ols core and returns its full result. fit_glm supports four family and link combinations. An empty link selects the canonical link.

Other families or links raise 1007. A binomial response can be 0/1 numbers, booleans, a factor or character column (the first level is failure and all other levels are success, as in R), proportions with integer trial weights, or cbind(successes, failures).

$weights is either an array aligned with the original rows or the name of a numeric column. Missing weights delete their rows. Weights must be positive.

Formula fits add these keys to the matrix result: formula, response, columns (aligned with coefficients), terms, assign, intercept, and named (column => coefficient), and predvars. predvars holds the state that prediction must reuse: each factor’s levels after deletion and each orthogonal poly() basis’s recurrence coefficients (R’s xlevels and predvars). fit_glm also returns family and link. All list keys of the matrix functions are unchanged.

Random-effect terms in fit_lm or fit_glm raise 1007 and name fit_lmer and fit_mixed as the alternatives.

fit_lmer and fit_glmer

fit_lmer maps lme4-style bars onto fit_mixed. REML is the default, as in lme4::lmer(). Supported bars include (1 | g), (x | g), (0 + x | g), (x || g) (diagonal covariance), several terms such as (1 | g1) + (1 | g2), interaction groups (1 | g1:g2), and nesting (1 | g1/g2). Nesting expands to (1 | g1) + (1 | g2:g1), with lme4’s group name and g2level:g1level labels. The result is the fit_mixed contract. Each random entry also has label, group, columns, and correlated. Its effects_cond_var lists one k×k conditional covariance per group, in group_labels order. These are the blocks that ranef(fit, condVar = TRUE) attaches as postVar. For (x || g), lme4 reports only the diagonal of each block.

fit_glmer supports exactly one random intercept (1 | g) with the binomial family and logit link, over the fit_logistic_mixed core. quadraturePoints corresponds to lme4’s nAGQ. Responses follow fit_glm. random_effects_cond_var gives the conditional variance of each group’s mode, as lme4’s postVar.

Prediction: predict

predict($fit, $newdata) predicts from any formula fit, as R’s predict() does for lm, glm and merMod objects. It rebuilds the fit’s formula over $newdata, which needs the predictor columns but not the response.

$fit = Phit\fit_glm('am ~ wt + poly(hp, 2)', $cars, family: 'binomial');
$p = Phit\predict($fit, ['wt' => [2.2, 3.4], 'hp' => [95, 180]], type: 'response', seFit: true);
// $p['fit'] probabilities, $p['se_fit'] delta-method standard errors

• Factors are coded with the levels the fit saw. A level the fit never saw raises 1001 with R’s message (factor g has new level z). Logical and character columns work the same way.
• poly(x, k) evaluates the fitted orthogonal basis at the new values, as R’s predict.poly() does. It is not refitted to the new data.
• offset() terms are evaluated over $newdata and added.
• A row with a missing predictor predicts null and is listed in missing_rows. Every list has one entry per $newdata row.
• type: 'link' (the default) returns the linear predictor. type: 'response' applies the inverse link, and seFit then uses the delta method, as R does.
• interval: 'confidence' or 'prediction' adds lower and upper for fit_lm fits, using t with df_residual at level. A prediction interval adds the residual variance; as in R, it assumes constant variance for a weighted fit.
• For fit_lmer and fit_glmer fits, predictions add each known group’s random effects (lme4’s default re.form = NULL). random: false gives population predictions (re.form = NA), and the grouping columns can then be absent. An unseen group raises 1001 unless allowNewLevels: true, which predicts it at the population level. Mixed predictions have no standard errors, as in lme4.

GLM standard errors use the covariance at the final coefficients. R’s predict.glm() uses vcov(), which reuses the weights of the iteration before the final one. The two differ by about the size of R’s last IRLS step: 1e-7 relative for am ~ wt + hp on mtcars, and up to about 1e-4 near separation. See R compatibility.

Differences from R

• Character levels use byte order. R’s default locale collation can order mixed-case or non-ASCII strings differently.
• Rank-deficient designs raise 1004. R reports NA (aliased) coefficients instead.
• Weights must be positive. R accepts zero weights. cbind() rows with zero trials are rejected; R gives them zero weight.
• A response that also appears on the right-hand side raises 1001. R drops the term with a warning.
• Only the functions listed above can be evaluated. factor() takes one argument, and poly() takes one variable.
• x %in% a is labelled x:a. R labels the term x %in% a; the columns are identical.
• model_matrix returns factor responses as 0-based codes with response_levels.
• fit_glm(family: 'gaussian') uses t-based intervals, as fit_ols does. R’s confint.default() is z-based and confint() profiles the likelihood.
• fit_lmer orders groups and effects by first appearance after deletion and keeps random terms in formula order. lme4 sorts levels and orders terms by their number of levels. Estimates are invariant to these orderings.
• Duplicate column names collapse in the named map; columns keeps every name.

Verification

tests/generate-formula-fixtures.R saves 62 R references under LC_COLLATE=C: 33 model.frame()/model.matrix() designs, 14 lm() fits, and 15 glm() fits. They use the repository’s housing, schools, learning, clinic, and cbpp data and a seeded synthetic data frame with factors, unused and reordered levels, ordered factors, logicals, character columns, missing values, weights, offsets, transforms, and no-intercept models. glm() references are restarted from R’s own solution until the coefficients stop changing, because the default deviance criterion fixes coefficients only to about the square root of its tolerance. tests/formula.php compares design columns, terms, and rows exactly, design values within 1e-12 relative, and fit quantities within the matrix suites’ 1e-9 + 1e-7*abs(reference).

tests/generate-formula-lme4-fixtures.R saves 10 lme4 references (7 lmer() and 3 glmer() fits), including crossed, nested, correlated, and uncorrelated terms and adaptive quadrature. The references include the ranef(condVar = TRUE) conditional variances, which are also compared for fit_mixed_intercept on the ML schools case. tests/formula-mixed.php uses the tolerances of the existing lme4 suites.

tests/generate-predict-fixtures.R saves 25 predict() references: lm() with orthogonal poly(), factors, logicals, ordered interactions, offsets, weights, no intercept, standard errors, and confidence and prediction intervals; glm() on the link and response scales for every family and link (including cbind() responses and offsets); and lmer()/glmer() conditional, population and response-scale predictions with an unseen group. The new data include values outside the fitted range, missing rows and a factor whose level set differs from the fit. GLM standard errors are compared tightly with the final-coefficient value recomputed in R and loosely (1e-3 relative) with R’s own se.fit. tests/predict.php also checks that predicting the training data reproduces the fitted values, and the error messages. All fixtures are indexed in the reference manifest.