Struct Fit 

pub struct Fit<'a, M, F>
where
    M: Glm,
    F: Float,
{
    pub result: Array1<F>,
    pub options: FitOptions<F>,
    pub model_like: F,
    pub n_iter: usize,
    pub history: Vec<IrlsStep<F>>,
    /* private fields */
}

the result of a successful GLM fit

Fields

result: Array1<F>

The parameter values that maximize the likelihood as given by the IRLS regression. If the dataset was internally standardized, this is transformed back

options: FitOptions<F>

The options used for this fit.

model_like: F

The value of the likelihood function for the fit result.

n_iter: usize

The number of overall iterations taken in the IRLS.

history: Vec<IrlsStep<F>>

The history of guesses and likelihoods over the IRLS iterations.

Implementations

impl<'a, M, F> Fit<'a, M, F>

where M: Glm, F: 'static + Float,

pub fn aic(&self) -> F

Returns the Akaike information criterion for the model fit.

\text{AIC} = D + 2K - 2\sum_{i} \ln w_{i}

where $D$ is the deviance, $K$ is the number of parameters, and $w_i$ are the variance weights. This is unique only to an additive constant, so only differences in AIC are meaningful.

pub fn bic(&self) -> F

Returns the Bayesian information criterion for the model fit.

\text{BIC} = K \ln(n) - 2l - 2\sum_{i} \ln w_{i}

where $K$ is the number of parameters, $n$ is the number of observations, and $l$ is the log-likelihood (including the variance weight normalization terms).

pub fn cooks(&self) -> RegressionResult<Array1<F>, F>

The Cook’s distance for each observation, which measures how much the predicted values change when leaving out each observation.

C_i = \frac{r_i^2 \, h_i}{K \, \hat\phi \, (1 - h_i)^2}

where $r_i$ is the Pearson residual, $h_i$ is the leverage, $K$ is the rank (number of parameters), and $\hat\phi$ is the estimated dispersion.

pub fn covariance(&self) -> RegressionResult<&Array2<F>, F>

The covariance matrix of the parameter estimates. When no regularization is used, this is:

\text{Cov}[\hat{\boldsymbol\beta}] = \hat\phi \, (\mathbf{X}^\mathsf{T}\mathbf{WSX})^{-1}

When regularization is active, the sandwich form is used to correctly account for the bias introduced by the penalty:

\text{Cov}[\hat{\boldsymbol\beta}] = \hat\phi \, \mathcal{I}_\text{reg}^{-1} \, \mathcal{I}_\text{data} \, \mathcal{I}_\text{reg}^{-1}

where $\mathcal{I}_\text{reg}$ is the regularized Fisher information and $\mathcal{I}_\text{data}$ is the unregularized (data-only) Fisher information. When unregularized, $\mathcal{I}_\text{reg} = \mathcal{I}_\text{data}$ and this reduces to the standard form. The result is cached on first access.

pub fn deviance(&self) -> F

Returns the deviance of the fit:

D = -2 \left[ l(\hat{\boldsymbol\beta}) - l_\text{sat} \right]

Asymptotically $\chi^2$-distributed with ndf() degrees of freedom. The unregularized likelihood is used.

pub fn dispersion(&self) -> F

The dispersion parameter $\hat\phi$ relating the variance to the variance function: $\text{Var}[y] = \phi \, V(\mu)$.

Identically one for logistic, binomial, and Poisson regression. For families with a free dispersion (linear, gamma), estimated as:

\hat\phi = \frac{D}{\left(1 - \frac{K}{n_\text{eff}}\right) \sum_i w_i}

which reduces to $D / (N - K)$ without variance weights.

pub fn fisher(&self, params: &Array1<F>) -> Array2<F>

Returns the Fisher information (the negative Hessian of the log-likelihood) at the parameter values given:

\mathcal{I}(\boldsymbol\beta) = \mathbf{X}^\mathsf{T}\mathbf{W}\eta'^2\mathbf{S}\mathbf{X}

where $\mathbf{S} = \text{diag}(V(\mu_i))$ and $\eta'$ is the derivative of the natural parameter in terms of the linear predictor ($\eta(\omega) = g_0(g^{-1}(\omega))$ where $g_0$ is the canonical link function). The regularization is included.

