Skip to content

ALM

Jump to bottom
Ivan Svetunkov edited this page Oct 22, 2026 · 9 revisions

ALM — Augmented Linear Model

The Augmented Linear Model (alm() in R, ALM in Python) is the core estimator of greybox. It is a regression whose parameters are estimated by maximising the likelihood of a chosen distribution or by minimising another loss, which covers continuous, positive, bounded, count and binary responses. It can also include an occurrence part (for responses with zeros), a model for the scale, and autoregressive and integrated terms of the response.

R uses the formula interface: alm(y ~ x1 + x2, data, distribution="dnorm"). Python follows scikit-learn: ALM(distribution="dnorm").fit(X, y), with formula() turning an R-style formula and a DataFrame into (y, X).

Function signatures

R

alm(formula, data, subset, na.action,
    distribution=c("dnorm","dlaplace","ds","dgnorm","dlogis","dt","dalaplace",
                   "dlnorm","dllaplace","dls","dlgnorm","dbcnorm",
                   "dinvgauss","dgamma","dexp",
                   "dfnorm","drectnorm",
                   "dpois","dnbinom","dbinom","dgeom",
                   "dbeta","dlogitnorm",
                   "plogis","pnorm"),
    loss=c("likelihood","MSE","MAE","HAM","LASSO","RIDGE","ROLE","QUALE"),
    occurrence=c("none","plogis","pnorm"),
    scale=NULL,
    orders=c(0,0,0),
    parameters=NULL, fast=FALSE, ...)

predict.alm(object, newdata=NULL, interval=c("none", "confidence", "prediction"),
            level=0.95, side=c("both","upper","lower"),
            occurrence=NULL, ...)

forecast.alm(object, newdata=NULL, h=NULL, ...)

Python

class ALM:
    def __init__(
        self,
        distribution="dnorm",
        loss="likelihood",
        occurrence="none",
        scale_formula=None,
        orders=(0, 0, 0),
        alpha=None,
        shape=None,
        lambda_bc=None,
        size=None,
        nu=None,
        trim=None,
        lambda_l1=None,
        lambda_l2=None,
        nlopt_kargs=None,
        verbose=0,
    ): ...

    def fit(self, X, y, formula=None, feature_names=None): ...

    def predict(
        self,
        X: np.ndarray | None = None,
        interval: Literal["none", "confidence", "prediction"] = "none",
        level: float | list[float] = 0.95,
        side: Literal["both", "upper", "lower"] = "both",
        occurrence: ArrayLike | None = None,
    ) -> PredictionResult: ...

Parameters

Parameter R Python Default Description
Formula formula — — y ~ x1 + x2. Python builds (y, X) with formula() and passes them to fit(); fit(formula=) stores the formula string.
Data data X, y of fit() — R: a data frame or matrix. Python: the design matrix and the response.
Feature names — feature_names of fit() None Names of the columns of X when it is an array.
Subset / NA action subset, na.action — — R only; subset the data before calling fit() in Python.
Distribution distribution distribution "dnorm" See Distributions.
Loss loss loss "likelihood" See Loss functions.
Occurrence occurrence occurrence "none" "none", "plogis", "pnorm", or a fitted occurrence model (alm / ALM). R also accepts a formula. With an occurrence part, the model is estimated in two steps: occurrence, then the sizes on the non-zero observations.
Scale model scale scale_formula NULL / None Regressors of the scale, with a log link. R: a formula. Python: a one-sided formula string ("~ x1 + x2"), the design matrix of the scale, or a ScaleModel from sm(). Needs loss="likelihood".
ARIMA orders orders orders c(0,0,0) / (0, 0, 0) (p, d, q) of the response. Only AR and I are supported; q must be 0.
Parameters parameters — NULL R only: a vector of parameters to use instead of estimating them.
Fast mode fast — FALSE R only: skips the checks of variability and multicollinearity of the regressors.
ALaplace level alpha (...) alpha estimated Quantile level of dalaplace.
Shape shape (...) shape estimated Shape of dgnorm and dlgnorm.
Box-Cox lambda lambdaBC (...) lambda_bc estimated Box-Cox parameter of dbcnorm.
Size size (...) size estimated Size of the Negative Binomial. The size of the Binomial is always the number of distinct values of y minus one.
Degrees of freedom nu (...) nu estimated Degrees of freedom of dt (and dchisq in Python).
Trim trim (...) trim 0.05 / None (= 0.05) Share of the sample trimmed by loss="ROLE".
Shrinkage lambda (...) lambda_l1 (LASSO), lambda_l2 (RIDGE) 0 with a warning Strength of shrinkage between 0 and 1. R uses one lambda for LASSO, RIDGE and QUALE.
Fisher information FI (...) — FALSE R: return $FI.
Optimiser settings B, algorithm, maxeval, … (...) nlopt_kargs (dict) see Optimiser
Verbosity print_level (...) verbose 0

