This is a modification of the lm function, fitting (multivariate) linear models with order constraints on the model coefficients.
orlm(formula, data, constr, rhs, nec, control = orlmcontrol()) # S3 method for formula orlm(formula, data, constr, rhs, nec, control = orlmcontrol())
| 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 \; \theta \geq 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 |
| control | a list of control arguments; see |
an object of class orlm
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] \theta = rhs[1:nec]\); and the remaing ones are the inequality contraints: \(Constr[nec+1:c_m, 1:tk] \theta \geq rhs[nec+1:c_m]\). Two requirements should be met:
The first \(nec\) constraints must be the equality contraints (i.e., \(Constr[1:nec, 1:tk] \theta = rhs[1:nec]\)) and the remaining ones the inequality contraints (i.e., \(Constr[nec+1:c_m, 1:tk] \theta \geq rhs[nec+1:c_m]\)).
When \(rhs\) is not zero, \(Constr\) should be of full rank (after discarding redundant restrictions).
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.
solve.QP, goric
######################## ## 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)))#> #> Call: #> orlm.formula(formula = cbind(SDH, SGOT, SGPT) ~ dose - 1, data = vinylidene, #> constr = constr, rhs = rep(0, nrow(constr)), nec = nrow(constr)) #> #> Coefficients: #> SDH SGOT SGPT #> dosed1 24.12 105.38 59.70 #> dosed2 24.12 105.38 59.70 #> dosed3 24.12 105.38 59.70 #> dosed4 24.12 105.37 59.70 #># 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)#> #> Call: #> orlm.formula(formula = cbind(SDH, SGOT, SGPT) ~ dose - 1, data = vinylidene, #> constr = constr, rhs = rep(0, nrow(constr)), nec = 0) #> #> logLik: -396.8 #> #> Coefficients: #> SDH SGOT SGPT #> dosed1 24.12 105.37 63.00 #> dosed2 24.12 105.37 63.00 #> dosed3 24.12 105.37 60.64 #> dosed4 24.12 105.37 52.16 #> #> Unconstrained solution: #> SDH SGOT SGPT #> dosed1 22.7 99.3 61.9 #> dosed2 22.8 108.4 63.8 #> dosed3 23.7 100.9 60.2 #> dosed4 27.3 112.9 52.9 #> #> Active constraints: 6 4 3 5 2 1 7# unconstrained model (fmunc <- orlm(cbind(SDH, SGOT, SGPT) ~ dose-1, data=vinylidene, constr=matrix(0, nrow=1, ncol=12), rhs=0, nec=0))#> #> Call: #> orlm.formula(formula = cbind(SDH, SGOT, SGPT) ~ dose - 1, data = vinylidene, #> constr = matrix(0, nrow = 1, ncol = 12), rhs = 0, nec = 0) #> #> Coefficients: #> SDH SGOT SGPT #> dosed1 22.7 99.3 61.9 #> dosed2 22.8 108.4 63.8 #> dosed3 23.7 100.9 60.2 #> dosed4 27.3 112.9 52.9 #>