This site is version 0.0.1. Most of the code on this site is drafted by an agent and needs review. I apologize for the AI-speak, we are working on rewriting it!
viv-lang.org

vp.stats.glm

Generalized linear models by R's IRLS: binomial, Poisson, Gamma and Gaussian, with R's warnings.

Part of the distribution vipyr; it runs from a checkout of the repository. Planned on PyPI as pip install vipyr.

vp.stats.glm(y, X, family=..., link=...) is R's glm.fit step for step (lib/stats/glm.viv): R's starting means, the working response and weights, each weighted least-squares step by the pivoted QR of vp.stats.ols, R's clamps in the link functions, and R's stopping rule |dev - dev_old| / (|dev| + 0.1) < epsilon. The result carries R's summary.glm numbers: dispersion 1 and z tests for binomial and Poisson, the Pearson estimate and t tests for Gamma and Gaussian, the deviance, the null deviance and R's AIC. When the fitted probabilities reach 0 or 1 (a separated binomial fit) or the iterations run out, the result's warnings hold R's text and a GLMWarning is raised as a Python warning.

Two named slots. fit="irls" (the default) is Fisher scoring; fit="newton" takes Newton steps with the observed information. cov= chooses the information the covariance inverts: "expected" is R's (X^T W X)^-1; "observed" is the negative Hessian of the log-likelihood, whose per-observation second derivative is taken by deriv (forward-mode, exact) of the score. For the canonical links the two are the same; for probit and the Gamma log link they differ, and the observed information is what a likelihood-based standard error uses.

Signature

vp.stats.glm(y, X, *, family='binomial', link=None, weights=None, offset=None, fit='irls', cov='expected', epsilon=1e-08, maxit=25, names=None)

Parameters

y : ndarray[float64], required
The response: 0/1 or proportions (binomial), counts (poisson), positive values (gamma).
X : ndarray[float64], required
The n by k design matrix.
family : str, default 'binomial'
'binomial', 'poisson', 'gamma' or 'gaussian'.
link : str, default None
The link: logit or probit (binomial), log (poisson), inverse or log (gamma), identity (gaussian); None: the canonical one.
weights : ndarray[float64], default None
Prior weights (binomial: the numbers of trials).
offset : ndarray[float64], default None
A known term added to the linear predictor.
fit : str, default 'irls'
The fitting slot: 'irls' (R's glm.fit, Fisher scoring) or 'newton' (Newton's method with the observed information).
cov : str, default 'expected'
'expected': R's (X^T W X)^-1; 'observed': the inverse of the observed information, its Hessian derivative by deriv (differs for probit and gamma-log).
epsilon : float, default 1e-08
R's convergence tolerance on the relative change of the deviance.
maxit : int, default 25
The iterations allowed (R's default).
names : list, default None
Names of the columns.

Returns

A frozen GLMResult. Iterating it gives the fields below in order; it also has params, tests, log and check().

coef : ndarray[float64]
The coefficients.
bse : ndarray[float64]
Their standard errors.
tvalues : ndarray[float64]
Coef / bse (z for binomial and poisson, t for gamma and gaussian).
pvalues : ndarray[float64]
Two-sided p-values, as R's summary.glm.
conf_int : ndarray[float64]
95% Wald intervals.
deviance : float
The residual deviance.
null_deviance : float
The deviance of the intercept-only model.
df_resid : int
N minus rank.
df_null : int
N minus the intercept.
dispersion : float
1 (binomial, poisson) or the Pearson chi-squared over df_resid.
aic : float
R's AIC.
iterations : int
The iterations taken.
converged : bool
Whether R's rule stopped it.
mu : ndarray[float64]
The fitted means.
eta : ndarray[float64]
The linear predictor.
resid_deviance : ndarray[float64]
The deviance residuals.
resid_pearson : ndarray[float64]
The Pearson residuals.
cov_params : ndarray[float64]
Dispersion times the unscaled covariance.
warnings : tuple
R's warnings, as text.
family : str
The family.
link : str
The link.
fit : str
The fitting method.
cov : str
Which information the covariance inverts.
names : list
The column names.
nobs : int
N.
rank : int
The rank of X.

Raises

ValueError
A family, link, fitting method or cov that is not one of the registered.
vipyr._stats.GLMWarning
A warning, with R's text: the iteration did not converge, or fitted probabilities or rates numerically 0 or 1 (separation).
TypeError
An argument of the wrong type, or an unknown keyword.
vipyr.VivBuildError
An ImportError: no C compiler, or the program's C does not build.

Examples

import numpy as np
import vipyr as vp
x = np.arange(1.0, 11.0)
y = np.array([0, 0, 1, 0, 0, 1, 1, 0, 1, 1.0])
X = np.column_stack([np.ones(10), x])
r = vp.stats.glm(y, X, family="binomial", link="probit", names=["const", "x"])
print(r.summary())
assert r.check()
Generalized linear model (vp.stats.glm: binomial, link probit, fit irls, cov expected)
==============================================================================
No. Observations:       10     Deviance:            10.8231   Null deviance:    13.8629
Df Residuals:            8     Dispersion:                1   AIC:              14.8231
Iterations:              4     Converged:              True
==============================================================================
                  coef    std err         z     P>|z|      [0.025      0.975]
-----------------------------------------------------------------------------
const           -1.505      1.024    -1.469     0.142      -3.513      0.5025 
x               0.2737     0.1682     1.627     0.104    -0.05592      0.6033

Each example ran as written when these pages were generated, and the output shown is what it printed.

Checked against

Reference
R 4.3.3's glm on its own datasets (esoph-style binomial with trials, mtcars logit and probit, warpbreaks Poisson, a Gamma fit with both links) and a separated logistic fit, to 10−10 at R's epsilon 10−14 (lib/stats/glm_check.py, python/tests/test_stats.py).

Source