API

Types defined in the package

RobustModels.AbstractEstimatorType

A robust estimator is a location or scale estimator associated to one (or more) loss function and solved using Iteratively Reweighted Least Square.

source
RobustModels.LossFunctionType

An estimator needs a cost/loss function for the modified (weighted) least squares problems of the form:

\[\min \sum_i \rho\left(\dfrac{r_i}{\hat{\sigma}} ight)\]

source
RobustModels.DensePredCGType
DensePredCG

A LinPred type with Conjugate Gradient and a dense X

Members

  • X: Model matrix of size n × p with n ≥ p. Should be full column rank.
  • beta0: base coefficient vector of length p
  • delbeta: increment to coefficient vector, also of length p
  • scratchbeta: scratch vector of length p, used in linpred! method
source
RobustModels.SparsePredCGType
SparsePredCG

A LinPred type with Conjugate Gradient and a sparse X

Members

  • X: Model matrix of size n × p with n ≥ p. Should be full column rank.
  • beta0: base coefficient vector of length p
  • delbeta: increment to coefficient vector, also of length p
  • scratchbeta: scratch vector of length p, used in linpred! method
source
GLM.DensePredCholType
DensePredChol{T}

A LinPred type with a dense Cholesky factorization of X'X

Members

  • X: model matrix of size n × p with n ≥ p. Should be full column rank.
  • beta0: base coefficient vector of length p
  • delbeta: increment to coefficient vector, also of length p
  • scratchbeta: scratch vector of length p, used in linpred! method
  • chol: a Cholesky object created from X'X, possibly using row weights.
  • scratchm1: scratch Matrix{T} of the same size as X
  • scratchm2: scratch Matrix{T} os the same size as X'X
GLM.SparsePredCholType
SparsePredChol{T<:BlasReal} <: LinPred

A LinPred type with a sparse Cholesky factorization of X'X

Members

  • X: model matrix of size n×p with n ≥ p. Should be full column rank.
  • beta0: base coefficient vector of length p
  • delbeta: increment to coefficient vector, also of length p
  • scratchbeta: scratch vector of length p, used in linpred! method
  • chol: a Cholesky object created from X'X, possibly using row weights.
source
RobustModels.RidgePredType
RidgePred

Regularized predictor using ridge regression on the p features.

Members

  • X: model matrix
  • λ: shrinkage parameter of the regularizer
  • G: regularizer matrix of size p×p.
  • βprior: regularizer prior of the coefficient values. Default to zeros(p).
  • pred: the non-regularized predictor using an extended model matrix.
  • pivot: for DensePredChol, if the decomposition was pivoted.
  • scratchbeta: scratch vector of length p, used in linpred! method
source
RobustModels.RobustLinRespType
RobustLinResp

Robust linear response structure.

Solve the following minimization problem:

\[\min \sum_i \rho\left(\dfrac{r_i}{\hat{\sigma}} ight)\]

Fields

  • est: estimator used for the model
  • y: response vector
  • μ: mean response vector
  • offset: offset added to to form μ. Can be of length 0
  • wts: prior case weights. Can be of length 0.
  • σ: current estimate of the scale or dispersion
  • devresid: the deviance residuals
  • wrkwt: working case weights for the Iteratively Reweighted Least Squares (IRLS) algorithm
  • wrkres: working residuals for IRLS
  • wrkscaledres: scaled residuals for IRLS
source
RobustModels.QuantileRegressionType
QuantileRegression

Quantile regression representation

Fields

  • τ: the quantile value
  • X: the model matrix
  • β: the coefficients
  • y: the response vector
  • wts: the weights
  • wrkres: the working residuals
  • fitdispersion: if true, the dispersion is estimated otherwise it is kept fixed
  • fitted: if true, the model was already fitted
source

Constructors for models

StatsBase.fitMethod
fit(::Type{M},
    X::Union{AbstractMatrix{T},SparseMatrixCSC{T}},
    y::AbstractVector{T},
    est::Estimator;
    method::Symbol       = :chol, # :cg
    dofit::Bool          = true,
    wts::FPVector        = similar(y, 0),
    offset::FPVector     = similar(y, 0),
    fitdispersion::Bool  = false,
    ridgeλ::Real         = 0,
    ridgeG::Union{UniformScaling, AbstractArray} = I,
    βprior::AbstractVector = [],
    quantile::Union{Nothing, AbstractFloat} = nothing,
    initial_scale::Union{Symbol, Real}=:mad,
    σ0::Union{Nothing, Symbol, Real}=initial_scale,
    initial_coef::AbstractVector=[], 
    β0::AbstractVector=initial_coef, 
    correct_leverage::Bool=false
    fitargs...) where {M<:RobustLinearModel, T<:AbstractFloat}

