orgls: Fitting generalized least squares regression models with...

In goric: Generalized Order-Restricted Information Criterion

Description

orgls is used to fit generalised least square models analogously to the function gls in package nlme but with order restrictions on the parameters.

Usage

orgls(formula, data, constr, rhs, nec, weights = NULL, correlation = NULL,
  control = orlmcontrol())

## S3 method for class 'formula'
orgls(formula, data, constr, rhs, nec, weights = NULL,
  correlation = NULL, 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 orgls 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.
weights	a varClasses object; more details are provided on the help pages in R package nlme
correlation	a corClasses object; more details are provided on the help pages in R package nlme
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).

Value

an object of class orgls

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

# generating example data
library(mvtnorm)
# group means
m <- c(0,5,5,7)
# compound symmetry structure of residuals
# (10 individuals per group, rho=0.7) 
cormat <- kronecker(diag(length(m)*10), matrix(0.7, nrow=length(m), ncol=length(m)))
diag(cormat) <- 1
# different variances per group
sds <- rep(c(1,2,0.5,1), times=10*length(m))
sigma <- crossprod(diag(sds), crossprod(cormat, diag(sds)))
response <- as.vector(rmvnorm(1, rep(m, times=10*length(m)), sigma=sigma))
dat <- data.frame(response,
                  grp=rep(LETTERS[1:length(m)], times=10*length(m)), 
                  ID=as.factor(rep(1:(10*length(m)), each=length(m))))
                  
## set of gls models:
# unconstrained model
m1 <- orgls(response ~ grp-1, data = dat,
            constr=rbind(c(0,0,0,0)), rhs=0, nec=0,
            weights=varIdent(form=~1|grp),
            correlation=corCompSymm(form=~1|ID))

# simple order
m2 <- orgls(response ~ grp-1, data = dat,
            constr=rbind(c(-1,1,0,0),c(0,-1,1,0),c(0,0,-1,1)), rhs=c(0,0,0), nec=0,
            weights=varIdent(form=~1|grp),
            correlation=corCompSymm(form=~1|ID))

# equality constraints
m3 <- orgls(response ~ grp-1, data = dat,
            constr=rbind(c(-1,1,0,0),c(0,-1,1,0),c(0,0,-1,1)), rhs=c(0,0,0), nec=3,
            weights=varIdent(form=~1|grp),
            correlation=corCompSymm(form=~1|ID))

goric documentation built on April 19, 2021, 5:06 p.m.