Statistics

Requirements and fit keys

The helpers need no extension to plot a fit array. Two operations build design rows from the fit’s formula with Phit\model_matrix(), which comes from the phit extension or its polyfill:

• predictions for formula fits in effect_plot(), predictions(), and marginal_means(), unless design: is given
• random-slope design columns in random_effects_plot(), when it rebuilds conditional variances from the data

Without either, these calls throw an InvalidArgumentException. Each helper reads these keys of the fit array:

Coefficient plots

coefplot() draws estimates with confidence intervals and a reference line at no effect. It drops the intercept unless intercept: true. It uses the fit’s conf_int by default, and level: recomputes the intervals from the standard errors with t or normal quantiles to match the fit.

use function Phiz\coefplot;

// The fit as phit returns it: Phit\fit_lm('mpg ~ wt + hp + am', $mtcars).
$fit = json_decode(file_get_contents(__DIR__.'/fits/mtcars-lm.json'), true);

$plot = coefplot($fit)->labs(title: 'mpg ~ wt + hp + am');

Several fits keyed by label are compared side by side, with one color per model. sort: true orders terms by estimate, and terms: selects and orders them.

use function Phiz\coefplot;

$full = json_decode(file_get_contents(__DIR__.'/fits/mtcars-lm.json'), true);
$reduced = json_decode(file_get_contents(__DIR__.'/fits/mtcars-lm-reduced.json'), true);

$plot = coefplot(['mpg ~ wt + hp + am' => $full, 'mpg ~ wt + hp' => $reduced]);

exponentiate: true shows odds or rate ratios from a GLM on a log scale with the reference at 1. The axis is labeled by the fit’s link: log odds or odds ratios for logit models, and log rates or rate ratios for log links.

use function Phiz\coefplot;

// Phit\fit_glm('am ~ wt + hp', $mtcars, 'binomial')
$glm = json_decode(file_get_contents(__DIR__.'/fits/mtcars-glm.json'), true);

$plot = coefplot($glm, exponentiate: true)->labs(x: 'Odds ratio (95% CI)');

Regression diagnostics

residual_plot(), qq_plot(), scale_location_plot(), and fitted_plot() draw the standard diagnostic panels. diagnostic_plots() arranges all four in a grid.

use function Phiz\diagnostic_plots;

// The fit as phit returns it: Phit\fit_lm('mpg ~ wt + hp + am', $mtcars).
$fit = json_decode(file_get_contents(__DIR__.'/fits/mtcars-lm.json'), true);

$plot = diagnostic_plots($fit, width: 240, height: 180);

qq_plot() also accepts a plain list of values. The QQ geoms place the sorted sample at the normal quantiles of R’s ppoints(n), per group, as qqnorm() and ggplot2’s stat_qq() do. The reference line passes through the sample quartiles, as qqline() and stat_qq_line() draw it.

Predicted values

effect_plot() draws predictions with confidence bands as one variable varies across its observed range. by: draws one line per level of a grouping variable. Other variables are held at their mean, or at their most common value when they are categorical.

For formula fits, phiz builds the prediction design with Phit\model_matrix() on the fitting data and the grid together, so factor levels and reference categories match the fit. For matrix fits, pass design:, a function from a grid row to a design row. GLM intervals are computed on the link scale and transformed to the response scale. The predictions match R’s predict() with interval = "confidence".

use Phiz\Datasets;

use function Phiz\effect_plot;

// The fit as phit returns it: Phit\fit_lm('mpg ~ wt + hp + am', $mtcars).
$fit = json_decode(file_get_contents(__DIR__.'/fits/mtcars-lm.json'), true);

// Formula fits build the prediction design with Phit\model_matrix(); here we pass it directly.
$plot = effect_plot($fit, Datasets::mtcars(), 'wt', by: 'am', design: fn (array $row) => [1, $row['wt'], $row['hp'], $row['am']])
    ->labs(x: 'Weight (1000 lbs)', y: 'mpg', color: 'Manual', fill: 'Manual');

use Phiz\Datasets;

use function Phiz\effect_plot;

// Phit\fit_glm('am ~ wt + hp', $mtcars, 'binomial'): predictions on the probability scale.
$glm = json_decode(file_get_contents(__DIR__.'/fits/mtcars-glm.json'), true);

$plot = effect_plot($glm, Datasets::mtcars(), 'wt', design: fn (array $row) => [1, $row['wt'], $row['hp']])
    ->labs(x: 'Weight (1000 lbs)', y: 'P(manual transmission)');