Create a robust model with the model matrix (or formula) X and response vector (or dataframe) y, using a robust estimator.

Arguments

  • X: the model matrix (it can be dense or sparse) or a formula
  • y: the response vector or a dataframe.
  • est: a robust estimator

Keywords

  • method::Symbol = :chol: the method to use for solving the weighted linear system, chol (default) or cg;
  • dofit::Bool = true: if false, return the model object without fitting;
  • wts::Vector = []: a weight vector, should be empty if no weights are used;
  • offset::Vector = []: an offset vector, should be empty if no offset is used;
  • fitdispersion::Bool = false: reevaluate the dispersion;
  • ridgeλ::Real = 0: if positive, perform a robust ridge regression with shrinkage parameter ridgeλ. RidgePred object will be used;
  • ridgeG::Union{UniformScaling, AbstractArray} = I: define a custom regularization matrix. Default to unity matrix (with 0 for the intercept);
  • βprior::AbstractVector = []: define a custom prior for the coefficients for ridge regression. Default to zeros(p);
  • quantile::Union{Nothing, AbstractFloat} = nothing: only for GeneralizedQuantileEstimator, define the quantile to estimate;
  • initial_scale::Union{Symbol, Real}=:mad: the initial scale estimate, for non-convex estimator it helps to find the global minimum. Automatic computation using :mad, L1 or extrema (non-robust).
  • σ0::Union{Nothing, Symbol, Real}=initial_scale: alias of initial_scale;
  • initial_coef::AbstractVector=[]: the initial coefficients estimate, for non-convex estimator it helps to find the global minimum.
  • β0::AbstractVector=initial_coef: alias of initial_coef;
  • correct_leverage::Bool=false: apply the leverage correction weights with leverage_weights.
  • fitargs...: other keyword arguments used to control the convergence of the IRLS algorithm (see pirls!).

Output

the RobustLinearModel object.

source
StatsBase.fitMethod
fit(::Type{M}, X::Union{AbstractMatrix{T},SparseMatrixCSC{T}},
    y::AbstractVector{T}; quantile::AbstractFloat=0.5,
    dofit::Bool          = true,
    wts::FPVector        = similar(y, 0),
    fitdispersion::Bool  = false,
    fitargs...) where {M<:QuantileRegression, T<:AbstractFloat}

Fit a quantile regression model with the model matrix (or formula) X and response vector (or dataframe) y.

It is solved using the exact interior method.

Arguments

  • X: the model matrix (it can be dense or sparse) or a formula
  • y: the response vector or a dataframe.

Keywords

  • quantile::AbstractFloat=0.5: the quantile value for the regression, between 0 and 1.
  • dofit::Bool = true: if false, return the model object without fitting;
  • wts::Vector = []: a weight vector, should be empty if no weights are used;
  • fitdispersion::Bool = false: reevaluate the dispersion;
  • fitargs...: other keyword arguments like verbose to print iteration details.

Output

the RobustLinearModel object.

source
RobustModels.rlmFunction
rlm(X, y, args...; kwargs...)

An alias for fit(RobustLinearModel, X, y, est; kwargs...).

The arguments X and y can be a Matrix and a Vector or a Formula and a DataFrame.

source
RobustModels.quantregFunction
quantreg(X, y, args...; kwargs...)

An alias for fit(QuantileRegression, X, y; kwargs...).

The arguments X and y can be a Matrix and a Vector or a Formula and a DataFrame.

source
StatsBase.fit!Function

Fit a statistical model in-place.

fit!(m::RobustLinearModel; initial_scale::Union{Symbol, Real}=:mad,
          σ0::Union{Nothing, Symbol, Real}=initial_scale,
          initial_coef::AbstractVector=[], 
          β0::AbstractVector=initial_coef, 
          correct_leverage::Bool=false, kwargs...)

Optimize the objective of a RobustLinearModel. When verbose is true the values of the objective and the parameters are printed on stdout at each iteration.

This function assumes that m was correctly initialized.

This function returns early if the model was already fitted, instead call refit!.

source
fit!(m::QuantileRegression;
     verbose::Bool=false,
     quantile::Union{Nothing, AbstractFloat}=nothing,
     correct_leverage::Bool=false,
     kwargs...)

Optimize the objective of a QuantileRegression. When verbose is true the values of the objective and the parameters are printed on stdout at each function evaluation.

This function assumes that m was correctly initialized. This function returns early if the model was already fitted, instead call refit!.

