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 |
|
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 |
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.
-
"osqp"
(the default) solves the quadratic program via the Alternating Direction Method of Multipliers (ADMM). Internally it calls the function solve_osqp. -
"quadprog"
performs a dual algorithm to solve the quadratic program. It relies on the function solve.QP. -
"coneproj"
solves the quadratic program by a cone projection method. It relies on the function qprog.
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)