The goric function calculates the order-restricted log likelihood, the penalty of the generalised order restricted information criterion (GORIC), the GORIC values, differences to the minimum GORIC value, and the GORIC weights for a set of hypotheses, where the penalty is based on \(iter\) iterations. The hypothesis with the lowest GORIC value is the preferred one. The GORIC weights reflect the support of each hypothesis in the set. To compare two hypotheses (and not one to the whole set), one should examine the ratio of the two corresponding GORIC weights. To safequard for weak hypotheses (i.e., hypotheses not supported by the data), one should include a model with no constraints (the so-called unconstrained model).
goric(object, ..., iter = 1e+05, type = "GORIC", dispersion = 1, mc.cores = 1) # S3 method for orlm goric(object, ..., iter = 1e+05, type = "GORIC", mc.cores = 1) # S3 method for orgls goric(object, ..., iter = 1e+05, type = "GORIC", mc.cores = 1) # S3 method for list goric(object, ..., iter = 1e+05, type = "GORIC", dispersion = 1, mc.cores = 1) # S3 method for orglm goric(object, ..., iter = 1e+05, type = "GORIC", dispersion = 1, mc.cores = 1)
| object | an object of class orlm, orgls, orglm, or a list of these objects |
|---|---|
| ... | further objects of class orlm, orgls, or orglm |
| iter | number of iterations to calculate GORIC penalty terms |
| type | if |
| dispersion | dispersion parameter to scale GORIC analogously to QAIC in generalized linear models |
| mc.cores | number of cores using a socket cluster implemented in package |
a data.frame with the information criteria or a single penalty term
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.
## 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) constr#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] #> [1,] 1 -1 0 0 0 0 0 0 0 0 0 0 #> [2,] 0 1 -1 0 0 0 0 0 0 0 0 0 #> [3,] 0 0 1 -1 0 0 0 0 0 0 0 0 #> [4,] 0 0 0 0 1 -1 0 0 0 0 0 0 #> [5,] 0 0 0 0 0 1 -1 0 0 0 0 0 #> [6,] 0 0 0 0 0 0 1 -1 0 0 0 0 #> [7,] 0 0 0 0 0 0 0 0 1 -1 0 0 #> [8,] 0 0 0 0 0 0 0 0 0 1 -1 0 #> [9,] 0 0 0 0 0 0 0 0 0 0 1 -1# 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 #># calculate GORIC # (only small number of iterations to decrease computation time, default: iter=100000) goric(fm0, fm1, fmunc, iter=1000)#> loglik penalty goric goric_weights #> fm0 -406.5447 4.000 821.0895 0.000 #> fm1 -396.8473 7.514 808.7225 0.072 #> fmunc -388.8036 13.000 803.6072 0.928