,
Antonio Gasparrini [aut] | Title: | Constrained Iteratively Reweighted Least Squares |
|---|---|
| Description: | Routines to fit generalized linear models with constrained coefficients, along with inference on the coefficients. Designed to be used in conjunction with the base glm() function. |
| Authors: | Pierre Masselot [aut, cre, cph]
,
Antonio Gasparrini [aut] |
| Maintainer: | Pierre Masselot <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 0.3.1.9002 |
| Built: | 2025-03-07 10:33:02 UTC |
| Source: | https://github.com/pierremasselot/cirls |
Checks a constraint matrix does not contains redundant rows
checkCmat(Cmat)checkCmat(Cmat)
Cmat |
A constraint matrix as passed to |
The user typically doesn't need to use checkCmat as it is internally called by cirls.control(). However, it might be useful to undertsand if Cmat can be reduced for inference purpose. See the note in confint.cirls().
A constraint matrix is irreducible if no row can be expressed as a positive linear combination of the other rows. When it happens, it means the constraint is actually implicitly included in other constraints in the matrix and can be dropped. Note that this a less restrictive condition than the constraint matrix having full row rank (see some examples).
The function starts by checking if some constraints are redundant and, if so, checks if they underline equality constraints. In the latter case, the constraint matrix can be reduced by expressing these constraints as a single equality constraint with identical lower and upper bounds (see cirls.fit()).
The function also checks whether there are "zero constraints" i.e. constraints with only zeros in Cmat in which case they will be labelled as redundant.
A list with two elements:
redundant |
Vector of indices of redundant constraints |
equality |
Indicates which constraints are part of an underlying equality constraint |
Meyer, M.C., 1999. An extension of the mixed primal–dual bases algorithm to the case of more constraints than dimensions. Journal of Statistical Planning and Inference 81, 13–31. doi:10.1016/S0378-3758(99)00025-7
################################################### # Example of reducible matrix # Constraints: successive coefficients should increase and be convex p <- 5 cmatic <- rbind(diff(diag(p)), diff(diag(p), diff = 2)) # Checking indicates that constraints 2 to 4 are redundant. # Intuitively, if the first two coefficients increase, # then convexity forces the rest to increase checkCmat(cmatic) # Check without contraints checkCmat(cmatic[-(2:4),]) ################################################### # Example of irreducible matrix # Constraints: coefficients form an S-shape p <- 4 cmats <- rbind( diag(p)[1,], # positive diff(diag(p))[c(1, p - 1),], # Increasing at both end diff(diag(p), diff = 2)[1:(p/2 - 1),], # First half convex -diff(diag(p), diff = 2)[(p/2):(p-2),] # second half concave ) # Note, this matrix is not of full row rank qr(t(cmats))$rank all.equal(cmats[2,] + cmats[4,] - cmats[5,], cmats[3,]) # However, it is irreducible: all constraints are necessary checkCmat(cmats) ################################################### # Example of underlying equality constraint # Contraint: Parameters sum is >= 0 and sum is <= 0 cmateq <- rbind(rep(1, 3), rep(-1, 3)) # Checking indicates that both constraints imply equality constraint (sum == 0) checkCmat(cmateq)################################################### # Example of reducible matrix # Constraints: successive coefficients should increase and be convex p <- 5 cmatic <- rbind(diff(diag(p)), diff(diag(p), diff = 2)) # Checking indicates that constraints 2 to 4 are redundant. # Intuitively, if the first two coefficients increase, # then convexity forces the rest to increase checkCmat(cmatic) # Check without contraints checkCmat(cmatic[-(2:4),]) ################################################### # Example of irreducible matrix # Constraints: coefficients form an S-shape p <- 4 cmats <- rbind( diag(p)[1,], # positive diff(diag(p))[c(1, p - 1),], # Increasing at both end diff(diag(p), diff = 2)[1:(p/2 - 1),], # First half convex -diff(diag(p), diff = 2)[(p/2):(p-2),] # second half concave ) # Note, this matrix is not of full row rank qr(t(cmats))$rank all.equal(cmats[2,] + cmats[4,] - cmats[5,], cmats[3,]) # However, it is irreducible: all constraints are necessary checkCmat(cmats) ################################################### # Example of underlying equality constraint # Contraint: Parameters sum is >= 0 and sum is <= 0 cmateq <- rbind(rep(1, 3), rep(-1, 3)) # Checking indicates that both constraints imply equality constraint (sum == 0) checkCmat(cmateq)
Internal function controlling the glm fit with linear constraints. Typically only used internally by cirls.fit, but may be used to construct a control argument.
cirls.control(epsilon = 1e-08, maxit = 25, trace = FALSE, Cmat = NULL, lb = 0L, ub = Inf, qp_solver = "quadprog", qp_pars = list())cirls.control(epsilon = 1e-08, maxit = 25, trace = FALSE, Cmat = NULL, lb = 0L, ub = Inf, qp_solver = "quadprog", qp_pars = list())
epsilon |
Positive convergence tolerance. The algorithm converges when the relative change in deviance is smaller than |
maxit |
Integer giving the maximal number of CIRLS iterations. |
trace |
Logical indicating if output should be produced for each iteration. |
Cmat |
Constraint matrix specifying the linear constraints applied to coefficients. Can also be provided as a list of matrices for specific terms. |
lb, ub
|
Lower and upper bound vectors for the linear constraints. Identical values in |
qp_solver |
The quadratic programming solver. One of |
qp_pars |
List of parameters specific to the quadratic programming solver. See respective packages help. |
The control argument of glm is by default passed to the control argument of cirls.fit, which uses its elements as arguments for cirls.control: the latter provides defaults and sanity checking. The control parameters can alternatively be passed through the ... argument of glm. See glm.control for details on general GLM fitting control, and cirls.fit for details on arguments specific to constrained GLMs.
A named list containing arguments to be used in cirls.fit.
the main function cirls.fit, and glm.control.
# 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))# 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))
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.
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)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)
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 |
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 but not both. If any parameter is passed through control, then ... will be ignored.
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.
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.
The function cirls.fit relies on a quadratic programming solver. Several solver are currently available.
"quadprog" (the default) performs a dual algorithm to solve the quadratic program. It relies on the function solve.QP.
"osqp" solves the quadratic program via the Alternating Direction Method of Multipliers (ADMM). Internally it calls the function solve_osqp.
"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. "quadprog" is set as the default being generally more reliable than the other solvers. "osqp" is faster but can be less accurate, in which case it is recommended to increase convergence tolerance at the cost of speed.
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 constraint matrix, and lower and upper bound vectors. If provided as lists, the full expanded matrix and vectors are returned. |
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
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.
#################################################### # 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)#################################################### # 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)
Computes a derivative matrix for B-splines that can then be used for shape-constraints. It is internally called by shapeConstr and should not be used directly.
dmat(d, s, knots, ord)dmat(d, s, knots, ord)
d |
Non-negative integer giving the order of derivation. Should be between 0 and |
s |
Sign of the derivative. |
knots |
Vector of ordered knots from the spline bases. |
ord |
Non-negative integer giving the order of the spline. |
Does the heavy lifting in shapeConstr to create a constraint matrix for shape-constrained B-splines. Only useful for advanced users to create constraint matrices without passing an object to one of the shapeConstr methods.
A matrix of weighted differences that can be used to constrain B-spline bases.
dmat doesn't perform any checks of the parameters so use carefully. In normal usage, checks are done by shapeConstr methods.
# A second derivative matrix for cubic B-Splines with regularly spaced knots # Can be used to enforce convexity cirls:::dmat(2, 1, 1:15, 4)# A second derivative matrix for cubic B-Splines with regularly spaced knots # Can be used to enforce convexity cirls:::dmat(2, 1, 1:15, 4)
Creates a constraint matrix to shape-constrain a set of coefficients. Mainly intended for splines but can constrain various bases or set of variables. Will typically be called from within cirls.fit but can be used to generate constraint matrices.
shapeConstr(x, shape, ...) ## Default S3 method: shapeConstr(x, shape, intercept = FALSE, ...)shapeConstr(x, shape, ...) ## Default S3 method: shapeConstr(x, shape, intercept = FALSE, ...)
x |
An object representing a design matrix of predictor variables, typically basis functions. See details for supported objects. |
shape |
A character vector indicating one or several shape-constraints. See details for supported shapes. |
... |
Additional parameters passed to or from other methods. |
intercept |
For the default method, a logical value indicating if the design matrix includes an intercept. In most cases will be automatically extracted from |
The recommended usage is to directly specify the shape constraint through the shape argument in the call to glm with cirls.fit. This method is then called internally to create the constraint matrix. However, shapeConstr can nonetheless be called directly to manually build or inspect the constraint matrix for a given shape and design matrix.
The parameters necessary to build the constraint matrix (e.g. knots and ord for splines) are typically extracted from the x object. This is also true for the intercept for most of the object, except for the default method for which it can be useful to explicitly provide it. In a typical usage in which shapeConstr would only be called within cirls.fit, intercept is automatically determined from the glm formula.
The shape argument allows to define a specific shape for the association between the expanded term in x and the response of the regression model. This shape can describe the relation between coefficients for the default method, or the shape of the smooth term for spline bases. At the moment, six different shapes are supported, with up to three allowed simultaneously (one from each category):
"pos" or "neg": Positive/Negative. Applies to the full association.
"inc" or "dec": Monotonically Increasing/Decreasing.
"cvx" or "ccv": Convex/Concave.
In addition to the default method, shapeConstr currently supports several objects, creating an appropriate shape-constraint matrix depending on the object. The full list can be obtained by methods(shapeConstr).
onebasis: General method for basis functions generated in the package.
ps: Penalised splines (P-Splines).
A constraint matrix to be passed to Cmat in cirls.fit.
Zhou, S. & Wolfe, D. A., 2000, On derivative estimation in spline regression. Statistica Sinica 10, 93–108.
cirls.fit() which typically calls shapeConstr internally.
# example code# example code
cirls object.Simulates coefficients for a fitted cirls object. confint and vcov compute confidence intervals and the Variance-Covariance matrix for coefficients from a fitted cirls object. These methods supersede the default methods for cirls objects.
simulCoef(object, nsim = 1, seed = NULL, complete = TRUE) ## S3 method for class 'cirls' confint(object, parm, level = 0.95, nsim = 1000, complete = TRUE, ...) ## S3 method for class 'cirls' vcov(object, complete = TRUE, nsim = 1000, constrained = TRUE, ...)simulCoef(object, nsim = 1, seed = NULL, complete = TRUE) ## S3 method for class 'cirls' confint(object, parm, level = 0.95, nsim = 1000, complete = TRUE, ...) ## S3 method for class 'cirls' vcov(object, complete = TRUE, nsim = 1000, constrained = TRUE, ...)
object |
A fitted |
nsim |
The number of simulations to perform. |
seed |
Either NULL or an integer that will be used in a call to |
complete |
If FALSE, doesn't return inference for undetermined coefficients in case of an over-determined model. |
parm |
A specification of which parameters to compute the confidence intervals for. Either a vector of numbers or a vector of names. If missing, all parameters are considered. |
level |
The confidence level required. |
... |
Further arguments passed to or from other methods. For |
constrained |
If set to |
confint and vcov are custom methods for cirls objects to supersede the default methods used for glm objects. Internally, they both call simulCoef to generate coefficient vectors from a Truncated Multivariate Normal Distribution using the TruncatedNormal::rtmvnorm() function. This distribution accounts for truncation by constraints, ensuring all coefficients are feasible with respect to the constraint matrix. simulCoef typically doesn't need to be used directly for confidence intervals and variance-covariance matrices, but it can be used to check other summaries of the coefficients distribution.
These methods only work when Cmat is of full row rank, i.e. if there are less constraints than variables in object.
For simulCoef, a matrix with nsim rows containing simulated coefficients.
For confint, a two-column matrix with columns giving lower and upper confidence limits for each parameter.
For vcov, a matrix of the estimated covariances between the parameter estimates of the model.
By default, the Variance-Covariance matrix generated by vcov is different than the one returned by summary(obj)$cov.scaled. The former accounts for the reduction in degrees of freedom resulting from the constraints, while the latter is the unconstrained GLM Variance-Covariance. Note that the unconstrained one can be obtained from vcov by setting constrained = FALSE.
Geweke, J.F., 1996. Bayesian Inference for Linear Models Subject to Linear Inequality Constraints, in: Lee, J.C., Johnson, W.O., Zellner, A. (Eds.), Modelling and Prediction Honoring Seymour Geisser. Springer, New York, NY, pp. 248–263. doi:10.1007/978-1-4612-2414-3_15
Botev, Z.I., 2017, The normal law under linear restrictions: simulation and estimation via minimax tilting, Journal of the Royal Statistical Society, Series B, 79 (1), pp. 1–24.
rtmvnorm for the underlying routine to simulate from a TMVN. checkCmat() to check if the contraint matrix can be reduced.
#################################################### # Isotonic regression #----- Perform isotonic regression # Generate data set.seed(222) p1 <- 5; p2 <- 3 x1 <- matrix(rnorm(100 * p1), 100, p1) x2 <- matrix(rnorm(100 * p2), 100, p2) b1 <- runif(p1) |> sort() b2 <- runif(p2) y <- x1 %*% b1 + x2 %*% b2 + rnorm(100, sd = 2) # Fit model Ciso <- diff(diag(p1)) resiso <- glm(y ~ x1 + x2, method = cirls.fit, Cmat = list(x1 = Ciso)) #----- Extract uncertainty # Extract variance covariance vcov(resiso) # Extract confidence intervals confint(resiso) # We can extract the usual unconstrained matrix vcov(resiso, constrained = FALSE) all.equal(vcov(resiso, constrained = FALSE), summary(resiso)$cov.scaled) # Simulate from the distribution of coefficients sims <- simulCoef(resiso, nsim = 10) # Check that all simulated coefficient vectors are feasible apply(resiso$Cmat %*% t(sims) >= resiso$lb, 2, all)#################################################### # Isotonic regression #----- Perform isotonic regression # Generate data set.seed(222) p1 <- 5; p2 <- 3 x1 <- matrix(rnorm(100 * p1), 100, p1) x2 <- matrix(rnorm(100 * p2), 100, p2) b1 <- runif(p1) |> sort() b2 <- runif(p2) y <- x1 %*% b1 + x2 %*% b2 + rnorm(100, sd = 2) # Fit model Ciso <- diff(diag(p1)) resiso <- glm(y ~ x1 + x2, method = cirls.fit, Cmat = list(x1 = Ciso)) #----- Extract uncertainty # Extract variance covariance vcov(resiso) # Extract confidence intervals confint(resiso) # We can extract the usual unconstrained matrix vcov(resiso, constrained = FALSE) all.equal(vcov(resiso, constrained = FALSE), summary(resiso)$cov.scaled) # Simulate from the distribution of coefficients sims <- simulCoef(resiso, nsim = 10) # Check that all simulated coefficient vectors are feasible apply(resiso$Cmat %*% t(sims) >= resiso$lb, 2, all)