pub fn hat(&self) -> RegressionResult<&Array2<F>, F>

Returns the hat matrix, also known as the “projection” or “influence” matrix:

P_{ij} = \frac{\partial \hat{y}_i}{\partial y_j}

Orthogonal to the response residuals at the fit result: $\mathbf{P}(\mathbf{y} - \hat{\mathbf{y}}) = 0$. This version is not symmetric, but the diagonal is invariant to this choice of convention.

pub fn infl_coef(&self) -> RegressionResult<Array2<F>, F>

The one-step approximation to the change in coefficients from excluding each observation:

\Delta\boldsymbol\beta^{(-i)} \approx \frac{1}{1-h_i}\,\mathcal{I}^{-1}\mathbf{x}^{(i)} w_i \eta'_i e_i

Each row $i$ should be subtracted from $\hat{\boldsymbol\beta}$ to approximate the coefficients that would result from excluding observation $i$. Exact for linear models; a one-step approximation for nonlinear models.

pub fn leverage(&self) -> RegressionResult<Array1<F>, F>

Returns the leverage $h_i = P_{ii}$ for each observation: the diagonal of the hat matrix. Indicates the sensitivity of each prediction to its corresponding observation.

pub fn loo_exact(&self) -> RegressionResult<Array2<F>, F>

Returns exact coefficients from leaving each observation out, one-at-a-time. This is a much more expensive operation than the original regression because a new one is performed for each observation.

pub fn lr_test(&self) -> F

Perform a likelihood-ratio test, returning the statistic:

\Lambda = -2 \ln \frac{L_0}{L} = -2(l_0 - l)

where $L_0$ is the null model likelihood (intercept only) and $L$ is the fit likelihood. By Wilks’ theorem, asymptotically $\chi^2$-distributed with test_ndf() degrees of freedom.

pub fn lr_test_against(&self, alternative: &Array1<F>) -> F

Perform a likelihood-ratio test against a general alternative model, not necessarily a null model. The alternative model is regularized the same way that the regression resulting in this fit was. The degrees of freedom cannot be generally inferred.

pub fn ndf(&self) -> F

Returns the residual degrees of freedom in the model, i.e. the number of data points minus the number of parameters. Not to be confused with test_ndf(), the degrees of freedom in the statistical tests of the fit parameters.

pub fn null_like(&self) -> F

Returns the likelihood given the null model, which fixes all parameters to zero except the intercept (if it is used). A total of test_ndf() parameters are constrained.

pub fn predict<S>( &self, data_x: &ArrayBase<S, Ix2>, lin_off: Option<&Array1<F>>, ) -> Array1<F>

where S: Data<Elem = F>,

Returns $\hat{\mathbf{y}} = g^{-1}(\mathbf{X}\hat{\boldsymbol\beta} + \boldsymbol\omega_0)$ given input data $\mathbf{X}$ and an optional linear offset $\boldsymbol\omega_0$. The data need not be the training data.

pub fn resid_dev(&self) -> Array1<F>

Return the deviance residuals for each point in the training data:

d_i = \text{sign}(y_i - \hat\mu_i)\sqrt{D_i}

where $D_i$ is the per-observation deviance contribution. This is usually a better choice than Pearson residuals for non-linear models.

pub fn resid_dev_std(&self) -> RegressionResult<Array1<F>, F>

Return the standardized deviance residuals (internally studentized):

d_i^* = \frac{d_i}{\sqrt{\hat\phi(1 - h_i)}}

where $d_i$ is the deviance residual, $\hat\phi$ is the dispersion, and $h_i$ is the leverage. Generally applicable for outlier detection.

pub fn resid_part(&self) -> Array2<F>

Return the partial residuals as an n_obs × n_predictors matrix. Each column $j$ contains working residuals plus the centered contribution of predictor $j$:

r^{(j)}_i = r^w_i + (x_{ij} - \bar x_j) \beta_j

