Description Usage Arguments Details Value Author(s) References See Also Examples
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.
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, ...)

object 
a fitted linear model object of class "lm", "mlm",
"rlm" or "glm". For class "rlm" only the loss function 
constraints 
there are two ways to constrain parameters.
First, the constraint syntax consists of one or more textbased
descriptions, where the syntax can be specified as a literal
string enclosed by single quotes. Only the names of Second, the constraint syntax consists of a matrix R (or a vector in case of one constraint) and defines the lefthand 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 " 
B 
integer; number of bootstrap draws for 
rhs 
vector on the righthand 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 
mix.weights 
if 
mix.bootstrap 
integer; number of bootstrap draws for

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 
seed 
seed value. 
control 
a list of control arguments:

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. 
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 textbased descriptions in the constraints syntax:
Equality constraints: The "==
" operator can be
used to define equality constraints (e.g., x1 == 1
or
x1 == x2
).
Inequality constraints: The "<
" or ">
"
operator can be used to define inequality constraints
(e.g., x1 > 1
or x1 < x2
).
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 socalled Delta method.
Variable names of interaction effects in objects of class lm,
rlm and glm contain a semicolon (:) between the variables. To impose
constraints on parameters of interaction effects, the semicolon
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.
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 Rsquared. 
R2.reduced 
restricted Rsquared. 
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 loglikelihood. 
Sigma 
variancecovariance matrix of unrestricted model. 
constraints 
matrix with restrictions. 
rhs 
vector of righthand side elements. 
neq 
number of equality restrictions. 
wt.bar 
chibarsquare 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. 
Leonard Vanbrabant and Yves Rosseel
Schoenberg, R. (1997). Constrained Maximum Likelihood. Computational Economics, 10, 251–266.
Shapiro, A. (1988). Towards a unified theory of inequalityconstrained testing in multivariate analysis. International Statistical Review 56, 49–62.
Silvapulle, M.J. and Sen, P.K. (2005). Constrained Statistical Inference. Wiley, New York
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
# HuberWhite (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 modelbased 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)

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.