Fitting multivariate regression models with order restrictions

Share:

Description

It determines the order-restricted maximum likelihood estimates and the corresponding log likelihood for the hypothesis of interest. Additionally it gives the (unconstrained) maximum likelihood estimates and the active contraints.

Usage

1
2
3
orlm(formula, data, constr, rhs, nec, control = orlmcontrol())
## S3 method for class 'formula'
orlm(formula, data, constr, rhs, nec, control = orlmcontrol())

Arguments

formula

an object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted.

data

an optional data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which lm is called.

constr

matrix with constraints; with rows as constraint definition, columns should be in line with the parameters of the model

rhs

vector of right hand side elements; Constr \; θ ≥q rhs; number should equal the number of rows of the constr matrix

nec

number of equality constraints; a numeric value treating the first nec constr rows as equality constraints, or a logical vector with TRUE for equality- and FALSE for inequality constraints.

control

a list of control arguments; see orlmcontrol for details.

Details

The contraints in the hypothesis of interest are defined by Constr, rhs, and nec. The first nec constraints are the equality contraints: Constr[1:nec, 1:tk] θ = rhs[1:nec]; and the remaing ones are the inequality contraints: Constr[nec+1:c_m, 1:tk] θ ≥q rhs[nec+1:c_m].

Two requirements should be met:

  1. The first nec constraints must be the equality contraints (i.e., Constr[1:nec, 1:tk] θ = rhs[1:nec]) and the remaining ones the inequality contraints (i.e., Constr[nec+1:c_m, 1:tk] θ ≥q rhs[nec+1:c_m]).

  2. When rhs is not zero, Constr should be of full rank (after discarding redundant restrictions). More information can be obtained from Kuiper, Hoijtink, and Silvapulle (2011) and Kuiper, Hoijtink, and Silvapulle (unpublished).

Value

an object of class orlm

Author(s)

Daniel Gerhard and Rebecca M. Kuiper

References

Kuiper R.M., Hoijtink H., Silvapulle M.J. (2011). An Akaike-type Information Criterion for Model Selection Under Inequality Constraints. Biometrika, 98, 495–501.

Kuiper R.M., Hoijtink H., Silvapulle M.J. (2012). Generalization of the Order-Restricted Information Criterion for Multivariate Normal Linear Models. Journal of Statistical Planning and Inference, 142, 2454-2463. doi:10.1016/j.jspi.2012.03.007.

Kuiper R.M. and Hoijtink H. (submitted). A Fortran 90 Program for the Generalization of the Order-Restricted Information Criterion. Journal of Statictical Software.

See Also

solve.QP, 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
########################
## Artificial example ##
########################
n <- 10
m <- c(1,2,1,5)
nm <- length(m)
dat <- data.frame(grp=as.factor(rep(1:nm, each=n)),
                  y=rnorm(n*nm, rep(m, each=n), 1))

# unrestricted linear model
cm1 <- matrix(0, nrow=1, ncol=4)
fm1 <- orlm(y ~ grp-1, data=dat, constr=cm1, rhs=0, nec=0)

# order restriction (increasing means)
cm2 <- rbind(c(-1,1,0,0),
             c(0,-1,1,0),
             c(0,0,-1,1))
fm2 <- orlm(y ~ grp-1, data=dat, constr=cm2,
            rhs=rep(0,nrow(cm2)), nec=0)

# order restriction (increasing at least by delta=1)
fm3 <- orlm(y ~ grp-1, data=dat, constr=cm2,
            rhs=rep(1,nrow(cm2)), nec=0)

# larger than average of the neighboring first 2 parameters
cm4 <- rbind(c(-0.5,-0.5,1,0),
             c(0,-0.5,-0.5,1))
fm4 <- orlm(y ~ grp-1, data=dat, constr=cm4,
            rhs=rep(0,nrow(cm4)), nec=0)

# equality constraints (all parameters equal)
fm5 <- orlm(y ~ grp-1, data=dat, constr=cm2,
            rhs=rep(0,nrow(cm2)), nec=nrow(cm2))
# alternatively
fm5 <- orlm(y ~ grp-1, data=dat, constr=cm2,
            rhs=rep(0,nrow(cm2)), nec=c(TRUE,TRUE,TRUE))

# constraining the 1st and the 4th parameter
# to their true values, and the 2nd and 3rd between them
cm6 <- rbind(c( 1,0,0,0),
             c(-1,1,0,0),
             c(0,-1,0,1),
             c(-1,0,1,0),
             c(0,0,-1,1),
             c(0,0, 0,1))
fm6 <- orlm(y ~ grp-1, data=dat, constr=cm6,
            rhs=c(1,rep(0,4),5), nec=c(TRUE,rep(FALSE,4),TRUE))


###############################################################
## Example from Kuiper, R.M. and Hoijtink, H. (Unpublished). ##
## A Fortran 90 program for the generalization of the        ##
## order restricted information criterion.                   ##
###############################################################

# constraint definition
cmat <- cbind(diag(3), 0) + cbind(0, -diag(3))
constr <- kronecker(diag(3), cmat)

# no effect model
fm0 <- orlm(cbind(SDH, SGOT, SGPT) ~ dose-1, data=vinylidene,
            constr=constr, rhs=rep(0, nrow(constr)), nec=nrow(constr))
fm0

# order constrained model (increasing serum levels with increasing doses)
fm1 <- orlm(cbind(SDH, SGOT, SGPT) ~ dose-1, data=vinylidene,
            constr=constr, rhs=rep(0, nrow(constr)), nec=0)
summary(fm1)

# unconstrained model
fmunc <- orlm(cbind(SDH, SGOT, SGPT) ~ dose-1, data=vinylidene,
              constr=matrix(0, nrow=1, ncol=12), rhs=0, nec=0)
fmunc