source
RobustModels.refit!Function
refit!(m::RobustLinearModel, [y::FPVector];
                             wts::Union{Nothing, FPVector} = nothing,
                             offset::Union{Nothing, FPVector} = nothing,
                             quantile::Union{Nothing, AbstractFloat} = nothing,
                             ridgeλ::Union{Nothing, Real} = nothing,
                             kwargs...)

Refit the RobustLinearModel.

This function assumes that m was correctly initialized and the model is refitted with the new values for the response, weights, offset, quantile and ridge shrinkage.

Defining a new quantile is only possible for GeneralizedQuantileEstimator.

Defining a new ridgeλ is only possible for RidgePred objects.

source
refit!(m::QuantileRegression, [y::FPVector ; verbose::Bool=false, quantile::Union{Nothing, AbstractFloat}=nothing])

Refit the QuantileRegression model with the new values for the response, weights and quantile. This function assumes that m was correctly initialized.

source

Model methods

StatsBase.coefFunction
coef(model::StatisticalModel)

Return the coefficients of the model.

StatsBase.coeftableFunction
coeftable(model::StatisticalModel; level::Real=0.95)

Return a table with coefficients and related statistics of the model. level determines the level for confidence intervals (by default, 95%).

The returned CoefTable object implements the Tables.jl interface, and can be converted e.g. to a DataFrame via using DataFrames; DataFrame(coeftable(model)).

StatsBase.confintFunction
confint(model::StatisticalModel; level::Real=0.95)

Compute confidence intervals for coefficients, with confidence level level (by default 95%).

StatsBase.devianceFunction
deviance(model::StatisticalModel)

Return the deviance of the model relative to a reference, which is usually when applicable the saturated model. It is equal, up to a constant, to $-2 \log L$, with $L$ the likelihood of the model.

deviance(m::RobustLinearModel)

The sum of twice the loss/objective applied to the scaled residuals.

It is consistent with the definition of the deviance for OLS.

source
deviance(obj::LinearModel)

For linear models, the deviance is equal to the residual sum of squares (RSS).

StatsBase.nulldevianceFunction
nulldeviance(model::StatisticalModel)

Return the deviance of the null model, that is the one including only the intercept.

nulldeviance(obj::LinearModel)

For linear models, the deviance of the null model is equal to the total sum of squares (TSS).

StatsBase.dofFunction
dof(model::StatisticalModel)

Return the number of degrees of freedom consumed in the model, including when applicable the intercept and the distribution's dispersion parameter.

StatsBase.dof_residualFunction
dof_residual(model::RegressionModel)

Return the residual degrees of freedom of the model.

StatsBase.nobsMethod
nobs(model::StatisticalModel)

Return the number of independent observations on which the model was fitted. Be careful when using this information, as the definition of an independent observation may vary depending on the model, on the format used to pass the data, on the sampling plan (if specified), etc.

StatsBase.isfittedFunction
isfitted(model::StatisticalModel)

Indicate whether the model has been fitted.

StatsBase.islinearFunction
islinear(model::StatisticalModel)

Indicate whether the model is linear.

StatsBase.loglikelihoodFunction
loglikelihood(model::StatisticalModel)

Return the log-likelihood of the model.

loglikelihood(model::StatisticalModel, ::Colon)

Return a vector of each observation's contribution to the log-likelihood of the model. In other words, this is the vector of the pointwise log-likelihood contributions.

In general, sum(loglikehood(model, :)) == loglikelihood(model).

loglikelihood(model::StatisticalModel, observation)

Return the contribution of observation to the log-likelihood of model.

StatsBase.nullloglikelihoodFunction
loglikelihood(model::StatisticalModel)

Return the log-likelihood of the null model corresponding to model. This is usually the model containing only the intercept.

StatsBase.stderrorFunction
stderror(model::StatisticalModel)

Return the standard errors for the coefficients of the model.

StatsBase.vcovFunction
vcov(model::StatisticalModel)

Return the variance-covariance matrix for the coefficients of the model.

StatsBase.weightsFunction
weights(vs)

Construct a Weights vector from array vs. See the documentation for Weights for more details.

weights(model::StatisticalModel)

Return the weights used in the model.

RobustModels.workingweightsFunction
workingweights(m::RobustLinearModel)

The robust weights computed by the model.

This can be used to detect outliers, as outliers weights are lower than the weights of valid data points.

source
StatsBase.fittedFunction
fitted(model::RegressionModel)

Return the fitted values of the model.

StatsBase.predictFunction
predict(model::RegressionModel, [newX])