The bar denotes the fully weighted column mean (combining both variance and frequency weights), which is the WLS-consistent choice: it ensures that $\sum_i w_i r^{(j)}_i = 0$ for each predictor. R’s residuals(model, type = "partial") uses only frequency weights in its centering (excluding variance weights), so results will differ when variance weights are present. For models without an intercept, no centering is applied. The intercept term is always excluded from the output.

pub fn resid_pear(&self) -> &Array1<F>

Return the Pearson residuals for each point in the training data:

r_i = \sqrt{w_i} \, \frac{y_i - \hat\mu_i}{\sqrt{V(\hat\mu_i)}}

where $V$ is the variance function and $w_i$ are the variance weights. Not scaled by the dispersion for families with a free dispersion parameter.

pub fn resid_pear_std(&self) -> RegressionResult<Array1<F>, F>

Return the standardized Pearson residuals (internally studentized):

r_i^* = \frac{r_i}{\sqrt{\hat\phi(1 - h_i)}}

where $r_i$ is the Pearson residual and $h_i$ is the leverage. These are expected to have unit variance.

pub fn resid_resp(&self) -> Array1<F>

Return the response residuals: $e_i^\text{resp} = y_i - \hat\mu_i$.

pub fn resid_student(&self) -> RegressionResult<Array1<F>, F>

Return the externally studentized residuals:

\tilde{t}_i = \text{sign}(e_i) \sqrt{\frac{D_i^{(-i)}}{\hat\phi^{(-i)}}}

where $D_i^{(-i)}$ and $\hat\phi^{(-i)}$ are the LOO deviance and dispersion approximated via one-step deletion. Under normality, $t$-distributed with $N - K - 1$ degrees of freedom. This is a robust and general method for outlier detection.

pub fn resid_work(&self) -> Array1<F>

Returns the working residuals:

e_i^\text{work} = g'(\hat\mu_i)\,(y_i - \hat\mu_i) = \frac{y_i - \hat\mu_i}{\eta'(\omega_i)\,V(\hat\mu_i)}

where $g'(\mu)$ is the derivative of the link function and $\eta'(\omega)$ is the derivative of the natural parameter with respect to the linear predictor. For canonical links $\eta'(\omega) = 1$, reducing this to $(y_i - \hat\mu_i)/V(\hat\mu_i)$.

These can be interpreted as the residual differences mapped into the linear predictor space of $\omega = \mathbf{x}\cdot\boldsymbol{\beta}$.

pub fn score(&self, params: Array1<F>) -> Array1<F>

Returns the score function $\nabla_{\boldsymbol\beta} l$ (the gradient of the regularized log-likelihood) at the parameter values given. Should be zero at the MLE. The input and output are in the external (unstandardized) parameter space.

pub fn score_test(&self) -> RegressionResult<F, F>

Returns the score test statistic:

S = \mathbf{J}(\boldsymbol\beta_0)^\mathsf{T} \, \mathcal{I}(\boldsymbol\beta_0)^{-1} \, \mathbf{J}(\boldsymbol\beta_0)

where $\mathbf{J}$ is the score and $\mathcal{I}$ is the Fisher information, both evaluated at the null parameters. Asymptotically $\chi^2$-distributed with test_ndf() degrees of freedom.

pub fn score_test_against( &self, alternative: Array1<F>, ) -> RegressionResult<F, F>

Returns the score test statistic compared to another set of model parameters, not necessarily a null model. The degrees of freedom cannot be generally inferred.

pub fn test_ndf(&self) -> usize

The degrees of freedom for the likelihood ratio test, the score test, and the Wald test. Not to be confused with ndf(), the degrees of freedom in the model fit.

pub fn wald_test(&self) -> F

Returns the Wald test statistic:

W = (\hat{\boldsymbol\beta} - \boldsymbol\beta_0)^\mathsf{T} \, \mathcal{I}(\boldsymbol\beta_0) \, (\hat{\boldsymbol\beta} - \boldsymbol\beta_0)

Compared to a null model with only an intercept (if one is used). Asymptotically $\chi^2$-distributed with test_ndf() degrees of freedom.

pub fn wald_test_against(&self, alternative: &Array1<F>) -> F

Returns the Wald test statistic compared to another specified model fit instead of the null model. The degrees of freedom cannot be generally inferred.

