cirls.fit {cirls}R Documentation

Constrained Iteratively Reweighted Least-Squares

Description

Fits a generalized linear model with linear constraints on the coefficients through a Constrained Iteratively Reweighted Least-Squares (CIRLS) algorithm. This function is the constrained counterpart to glm.fit and is meant to be called by glm through its method argument. See details for the main differences.

Usage

cirls.fit(x, y, weights = rep.int(1, nobs), start = NULL,
  etastart = NULL, mustart = NULL, offset = rep.int(0, nobs),
  family = stats::gaussian(), control = list(), intercept = TRUE,
  singular.ok = TRUE)

Arguments

x, y

x is a design matrix and y is a vector of response observations. Usually internally computed by glm.

weights

An optional vector of observation weights.

start

Starting values for the parameters in the linear predictor.

etastart

Starting values for the linear predictor.

mustart

Starting values for the vector or means.

offset

An optional vector specifying a known component in the model. See model.offset.

family

The result of a call to a family function, describing the error distribution and link function of the model. See family for details of available family functions.

control

A list of parameters controlling the fitting process. See details and cirls.control.

intercept

Logical. Should an intercept be included in the null model?

singular.ok

Logical. If FALSE, the function returns an error for singular fits.

Details

This function is a plug-in for glm and works similarly to glm.fit. In addition to the parameters already available in glm.fit, cirls.fit allows the specification of a constraint matrix Cmat with bound vectors lb and ub on the regression coefficients. These additional parameters can be passed through the control list or through ... in glm.

The CIRLS algorithm is a modification of the classical IRLS algorithm in which each update of the regression coefficients is performed by a quadratic program (QP), ensuring the update stays within the feasible region defined by Cmat, lb and ub. More specifically, this feasible region is defined as

⁠lb <= Cmat %*% coefficients <= ub⁠

where coefficients is the coefficient vector returned by the model. This specification allows for any linear constraint, including equality ones.

Specifying Cmat, lb and ub

Cmat is a matrix that defines the linear constraints. If provided directly as a matrix, the number of columns in Cmat must match the number of coefficients estimated by glm. This includes all variables that are not involved in any constraint potential expansion such as factors or splines for instance, as well as the intercept. Columns not involved in any constraint will be filled by 0s.

Alternatively, it may be more convenient to pass Cmat as a list of constraint matrices for specific terms. This is advantageous if a single term should be constrained in a model containing many terms. If provided as a list, Cmat is internally expanded to create the full constraint matrix. See examples of constraint matrices below.

lb and ub are vectors defining the bounds of the constraints. By default they are set to 0 and Inf, meaning that the linear combinations defined by Cmat should be positive, but any bounds are possible. When some elements of lb and ub are identical, they define equality constraints. Setting lb = -Inf and ub = Inf disable the constraints.

Quadratic programming solvers

The function cirls.fit relies on a quadratic programming solver. Several solver are currently available.

Each solver has specific parameters that can be controlled through the argument qp_pars. Sensible defaults are set within cirls.control and the user typically doesn't need to provide custom parameters.

Value

A cirls object inheriting from the class glm. At the moment, two non-standard methods specific to cirls objects are available: vcov.cirls to obtain the coefficients variance-covariance matrix and confint.cirls to obtain confidence intervals. These custom methods account for the reduced degrees of freedom resulting from the constraints, see vcov.cirls and confint.cirls. Any method for glm objects can be used, including the generic coef or summary for instance.

An object of class cirls includes all components from glm objects, with the addition of:

active.cons

vector of indices of the active constraints in the fitted model.

inner.iter

number of iterations performed by the last call to the QP solver.

Cmat, lb, ub

the (expanded) constraint matrix, and lower and upper bound vectors.

References

Goldfarb, D., Idnani, A., 1983. A numerically stable dual method for solving strictly convex quadratic programs. Mathematical Programming 27, 1–33. doi:10.1007/BF02591962

Meyer, M.C., 2013. A Simple New Algorithm for Quadratic Programming with Applications in Statistics. Communications in Statistics - Simulation and Computation 42, 1126–1139. doi:10.1080/03610918.2012.659820

Stellato, B., Banjac, G., Goulart, P., Bemporad, A., Boyd, S., 2020. OSQP: an operator splitting solver for quadratic programs. Math. Prog. Comp. 12, 637–672. doi:10.1007/s12532-020-00179-2

See Also

vcov.cirls, confint.cirls for methods specific to cirls objects. cirls.control for fitting parameters specific to cirls.fit. glm for details on glm objects.

Examples

####################################################
# Simple non-negative least squares

# Simulate predictors and response with some negative coefficients
set.seed(111)
n <- 100
p <- 10
betas <- rep_len(c(1, -1), p)
x <- matrix(rnorm(n * p), nrow = n)
y <- x %*% betas + rnorm(n)

# Define constraint matrix (includes intercept)
# By default, bounds are 0 and +Inf
Cmat <- cbind(0, diag(p))

# Fit GLM by CIRLS
res1 <- glm(y ~ x, method = cirls.fit, Cmat = Cmat)
coef(res1)

# Same as passing Cmat through the control argument
res2 <- glm(y ~ x, method = cirls.fit, control = list(Cmat = Cmat))
identical(coef(res1), coef(res2))

####################################################
# Increasing coefficients

# Generate two group of variables: an isotonic one and an unconstrained one
set.seed(222)
p1 <- 5; p2 <- 3
x1 <- matrix(rnorm(100 * p1), 100, p1)
x2 <- matrix(rnorm(100 * p2), 100, p2)

# Generate coefficients: those in b1 should be increasing
b1 <- runif(p1) |> sort()
b2 <- runif(p2)

# Generate full data
y <- x1 %*% b1 + x2 %*% b2 + rnorm(100, sd = 2)

#----- Fit model

# Create constraint matrix and expand for intercept and unconstrained variables
Ciso <- diff(diag(p1))
Cmat <- cbind(0, Ciso, matrix(0, nrow(Ciso), p2))

# Fit model
resiso <- glm(y ~ x1 + x2, method = cirls.fit, Cmat = Cmat)
coef(resiso)

# Compare with unconstrained
plot(c(0, b1, b2), pch = 16)
points(coef(resiso), pch = 16, col = 3)
points(coef(glm(y ~ x1 + x2)), col = 2)

#----- More convenient specification

# Cmat can be provided as a list
resiso2 <- glm(y ~ x1 + x2, method = cirls.fit, Cmat = list(x1 = Ciso))

# Internally Cmat is expanded and we obtain the same result
identical(resiso$Cmat, resiso2$Cmat)
identical(coef(resiso), coef(resiso2))

#----- Adding bounds to the constraints
# Difference between coefficients must be above a lower bound and below 1
lb <- 1 / (p1 * 2)
ub <- 1

# Re-fit the model
resiso3 <- glm(y ~ x1 + x2, method = cirls.fit, Cmat = list(x1 = Ciso),
  lb = lb, ub = ub)

# Compare the fit
plot(c(0, b1, b2), pch = 16)
points(coef(resiso), pch = 16, col = 3)
points(coef(glm(y ~ x1 + x2)), col = 2)
points(coef(resiso3), pch = 16, col = 4)

[Package cirls version 0.3.1 Index]