Form the predicted response of model. An object with new covariate values newX can be supplied, which should have the same type and structure as that used to fit model; e.g. for a GLM it would generally be a DataFrame with the same variable names as the original predictors.

StatsBase.leverageFunction
leverage(model::RegressionModel)

Return the diagonal of the projection matrix of the model.

StatsBase.modelmatrixFunction
modelmatrix(model::RegressionModel)

Return the model matrix (a.k.a. the design matrix).

GLM.dispersionMethod
dispersion(m::RobustLinearModel, sqr::Bool=false)

The dispersion is the (weighted) sum of robust residuals. If sqr is true, return the squared dispersion.

source
StatsBase.responseFunction
response(model::RegressionModel)

Return the model response (a.k.a. the dependent variable).

StatsBase.residualsFunction
residuals(model::RegressionModel)

Return the residuals of the model.

RobustModels.scaleFunction
scale(m::RobustLinearModel, sqr::Bool=false)

The robust scale estimate used for the robust estimation.

If sqr is true, the square of the scale is returned.

source
RobustModels.tauscaleFunction
tauscale(m::RobustLinearModel, sqr::Bool=false; kwargs...)

The robust τ-scale that is minimized in τ-estimation.

If sqr is true, the square of the τ-scale is returned.

source
RobustModels.location_varianceFunction
location_variance(r::RobustLinResp, sqr::Bool=false)

Compute the part of the variance of the coefficients β that is due to the encertainty from the location.

If sqr is false, return the standard deviation instead.

From Maronna et al., Robust Statistics: Theory and Methods, Equation 4.49

source

Estimators

RobustModels.MEstimatorType
MEstimator{L<:LossFunction} <: AbstractEstimator

M-estimator for a given loss function.

The M-estimator is obtained by minimizing the loss function:

\[\hat{\mathbf{\beta}} = \underset{\mathbf{\beta}}{\textrm{argmin}} \sum_{i=1}^n \rho\left(\dfrac{r_i}{\hat{\sigma}}\right)\]

with the residuals $\mathbf{r} = \mathbf{y} - \mathbf{X} \mathbf{\beta}$ , and a robust scale estimate $\hat{\sigma}$.

Fields

source
RobustModels.SEstimatorType
SEstimator{L<:BoundedLossFunction} <: AbstractEstimator

S-estimator for a given bounded loss function.

The S-estimator is obtained by minimizing the scale estimate:

\[\hat{\mathbf{\beta}} = \underset{\mathbf{\beta}}{\textrm{argmin }} \hat{\sigma}^2\]

where the robust scale estimate $\hat{\sigma}}$ is solution of:

\[\dfrac{1}{n} \sum_{i=1}^n \rho\left(\dfrac{r_i}{\hat{\sigma}}\right) = \delta\]

with the residuals $\mathbf{r} = \mathbf{y} - \mathbf{X} \mathbf{\beta}$ , $\rho$ is a bounded loss function with $\underset{r \to \infty}{\lim} \rho(r) = 1$ and $\delta$ is the finite breakdown point, usually 0.5.

Fields

source
RobustModels.MMEstimatorType
MMEstimator{L1<:BoundedLossFunction, L2<:LossFunction} <: AbstractEstimator

MM-estimator for the given loss functions.

The MM-estimator is obtained using a two-step process:

  1. compute a robust scale estimate with a high breakdown point using a S-estimate and the loss function L1.
  2. compute an efficient estimate using a M-estimate with the loss function L2.

Fields

  • loss1: the BoundedLossFunction used for the high breakdown point S-estimation.
  • loss2: the LossFunction used for the efficient M-estimation.
  • scaleest: boolean specifying the if the estimation is in the S-estimation step (true)

or the M-estimation step (false).

source
RobustModels.TauEstimatorType
TauEstimator{L1<:BoundedLossFunction, L2<:BoundedLossFunction} <: AbstractEstimator

τ-estimator for the given loss functions.

The τ-estimator corresponds to a M-estimation, where the loss function is a weighted sum of a high breakdown point loss and an efficient loss. The weight is recomputed at every step of the Iteratively Reweighted Least Square, so the estimate is both robust (high breakdown point) and efficient.

Fields

source
RobustModels.GeneralizedQuantileEstimatorType
GeneralizedQuantileEstimator{L<:LossFunction} <: AbstractQuantileEstimator

Generalized Quantile Estimator is an M-Estimator with asymmetric loss function.

For L1Loss, this corresponds to quantile regression (although it is better to use quantreg for quantile regression because it gives the exact solution).

For L2Loss, this corresponds to Expectile regression (see ExpectileEstimator).