scale: 'link' draws predictions on the link scale instead, where a GLM is linear in its predictors:

use Phiz\Datasets;

use function Phiz\effect_plot;

// The same logistic fit, on the log-odds scale, where the model is linear.
$glm = json_decode(file_get_contents(__DIR__.'/fits/mtcars-glm.json'), true);

$plot = effect_plot($glm, Datasets::mtcars(), 'wt', design: fn (array $row) => [1, $row['wt'], $row['hp']], scale: 'link')
    ->labs(x: 'Weight (1000 lbs)');

predictions($fit, $grid) returns the predicted rows without plotting them, and takes the same scale: option.

Marginal means

marginal_means() estimates the mean response at each level of a factor, as R’s emmeans package does. For each level it averages the predictions with equal weight over the levels of the model’s other categorical variables, with numeric variables at their means. The interval comes from the averaged design row. by: estimates the means within each level of a second factor. For a GLM, the means and intervals are computed on the link scale and transformed to the response scale, as emmeans’ type = "response" does.

The model’s variables come from its formula: names inside factor() and non-numeric columns are categorical. For a matrix fit, pass design: and variables: ['cyl' => 'categorical', 'wt' => 'numeric'].

use Phiz\Datasets;

use function Phiz\marginal_means;

// Phit\fit_lm('mpg ~ factor(cyl) + factor(am) + wt', $mtcars)
$fit = json_decode(file_get_contents(__DIR__.'/fits/mtcars-lm-factors.json'), true);

$plot = marginal_means($fit, Datasets::mtcars(), 'cyl', by: 'am')
    ->labs(x: 'Cylinders', color: 'Manual');

use Phiz\Datasets;

use function Phiz\marginal_means;

// Phit\fit_glm('am ~ factor(cyl) + wt', $mtcars, 'binomial'): means on the probability scale.
$glm = json_decode(file_get_contents(__DIR__.'/fits/mtcars-glm-factor.json'), true);

$plot = marginal_means($glm, Datasets::mtcars(), 'cyl')
    ->labs(x: 'Cylinders', y: 'P(manual), marginal mean (95% CI)');

The linear-model means match a direct computation in R to 10 digits. GLM means match too. Their intervals differ from R’s by about 1e-5. The difference comes from the fit’s covariance matrix. For am ~ wt + hp, phit’s coefficients agree with R’s glm() to 12 digits, while the variance of the wt coefficient is 9.4167679 in phit and 9.4167669 in R’s vcov().

Random effects

random_effects_plot() draws a caterpillar plot of a mixed model’s random effects. It works with fits from fit_lmer and fit_glmer. Current phit fits report the conditional variances themselves (effects_cond_var, and random_effects_cond_var for GLMMs), and the plot uses them directly. For older fits, pass the data the model was fitted on, and phiz rebuilds the variances from it. Each group’s conditional mode is shown with an interval from its conditional variance. Groups are sorted by their first effect, and random slopes get their own panels. The conditional covariance for group j is (Z_j’W_jZ_j / σ² + G⁻¹)⁻¹, treating the fixed effects and variance components as known, as lme4’s ranef(fit, condVar = TRUE) does. W holds the working weights of binomial and Poisson GLMMs. For intercepts, the standard deviations match lme4 to 1e-8. group: selects a term when the model has several. When the fit is singular, with a variance estimated at zero or random effects perfectly correlated, the subtitle says so. The intervals still follow lme4, which computes them without inverting G.

use Phiz\Datasets;

use function Phiz\random_effects_plot;

// Phit\fit_lmer('hwy ~ displ + (1 + displ | manufacturer)', $mpg)
$mixed = json_decode(file_get_contents(__DIR__.'/fits/mpg-lmer-slopes.json'), true);

$plot = random_effects_plot($mixed, Datasets::mpg());

Posterior intervals

interval_plot() draws posterior intervals by variable, as bayesplot’s mcmc_intervals() does. It shows a thick inner interval (50% by default), a thin outer interval (90%), and the median. It takes draws in the shape CmdStanPHP’s draws_columns() returns, or a summary from summary_columns(). Selecting a base name such as theta selects all of its elements. These draws come from the eight schools model, sampled with CmdStanPHP.

use function Phiz\interval_plot;

// Draws as CmdStanPHP's $fit->draws_columns() returns them: variable => draws, plus .chain, .iteration, .draw.
$draws = json_decode(file_get_contents(__DIR__.'/fits/eight-schools-draws.json'), true);

$plot = interval_plot($draws, variables: ['mu', 'tau', 'theta']);