restriktor: Estimating linear regression models with (in)equality...

Description Usage Arguments Details Value Author(s) References See Also Examples

View source: R/restriktor.R

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.