Distributions

R alm() accepts 25 distributions. Python ALM accepts the same 25 plus dchisq. The location mu is linked to the regressors X B as shown.

Code Name Location Extra parameter R Python
dnorm Normal X B — Yes Yes
dlaplace Laplace X B — Yes Yes
ds S X B — Yes Yes
dgnorm Generalised Normal X B shape Yes Yes
dlogis Logistic X B — Yes Yes
dt Student's t X B nu Yes Yes
dalaplace Asymmetric Laplace X B alpha Yes Yes
dlnorm Log-Normal X B on log(y) — Yes Yes
dllaplace Log-Laplace X B on log(y) — Yes Yes
dls Log-S X B on log(y) — Yes Yes
dlgnorm Log-Generalised Normal X B on log(y) shape Yes Yes
dbcnorm Box-Cox Normal X B on the Box-Cox transform lambdaBC / lambda_bc Yes Yes
dinvgauss Inverse Gaussian exp(X B) — Yes Yes
dgamma Gamma exp(X B) — Yes Yes
dexp Exponential exp(X B) — Yes Yes
dchisq Chi-squared (X B)^2 nu — Yes
dfnorm Folded Normal X B — Yes Yes
drectnorm Rectified Normal X B — Yes Yes
dpois Poisson exp(X B) — Yes Yes
dnbinom Negative Binomial exp(X B) size Yes Yes
dbinom Binomial exp(X B) size (from the data) Yes Yes
dgeom Geometric exp(X B) — Yes Yes
dbeta Beta exp(X B) for each of the two shapes — Yes Yes
dlogitnorm Logit-Normal X B on logit(y) — Yes Yes
plogis Logistic CDF (binary response) X B — Yes Yes
pnorm Normal CDF (binary response) X B — Yes Yes

dbeta estimates two sets of parameters, one for each shape. The densities themselves are described on distributions.

Loss functions

Loss Description R Python
likelihood Maximum likelihood of distribution (default). Yes Yes
MSE Mean Squared Error. Yes Yes
MAE Mean Absolute Error. Yes Yes
HAM Half Absolute Moment. Yes Yes
LASSO MSE with an L1 penalty on the parameters, divided by the standard deviations of the regressors. lambda lambda_l1
RIDGE MSE with an L2 penalty, likewise. lambda lambda_l2
ROLE RObust Likelihood Estimator: the likelihood with a trimmed mean of its terms. trim trim
QUALE QUAntile Likelihood Estimator: minimises a quantile of the likelihood terms. lambda —
custom A function loss(actual, fitted, B, xreg) returning the value to minimise. Yes Yes

QUALE is R only (see Roadmap).

Usage

R

library(greybox)

model <- alm(mpg ~ ., mtcars, distribution="dlnorm", subset=c(1:30))
summary(model)

predict(model, newdata=tail(mtcars, 2), interval="prediction", level=c(0.9,0.95))
#          Mean Lower bound (5%) Lower bound (2.5%) Upper bound (95%)
# [1,] 12.72061         9.225253           8.674573          17.54032
# [2,] 23.72897        19.264306          18.510198          29.22835
#      Upper bound (97.5%)
# [1,]            18.65381
# [2,]            30.41912