pub fn wald_z(&self) -> RegressionResult<Array1<F>, F>

Returns the per-parameter Wald $z$-statistic:

z_k = \frac{\hat\beta_k}{\sqrt{\text{Cov}[\hat{\boldsymbol\beta}]_{kk}}}

Since it does not account for covariance between parameters it may not be accurate.

impl<'a, F> Fit<'a, Linear, F>

where F: 'static + Float,

Specialized functions for OLS.

pub fn r_sq(&self) -> F

Returns the coefficient of determination $R^2$:

R^2 = 1 - \frac{D}{D_0}

where $D$ is the deviance of the fitted model and $D_0$ is the null deviance (deviance of the intercept-only model). For Gaussian with no weights this reduces to $1 - \text{RSS}/\text{TSS}$. Using the deviance ratio correctly handles variance and frequency weights.

pub fn rss(&self) -> F

Returns the weighted residual sum of squares (RSS): $\text{RSS} = \sum_i w_i (y_i - \hat\mu_i)^2$ where $w_i$ combines variance and frequency weights.

impl<'a, M, F> Fit<'a, M, F>

where M: Glm, F: 'static + Float,

Methods that require the stats feature for distribution CDF evaluation.

pub fn pvalue_lr_test(&self) -> F

Returns the p-value for the omnibus likelihood-ratio test (full model vs. null/intercept-only model).

The LR statistic is asymptotically $\chi^2$-distributed with test_ndf() degrees of freedom, so the p-value is the upper-tail probability:

p = 1 - F_{\chi^2}(\Lambda;\, \text{test\_ndf})

pub fn pvalue_wald(&self) -> RegressionResult<Array1<F>, F>

Returns per-parameter p-values from the Wald $z$-statistics.

The reference distribution depends on whether the family has a free dispersion parameter:

• No dispersion (logistic, Poisson): standard normal — two-tailed $p_k = 2\bigl[1 - \Phi(|z_k|)\bigr]$
• Free dispersion (linear): Student-$t$ with ndf() degrees of freedom — two-tailed $p_k = 2\bigl[1 - F_t(|z_k|;\, \text{ndf})\bigr]$

IMPORTANT: Note that this test is an approximation for non-linear models and is known to sometimes yield misleading values compared to an exact test. It is not hard to find it give p-values that may imply significantly different conclusions for your analysis (e.g. p<0.07 vs. p<0.02 in one of our tests).

pub fn pvalue_exact(&self) -> RegressionResult<Array1<F>, F>

Returns per-parameter p-values from drop-one analysis of deviance (exact, expensive).

For each parameter $k$, a reduced model is fit with that parameter removed and the deviance difference $\Delta D = D_{\text{reduced}} - D_{\text{full}}$ is used:

• No dispersion: $p = 1 - F_{\chi^2}(\Delta D;\, 1)$
• Free dispersion: $F = \Delta D / \hat\phi$, $p = 1 - F_F(F;\, 1,\, \text{ndf})$

For the intercept (if present), the reduced model is fit without an intercept. For all other parameters, the reduced model is fit with that column removed from the design matrix.

impl<'a, M, F> Fit<'a, M, F>

where M: Glm + Response, F: 'static + Float,

pub fn resid_quantile(&self) -> Array1<F>

where <M as Response>::DistributionType: ContinuousCDF<f64, f64>,

Returns the quantile residuals, which are standard-normal distributed by construction if the residuals are distributed according to the GLM family’s response function.

q_i = \Phi^{-1}(F(y_i | \mathbf{x}_i \cdot \boldsymbol{\beta}))

where $\Phi$ is the standard normal quantile function and $F$ is the CDF of the response distribution, conditioned on the predicted mean.

Only implemented for families with continuous response distributions. The variance weights are used to scale the spread for families that use free dispersion. For the Linear model, this is an expensive way of getting $\frac{(y_i - \hat \mu_i)}{\sqrt{\hat\phi / w_i}} = \frac{(y_i - \hat \mu_i)}{\hat\sigma_i}$.

These residuals are not standardized/studentized in any way, meaning the impact of each observation is present in the corresponding response distribution for that point.