# restriktor: Estimating linear regression models with (in)equality... In restriktor: Restricted Statistical Estimation and Inference for Linear Models

## Description

Function `restriktor` estimates the parameters of an univariate and a multivariate linear model (`lm`), a robust estimation of the linear model (`rlm`) and a generalized linear model (`glm`) subject to linear equality and linear inequality restrictions. It is a convenience function. The real work horses are the `conLM`, `conMLM`, `conRLM` and the `conGLM` functions.

## Usage

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29``` ```restriktor(object, constraints = NULL, ...) ## S3 method for class 'lm' conLM(object, constraints = NULL, se = "standard", B = 999, rhs = NULL, neq = 0L, mix.weights = "pmvnorm", mix.bootstrap = 99999L, parallel = "no", ncpus = 1L, cl = NULL, seed = NULL, control = list(), verbose = FALSE, debug = FALSE, ...) ## S3 method for class 'rlm' conRLM(object, constraints = NULL, se = "standard", B = 999, rhs = NULL, neq = 0L, mix.weights = "pmvnorm", mix.bootstrap = 99999L, parallel = "no", ncpus = 1L, cl = NULL, seed = NULL, control = list(), verbose = FALSE, debug = FALSE, ...) ## S3 method for class 'glm' conGLM(object, constraints = NULL, se = "standard", B = 999, rhs = NULL, neq = 0L, mix.weights = "pmvnorm", mix.bootstrap = 99999L, parallel = "no", ncpus = 1L, cl = NULL, seed = NULL, control = list(), verbose = FALSE, debug = FALSE, ...) ## S3 method for class 'mlm' conMLM(object, constraints = NULL, se = "none", B = 999, rhs = NULL, neq = 0L, mix.weights = "pmvnorm", mix.bootstrap = 99999L, parallel = "no", ncpus = 1L, cl = NULL, seed = NULL, control = list(), verbose = FALSE, debug = FALSE, ...) ```

## Arguments

 `object` a fitted linear model object of class "lm", "mlm", "rlm" or "glm". For class "rlm" only the loss function `bisquare` is supported for now, otherwise the function gives an error. `constraints` there are two ways to constrain parameters. First, the constraint syntax consists of one or more text-based descriptions, where the syntax can be specified as a literal string enclosed by single quotes. Only the names of `coef(model)` can be used as names. See details for more information. Note that objects of class "mlm" do not (yet) support this method. Second, the constraint syntax consists of a matrix R (or a vector in case of one constraint) and defines the left-hand side of the constraint Rθ ≥ rhs, where each row represents one constraint. The number of columns needs to correspond to the number of parameters estimated (θ) by model. The rows should be linear independent, otherwise the function gives an error. For more information about constructing the matrix R and rhs see details. `se` if "`standard`" (default), conventional standard errors are computed based on inverting the observed augmented information matrix. If "const", homoskedastic standard errors are computed. If "`HC0`" or just "`HC`", heteroskedastic robust standard errors are computed (a.k.a Huber White). The options "`HC1`", "`HC2`", "`HC3`", "`HC4`", "`HC4m`", and "`HC5`" are refinements of "`HC0`". For more details about heteroskedastic robust standard errors see the sandwich package. If "`boot.standard`", bootstrapped standard errors are computed using standard bootstrapping. If "`boot.model.based`" or "`boot.residual`", bootstrapped standard errors are computed using model-based bootstrapping. If "`none`", no standard errors are computed. Note that for objects of class "mlm" no standard errors are available (yet). `B` integer; number of bootstrap draws for `se`. The default value is set to 999. Parallel support is available. `rhs` vector on the right-hand side of the constraints; Rθ ≥ rhs. The length of this vector equals the number of rows of the constraints matrix R and consists of zeros by default. Note: only used if constraints input is a matrix or vector. `neq` integer (default = 0) treating the number of constraints rows as equality constraints instead of inequality constraints. For example, if `neq = 2`, this means that the first two rows of the constraints matrix R are treated as equality constraints. Note: only used if constraints input is a matrix or vector. `mix.weights` if `"pmvnorm"` (default), the chi-bar-square weights are computed based on the multivariate normal distribution function with additional Monte Carlo steps. If `"boot"`, the chi-bar-square weights are computed using parametric bootstrapping. If `"none"`, no chi-bar-square weights are computed. The weights are necessary in the `restriktor.summary` function for computing the GORIC. Moreover, the weights are re-used in the `iht` function for computing the p-value for the test-statistic, unless the p-value is computed directly via bootstrapping. `mix.bootstrap` integer; number of bootstrap draws for `mix.weights = "boot"`. The default value is set to 99999. Parallel support is available. `parallel` the type of parallel operation to be used (if any). If missing, the default is set "no". `ncpus` integer: number of processes to be used in parallel operation: typically one would chose this to the number of available CPUs. `cl` an optional parallel or snow cluster for use if parallel = "snow". If not supplied, a cluster on the local machine is created for the duration of the `restriktor` call. `seed` seed value. `control` a list of control arguments: `absval` tolerance criterion for convergence (default = sqrt(.Machine\$double.eps)). `maxit` the maximum number of iterations for the optimizer (default = 10000). `tol` numerical tolerance value. Estimates smaller than `tol` are set to 0. `verbose` logical; if TRUE, information is shown at each bootstrap draw. `debug` if TRUE, debugging information about the constraints is printed out. `...` no additional arguments for now.

