restriktor (version 0.4-501)

goric: Generalized Order-Restricted Information Criterion (Approximation) Weights

Description

The goric function computes GORIC(A) weights, which are comparable to the Akaike weights.

Usage

goric(object, ...)
# S3 method for default
goric(object, ..., constraints = NULL,
      comparison = c("unconstrained", "complement", "none"), 
      VCOV = NULL, sample.nobs = NULL, type = "goric", 
      auto.bound = FALSE, debug = FALSE)
# S3 method for restriktor
goric(object, ..., constraints = NULL,
      comparison = "unconstrained", type = "goric", 
      auto.bound = NULL, debug = FALSE)
# S3 method for lm
goric(object, ..., constraints = NULL,
      comparison = "unconstrained", type = "goric",
      missing = "none", auxiliary = c(), emControl = list(),
      auto.bound = NULL, debug = FALSE)
# S3 method for numeric
goric(object, ..., constraints = NULL,
      VCOV = NULL, comparison = "unconstrained",
      type = "gorica", sample.nobs = NULL,
      auto.bound = NULL, debug = FALSE)
# S3 method for lavaan
goric(object, ..., constraints = NULL,
      comparison = "unconstrained", type = "gorica",
      standardized = FALSE, auto.bound = NULL, debug = FALSE)
# S3 method for CTmeta
goric(object, ..., constraints = NULL,
      comparison = "unconstrained", type = "gorica", 
      sample.nobs = NULL, auto.bound = NULL, debug = FALSE)
# S3 method for rma
goric(object, ..., constraints = NULL,
      VCOV = NULL, comparison = "unconstrained", type = "gorica", 
      sample.nobs = NULL, auto.bound = NULL, debug = FALSE)

# S3 method for con_goric
print(x, digits = max(3, getOption("digits") - 4), ...)
# S3 method for con_goric
summary(object, brief = TRUE, digits = max(3, getOption("digits") - 4), ...)
# S3 method for con_goric
coef(object, ...)

Value

The function returns a dataframe with the log-likelihood, penalty term, GORIC(A) values and the GORIC(A) weights. Furthermore, a dataframe is returned with the relative GORIC(A) weights.

Arguments

object

an object containing the outcome of an unconstrained statistical analysis. Currently, the following objects can be processed:

• a fitted object of class restriktor.

• a fitted unconstrained object of class lm, rlm or glm.

• a numeric vector containing the unconstrained estimates resulting from any statistical analysis.

• a fitted object of class lavaan.

• a fitted object of class CTmeta.

• a fitted object of class rma.

x

an object of class con_goric.

...

this depends on the class of the object. If object is of class restriktor, further objects of class restriktor can be passed. Note that, the objects have to be of the same class. If object is of class lavaan, the standardized or unstandardized vcov can be used, using setting standardized = TRUE. See details for more information.

constraints

list; 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) or names(vector) 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\theta \ge rhs\), where each row represents one constraint. The number of columns needs to correspond to the number of parameters estimated (\(\theta\)) 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.

comparison

if "unconstrained" (default) the unconstrained model is included in the set of models. If "complement" then the restricted object is compared against its complement. Note that the complement can only be computed for one model/hypothesis at a time (for now). If "none" the model is only compared against the models provided by the user.

VCOV

variance-coviance matrix. Only needed if object is of class numeric and type = "gorica" or type = "goricac".

sample.nobs

the number of observations if type = "goricac". Note that, if type = "goricc", the number of observations are inherited from the fitted object.

type

if "goric", the generalized order-restricted information criterion value is computed. If "gorica" the log-likihood is computed using the multivariate normal distribution function.

missing

the default setting for objects of class "lm" is listwise: all cases with missing values are removed from the data before the analysis. This is only valid if the data are missing completely at random (MCAR). Another option is to use "two.stage". In this approach, missing data are imputed using an EM algorithm. However, we cannot use the complete data as input for futher analyses, because the resulting complete data variance-covariance matrix will not be correct. Therefore, we compute the correct aymptotic covariance (Savalei and Bentler, 2009) and use it as input for the goric.numeric function to compute a GORICA(C) value. Note that, the parameter estimates are also recomputed using the complete data.

auxiliary

Vector. The inclusion of auxiliary variables can improve the imputation model. These auxiliary variables are not part of the target model.

emControl

a list of control arguments for the emnorm function from the norm2 package.

standardized

if TRUE, standardized parameter estimates are used.

auto.bound

not used yet.

digits

the number of significant digits to use when printing.

debug

if TRUE, debugging information is printed out.

brief

if TRUE, a short overview is printed.

Author

Leonard Vanbrabant and Rebecca Kuiper

Details

The GORIC(A) values themselves are not interpretable and the interest lie in their differences. The GORIC(A) weights reflect the support of each hypothesis in the set. To compare two hypotheses (and not one to the whole set), one can examine the ratio of the two corresponding GORIC(A) weights. To avoid selecting a weakly supported hypothesis as the best one, the unconstrained hypothesis is usually included as safeguard.

In case of one order-constrained hypothesis, say H1, the complement Hc can be computed as competing hypothesis. The complement is defined as Hc = not H1.