Python

from greybox import ALM, formula, mtcars

y_vec, X_mat = formula("mpg ~ .", mtcars.head(30))
model = ALM(distribution="dlnorm")
model.fit(X_mat, y_vec)
print(model.summary())

_, X_new = formula("mpg ~ .", mtcars.tail(2))
model.predict(X_new, interval="prediction", level=[0.9, 0.95])
#         mean    lower_0    lower_1    upper_0    upper_1
# 0  12.720608   9.225253   8.674573  17.540319  18.653815
# 1  23.728968  19.264306  18.510198  29.228353  30.419119

The estimates, standard errors and information criteria of the two are the same (AICc 165.5955).

Fitted model

Elements

Element R Python Description
Parameters coefficients, coef() coefficients All the estimated parameters, intercept first.
Slopes / intercept — coef, intercept_ Python only: coef excludes the intercept, which is in intercept_.
All optimised values B B_ The optimised vector, the extra parameter of the distribution first.
Fitted values fitted, fitted() fitted, fitted_values_
Residuals residuals, residuals() residuals, residuals_
Location mu mu_ The location parameter of the distribution.
Scale scale scale The variance σ² for dnorm, dlnorm (of log(y)), dbcnorm, dlogitnorm, dfnorm and drectnorm; the scale otherwise. A scale model when one was estimated.
Scale model scale (class "scale") scale_model
Occurrence model occurrence occurrence_ The fitted occurrence model, if any.
Actuals data[,1], actuals() actuals, data
Distribution / loss distribution, loss distribution_, loss_
Loss value lossValue loss_value
Log-likelihood logLik loglik, log_lik
Degrees of freedom df, df.residual, rank df_, df_residual_, rank_
Extra parameters other other_ alpha, shape, nu, size, …
Fisher information FI (with FI=TRUE) —
Optimiser output res nlopt_result_, n_iter_
Time elapsed timeElapsed time_elapsed
ARIMA — arima_string_, arima_polynomial_ Python only: e.g. "ARIMA(1,0,0)" and the AR+I polynomial keyed by the lag names.
Call, terms call, terms, xlevels, subset formula

Methods

Method R Python
Summary summary(model, level=0.95, bootstrap=FALSE) model.summary(level=0.95, bootstrap=False)
Prediction predict(model, newdata, interval, level, side, occurrence) model.predict(X, interval, level, side, occurrence)
Forecast forecast(model, newdata=NULL, h=NULL) model.predict(X)
Confidence intervals confint(model, parm, level=0.95, bootstrap=FALSE) model.confint(parm, level=0.95, bootstrap=False)
Covariance matrix vcov(model, bootstrap=FALSE) model.vcov(bootstrap=False)
Information criteria AIC(), AICc(), BIC(), BICc() aic, aicc, bic, bicc
Log-likelihood logLik(model) model.loglik
Number of observations / parameters nobs(), nparam() nobs, nparam, n_param
Scale and sigma extractScale(), extractSigma(), sigma() extract_scale(), extract_sigma(), sigma
Diagnostics hatvalues(), rstandard(), rstudent(), cooks.distance() hatvalues(), rstandard(), rstudent()
Diagnostic plots plot(model, which=c(1,2,4,6)) model.plot(which=(1, 2, 4, 6))
Dynamic multipliers multipliers(model, parm, h=10) multipliers(model, parm, h=10), model.multipliers(parm, h=10)
Point likelihood pointLik(model) point_lik(model)
Bootstrap of parameters coefbootstrap(model) coefbootstrap(model)
Scale model sm(model, formula), implant(model, scale) sm(model, formula), implant(model, scale)
Outliers outlierdummy(model) outlier_dummy(model)
sklearn interface — get_params(), set_params(), score(X, y, metric)

Predictions

R predict() returns a list of class "predict.greybox" with model, mean, lower, upper, level, newdata, variances, newdataProvided, location, scale and distribution.