## Details

The constraint syntax can be specified in two ways. First as a literal string enclosed by single quotes as shown below:

 ```1 2 3 4 5 6 7 8``` ```myConstraints <- ' # 1. inequality constraints x1 > 0 x1 < x2 ! 2. equality constraints x3 == x4; x4 == x5 ' ```

The variable names x1 to x5 refer to the corresponding regression coefficient. Thus, constraints are impose on regression coefficients and not on the data.

Second, the above constraints syntax can also be written in matrix/vector notation as:

(The first column refers to the intercept, the remaining five columns refer to the regression coefficients x1 to x5.)

 ``` 1 2 3 4 5 6 7 8 9 10 11``` ```myConstraints <- rbind(c(0, 0, 0,-1, 1, 0), #equality constraint x3 == x4 c(0, 0, 0, 0,-1, 1), #equality constraint x4 == x5 c(0, 1, 0, 0, 0, 0), #inequality constraint x1 > rhs c(0,-1, 1, 0, 0, 0)) #inequality constraint x1 < x2 # the length of rhs is equal to the number of myConstraints rows. myRhs <- c(0,0,0,0) # the first two rows should be considered as equality constraints myNeq <- 2 ```

Blank lines and comments can be used in between the constraints, and constraints can be split over multiple lines. Both the hashtag (#) and the exclamation (!) characters can be used to start a comment. Multiple constraints can be placed on a single line if they are separated by a semicolon (;).

There can be three types of text-based descriptions in the constraints syntax:

1. Equality constraints: The "`==`" operator can be used to define equality constraints (e.g., `x1 == 1` or `x1 == x2`).

2. Inequality constraints: The "`<`" or "`>`" operator can be used to define inequality constraints (e.g., `x1 > 1` or `x1 < x2`).

3. Newly defined parameters: The "`:=`" operator can be used to define new parameters, which take on values that are an arbitrary function of the original model parameters. The function must be specified in terms of the parameter names in `coef(model)` (e.g., `new := x1 + 2*x2`). By default, the standard errors for these defined parameters are computed by using the so-called Delta method.

Variable names of interaction effects in objects of class lm, rlm and glm contain a semi-colon (:) between the variables. To impose constraints on parameters of interaction effects, the semi-colon must be replaced by a dot (.) (e.g., `x3:x4` becomes `x3.x4`). In addition, the intercept variable names is shown as "`(Intercept)`". To impose restrictions on the intercept both parentheses must be replaced by a dot "`.Intercept.`" (e.g.,`.Intercept. > 10`). Note: in most practical situations we do not impose restrictions on the intercept because we do not have prior knowledge about the intercept. Moreover, the sign of the intercept can be changed arbitrarily by shifting the response variable y.

Each element can be modified using arithmetic operators. For example, if `x2` is expected to be twice as large as `x1`, then "`2*x2 == x1`".

If `constraints = NULL`, the unrestricted model is fitted.

## Value

An object of class restriktor, for which a print and a summary method are available. More specifically, it is a list with the following items:

 `CON` a list with useful information about the restrictions. `call` the matched call. `timing` how much time several tasks take. `parTable` a parameter table with information about the observed variables in the model and the imposed restrictions. `b.unrestr` unrestricted regression coefficients. `b.restr` restricted regression coefficients. `residuals` restricted residuals. `wresid` a working residual, weighted for "inv.var" weights only (rlm only) `fitted` restricted fitted mean values. `weights` (only for weighted fits) the specified weights. `wgt` the weights used in the IWLS process (rlm only). `scale` the robust scale estimate used (rlm only). `stddev` a scale estimate used for the standard errors. `R2.org` unrestricted R-squared. `R2.reduced` restricted R-squared. `df.residual` the residual degrees of freedom `s2.unrestr` mean squared error of unrestricted model. `s2.restr` mean squared error of restricted model. `loglik` restricted log-likelihood. `Sigma` variance-covariance matrix of unrestricted model. `constraints` matrix with restrictions. `rhs` vector of right-hand side elements. `neq` number of equality restrictions. `wt.bar` chi-bar-square mixing weights or a.k.a. level probabilities. `iact` active restrictions. `converged` did the IWLS converge (rlm only)? `iter` number of iteration needed for convergence (rlm only). `bootout` object of class boot. Only available if bootstrapped standard errors are requested, else bootout = NULL. `control` list with control options. `model.org` original model. `se` as input. This information is needed in the summary function. `information` observed information matrix with the inverted information matrix and the augmented information matrix as attributes.

## Author(s)

Leonard Vanbrabant and Yves Rosseel

## References

Schoenberg, R. (1997). Constrained Maximum Likelihood. Computational Economics, 10, 251–266.

Shapiro, A. (1988). Towards a unified theory of inequality-constrained testing in multivariate analysis. International Statistical Review 56, 49–62.

Silvapulle, M.J. and Sen, P.K. (2005). Constrained Statistical Inference. Wiley, New York

## See Also

`iht`, `goric`

## Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138``` ```## lm ## unrestricted linear model for ages (in months) at which an ## infant starts to walk alone. # prepare data DATA1 <- subset(ZelazoKolb1972, Group != "Control") # fit unrestricted linear model fit1.lm <- lm(Age ~ -1 + Group, data = DATA1) # the variable names can be used to impose restrictions on # the corresponding regression parameters. coef(fit1.lm) # restricted linear model with restrictions that the walking # exercises would not have a negative effect of increasing the # mean age at which a child starts to walk. fit1.con <- restriktor(fit1.lm, constraints = ' GroupActive < GroupPassive; GroupPassive < GroupNo ') summary(fit1.con) ## Not run: # Or in matrix notation. myConstraints1 <- rbind(c(-1, 1, 0), c( 0,-1, 1)) myRhs1 <- rep(0L, nrow(R1)) myNeq1 <- 0 fit1.con <- restriktor(fit1.lm, constraints = myConstraints1, rhs = myRhs1, neq = myNeq1) summary(fit1.con) ## End(Not run) ######################### ## Artificial examples ## ######################### library(MASS) ## mlm # generate data n <- 30 mu <- rep(0, 4) Sigma <- matrix(5,4,4) diag(Sigma) <- c(10,10,10,10) # 4 Y's. Y <- mvrnorm(n, mu, Sigma) # fit unrestricted multivariate linear model fit.mlm <- lm(Y ~ 1) # constraints myConstraints2 <- diag(0,4) diag(myConstraints2) <- 1 # fit restricted multivariate linear model fit2.con <- restriktor(fit.mlm, constraints = myConstraints2) summary(fit2.con) ## rlm # generate data n <- 10 means <- c(1,2,1,3) nm <- length(means) group <- as.factor(rep(1:nm, each = n)) y <- rnorm(n * nm, rep(means, each = n)) DATA2 <- data.frame(y, group) # fit unrestricted robust linear model fit3.rlm <- rlm(y ~ -1 + group, data = DATA2, method = "MM") coef(fit3.rlm) ## increasing means myConstraints3 <- ' group1 < group2 group2 < group3 group3 < group4 ' # fit restricted robust linear model and compute # Huber-White (robust) standard errors. fit3.con <- restriktor(fit3.rlm, constraints = myConstraints2, se = "HC0") summary(fit3.con) ## Not run: ## increasing means in matrix notation. myConstraints3 <- rbind(c(-1, 1, 0, 0), c( 0,-1, 1, 0), c( 0, 0,-1, 1)) myRhs3 <- rep(0L, nrow(myConstraints3)) myNeq2 <- 0 fit3.con <- restriktor(fit3.rlm, constraints = myConstraints3, rhs = myRhs2, neq = myNeq2, se = "HC0") summary(fit3.con) ## End(Not run) ## equality restrictions only. myConstraints4 <- ' group1 == group2 group2 == group3 group3 == group4 ' fit4.con <- restriktor(fit3.rlm, constraints = myConstraints4) summary(fit4.con) ## combination of equality and inequality restrictions. myConstraints5 <- ' group1 == group2 group3 < group4 ' # fit restricted model and compute model-based bootstrapped # standard errors. We only generate 9 bootstrap samples in this # example; in practice you may wish to use a much higher number. fit5.con <- restriktor(fit3.rlm, constraints = myConstraints4, se = "boot.model.based", B = 9) # an error is probably thrown, due to a too low number of bootstrap draws. summary(fit5.con) # restriktor can also be used to define effects using the := operator # and impose restrictions on them. For example, compute the average # effect (AVE) and impose the restriction AVE > 0. # generate data n <- 30 b0 <- 10; b1 = 0.5; b2 = 1; b3 = 1.5 X <- c(rep(c(0), n/2), rep(c(1), n/2)) set.seed(90) Z <- rnorm(n, 16, 5) y <- b0 + b1*X + b2*Z + b3*X*Z + rnorm(n, 0, sd = 10) DATA3 = data.frame(cbind(y, X, Z)) # fit linear model with interaction fit6.lm <- lm(y ~ X*Z, data = DATA3) fit6.con <- restriktor(fit6.lm, constraints = ' AVE := X + 16.86137*X.Z; AVE > 0 ') summary(fit6.con) ```

restriktor documentation built on Feb. 25, 2020, 5:08 p.m.