If the object(s) is of class restriktor the constraints are inherited. Otherwise, the constraint syntax can be parsed via the constraint argument. If the object is an unconstrained model of class lm, rlm or glm, then the constraints can be specified in two ways, see restriktor. Note that if the constraints are written in matrix notation, then the constraints for each model/hypothesis is put in a named list with specific names amat, rhs, and neq. For example with three parameters x1, x2, x3, and x1 > 0: list(model1 = list(amat = rbind(c(1, 0, 0)), rhs = 0, neq = 0))). The rhs and neq are not required if they are equal to 0. If type = "gorica", then the object might be a (named) numeric vector. The constraints can again be specified in two ways, see restriktor. For examples, see below.

To determine the penalty term values, the chi-bar-square weights (a.k.a. level probabilities) must be computed. If "mix.weights = "pmvnorm" " (default), the chi-bar-square weights are computed based on the multivariate normal distribution function with additional Monte Carlo steps. If "mix.weights = "boot" ", the chi-bar-square weights are computed using parametric bootstrapping (see restriktor).

The "two.stage" approach for missing data uses the EM algorithm from the norm2 package. The response variables are assumed to be jointly normal. In practice, missing-data procedures designed for variables that are normal are sometimes applied to variables that are not. Binary and ordinal variables are sometimes imputed under a normal model, and the imputed values may be classified or rounded. This is also how restriktor handles (ordered) factors for now.

A better strategy (not implemented yet) would be to convert (ordered) factors into a pair of dummy variables. If the (ordered) factors have missing values, the dummy variables could be included as columns of Y and imputed, but then you have to decide how to convert the continuously distributed imputed values for these dummy codes back into categories.

References

Kuiper, R.M., Hoijtink, H., and Silvapulle, M.J. (2011). An Akaike-type information criterion for model selection under inequality constraints. Biometrika, 98, 2, 495--501.

Vanbrabant, L. and Kuiper, R. (2020). Evaluating a theory-based hypothesis against its complement using an AIC-type information criterion with an application to facial burn injury. Psychological Methods.

Victoria Savalei and Peter M. Bentler (2009) A Two-Stage Approach to Missing Data: Theory and Application to Auxiliary Variables, Structural Equation Modeling: A Multidisciplinary Journal, 16:3, 477-497, DOI: 10.1080/10705510903008238

Examples

library(MASS)
## lm
## unrestricted linear model for ages (in months) at which an 
## infant starts to walk alone.

# prepare data
DATA <- subset(ZelazoKolb1972, Group != "Control")
  
# fit unrestrikted linear model
fit1.lm <- lm(Age ~ Group, data = DATA)

# some artificial restrictions
H1 <- "GroupPassive > 0; GroupPassive < GroupNo"
H2 <- "GroupPassive > 0; GroupPassive > GroupNo"
H3 <- "GroupPassive = 0; GroupPassive < GroupNo"


# object is of class lm
goric(fit1.lm, constraints = list(H1 = H1, H2 = H2, H3 = H3))

# same result, but using objects of class restriktor
fit1.con <- restriktor(fit1.lm, constraints = H1)
fit2.con <- restriktor(fit1.lm, constraints = H2)
fit3.con <- restriktor(fit1.lm, constraints = H3)

goric(fit1.con, fit2.con, fit3.con)

# same result, but using the parameter estimates and covariance matrix as input
# Note, that in case of a numeric input only the gorica(c) can be computed. 
goric(coef(fit1.lm), VCOV = vcov(fit1.lm), constraints = list(H1 = H1, H2 = H2, H3 = H3))


# hypothesis H1 versus the complement (i.e., not H1)
goric(fit1.lm, constraints = list(H1 = H1), comparison = "complement")


## GORICA
# generate data
n <- 10
x1 <- rnorm(n)
x2 <- rnorm(n)
y <- 1 + x1 + x2 + rnorm(n)
# fit unconstrained linear model
fit.lm <- lm(y ~ x1 + x2)

# extract unconstrained estimates
est <- coef(fit.lm)
# unconstrained variance-covariance matrix
VCOV <- vcov(fit.lm)

## constraint syntax (character)
h1 <- "x1 > 0"
h2 <- "x1 > 0; x2 > 0"
# use fitted unconstrained linear model
goric(fit.lm, constraints = list(h1, h2), type = "gorica")
# use unconstrained estimates
goric(est, VCOV = VCOV, constraints = list(h1, h2), type = "gorica")

## constraint syntax (matrix notation)
h1 <- list(amat = c(0,1,0))
h2 <- list(amat = rbind(c(0,1,0), c(0,0,1)), rhs = c(0.5, 1), neq = 0)
goric(fit.lm, constraints = list(h1, h2), type = "gorica")
goric(est, VCOV = VCOV, constraints = list(h1, h2), type = "gorica")


## use parallel processing for computing the level probabilities (i.e., mixing weights)
if (FALSE) {
goric(fit.lm, constraints = list(H1), comparison = "complement",
      mix.weights = "boot", parallel = "snow", ncpus = 8, mix.bootstrap = 99999)
}