Nothing
lav_constraints_parse <- function(partable = NULL, constraints = NULL,
theta = NULL,
debug = FALSE) {
if (!missing(debug)) {
current.debug <- lav_debug()
if (lav_debug(debug))
on.exit(lav_debug(current.debug), TRUE)
}
# just in case we do not have a $free column in partable
if (is.null(partable$free)) {
partable$free <- seq_len(length(partable$lhs))
}
# from the partable: free parameters
if (!is.null(theta)) {
# nothing to do
} else if (!is.null(partable$est)) {
theta <- partable$est[partable$free > 0L]
} else if (!is.null(partable$start)) {
theta <- partable$start[partable$free > 0L]
} else {
theta <- rep(0, length(partable$lhs))
}
# number of free (but possibliy constrained) parameters
npar <- length(theta)
# parse the constraints
if (is.null(constraints)) {
LIST <- NULL
} else if (!is.character(constraints)) {
lav_msg_stop(gettext("constraints should be a string"))
} else {
FLAT <- lavParseModelString(constraints)
CON <- attr(FLAT, "constraints")
LIST <- list()
if (length(CON) > 0L) {
lhs <- unlist(lapply(CON, "[[", "lhs"))
op <- unlist(lapply(CON, "[[", "op"))
rhs <- unlist(lapply(CON, "[[", "rhs"))
LIST$lhs <- c(LIST$lhs, lhs)
LIST$op <- c(LIST$op, op)
LIST$rhs <- c(LIST$rhs, rhs)
} else {
lav_msg_stop(gettext("no constraints found in constraints argument"))
}
}
# simple equality constraints?
ceq.simple <- FALSE
if (!is.null(partable$unco)) {
ceq.simple <- TRUE
}
# variable definitions
def.function <- lav_partable_constraints_def(partable,
con = LIST,
debug = debug
)
# construct ceq/ciq functions
ceq.function <- lav_partable_constraints_ceq(partable,
con = LIST,
debug = debug
)
# linear or nonlinear?
ceq.linear.idx <- lav_constraints_linear_idx(
func = ceq.function,
npar = npar
)
ceq.nonlinear.idx <- lav_constraints_nonlinear_idx(
func = ceq.function,
npar = npar
)
# inequalities
cin.function <- lav_partable_constraints_ciq(partable,
con = LIST,
debug = debug
)
# linear or nonlinear?
cin.linear.idx <- lav_constraints_linear_idx(
func = cin.function,
npar = npar
)
cin.nonlinear.idx <- lav_constraints_nonlinear_idx(
func = cin.function,
npar = npar
)
# Jacobians
if (!is.null(body(ceq.function))) {
ceq.JAC <- try(lav_func_jacobian_complex(
func = ceq.function,
x = theta
), silent = TRUE)
if (inherits(ceq.JAC, "try-error")) { # eg. pnorm()
ceq.JAC <- lav_func_jacobian_simple(func = ceq.function, x = theta)
}
# constants
# do we have a non-zero 'rhs' elements? FIXME!!! is this reliable??
ceq.rhs <- -1 * ceq.function(numeric(npar))
# evaluate constraints
ceq.theta <- ceq.function(theta)
} else {
ceq.JAC <- matrix(0, nrow = 0L, ncol = npar)
ceq.rhs <- numeric(0L)
ceq.theta <- numeric(0L)
}
if (!is.null(body(cin.function))) {
cin.JAC <- try(lav_func_jacobian_complex(
func = cin.function,
x = theta
), silent = TRUE)
if (inherits(cin.JAC, "try-error")) { # eg. pnorm()
cin.JAC <- lav_func_jacobian_simple(func = cin.function, x = theta)
}
# constants
# do we have a non-zero 'rhs' elements? FIXME!!! is this reliable??
cin.rhs <- -1 * cin.function(numeric(npar))
# evaluate constraints
cin.theta <- cin.function(theta)
} else {
cin.JAC <- matrix(0, nrow = 0L, ncol = npar)
cin.rhs <- numeric(0L)
cin.theta <- numeric(0L)
}
# shortcut flags
ceq.linear.flag <- length(ceq.linear.idx) > 0L
ceq.nonlinear.flag <- length(ceq.nonlinear.idx) > 0L
ceq.flag <- ceq.linear.flag || ceq.nonlinear.flag
cin.linear.flag <- length(cin.linear.idx) > 0L
cin.nonlinear.flag <- length(cin.nonlinear.idx) > 0L
cin.flag <- cin.linear.flag || cin.nonlinear.flag
ceq.only.flag <- ceq.flag && !cin.flag
cin.only.flag <- cin.flag && !ceq.flag
ceq.linear.only.flag <- (ceq.linear.flag &&
!ceq.nonlinear.flag &&
!cin.flag)
ceq.simple.only <- ceq.simple && !ceq.flag && !cin.flag
# additional info if ceq.linear.flag
if (ceq.linear.flag) {
## NEW: 18 nov 2014: handle general *linear* constraints
##
## see Nocedal & Wright (2006) 15.3
## - from x to x.red:
## x.red <- MASS::ginv(Q2) %*% (x - Q1 %*% solve(t(R)) %*% b)
## or
## x.red <- as.numeric((x - b %*% qr.coef(QR,diag(npar))) %*% Q2)
##
## - from x.red to x
## x <- as.numeric(Q1 %*% solve(t(R)) %*% b + Q2 %*% x.red)
## or
## x <- as.numeric(b %*% qr.coef(QR, diag(npar))) +
## as.numeric(Q2 %*% x.red)
##
## we write eq.constraints.K = Q2
## eq.constraints.k0 = b %*% qr.coef(QR, diag(npar)))
# compute range+null space of the jacobion (JAC) of the constraint
# matrix
# JAC <- lav_func_jacobian_complex(func = ceq.function,
# x = lavpartable$start[lavpartable$free > 0L]
QR <- qr(t(ceq.JAC))
ranK <- QR$rank
Q <- qr.Q(QR, complete = TRUE)
# Q1 <- Q[,1:ranK, drop = FALSE] # range space
# Q2 <- Q[,-seq_len(ranK), drop = FALSE] # null space
# R <- qr.R(QR)
ceq.JAC.NULL <- Q[, -seq_len(ranK), drop = FALSE]
if (all(ceq.rhs == 0)) {
ceq.rhs.NULL <- numeric(npar)
} else {
tmp <- qr.coef(QR, diag(npar))
NA.idx <- which(is.na(rowSums(tmp))) # catch NAs
if (length(NA.idx) > 0L) {
tmp[NA.idx, ] <- 0
}
ceq.rhs.NULL <- as.numeric(ceq.rhs %*% tmp)
}
} else {
ceq.JAC.NULL <- matrix(0, 0L, 0L)
ceq.rhs.NULL <- numeric(0L)
}
# if simple equalities only, create 'K' matrix
ceq.simple.K <- matrix(0, 0, 0)
if (ceq.simple.only) {
n.unco <- max(partable$unco)
n.free <- max(partable$free)
ceq.simple.K <- matrix(0, nrow = n.unco, ncol = n.free)
#####
##### FIXME !
#####
idx.free <- partable$free[partable$free > 0]
for (k in 1:n.unco) {
c <- idx.free[k]
ceq.simple.K[k, c] <- 1
}
}
# dummy jacobian 'function'
ceq.jacobian <- function() NULL
cin.jacobian <- function() NULL
OUT <- list(
def.function = def.function,
ceq.function = ceq.function,
ceq.JAC = ceq.JAC,
ceq.jacobian = ceq.jacobian,
ceq.rhs = ceq.rhs,
ceq.theta = ceq.theta,
ceq.linear.idx = ceq.linear.idx,
ceq.nonlinear.idx = ceq.nonlinear.idx,
ceq.linear.flag = ceq.linear.flag,
ceq.nonlinear.flag = ceq.nonlinear.flag,
ceq.flag = ceq.flag,
ceq.linear.only.flag = ceq.linear.only.flag,
ceq.JAC.NULL = ceq.JAC.NULL,
ceq.rhs.NULL = ceq.rhs.NULL,
ceq.simple.only = ceq.simple.only,
ceq.simple.K = ceq.simple.K,
cin.function = cin.function,
cin.JAC = cin.JAC,
cin.jacobian = cin.jacobian,
cin.rhs = cin.rhs,
cin.theta = cin.theta,
cin.linear.idx = cin.linear.idx,
cin.nonlinear.idx = cin.nonlinear.idx,
cin.linear.flag = cin.linear.flag,
cin.nonlinear.flag = cin.nonlinear.flag,
cin.flag = cin.flag,
cin.only.flag = cin.only.flag
)
OUT
}
lav_constraints_linear_idx <- function(func = NULL, npar = NULL) {
if (is.null(func) || is.null(body(func))) {
return(integer(0L))
}
# seed 1: rnorm
A0 <- lav_func_jacobian_complex(func = func, x = rnorm(npar))
# seed 2: rnorm
A1 <- lav_func_jacobian_complex(func = func, x = rnorm(npar))
A0minA1 <- A0 - A1
linear <- apply(A0minA1, 1, function(x) all(x == 0))
which(linear)
}
lav_constraints_nonlinear_idx <- function(func = NULL, npar = NULL) {
if (is.null(func) || is.null(body(func))) {
return(integer(0L))
}
# seed 1: rnorm
A0 <- lav_func_jacobian_complex(func = func, x = rnorm(npar))
# seed 2: rnorm
A1 <- lav_func_jacobian_complex(func = func, x = rnorm(npar))
A0minA1 <- A0 - A1
linear <- apply(A0minA1, 1, function(x) all(x == 0))
which(!linear)
}
# FIXME: is there a more elegant/robust way to do this??
lav_constraints_check_linear <- function(model) {
# seed 1: rnorm
A.ceq <- A.cin <- matrix(0, model@nx.free, 0)
if (!is.null(body(model@ceq.function))) {
A.ceq <- t(lav_func_jacobian_complex(func = model@ceq.function, x = rnorm(model@nx.free)))
}
if (!is.null(body(model@cin.function))) {
A.cin <- t(lav_func_jacobian_complex(func = model@cin.function, x = rnorm(model@nx.free)))
}
A0 <- cbind(A.ceq, A.cin)
# seed 2: rnorm
A.ceq <- A.cin <- matrix(0, model@nx.free, 0)
if (!is.null(body(model@ceq.function))) {
A.ceq <- t(lav_func_jacobian_complex(func = model@ceq.function, x = rnorm(model@nx.free)))
}
if (!is.null(body(model@cin.function))) {
A.cin <- t(lav_func_jacobian_complex(func = model@cin.function, x = rnorm(model@nx.free)))
}
A1 <- cbind(A.ceq, A.cin)
A0minA1 <- all.equal(A0, A1)
if (is.logical(A0minA1) && A0minA1 == TRUE) {
return(TRUE)
} else {
return(FALSE)
}
}
# check if the equality constraints are 'simple' (a == b)
lav_constraints_check_simple <- function(lavmodel = NULL) {
ones <- (lavmodel@ceq.JAC == 1 | lavmodel@ceq.JAC == -1)
simple <- all(lavmodel@ceq.rhs == 0) &&
all(apply(lavmodel@ceq.JAC != 0, 1, sum) == 2) &&
all(apply(ones, 1, sum) == 2) &&
length(lavmodel@ceq.nonlinear.idx) == 0
# TRUE or FALSE
simple
}
lav_constraints_R2K <- function(lavmodel = NULL, R = NULL) {
# constraint matrix
if (!is.null(lavmodel)) {
R <- lavmodel@ceq.JAC
}
stopifnot(!is.null(R))
npar.full <- NCOL(R)
npar.red <- npar.full - NROW(R)
K <- diag(npar.full)
for (i in 1:NROW(R)) {
idx1 <- which(R[i, ] == 1)
idx2 <- which(R[i, ] == -1)
K[idx2, idx1] <- 1
}
# remove redundant columns
neg.idx <- which(colSums(R) < 0)
K <- K[, -neg.idx]
K
}
lav_constraints_lambda_pre <- function(lavobject = NULL, method = "Don") {
# compute factor 'pre' so that pre %*% g = lambda
method <- tolower(method)
R <- lavobject@Model@con.jac[, ]
if (is.null(R) || length(R) == 0L) {
return(numeric(0L))
}
INFO <- lavTech(lavobject, "information.first.order")
npar <- nrow(INFO)
# Don 1985
if (method == "don") {
R.plus <- MASS::ginv(R)
# construct augmented matrix
Z <- rbind(
cbind(INFO, t(R)),
cbind(R, matrix(0, nrow = nrow(R), ncol = nrow(R)))
)
Z.plus <- MASS::ginv(Z)
P.star <- Z.plus[1:npar, 1:npar]
PRE <- t(R.plus) %*% (diag(npar) - INFO %*% P.star)
# Bentler EQS manual
} else if (method == "bentler") {
INFO.inv <- solve(INFO)
PRE <- solve(R %*% INFO.inv %*% t(R)) %*% R %*% INFO.inv
}
PRE
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.