Python predict() returns a PredictionResult with the same slots (newdata_provided for newdataProvided), plus side, interval and occurrence, its plot() method (R's plot.predict.greybox()), and to_dataframe(). It prints as a DataFrame: mean, then lower/upper, or lower_0, upper_0, … for several levels. columns, shape, index, values, len() and pred["mean"] read the DataFrame.

The bounds come from the quantile functions of the distribution: dlnorm intervals are asymmetric and positive, dpois/dnbinom intervals are integers, plogis/pnorm bounds lie in (0, 1).

Examples

Normal regression and prediction intervals

R

model <- alm(mpg ~ wt + hp, mtcars, distribution="dnorm")
summary(model)
# Response variable: mpg
# Distribution used in the estimation: Normal
# Loss function used in estimation: likelihood
# Coefficients:
#             Estimate Std. Error Lower 2.5% Upper 97.5%
# (Intercept)  37.2273     1.6271    33.8943     40.5602 *
# wt           -3.8778     0.6439    -5.1969     -2.5588 *
# hp           -0.0318     0.0092    -0.0506     -0.0129 *
#
# Error standard deviation: 2.6393
# Sample size: 32
# Number of estimated parameters: 4
# Number of degrees of freedom: 28
# Information criteria:
#      AIC     AICc      BIC     BICc
# 156.6523 158.1338 162.5153 165.0825

newdata <- data.frame(wt=c(2.5,3.5), hp=c(100,150))
predict(model, newdata, interval="prediction", level=0.95)
#          Mean Lower bound (2.5%) Upper bound (97.5%)
# [1,] 24.35540           18.81333            29.89747
# [2,] 18.88892           13.38848            24.38936

predict(model, newdata, interval="prediction", level=0.95, side="upper")
#          Mean Lower 0% Upper bound (95%)
# [1,] 24.35540     -Inf          28.95789
# [2,] 18.88892     -Inf          23.45685

Python

import pandas as pd
from greybox import ALM, formula, mtcars

y, X = formula("mpg ~ wt + hp", mtcars)
model = ALM(distribution="dnorm").fit(X, y)
print(model.summary())
# Time elapsed: 0.00 seconds
#
# Response variable: mpg
# Distribution used in the estimation: Normal
# Loss function used in estimation: likelihood
#
# Coefficients:
#                  Estimate   Std. Error   Lower 2.5%  Upper 97.5%
# ------------------------------------------------------------
#  (Intercept)      37.2273       1.6271      33.8943      40.5602 *
#           wt      -3.8778       0.6439      -5.1969      -2.5588 *
#           hp      -0.0318       0.0092      -0.0506      -0.0129 *
#
# Error standard deviation: 2.6393
# Sample size: 32
# Number of estimated parameters: 4
# Number of degrees of freedom: 28
#
# Information criteria:
#               AIC        AICc         BIC        BICc
#          156.6523    158.1338    162.5153    165.0825

model.coefficients   # [37.2273, -3.8778, -0.0318], as R's coef()
model.coef           # [-3.8778, -0.0318], without the intercept

newdata = pd.DataFrame({"mpg": [0, 0], "wt": [2.5, 3.5], "hp": [100, 150]})
_, X_new = formula("mpg ~ wt + hp", newdata)
model.predict(X_new, interval="prediction", level=0.95)
#         mean      lower      upper
# 0  24.355399  18.813332  29.897465
# 1  18.888920  13.388479  24.389362

model.predict(X_new, interval="prediction", level=0.95, side="upper")
#         mean      upper
# 0  24.355399  28.957892
# 1  18.888920  23.456846

Other distributions and losses

R

alm(mpg ~ wt + hp, mtcars, distribution="dlaplace")
alm(mpg ~ wt + hp + qsec, mtcars, loss="LASSO", lambda=0.5)
alm(mpg ~ wt + hp + qsec, mtcars, loss="ROLE")      # trim=0.05

lossFunction <- function(actual, fitted, B, xreg) return(mean(abs(actual-fitted)))
alm(mpg ~ wt + hp, mtcars, loss=lossFunction)

Python

import numpy as np

ALM(distribution="dlaplace").fit(X, y)
y3, X3 = formula("mpg ~ wt + hp + qsec", mtcars)
ALM(loss="LASSO", lambda_l1=0.5).fit(X3, y3)
ALM(loss="ROLE").fit(X3, y3)                         # trim=0.05

ALM(loss=lambda actual, fitted, B, xreg: np.mean(np.abs(actual - fitted))).fit(X, y)

The two languages give the same parameters for each of these, e.g. 33.2284, -3.6756, -0.0239, 0.1227 for LASSO with lambda=0.5 and 27.3859, -4.0871, -0.0147, 0.4285 for ROLE.

AR(1) model

orders=c(p, d, 0) adds p lags of the response (and d differences) as regressors, named yLag1, yLag2, … after the response. Without new regressors, R's forecast() produces h steps ahead (predict() without newdata returns the in-sample values); Python's predict() takes a matrix with the intercept column for each step.

R

ts_data <- c(2.1, 3.4, 3.0, 4.2, 5.1, 4.8, 5.5, 6.0, 5.8, 6.5,
             7.1, 6.8, 7.5, 8.0, 7.8)
model_ar <- alm(y ~ 1, data.frame(y=ts_data), distribution="dnorm", orders=c(1,0,0))
summary(model_ar)
# Response variable: y
# Distribution used in the estimation: Normal
# Loss function used in estimation: likelihood
# ARIMA(1,0,0) components were included in the model
# Coefficients:
#             Estimate Std. Error Lower 2.5% Upper 97.5%
# (Intercept)   0.8964     0.4419    -0.0664      1.8593
# yLag1         0.9018     0.0803     0.7269      1.0767 *
# ...

forecast(model_ar, h=3, interval="prediction", level=0.9)
# Warning: No newdata provided, the values will be forecasted
#          Mean Lower bound (5%) Upper bound (95%)
# [1,] 7.930359         6.808784          9.051935
# [2,] 8.047916         6.590294          9.505537
# [3,] 8.153926         6.471261          9.836591

Python

import numpy as np
from greybox import ALM, formula

ts = np.array([2.1, 3.4, 3.0, 4.2, 5.1, 4.8, 5.5, 6.0, 5.8, 6.5,
               7.1, 6.8, 7.5, 8.0, 7.8])
y, X = formula("y ~ 1", {"y": ts})

model_ar = ALM(distribution="dnorm", orders=(1, 0, 0))
model_ar.fit(X, y)
print(model_ar.summary())   # the same estimates as R: 0.8964 and 0.9018

model_ar.arima_string_       # 'ARIMA(1,0,0)'
model_ar.arima_polynomial_   # {'yLag1': 0.9017846462831005}

model_ar.predict(np.ones((3, 1)), interval="prediction", level=0.9)
#        mean     lower     upper
# 0  7.930359  6.808784  9.051935
# 1  8.047916  6.590294  9.505537
# 2  8.153926  6.471261  9.836591

Distributed lags

An autoregressive distributed lag (ARDL) model combines the AR terms of orders with lags of a regressor made by the backshift function B() inside the formula. multipliers() then gives the dynamic multipliers of that regressor. Both functions, and the related xregExpander(), are documented on manipulations.

R

model <- alm(drivers ~ kms + B(kms, 1) + B(kms, 2),
             Seatbelts, distribution="dnorm", orders=c(1, 0, 0))
summary(model)
# ...
#             Estimate Std. Error Lower 2.5% Upper 97.5%
# (Intercept) 673.4045   179.3085   319.6646   1027.1443 *
# kms          -0.0213     0.0140    -0.0489      0.0062
# B(kms, 1)    -0.0194     0.0181    -0.0550      0.0163
# B(kms, 2)     0.0390     0.0124     0.0146      0.0634 *
# driversLag1   0.6137     0.0664     0.4828      0.7446 *
# ...

round(multipliers(model, "kms", h=10), 4)
#      h1      h2      h3      h4      h5      h6      h7      h8      h9     h10
# -0.0213 -0.0324  0.0191  0.0117  0.0072  0.0044  0.0027  0.0017  0.0010  0.0006

Python

import numpy as np
from greybox import ALM, formula, B, multipliers

rng = np.random.default_rng(42)
n = 60
x = rng.normal(size=n).cumsum()
y = 0.5 * x + 0.3 * B(x, 1) + 0.2 * B(x, 2) + rng.normal(size=n) * 0.5

y_vec, X_mat = formula("y ~ x + B(x, 1) + B(x, 2)", {"y": y, "x": x})
model = ALM(distribution="dnorm", orders=(1, 0, 0))
model.fit(X_mat, y_vec)
print(model.summary())
# ...
#                  Estimate   Std. Error   Lower 2.5%  Upper 97.5%
# ------------------------------------------------------------
#  (Intercept)      -0.0769       0.0574      -0.1921       0.0382
#            x       0.4669       0.0668       0.3330       0.6007 *
#       B(x,1)       0.3223       0.1084       0.1049       0.5397 *
#       B(x,2)       0.1182       0.0941      -0.0704       0.3068
#        yLag1       0.0861       0.1260      -0.1665       0.3387
# ...

multipliers(model, "x", h=5)
# {'h1': 0.4669, 'h2': 0.3625, 'h3': 0.1494, 'h4': 0.0129, 'h5': 0.0011}   (rounded)

The multipliers are the impulse responses of y to a unit change of the regressor: h1 is the coefficient of x, and h2 = coefficient of B(x,1) + φ₁·h1, and so on. They are not cumulative. Python names the lagged columns B(x,1) (without the space R puts in B(kms, 1)).

Optimiser

Both languages optimise with NLopt (R through nloptr, Python through nlopt). R takes the settings in ..., Python in the nlopt_kargs dict, with the same names and defaults.

Setting Default Description
algorithm "NLOPT_LN_NELDERMEAD" NLopt algorithm.
maxeval 40 per parameter; 200 per parameter for LASSO / RIDGE and the ARI models Maximum number of evaluations. dnorm, dlnorm and dlogitnorm with likelihood or MSE and no AR terms or scale model are solved in closed form (maxeval=1), as is RIDGE.
maxtime no limit Maximum time in seconds.
xtol_rel 1e-6 Relative tolerance on the parameters.
xtol_abs 1e-8 Absolute tolerance on the parameters.
ftol_rel 1e-4; 1e-8 for LASSO / RIDGE and the ARI models Relative tolerance on the loss.
ftol_abs 0 Absolute tolerance on the loss.
print_level 0 Optimiser output.
B computed R only: starting values of the parameters (location, then scale).

Python accepts these algorithm names: NLOPT_LN_NELDERMEAD, NLOPT_LN_SBPLX, NLOPT_LN_COBYLA, NLOPT_LN_BOBYQA, NLOPT_LN_NEWUOA, NLOPT_LN_PRAXIS, NLOPT_LD_LBFGS, NLOPT_LD_SLSQP, NLOPT_LD_MMA.

R

alm(mpg ~ wt + hp, mtcars, distribution="dlaplace",
    algorithm="NLOPT_LN_SBPLX", maxeval=1000, xtol_rel=1e-8)

Python

ALM(distribution="dlaplace",
    nlopt_kargs={"algorithm": "NLOPT_LN_SBPLX", "maxeval": 1000, "xtol_rel": 1e-8}).fit(X, y)

R vs Python differences

Aspect R Python
Interface alm(formula, data) ALM().fit(X, y), with formula() for (y, X)
dchisq — Accepted
QUALE loss Accepted —
fast, parameters, subset, na.action Accepted —
Shrinkage parameter lambda lambda_l1 / lambda_l2
coef All parameters coefficients holds all parameters; coef omits the intercept
Forecasting an AR model forecast(model, h) predict(np.ones((h, 1)))
sklearn interface — get_params(), set_params(), score()

Numerical parity is described on R-Python-differences.

References

Clone this wiki locally