Fields

  • loss: the LossFunction.
  • τ: the quantile value to estimate, between 0 and 1.

Properties

  • tau, q, quantile are aliases for τ.
source
RobustModels.ExpectileEstimatorType

The expectile estimator is a generalization of the L2 estimator, for other quantile τ ∈ [0,1].

[1] Schnabel, Eilers - Computational Statistics and Data Analysis 53 (2009) 4168–4177 - Optimal expectile smoothing doi:10.1016/j.csda.2009.05.002

source

Loss functions

RobustModels.L1LossType

The standard L1 loss function takes the absolute value of the residual, and is convex but non-smooth. It is not a real L1 loss but a Huber loss with very small tuning constant.

source
RobustModels.HuberLossType

The convex Huber loss function switches from between quadratic and linear cost/loss function at a certain cutoff.

source
RobustModels.L1L2LossType

The convex L1-L2 loss interpolates smoothly between L2 behaviour for small residuals and L1 for outliers.

source
RobustModels.FairLossType

The (convex) "fair" loss switches from between quadratic and linear cost/loss function at a certain cutoff, and is C3 but non-analytic.

source
RobustModels.CauchyLossType

The non-convex Cauchy loss function switches from between quadratic behaviour to logarithmic tails. This rejects outliers but may result in multiple minima. For scale estimate, r.ψ(r) is used as a loss, which is the same as for Geman loss.

source
RobustModels.GemanLossType

The non-convex Geman-McClure for strong supression of outliers and does not guarantee a unique solution. For the S-Estimator, it is equivalent to the Cauchy loss.

source
RobustModels.TukeyLossType

The non-convex Tukey biweight estimator which completely suppresses the outliers, and does not guaranty a unique solution

source
RobustModels.YohaiZamarLossType

The non-convex (and bounded) optimal Yohai-Zamar loss function that minimizes the estimator bias. It was originally introduced in Optimal locally robust M-estimates of regression (1997) by Yohai and Zamar with a slightly different formula.

source

Estimator and Loss functions methods

RobustModels.psiFunction

The influence function ψ is the derivative of the loss function for the M-estimator, multiplied by the square of the tuning constant.

source
RobustModels.efficiency_tuning_constantFunction

The tuning constant c is computed so the efficiency for Normal distributed residuals is 0.95. The efficiency of the mean estimate μ is defined by:

eff_μ = (E[ψ'])²/E[ψ²]

source
RobustModels.mscale_lossFunction
mscale_loss(loss::L, x)

The rho-function that is used for M-scale estimation.

For monotone (convex) functions, χ(r) = r.ψ(r).

For bounded functions, χ(r) = ρ(r)/ρ(∞) so χ(∞) = 1.

source
RobustModels.breakdown_point_tuning_constantFunction

The M-estimate of scale is computed by solving:

\[\dfrac{1}{n} \sum_i \chi\left( \dfrac{r_i}{\hat{\sigma}}\right) = \delta\]

For monotone (convex) functions, χ(r) = r.ψ(r) and δ is defined as E[χ(r)] = δ for the Normal distribution N(0,1) For bounded functions, χ(r) = ρ(r)/ρ(∞) with χ(∞) = 1 and δ = E[χ]/χ(∞) with expectation w.r.t. Normal density.

The tuning constant c corresponding to a high breakdown point (0.5) is such that δ = 1/2, from 1/n Σ χ(r/ŝ) = δ

source
RobustModels.scale_estimateFunction
scale_estimate(loss, res; σ0=1.0, wts=[], verbose=false,
                         order=1, approx=false, nmax=30, 
                         rtol=1e-4, atol=0.1)

Compute the M-scale estimate from the loss function. If the loss is bounded, ρ is used as the function χ in the sum, otherwise r.ψ(r) is used if the loss is not bounded, to coincide with the Maximum Likelihood Estimator. Also, for bounded estimator, because f(s) = 1/(nδ) Σ ρ(ri/s) is decreasing the iteration step is not using the weights but is multiplicative.

source
RobustModels.update_weight!Function
update_weight!(E::TauEstimator, res::AbstractArray{T}; wts::AbstractArray{T}=T[])

Update the weight between the two estimators of a τ-estimator using the scaled residual.

source
RobustModels.tau_scale_estimateFunction
tau_scale_estimate!(E::TauEstimator, res::AbstractArray{T}, σ::Real, sqr::Bool=false;
                    wts::AbstractArray{T}=T[], bound::AbstractFloat=0.5) where {T<:AbstractFloat}

The τ-scale estimate, where σ is the scale estimate from the robust M-scale. If sqr is true, return the squared value.

source