Nothing
R/lav_model_estimate.R
In lavaan: Latent Variable Analysis
Defines functions
lav_model_estimate
# model estimation
lav_model_estimate <- function(lavmodel = NULL,
lavpartable = NULL, # for parscale = "stand"
lavh1 = NULL, # for multilevel + parsc
lavsamplestats = NULL,
lavdata = NULL,
lavoptions = NULL,
lavcache = list(),
start = "model",
do.fit = TRUE) {
estimator <- lavoptions$estimator
verbose <- lavoptions$verbose
debug <- lavoptions$debug
ngroups <- lavsamplestats@ngroups
if(lavsamplestats@missing.flag || estimator == "PML") {
group.weight <- FALSE
} else {
group.weight <- TRUE
}
# backwards compatibility < 0.6-11
if(is.null(lavoptions$optim.partrace)) {
lavoptions$optim.partrace <- FALSE
}
if(lavoptions$optim.partrace) {
# fx + parameter values
PENV <- new.env()
PENV$PARTRACE <- matrix(NA, nrow=0, ncol=lavmodel@nx.free + 1L)
}
# starting values (ignoring equality constraints)
x.unpack <- lav_model_get_parameters(lavmodel)
# override? use simple instead? (new in 0.6-7)
if(start == "simple") {
START <- numeric(length(lavpartable$lhs))
# set loadings to 0.7
loadings.idx <- which(lavpartable$free > 0L &
lavpartable$op == "=~")
if(length(loadings.idx) > 0L) {
START[loadings.idx] <- 0.7
}
# set (only) variances to 1
var.idx <- which(lavpartable$free > 0L &
lavpartable$op == "~~" &
lavpartable$lhs == lavpartable$rhs)
if(length(var.idx) > 0L) {
START[var.idx] <- 1
}
if(lavmodel@ceq.simple.only) {
x.unpack <- START[ lavpartable$free > 0L &
!duplicated(lavpartable$free) ]
} else {
x.unpack <- START[ lavpartable$free > 0L ]
}
}
# 1. parameter scaling (to handle data scaling, not parameter scaling)
parscale <- rep(1.0, length(x.unpack))
# for < 0.6 compatibility
if(is.null(lavoptions$optim.parscale)) {
lavoptions$optim.parscale <- "none"
}
if(lavoptions$optim.parscale == "none") {
# do nothing, but still set SCALE, as before
} else if(lavoptions$optim.parscale %in% c("stand", "st", "standardize",
"standarized", "stand.all")) {
# rescale parameters as if the data was standardized
# new in 0.6-2
#
# FIXME: this works well, as long as the variances of the
# latent variables (which we do not know) are more or less
# equal to 1.0 (eg std.lv = TRUE)
#
# Once we have better estimates of those variances, we could
# use them to set the scale
#
if(lavdata@nlevels > 1L) {
if(length(lavh1) > 0L) {
OV.VAR <- lapply(lavh1$implied$cov, diag)
} else {
OV.VAR <- lapply(do.call(c, lapply(lavdata@Lp, "[[", "ov.idx")),
function(x) rep(1, length(x) ))
}
} else {
if(lavoptions$conditional.x) {
OV.VAR <- lavsamplestats@res.var
} else {
OV.VAR <- lavsamplestats@var
}
}
if(lavoptions$std.lv) {
parscale <- lav_standardize_all(lavobject = NULL,
est = rep(1, length(lavpartable$lhs)),
est.std = rep(1, length(lavpartable$lhs)),
cov.std = FALSE, ov.var = OV.VAR,
lavmodel = lavmodel, lavpartable = lavpartable,
cov.x = lavsamplestats@cov.x)
} else {
# needs good estimates for lv variances!
# if there is a single 'marker' indicator, we could use
# its observed variance as an upper bound
# for the moment, set them to 1.0 (instead of 0.05)
# TODO: USE Bentler's 1982 approach to get an estimate of
# VETA; use those diagonal elements...
# but only if we have 'marker' indicators for each LV
LV.VAR <- vector("list", lavmodel@ngroups)
for(g in seq_len(lavmodel@ngroups)) {
mm.in.group <- 1:lavmodel@nmat[g] + cumsum(c(0,lavmodel@nmat))[g]
MLIST <- lavmodel@GLIST[ mm.in.group ]
LAMBDA <- MLIST$lambda
n.lv <- ncol(LAMBDA)
LV.VAR[[g]] <- rep(1.0, n.lv)
}
parscale <- lav_standardize_all(lavobject = NULL,
est = rep(1, length(lavpartable$lhs)),
#est.std = rep(1, length(lavpartable$lhs)),
# here, we use whatever the starting values are
# for the latent variances...
cov.std = FALSE, ov.var = OV.VAR,
lv.var = LV.VAR,
lavmodel = lavmodel, lavpartable = lavpartable,
cov.x = lavsamplestats@cov.x)
}
# in addition, take sqrt for variance parameters
var.idx <- which(lavpartable$op == "~~" &
lavpartable$lhs == lavpartable$rhs)
if(length(var.idx) > 0L) {
parscale[var.idx] <- sqrt(abs(parscale[var.idx]))
}
if(lavmodel@ceq.simple.only) {
parscale <- parscale[ lavpartable$free > 0 &
!duplicated(lavpartable$free) ]
} else {
parscale <- parscale[ lavpartable$free > 0 ]
}
}
# parscale should obey the equality constraints
if(lavmodel@eq.constraints && lavoptions$optim.parscale != "none") {
# pack
p.pack <- as.numeric( (parscale - lavmodel@eq.constraints.k0) %*%
lavmodel@eq.constraints.K )
# unpack
parscale <- as.numeric(lavmodel@eq.constraints.K %*% p.pack) +
lavmodel@eq.constraints.k0
}
if(debug) {
cat("parscale = ", parscale, "\n")
}
z.unpack <- x.unpack * parscale
# 2. pack (apply equality constraints)
if(lavmodel@eq.constraints) {
z.pack <- as.numeric( (z.unpack - lavmodel@eq.constraints.k0) %*%
lavmodel@eq.constraints.K )
} else {
z.pack <- z.unpack
}
# 3. transform (already constrained) variances to standard deviations?
# TODO
#if(lavoptions$optim.var.transform == "sqrt" &&
# length(lavmodel@x.free.var.idx) > 0L) {
# # transforming variances using atan (or another sigmoid function?)
# # FIXME: better approach?
# #start.x[lavmodel@x.free.var.idx] <-
# # atan(start.x[lavmodel@x.free.var.idx])
# start.x[lavmodel@x.free.var.idx] <-
# sqrt(start.x[lavmodel@x.free.var.idx]) # assuming positive var
#}
# final starting values for optimizer
start.x <- z.pack
if(debug) {
cat("start.x = ", start.x, "\n")
}
# user-specified bounds? (new in 0.6-2)
if(is.null(lavpartable$lower)) {
lower <- -Inf
} else {
if(lavmodel@ceq.simple.only) {
free.idx <- which(lavpartable$free > 0L &
!duplicated(lavpartable$free))
lower <- lavpartable$lower[free.idx]
} else if(lavmodel@eq.constraints) {
# bounds have no effect any longer....
warning("lavaan warning: bounds have no effect in the presence of linear equality constraints")
lower <- -Inf
} else {
lower <- lavpartable$lower[ lavpartable$free > 0L ]
}
}
if(is.null(lavpartable$upper)) {
upper <- +Inf
} else {
if(lavmodel@ceq.simple.only) {
free.idx <- which(lavpartable$free > 0L &
!duplicated(lavpartable$free))
upper <- lavpartable$upper[free.idx]
} else if(lavmodel@eq.constraints) {
# bounds have no effect any longer....
if(is.null(lavpartable$lower)) {
# bounds have no effect any longer....
warning("lavaan warning: bounds have no effect in the presence of linear equality constraints")
}
upper <- +Inf
} else {
upper <- lavpartable$upper[ lavpartable$free > 0L ]
}
}
# check for inconsistent lower/upper bounds
# this may happen if we have equality constraints; qr() may switch
# the sign...
bad.idx <- which(lower > upper)
if(length(bad.idx) > 0L) {
# switch
#tmp <- lower[bad.idx]
#lower[bad.idx] <- upper[bad.idx]
#upper[bad.idx] <- tmp
lower[bad.idx] <- -Inf
upper[bad.idx] <- +Inf
}
# function to be minimized
objective_function <- function(x, verbose = FALSE, infToMax = FALSE,
debug = FALSE) {
# 3. standard deviations to variances
# WARNING: x is still packed here!
#if(lavoptions$optim.var.transform == "sqrt" &&
# length(lavmodel@x.free.var.idx) > 0L) {
# #x[lavmodel@x.free.var.idx] <- tan(x[lavmodel@x.free.var.idx])
# x.var <- x[lavmodel@x.free.var.idx]
# x.var.sign <- sign(x.var)
# x[lavmodel@x.free.var.idx] <- x.var.sign * (x.var * x.var) # square!
#}
# 2. unpack
if(lavmodel@eq.constraints) {
x <- as.numeric(lavmodel@eq.constraints.K %*% x) +
lavmodel@eq.constraints.k0
}
# 1. unscale
x <- x / parscale
# update GLIST (change `state') and make a COPY!
GLIST <- lav_model_x2GLIST(lavmodel, x = x)
fx <- lav_model_objective(lavmodel = lavmodel,
GLIST = GLIST,
lavsamplestats = lavsamplestats,
lavdata = lavdata,
lavcache = lavcache,
verbose = verbose)
# only for PML: divide by N (to speed up convergence)
if(estimator == "PML") {
fx <- fx / lavsamplestats@ntotal
}
if(debug || verbose) {
cat(" objective function = ",
sprintf("%18.16f", fx), "\n", sep="")
}
if(debug) {
#cat("Current unconstrained parameter values =\n")
#tmp.x <- lav_model_get_parameters(lavmodel, GLIST=GLIST, type="unco")
#print(tmp.x); cat("\n")
cat("Current free parameter values =\n"); print(x); cat("\n")
}
if(lavoptions$optim.partrace) {
PENV$PARTRACE <- rbind(PENV$PARTRACE, c(fx, x))
}
# for L-BFGS-B
#if(infToMax && is.infinite(fx)) fx <- 1e20
if(!is.finite(fx)) {
fx.group <- attr(fx, "fx.group")
fx <- 1e20
attr(fx, "fx.group") <- fx.group # only for lav_model_fit()
}
fx
}
gradient_function <- function(x, verbose = FALSE, infToMax = FALSE,
debug = FALSE) {
# transform variances back
#if(lavoptions$optim.var.transform == "sqrt" &&
# length(lavmodel@x.free.var.idx) > 0L) {
# #x[lavmodel@x.free.var.idx] <- tan(x[lavmodel@x.free.var.idx])
# x.var <- x[lavmodel@x.free.var.idx]
# x.var.sign <- sign(x.var)
# x[lavmodel@x.free.var.idx] <- x.var.sign * (x.var * x.var) # square!
#}
# 2. unpack
if(lavmodel@eq.constraints) {
x <- as.numeric(lavmodel@eq.constraints.K %*% x) +
lavmodel@eq.constraints.k0
}
# 1. unscale
x <- x / parscale
# update GLIST (change `state') and make a COPY!
GLIST <- lav_model_x2GLIST(lavmodel, x = x)
dx <- lav_model_gradient(lavmodel = lavmodel,
GLIST = GLIST,
lavsamplestats = lavsamplestats,
lavdata = lavdata,
lavcache = lavcache,
type = "free",
group.weight = group.weight, ### check me!!
verbose = verbose,
ceq.simple = lavmodel@ceq.simple.only)
if(debug) {
cat("Gradient function (analytical) =\n"); print(dx); cat("\n")
}
# 1. scale (note: divide, not multiply!)
dx <- dx / parscale
# 2. pack
if(lavmodel@eq.constraints) {
dx <- as.numeric( dx %*% lavmodel@eq.constraints.K )
}
# 3. transform variances back
#if(lavoptions$optim.var.transform == "sqrt" &&
# length(lavmodel@x.free.var.idx) > 0L) {
# x.var <- x[lavmodel@x.free.var.idx] # here in 'var' metric
# x.var.sign <- sign(x.var)
# x.var <- abs(x.var)
# x.sd <- sqrt(x.var)
# dx[lavmodel@x.free.var.idx] <-
# ( 2 * x.var.sign * dx[lavmodel@x.free.var.idx] * x.sd )
#}
# only for PML: divide by N (to speed up convergence)
if(estimator == "PML") {
dx <- dx / lavsamplestats@ntotal
}
if(debug) {
cat("Gradient function (analytical, after eq.constraints.K) =\n"); print(dx); cat("\n")
}
dx
}
gradient_function_numerical <- function(x, verbose=FALSE, debug = FALSE) {
# NOTE: no need to 'tranform' anything here (var/eq)
# this is done anyway in objective_function
# numerical approximation using the Richardson method
npar <- length(x)
h <- 10e-6
dx <- numeric( npar )
## FIXME: call lav_model_objective directly!!
for(i in 1:npar) {
x.left <- x.left2 <- x.right <- x.right2 <- x
x.left[i] <- x[i] - h; x.left2[i] <- x[i] - 2*h
x.right[i] <- x[i] + h; x.right2[i] <- x[i] + 2*h
fx.left <- objective_function(x.left, verbose = FALSE, debug = FALSE)
fx.left2 <- objective_function(x.left2, verbose = FALSE, debug = FALSE)
fx.right <- objective_function(x.right, verbose = FALSE, debug = FALSE)
fx.right2 <- objective_function(x.right2, verbose = FALSE, debug = FALSE)
dx[i] <- (fx.left2 - 8*fx.left + 8*fx.right - fx.right2)/(12*h)
}
#dx <- lavGradientC(func=objective_function, x=x)
# does not work if pnorm is involved... (eg PML)
if(debug) {
cat("Gradient function (numerical) =\n"); print(dx); cat("\n")
}
dx
}
gradient_function_numerical_complex <- function(x, verbose=FALSE, debug = FALSE) {
dx <- Re(lav_func_gradient_complex(func = objective_function, x = x,
h = sqrt(.Machine$double.eps)))
# does not work if pnorm is involved... (eg PML)
if(debug) {
cat("Gradient function (numerical complex) =\n"); print(dx); cat("\n")
}
dx
}
# check if the initial values produce a positive definite Sigma
# to begin with -- but only for estimator="ML"
if(estimator %in% c("ML","FML","MML")) {
Sigma.hat <- computeSigmaHat(lavmodel, extra=TRUE, debug=lavoptions$debug)
for(g in 1:ngroups) {
if(!attr(Sigma.hat[[g]], "po")) {
group.txt <- ifelse(ngroups > 1,
paste(" in group ",g,".",sep=""), ".")
if(debug) {
print(Sigma.hat[[g]][,])
}
warning("lavaan ERROR: initial model-implied matrix (Sigma) is not positive definite;\n check your model and/or starting parameters", group.txt)
x <- start.x
fx <- as.numeric(NA)
attr(fx, "fx.group") <- rep(as.numeric(NA), ngroups)
attr(x, "converged") <- FALSE
attr(x, "iterations") <- 0L
attr(x, "control") <- lavoptions@control
attr(x, "fx") <- fx
return(x)
}
}
}
# parameter scaling
# FIXME: what is the best way to set the scale??
# current strategy: if startx > 1.0, we rescale by using
# 1/startx
SCALE <- rep(1.0, length(start.x))
if(lavoptions$optim.parscale == "none") {
idx <- which(abs(start.x) > 1.0)
if(length(idx) > 0L) {
SCALE[idx] <- abs(1.0/start.x[idx])
}
}
if(debug) {
cat("SCALE = ", SCALE, "\n")
}
# first try: check if starting values return a finite value
fx <- objective_function(start.x, verbose = verbose, debug = debug)
if(!is.finite(fx)) {
# emergency change of start.x
start.x <- start.x / 10
}
# first some nelder mead steps? (default = FALSE)
INIT_NELDER_MEAD <- lavoptions$optim.init_nelder_mead
# gradient: analytic, numerical or NULL?
if(is.character(lavoptions$optim.gradient)) {
if(lavoptions$optim.gradient %in% c("analytic","analytical")) {
GRADIENT <- gradient_function
} else if(lavoptions$optim.gradient %in% c("numerical", "numeric")) {
GRADIENT <- gradient_function_numerical
} else if(lavoptions$optim.gradient %in% c("numeric.complex", "complex")){
GRADIENT <- gradient_function_numerical_complex
} else if(lavoptions$optim.gradient %in% c("NULL", "null")) {
GRADIENT <- NULL
} else {
warning("lavaan WARNING: gradient should be analytic, numerical or NULL")
}
} else if(is.logical(lavoptions$optim.gradient)) {
if(lavoptions$optim.gradient) {
GRADIENT <- gradient_function
} else {
GRADIENT <- NULL
}
} else if(is.null(lavoptions$optim.gradient)) {
GRADIENT <- gradient_function
}
# default optimizer
if(length(lavmodel@ceq.nonlinear.idx) == 0L &&
length(lavmodel@cin.linear.idx) == 0L &&
length(lavmodel@cin.nonlinear.idx) == 0L) {
if(is.null(lavoptions$optim.method)) {
OPTIMIZER <- "NLMINB"
#OPTIMIZER <- "BFGS" # slightly slower, no bounds; better scaling!
#OPTIMIZER <- "L-BFGS-B" # trouble with Inf values for fx!
} else {
OPTIMIZER <- toupper(lavoptions$optim.method)
stopifnot(OPTIMIZER %in% c("NLMINB0", "NLMINB1", "NLMINB2",
"NLMINB", "BFGS", "L-BFGS-B", "NONE"))
if(OPTIMIZER == "NLMINB1") {
OPTIMIZER <- "NLMINB"
}
}
} else {
if(is.null(lavoptions$optim.method)) {
OPTIMIZER <- "NLMINB.CONSTR"
} else {
OPTIMIZER <- toupper(lavoptions$optim.method)
stopifnot(OPTIMIZER %in% c("NLMINB.CONSTR", "NLMINB", "NONE"))
}
if(OPTIMIZER == "NLMINB") {
OPTIMIZER <- "NLMINB.CONSTR"
}
}
if(INIT_NELDER_MEAD) {
if(verbose) cat(" initial Nelder-Mead step:\n")
trace <- 0L; if(verbose) trace <- 1L
optim.out <- optim(par=start.x,
fn=objective_function,
method="Nelder-Mead",
#control=list(maxit=10L,
# parscale=SCALE,
# trace=trace),
hessian=FALSE,
verbose=verbose, debug=debug)
cat("\n")
start.x <- optim.out$par
}
if(OPTIMIZER == "NLMINB0") {
if(verbose) cat(" quasi-Newton steps using NLMINB0 (no analytic gradient):\n")
#if(debug) control$trace <- 1L;
control.nlminb <- list(eval.max=20000L,
iter.max=10000L,
trace=0L,
#abs.tol=1e-20, ### important!! fx never negative
abs.tol=(.Machine$double.eps * 10),
rel.tol=1e-10,
#step.min=2.2e-14, # in =< 0.5-12
step.min=1.0, # 1.0 in < 0.5-21
step.max=1.0,
x.tol=1.5e-8,
xf.tol=2.2e-14)
control.nlminb <- modifyList(control.nlminb, lavoptions$control)
control <- control.nlminb[c("eval.max", "iter.max", "trace",
"step.min", "step.max",
"abs.tol", "rel.tol", "x.tol", "xf.tol")]
#cat("DEBUG: control = "); print(str(control.nlminb)); cat("\n")
optim.out <- nlminb(start=start.x,
objective=objective_function,
gradient=NULL,
lower=lower,
upper=upper,
control=control,
scale=SCALE,
verbose=verbose, debug=debug)
if(verbose) {
cat(" convergence status (0=ok): ", optim.out$convergence, "\n")
cat(" nlminb message says: ", optim.out$message, "\n")
cat(" number of iterations: ", optim.out$iterations, "\n")
cat(" number of function evaluations [objective, gradient]: ",
optim.out$evaluations, "\n")
}
# try again
if(optim.out$convergence != 0L) {
optim.out <- nlminb(start=start.x,
objective=objective_function,
gradient=NULL,
lower=lower,
upper=upper,
control=control,
scale=SCALE,
verbose=verbose, debug=debug)
}
iterations <- optim.out$iterations
x <- optim.out$par
if(optim.out$convergence == 0L) {
converged <- TRUE
} else {
converged <- FALSE
}
} else if(OPTIMIZER == "NLMINB") {
if(verbose) cat(" quasi-Newton steps using NLMINB:\n")
#if(debug) control$trace <- 1L;
control.nlminb <- list(eval.max=20000L,
iter.max=10000L,
trace=0L,
#abs.tol=1e-20, ### important!! fx never negative
abs.tol=(.Machine$double.eps * 10),
rel.tol=1e-10,
#step.min=2.2e-14, # in =< 0.5-12
step.min=1.0, # 1.0 in < 0.5-21
step.max=1.0,
x.tol=1.5e-8,
xf.tol=2.2e-14)
control.nlminb <- modifyList(control.nlminb, lavoptions$control)
control <- control.nlminb[c("eval.max", "iter.max", "trace",
"step.min", "step.max",
"abs.tol", "rel.tol", "x.tol", "xf.tol")]
#cat("DEBUG: control = "); print(str(control.nlminb)); cat("\n")
optim.out <- nlminb(start=start.x,
objective=objective_function,
gradient=GRADIENT,
lower=lower,
upper=upper,
control=control,
scale=SCALE,
verbose=verbose, debug=debug)
if(verbose) {
cat(" convergence status (0=ok): ", optim.out$convergence, "\n")
cat(" nlminb message says: ", optim.out$message, "\n")
cat(" number of iterations: ", optim.out$iterations, "\n")
cat(" number of function evaluations [objective, gradient]: ",
optim.out$evaluations, "\n")
}
iterations <- optim.out$iterations
x <- optim.out$par
if(optim.out$convergence == 0L) {
converged <- TRUE
} else {
converged <- FALSE
}
} else if(OPTIMIZER == "BFGS") {
# warning: Bollen example with estimator=GLS does NOT converge!
# (but WLS works!)
# - BB.ML works too
control.bfgs <- list(trace=0L, fnscale=1,
parscale=SCALE, ## or not?
ndeps=1e-3,
maxit=10000,
abstol=1e-20,
reltol=1e-10,
REPORT=1L)
control.bfgs <- modifyList(control.bfgs, lavoptions$control)
control <- control.bfgs[c("trace", "fnscale", "parscale", "ndeps",
"maxit", "abstol", "reltol", "REPORT")]
#trace <- 0L; if(verbose) trace <- 1L
optim.out <- optim(par=start.x,
fn=objective_function,
gr=GRADIENT,
method="BFGS",
control=control,
hessian=FALSE,
verbose=verbose, debug=debug)
if(verbose) {
cat(" convergence status (0=ok): ", optim.out$convergence, "\n")
cat(" optim BFGS message says: ", optim.out$message, "\n")
#cat("number of iterations: ", optim.out$iterations, "\n")
cat(" number of function evaluations [objective, gradient]: ",
optim.out$counts, "\n")
}
#iterations <- optim.out$iterations
iterations <- optim.out$counts[1]
x <- optim.out$par
if(optim.out$convergence == 0L) {
converged <- TRUE
} else {
converged <- FALSE
}
} else if(OPTIMIZER == "L-BFGS-B") {
# warning, does not cope with Inf values!!
control.lbfgsb <- list(trace=0L, fnscale=1,
parscale=SCALE, ## or not?
ndeps=1e-3,
maxit=10000,
REPORT=1L,
lmm=5L,
factr=1e7,
pgtol=0)
control.lbfgsb <- modifyList(control.lbfgsb, lavoptions$control)
control <- control.lbfgsb[c("trace", "fnscale", "parscale",
"ndeps", "maxit", "REPORT", "lmm",
"factr", "pgtol")]
optim.out <- optim(par=start.x,
fn=objective_function,
gr=GRADIENT,
method="L-BFGS-B",
lower=lower,
upper=upper,
control=control,
hessian=FALSE,
verbose=verbose, debug=debug,
infToMax=TRUE)
if(verbose) {
cat(" convergence status (0=ok): ", optim.out$convergence, "\n")
cat(" optim L-BFGS-B message says: ", optim.out$message, "\n")
#cat("number of iterations: ", optim.out$iterations, "\n")
cat(" number of function evaluations [objective, gradient]: ",
optim.out$counts, "\n")
}
#iterations <- optim.out$iterations
iterations <- optim.out$counts[1]
x <- optim.out$par
if(optim.out$convergence == 0L) {
converged <- TRUE
} else {
converged <- FALSE
}
} else if(OPTIMIZER == "NLMINB.CONSTR") {
ocontrol <- list(verbose=verbose)
if(!is.null(lavoptions$control$control.outer)) {
ocontrol <- c(lavoptions$control$control.outer, verbose=verbose)
}
control.nlminb <- list(eval.max=20000L,
iter.max=10000L,
trace=0L,
#abs.tol=1e-20,
abs.tol=(.Machine$double.eps * 10),
rel.tol=1e-9, # 1e-10 seems 'too strict'
step.min=1.0, # 1.0 in < 0.5-21
step.max=1.0,
x.tol=1.5e-8,
xf.tol=2.2e-14)
control.nlminb <- modifyList(control.nlminb, lavoptions$control)
control <- control.nlminb[c("eval.max", "iter.max", "trace",
"abs.tol", "rel.tol")]
cin <- cin.jac <- ceq <- ceq.jac <- NULL
if(!is.null(body(lavmodel@cin.function))) cin <- lavmodel@cin.function
if(!is.null(body(lavmodel@cin.jacobian))) cin.jac <- lavmodel@cin.jacobian
if(!is.null(body(lavmodel@ceq.function))) ceq <- lavmodel@ceq.function
if(!is.null(body(lavmodel@ceq.jacobian))) ceq.jac <- lavmodel@ceq.jacobian
trace <- FALSE; if(verbose) trace <- TRUE
optim.out <- nlminb.constr(start = start.x,
objective=objective_function,
gradient=GRADIENT,
control=control,
scale=SCALE,
verbose=verbose, debug=debug,
lower=lower,
upper=upper,
cin = cin, cin.jac = cin.jac,
ceq = ceq, ceq.jac = ceq.jac,
control.outer = ocontrol
)
if(verbose) {
cat(" convergence status (0=ok): ", optim.out$convergence, "\n")
cat(" nlminb.constr message says: ", optim.out$message, "\n")
cat(" number of outer iterations: ", optim.out$outer.iterations, "\n")
cat(" number of inner iterations: ", optim.out$iterations, "\n")
cat(" number of function evaluations [objective, gradient]: ",
optim.out$evaluations, "\n")
}
iterations <- optim.out$iterations
x <- optim.out$par
if(optim.out$convergence == 0) {
converged <- TRUE
} else {
converged <- FALSE
}
} else if(OPTIMIZER == "NONE") {
x <- start.x
iterations <- 0L
converged <- TRUE
control <- list()
# if inequality constraints, add con.jac/lambda
# needed for df!
if(length(lavmodel@ceq.nonlinear.idx) == 0L &&
length(lavmodel@cin.linear.idx) == 0L &&
length(lavmodel@cin.nonlinear.idx) == 0L) {
optim.out <- list()
} else {
# if inequality constraints, add con.jac/lambda
# needed for df!
optim.out <- list()
if(is.null(body(lavmodel@ceq.function))) {
ceq <- function(x, ...) { return( numeric(0) ) }
} else {
ceq <- lavmodel@ceq.function
}
if(is.null(body(lavmodel@cin.function))) {
cin <- function(x, ...) { return( numeric(0) ) }
} else {
cin <- lavmodel@cin.function
}
ceq0 <- ceq(start.x)
cin0 <- cin(start.x)
con0 <- c(ceq0, cin0)
JAC <- rbind(numDeriv::jacobian(ceq, x = start.x),
numDeriv::jacobian(cin, x = start.x))
nceq <- length( ceq(start.x) )
ncin <- length( cin(start.x) )
ncon <- nceq + ncin
ceq.idx <- cin.idx <- integer(0)
if(nceq > 0L) ceq.idx <- 1:nceq
if(ncin > 0L) cin.idx <- nceq + 1:ncin
cin.flag <- rep(FALSE, length(ncon))
if(ncin > 0L) cin.flag[cin.idx] <- TRUE
inactive.idx <- integer(0L)
cin.idx <- which(cin.flag)
if(ncin > 0L) {
slack <- 1e-05
inactive.idx <- which(cin.flag & con0 > slack)
}
attr(JAC, "inactive.idx") <- inactive.idx
attr(JAC, "cin.idx") <- cin.idx
attr(JAC, "ceq.idx") <- ceq.idx
optim.out$con.jac <- JAC
optim.out$lambda <- rep(0, ncon)
}
}
fx <- objective_function(x) # to get "fx.group" attribute
# check convergence
warn.txt <- ""
if(converged) {
# check.gradient
if(!is.null(GRADIENT) &&
OPTIMIZER %in% c("NLMINB", "BFGS", "L-BFGS-B")) {
# compute unscaled gradient
dx <- GRADIENT(x)
# NOTE: unscaled gradient!!!
if(converged && lavoptions$check.gradient &&
any(abs(dx) > lavoptions$optim.dx.tol)) {
# ok, identify the non-zero elements
non.zero <- which(abs(dx) > lavoptions$optim.dx.tol)
# which ones are 'boundary' points, defined by lower/upper?
bound.idx <- integer(0L)
if(!is.null(lavpartable$lower)) {
bound.idx <- c(bound.idx, which(lower == x))
}
if(!is.null(lavpartable$upper)) {
bound.idx <- c(bound.idx, which(upper == x))
}
if(length(bound.idx) > 0L) {
non.zero <- non.zero[- which(non.zero %in% bound.idx) ]
}
# this has many implications ... so should be careful to
# avoid false alarm
if(length(non.zero) > 0L) {
converged <- FALSE
warn.txt <- paste("the optimizer (", OPTIMIZER, ") ",
"claimed the model converged,\n",
" but not all elements of the gradient are (near) zero;\n",
" the optimizer may not have found a local solution\n",
" use check.gradient = FALSE to skip this check.",
sep = "")
}
}
} else {
dx <- numeric(0L)
}
} else {
dx <- numeric(0L)
warn.txt <- "the optimizer warns that a solution has NOT been found!"
}
# transform back
# 3.
#if(lavoptions$optim.var.transform == "sqrt" &&
# length(lavmodel@x.free.var.idx) > 0L) {
# #x[lavmodel@x.free.var.idx] <- tan(x[lavmodel@x.free.var.idx])
# x.var <- x[lavmodel@x.free.var.idx]
# x.var.sign <- sign(x.var)
# x[lavmodel@x.free.var.idx] <- x.var.sign * (x.var * x.var) # square!
#}
# 2. unpack
if(lavmodel@eq.constraints) {
x <- as.numeric(lavmodel@eq.constraints.K %*% x) +
lavmodel@eq.constraints.k0
}
# 1. unscale
x <- x / parscale
attr(x, "converged") <- converged
attr(x, "warn.txt") <- warn.txt
attr(x, "iterations") <- iterations
attr(x, "control") <- control
attr(x, "fx") <- fx
attr(x, "dx") <- dx
attr(x, "parscale") <- parscale
if(!is.null(optim.out$con.jac)) attr(x, "con.jac") <- optim.out$con.jac
if(!is.null(optim.out$lambda)) attr(x, "con.lambda") <- optim.out$lambda
if(lavoptions$optim.partrace) {
attr(x, "partrace") <- PENV$PARTRACE
}
x
}
# backwards compatibility
# estimateModel <- lav_model_estimate
Try the lavaan package in your browser
Any scripts or data that you put into this service are public.
lavaan documentation built on July 26, 2023, 5:08 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions lav_model_estimate
# model estimation
lav_model_estimate <- function(lavmodel = NULL,
lavpartable = NULL, # for parscale = "stand"
lavh1 = NULL, # for multilevel + parsc
lavsamplestats = NULL,
lavdata = NULL,
lavoptions = NULL,
lavcache = list(),
start = "model",
do.fit = TRUE) {
estimator <- lavoptions$estimator
verbose <- lavoptions$verbose
debug <- lavoptions$debug
ngroups <- lavsamplestats@ngroups
if(lavsamplestats@missing.flag || estimator == "PML") {
group.weight <- FALSE
} else {
group.weight <- TRUE
}
# backwards compatibility < 0.6-11
if(is.null(lavoptions$optim.partrace)) {
lavoptions$optim.partrace <- FALSE
}
if(lavoptions$optim.partrace) {
# fx + parameter values
PENV <- new.env()
PENV$PARTRACE <- matrix(NA, nrow=0, ncol=lavmodel@nx.free + 1L)
}
# starting values (ignoring equality constraints)
x.unpack <- lav_model_get_parameters(lavmodel)
# override? use simple instead? (new in 0.6-7)
if(start == "simple") {
START <- numeric(length(lavpartable$lhs))
# set loadings to 0.7
loadings.idx <- which(lavpartable$free > 0L &
lavpartable$op == "=~")
if(length(loadings.idx) > 0L) {
START[loadings.idx] <- 0.7
}
# set (only) variances to 1
var.idx <- which(lavpartable$free > 0L &
lavpartable$op == "~~" &
lavpartable$lhs == lavpartable$rhs)
if(length(var.idx) > 0L) {
START[var.idx] <- 1
}
if(lavmodel@ceq.simple.only) {
x.unpack <- START[ lavpartable$free > 0L &
!duplicated(lavpartable$free) ]
} else {
x.unpack <- START[ lavpartable$free > 0L ]
}
}
# 1. parameter scaling (to handle data scaling, not parameter scaling)
parscale <- rep(1.0, length(x.unpack))
# for < 0.6 compatibility
if(is.null(lavoptions$optim.parscale)) {
lavoptions$optim.parscale <- "none"
}
if(lavoptions$optim.parscale == "none") {
# do nothing, but still set SCALE, as before
} else if(lavoptions$optim.parscale %in% c("stand", "st", "standardize",
"standarized", "stand.all")) {
# rescale parameters as if the data was standardized
# new in 0.6-2
#
# FIXME: this works well, as long as the variances of the
# latent variables (which we do not know) are more or less
# equal to 1.0 (eg std.lv = TRUE)
#
# Once we have better estimates of those variances, we could
# use them to set the scale
#
if(lavdata@nlevels > 1L) {
if(length(lavh1) > 0L) {
OV.VAR <- lapply(lavh1$implied$cov, diag)
} else {
OV.VAR <- lapply(do.call(c, lapply(lavdata@Lp, "[[", "ov.idx")),
function(x) rep(1, length(x) ))
}
} else {
if(lavoptions$conditional.x) {
OV.VAR <- lavsamplestats@res.var
} else {
OV.VAR <- lavsamplestats@var
}
}
if(lavoptions$std.lv) {
parscale <- lav_standardize_all(lavobject = NULL,
est = rep(1, length(lavpartable$lhs)),
est.std = rep(1, length(lavpartable$lhs)),
cov.std = FALSE, ov.var = OV.VAR,
lavmodel = lavmodel, lavpartable = lavpartable,
cov.x = lavsamplestats@cov.x)
} else {
# needs good estimates for lv variances!
# if there is a single 'marker' indicator, we could use
# its observed variance as an upper bound
# for the moment, set them to 1.0 (instead of 0.05)
# TODO: USE Bentler's 1982 approach to get an estimate of
# VETA; use those diagonal elements...
# but only if we have 'marker' indicators for each LV
LV.VAR <- vector("list", lavmodel@ngroups)
for(g in seq_len(lavmodel@ngroups)) {
mm.in.group <- 1:lavmodel@nmat[g] + cumsum(c(0,lavmodel@nmat))[g]
MLIST <- lavmodel@GLIST[ mm.in.group ]
LAMBDA <- MLIST$lambda
n.lv <- ncol(LAMBDA)
LV.VAR[[g]] <- rep(1.0, n.lv)
}
parscale <- lav_standardize_all(lavobject = NULL,
est = rep(1, length(lavpartable$lhs)),
#est.std = rep(1, length(lavpartable$lhs)),
# here, we use whatever the starting values are
# for the latent variances...
cov.std = FALSE, ov.var = OV.VAR,
lv.var = LV.VAR,
lavmodel = lavmodel, lavpartable = lavpartable,
cov.x = lavsamplestats@cov.x)
}
# in addition, take sqrt for variance parameters
var.idx <- which(lavpartable$op == "~~" &
lavpartable$lhs == lavpartable$rhs)
if(length(var.idx) > 0L) {
parscale[var.idx] <- sqrt(abs(parscale[var.idx]))
}
if(lavmodel@ceq.simple.only) {
parscale <- parscale[ lavpartable$free > 0 &
!duplicated(lavpartable$free) ]
} else {
parscale <- parscale[ lavpartable$free > 0 ]
}
}
# parscale should obey the equality constraints
if(lavmodel@eq.constraints && lavoptions$optim.parscale != "none") {
# pack
p.pack <- as.numeric( (parscale - lavmodel@eq.constraints.k0) %*%
lavmodel@eq.constraints.K )
# unpack
parscale <- as.numeric(lavmodel@eq.constraints.K %*% p.pack) +
lavmodel@eq.constraints.k0
}
if(debug) {
cat("parscale = ", parscale, "\n")
}
z.unpack <- x.unpack * parscale
# 2. pack (apply equality constraints)
if(lavmodel@eq.constraints) {
z.pack <- as.numeric( (z.unpack - lavmodel@eq.constraints.k0) %*%
lavmodel@eq.constraints.K )
} else {
z.pack <- z.unpack
}
# 3. transform (already constrained) variances to standard deviations?
# TODO
#if(lavoptions$optim.var.transform == "sqrt" &&
# length(lavmodel@x.free.var.idx) > 0L) {
# # transforming variances using atan (or another sigmoid function?)
# # FIXME: better approach?
# #start.x[lavmodel@x.free.var.idx] <-
# # atan(start.x[lavmodel@x.free.var.idx])
# start.x[lavmodel@x.free.var.idx] <-
# sqrt(start.x[lavmodel@x.free.var.idx]) # assuming positive var
#}
# final starting values for optimizer
start.x <- z.pack
if(debug) {
cat("start.x = ", start.x, "\n")
}
# user-specified bounds? (new in 0.6-2)
if(is.null(lavpartable$lower)) {
lower <- -Inf
} else {
if(lavmodel@ceq.simple.only) {
free.idx <- which(lavpartable$free > 0L &
!duplicated(lavpartable$free))
lower <- lavpartable$lower[free.idx]
} else if(lavmodel@eq.constraints) {
# bounds have no effect any longer....
warning("lavaan warning: bounds have no effect in the presence of linear equality constraints")
lower <- -Inf
} else {
lower <- lavpartable$lower[ lavpartable$free > 0L ]
}
}
if(is.null(lavpartable$upper)) {
upper <- +Inf
} else {
if(lavmodel@ceq.simple.only) {
free.idx <- which(lavpartable$free > 0L &
!duplicated(lavpartable$free))
upper <- lavpartable$upper[free.idx]
} else if(lavmodel@eq.constraints) {
# bounds have no effect any longer....
if(is.null(lavpartable$lower)) {
# bounds have no effect any longer....
warning("lavaan warning: bounds have no effect in the presence of linear equality constraints")
}
upper <- +Inf
} else {
upper <- lavpartable$upper[ lavpartable$free > 0L ]
}
}
# check for inconsistent lower/upper bounds
# this may happen if we have equality constraints; qr() may switch
# the sign...
bad.idx <- which(lower > upper)
if(length(bad.idx) > 0L) {
# switch
#tmp <- lower[bad.idx]
#lower[bad.idx] <- upper[bad.idx]
#upper[bad.idx] <- tmp
lower[bad.idx] <- -Inf
upper[bad.idx] <- +Inf
}
# function to be minimized
objective_function <- function(x, verbose = FALSE, infToMax = FALSE,
debug = FALSE) {
# 3. standard deviations to variances
# WARNING: x is still packed here!
#if(lavoptions$optim.var.transform == "sqrt" &&
# length(lavmodel@x.free.var.idx) > 0L) {
# #x[lavmodel@x.free.var.idx] <- tan(x[lavmodel@x.free.var.idx])
# x.var <- x[lavmodel@x.free.var.idx]
# x.var.sign <- sign(x.var)
# x[lavmodel@x.free.var.idx] <- x.var.sign * (x.var * x.var) # square!
#}
# 2. unpack
if(lavmodel@eq.constraints) {
x <- as.numeric(lavmodel@eq.constraints.K %*% x) +
lavmodel@eq.constraints.k0
}
# 1. unscale
x <- x / parscale
# update GLIST (change `state') and make a COPY!
GLIST <- lav_model_x2GLIST(lavmodel, x = x)
fx <- lav_model_objective(lavmodel = lavmodel,
GLIST = GLIST,
lavsamplestats = lavsamplestats,
lavdata = lavdata,
lavcache = lavcache,
verbose = verbose)
# only for PML: divide by N (to speed up convergence)
if(estimator == "PML") {
fx <- fx / lavsamplestats@ntotal
}
if(debug || verbose) {
cat(" objective function = ",
sprintf("%18.16f", fx), "\n", sep="")
}
if(debug) {
#cat("Current unconstrained parameter values =\n")
#tmp.x <- lav_model_get_parameters(lavmodel, GLIST=GLIST, type="unco")
#print(tmp.x); cat("\n")
cat("Current free parameter values =\n"); print(x); cat("\n")
}
if(lavoptions$optim.partrace) {
PENV$PARTRACE <- rbind(PENV$PARTRACE, c(fx, x))
}
# for L-BFGS-B
#if(infToMax && is.infinite(fx)) fx <- 1e20
if(!is.finite(fx)) {
fx.group <- attr(fx, "fx.group")
fx <- 1e20
attr(fx, "fx.group") <- fx.group # only for lav_model_fit()
}
fx
}
gradient_function <- function(x, verbose = FALSE, infToMax = FALSE,
debug = FALSE) {
# transform variances back
#if(lavoptions$optim.var.transform == "sqrt" &&
# length(lavmodel@x.free.var.idx) > 0L) {
# #x[lavmodel@x.free.var.idx] <- tan(x[lavmodel@x.free.var.idx])
# x.var <- x[lavmodel@x.free.var.idx]
# x.var.sign <- sign(x.var)
# x[lavmodel@x.free.var.idx] <- x.var.sign * (x.var * x.var) # square!
#}
# 2. unpack
if(lavmodel@eq.constraints) {
x <- as.numeric(lavmodel@eq.constraints.K %*% x) +
lavmodel@eq.constraints.k0
}
# 1. unscale
x <- x / parscale
# update GLIST (change `state') and make a COPY!
GLIST <- lav_model_x2GLIST(lavmodel, x = x)
dx <- lav_model_gradient(lavmodel = lavmodel,
GLIST = GLIST,
lavsamplestats = lavsamplestats,
lavdata = lavdata,
lavcache = lavcache,
type = "free",
group.weight = group.weight, ### check me!!
verbose = verbose,
ceq.simple = lavmodel@ceq.simple.only)
if(debug) {
cat("Gradient function (analytical) =\n"); print(dx); cat("\n")
}
# 1. scale (note: divide, not multiply!)
dx <- dx / parscale
# 2. pack
if(lavmodel@eq.constraints) {
dx <- as.numeric( dx %*% lavmodel@eq.constraints.K )
}
# 3. transform variances back
#if(lavoptions$optim.var.transform == "sqrt" &&
# length(lavmodel@x.free.var.idx) > 0L) {
# x.var <- x[lavmodel@x.free.var.idx] # here in 'var' metric
# x.var.sign <- sign(x.var)
# x.var <- abs(x.var)
# x.sd <- sqrt(x.var)
# dx[lavmodel@x.free.var.idx] <-
# ( 2 * x.var.sign * dx[lavmodel@x.free.var.idx] * x.sd )
#}
# only for PML: divide by N (to speed up convergence)
if(estimator == "PML") {
dx <- dx / lavsamplestats@ntotal
}
if(debug) {
cat("Gradient function (analytical, after eq.constraints.K) =\n"); print(dx); cat("\n")
}
dx
}
gradient_function_numerical <- function(x, verbose=FALSE, debug = FALSE) {
# NOTE: no need to 'tranform' anything here (var/eq)
# this is done anyway in objective_function
# numerical approximation using the Richardson method
npar <- length(x)
h <- 10e-6
dx <- numeric( npar )
## FIXME: call lav_model_objective directly!!
for(i in 1:npar) {
x.left <- x.left2 <- x.right <- x.right2 <- x
x.left[i] <- x[i] - h; x.left2[i] <- x[i] - 2*h
x.right[i] <- x[i] + h; x.right2[i] <- x[i] + 2*h
fx.left <- objective_function(x.left, verbose = FALSE, debug = FALSE)
fx.left2 <- objective_function(x.left2, verbose = FALSE, debug = FALSE)
fx.right <- objective_function(x.right, verbose = FALSE, debug = FALSE)
fx.right2 <- objective_function(x.right2, verbose = FALSE, debug = FALSE)
dx[i] <- (fx.left2 - 8*fx.left + 8*fx.right - fx.right2)/(12*h)
}
#dx <- lavGradientC(func=objective_function, x=x)
# does not work if pnorm is involved... (eg PML)
if(debug) {
cat("Gradient function (numerical) =\n"); print(dx); cat("\n")
}
dx
}
gradient_function_numerical_complex <- function(x, verbose=FALSE, debug = FALSE) {
dx <- Re(lav_func_gradient_complex(func = objective_function, x = x,
h = sqrt(.Machine$double.eps)))
# does not work if pnorm is involved... (eg PML)
if(debug) {
cat("Gradient function (numerical complex) =\n"); print(dx); cat("\n")
}
dx
}
# check if the initial values produce a positive definite Sigma
# to begin with -- but only for estimator="ML"
if(estimator %in% c("ML","FML","MML")) {
Sigma.hat <- computeSigmaHat(lavmodel, extra=TRUE, debug=lavoptions$debug)
for(g in 1:ngroups) {
if(!attr(Sigma.hat[[g]], "po")) {
group.txt <- ifelse(ngroups > 1,
paste(" in group ",g,".",sep=""), ".")
if(debug) {
print(Sigma.hat[[g]][,])
}
warning("lavaan ERROR: initial model-implied matrix (Sigma) is not positive definite;\n check your model and/or starting parameters", group.txt)
x <- start.x
fx <- as.numeric(NA)
attr(fx, "fx.group") <- rep(as.numeric(NA), ngroups)
attr(x, "converged") <- FALSE
attr(x, "iterations") <- 0L
attr(x, "control") <- lavoptions@control
attr(x, "fx") <- fx
return(x)
}
}
}
# parameter scaling
# FIXME: what is the best way to set the scale??
# current strategy: if startx > 1.0, we rescale by using
# 1/startx
SCALE <- rep(1.0, length(start.x))
if(lavoptions$optim.parscale == "none") {
idx <- which(abs(start.x) > 1.0)
if(length(idx) > 0L) {
SCALE[idx] <- abs(1.0/start.x[idx])
}
}
if(debug) {
cat("SCALE = ", SCALE, "\n")
}
# first try: check if starting values return a finite value
fx <- objective_function(start.x, verbose = verbose, debug = debug)
if(!is.finite(fx)) {
# emergency change of start.x
start.x <- start.x / 10
}
# first some nelder mead steps? (default = FALSE)
INIT_NELDER_MEAD <- lavoptions$optim.init_nelder_mead
# gradient: analytic, numerical or NULL?
if(is.character(lavoptions$optim.gradient)) {
if(lavoptions$optim.gradient %in% c("analytic","analytical")) {
GRADIENT <- gradient_function
} else if(lavoptions$optim.gradient %in% c("numerical", "numeric")) {
GRADIENT <- gradient_function_numerical
} else if(lavoptions$optim.gradient %in% c("numeric.complex", "complex")){
GRADIENT <- gradient_function_numerical_complex
} else if(lavoptions$optim.gradient %in% c("NULL", "null")) {
GRADIENT <- NULL
} else {
warning("lavaan WARNING: gradient should be analytic, numerical or NULL")
}
} else if(is.logical(lavoptions$optim.gradient)) {
if(lavoptions$optim.gradient) {
GRADIENT <- gradient_function
} else {
GRADIENT <- NULL
}
} else if(is.null(lavoptions$optim.gradient)) {
GRADIENT <- gradient_function
}
# default optimizer
if(length(lavmodel@ceq.nonlinear.idx) == 0L &&
length(lavmodel@cin.linear.idx) == 0L &&
length(lavmodel@cin.nonlinear.idx) == 0L) {
if(is.null(lavoptions$optim.method)) {
OPTIMIZER <- "NLMINB"
#OPTIMIZER <- "BFGS" # slightly slower, no bounds; better scaling!
#OPTIMIZER <- "L-BFGS-B" # trouble with Inf values for fx!
} else {
OPTIMIZER <- toupper(lavoptions$optim.method)
stopifnot(OPTIMIZER %in% c("NLMINB0", "NLMINB1", "NLMINB2",
"NLMINB", "BFGS", "L-BFGS-B", "NONE"))
if(OPTIMIZER == "NLMINB1") {
OPTIMIZER <- "NLMINB"
}
}
} else {
if(is.null(lavoptions$optim.method)) {
OPTIMIZER <- "NLMINB.CONSTR"
} else {
OPTIMIZER <- toupper(lavoptions$optim.method)
stopifnot(OPTIMIZER %in% c("NLMINB.CONSTR", "NLMINB", "NONE"))
}
if(OPTIMIZER == "NLMINB") {
OPTIMIZER <- "NLMINB.CONSTR"
}
}
if(INIT_NELDER_MEAD) {
if(verbose) cat(" initial Nelder-Mead step:\n")
trace <- 0L; if(verbose) trace <- 1L
optim.out <- optim(par=start.x,
fn=objective_function,
method="Nelder-Mead",
#control=list(maxit=10L,
# parscale=SCALE,
# trace=trace),
hessian=FALSE,
verbose=verbose, debug=debug)
cat("\n")
start.x <- optim.out$par
}
if(OPTIMIZER == "NLMINB0") {
if(verbose) cat(" quasi-Newton steps using NLMINB0 (no analytic gradient):\n")
#if(debug) control$trace <- 1L;
control.nlminb <- list(eval.max=20000L,
iter.max=10000L,
trace=0L,
#abs.tol=1e-20, ### important!! fx never negative
abs.tol=(.Machine$double.eps * 10),
rel.tol=1e-10,
#step.min=2.2e-14, # in =< 0.5-12
step.min=1.0, # 1.0 in < 0.5-21
step.max=1.0,
x.tol=1.5e-8,
xf.tol=2.2e-14)
control.nlminb <- modifyList(control.nlminb, lavoptions$control)
control <- control.nlminb[c("eval.max", "iter.max", "trace",
"step.min", "step.max",
"abs.tol", "rel.tol", "x.tol", "xf.tol")]
#cat("DEBUG: control = "); print(str(control.nlminb)); cat("\n")
optim.out <- nlminb(start=start.x,
objective=objective_function,
gradient=NULL,
lower=lower,
upper=upper,
control=control,
scale=SCALE,
verbose=verbose, debug=debug)
if(verbose) {
cat(" convergence status (0=ok): ", optim.out$convergence, "\n")
cat(" nlminb message says: ", optim.out$message, "\n")
cat(" number of iterations: ", optim.out$iterations, "\n")
cat(" number of function evaluations [objective, gradient]: ",
optim.out$evaluations, "\n")
}
# try again
if(optim.out$convergence != 0L) {
optim.out <- nlminb(start=start.x,
objective=objective_function,
gradient=NULL,
lower=lower,
upper=upper,
control=control,
scale=SCALE,
verbose=verbose, debug=debug)
}
iterations <- optim.out$iterations
x <- optim.out$par
if(optim.out$convergence == 0L) {
converged <- TRUE
} else {
converged <- FALSE
}
} else if(OPTIMIZER == "NLMINB") {
if(verbose) cat(" quasi-Newton steps using NLMINB:\n")
#if(debug) control$trace <- 1L;
control.nlminb <- list(eval.max=20000L,
iter.max=10000L,
trace=0L,
#abs.tol=1e-20, ### important!! fx never negative
abs.tol=(.Machine$double.eps * 10),
rel.tol=1e-10,
#step.min=2.2e-14, # in =< 0.5-12
step.min=1.0, # 1.0 in < 0.5-21
step.max=1.0,
x.tol=1.5e-8,
xf.tol=2.2e-14)
control.nlminb <- modifyList(control.nlminb, lavoptions$control)
control <- control.nlminb[c("eval.max", "iter.max", "trace",
"step.min", "step.max",
"abs.tol", "rel.tol", "x.tol", "xf.tol")]
#cat("DEBUG: control = "); print(str(control.nlminb)); cat("\n")
optim.out <- nlminb(start=start.x,
objective=objective_function,
gradient=GRADIENT,
lower=lower,
upper=upper,
control=control,
scale=SCALE,
verbose=verbose, debug=debug)
if(verbose) {
cat(" convergence status (0=ok): ", optim.out$convergence, "\n")
cat(" nlminb message says: ", optim.out$message, "\n")
cat(" number of iterations: ", optim.out$iterations, "\n")
cat(" number of function evaluations [objective, gradient]: ",
optim.out$evaluations, "\n")
}
iterations <- optim.out$iterations
x <- optim.out$par
if(optim.out$convergence == 0L) {
converged <- TRUE
} else {
converged <- FALSE
}
} else if(OPTIMIZER == "BFGS") {
# warning: Bollen example with estimator=GLS does NOT converge!
# (but WLS works!)
# - BB.ML works too
control.bfgs <- list(trace=0L, fnscale=1,
parscale=SCALE, ## or not?
ndeps=1e-3,
maxit=10000,
abstol=1e-20,
reltol=1e-10,
REPORT=1L)
control.bfgs <- modifyList(control.bfgs, lavoptions$control)
control <- control.bfgs[c("trace", "fnscale", "parscale", "ndeps",
"maxit", "abstol", "reltol", "REPORT")]
#trace <- 0L; if(verbose) trace <- 1L
optim.out <- optim(par=start.x,
fn=objective_function,
gr=GRADIENT,
method="BFGS",
control=control,
hessian=FALSE,
verbose=verbose, debug=debug)
if(verbose) {
cat(" convergence status (0=ok): ", optim.out$convergence, "\n")
cat(" optim BFGS message says: ", optim.out$message, "\n")
#cat("number of iterations: ", optim.out$iterations, "\n")
cat(" number of function evaluations [objective, gradient]: ",
optim.out$counts, "\n")
}
#iterations <- optim.out$iterations
iterations <- optim.out$counts[1]
x <- optim.out$par
if(optim.out$convergence == 0L) {
converged <- TRUE
} else {
converged <- FALSE
}
} else if(OPTIMIZER == "L-BFGS-B") {
# warning, does not cope with Inf values!!
control.lbfgsb <- list(trace=0L, fnscale=1,
parscale=SCALE, ## or not?
ndeps=1e-3,
maxit=10000,
REPORT=1L,
lmm=5L,
factr=1e7,
pgtol=0)
control.lbfgsb <- modifyList(control.lbfgsb, lavoptions$control)
control <- control.lbfgsb[c("trace", "fnscale", "parscale",
"ndeps", "maxit", "REPORT", "lmm",
"factr", "pgtol")]
optim.out <- optim(par=start.x,
fn=objective_function,
gr=GRADIENT,
method="L-BFGS-B",
lower=lower,
upper=upper,
control=control,
hessian=FALSE,
verbose=verbose, debug=debug,
infToMax=TRUE)
if(verbose) {
cat(" convergence status (0=ok): ", optim.out$convergence, "\n")
cat(" optim L-BFGS-B message says: ", optim.out$message, "\n")
#cat("number of iterations: ", optim.out$iterations, "\n")
cat(" number of function evaluations [objective, gradient]: ",
optim.out$counts, "\n")
}
#iterations <- optim.out$iterations
iterations <- optim.out$counts[1]
x <- optim.out$par
if(optim.out$convergence == 0L) {
converged <- TRUE
} else {
converged <- FALSE
}
} else if(OPTIMIZER == "NLMINB.CONSTR") {
ocontrol <- list(verbose=verbose)
if(!is.null(lavoptions$control$control.outer)) {
ocontrol <- c(lavoptions$control$control.outer, verbose=verbose)
}
control.nlminb <- list(eval.max=20000L,
iter.max=10000L,
trace=0L,
#abs.tol=1e-20,
abs.tol=(.Machine$double.eps * 10),
rel.tol=1e-9, # 1e-10 seems 'too strict'
step.min=1.0, # 1.0 in < 0.5-21
step.max=1.0,
x.tol=1.5e-8,
xf.tol=2.2e-14)
control.nlminb <- modifyList(control.nlminb, lavoptions$control)
control <- control.nlminb[c("eval.max", "iter.max", "trace",
"abs.tol", "rel.tol")]
cin <- cin.jac <- ceq <- ceq.jac <- NULL
if(!is.null(body(lavmodel@cin.function))) cin <- lavmodel@cin.function
if(!is.null(body(lavmodel@cin.jacobian))) cin.jac <- lavmodel@cin.jacobian
if(!is.null(body(lavmodel@ceq.function))) ceq <- lavmodel@ceq.function
if(!is.null(body(lavmodel@ceq.jacobian))) ceq.jac <- lavmodel@ceq.jacobian
trace <- FALSE; if(verbose) trace <- TRUE
optim.out <- nlminb.constr(start = start.x,
objective=objective_function,
gradient=GRADIENT,
control=control,
scale=SCALE,
verbose=verbose, debug=debug,
lower=lower,
upper=upper,
cin = cin, cin.jac = cin.jac,
ceq = ceq, ceq.jac = ceq.jac,
control.outer = ocontrol
)
if(verbose) {
cat(" convergence status (0=ok): ", optim.out$convergence, "\n")
cat(" nlminb.constr message says: ", optim.out$message, "\n")
cat(" number of outer iterations: ", optim.out$outer.iterations, "\n")
cat(" number of inner iterations: ", optim.out$iterations, "\n")
cat(" number of function evaluations [objective, gradient]: ",
optim.out$evaluations, "\n")
}
iterations <- optim.out$iterations
x <- optim.out$par
if(optim.out$convergence == 0) {
converged <- TRUE
} else {
converged <- FALSE
}
} else if(OPTIMIZER == "NONE") {
x <- start.x
iterations <- 0L
converged <- TRUE
control <- list()
# if inequality constraints, add con.jac/lambda
# needed for df!
if(length(lavmodel@ceq.nonlinear.idx) == 0L &&
length(lavmodel@cin.linear.idx) == 0L &&
length(lavmodel@cin.nonlinear.idx) == 0L) {
optim.out <- list()
} else {
# if inequality constraints, add con.jac/lambda
# needed for df!
optim.out <- list()
if(is.null(body(lavmodel@ceq.function))) {
ceq <- function(x, ...) { return( numeric(0) ) }
} else {
ceq <- lavmodel@ceq.function
}
if(is.null(body(lavmodel@cin.function))) {
cin <- function(x, ...) { return( numeric(0) ) }
} else {
cin <- lavmodel@cin.function
}
ceq0 <- ceq(start.x)
cin0 <- cin(start.x)
con0 <- c(ceq0, cin0)
JAC <- rbind(numDeriv::jacobian(ceq, x = start.x),
numDeriv::jacobian(cin, x = start.x))
nceq <- length( ceq(start.x) )
ncin <- length( cin(start.x) )
ncon <- nceq + ncin
ceq.idx <- cin.idx <- integer(0)
if(nceq > 0L) ceq.idx <- 1:nceq
if(ncin > 0L) cin.idx <- nceq + 1:ncin
cin.flag <- rep(FALSE, length(ncon))
if(ncin > 0L) cin.flag[cin.idx] <- TRUE
inactive.idx <- integer(0L)
cin.idx <- which(cin.flag)
if(ncin > 0L) {
slack <- 1e-05
inactive.idx <- which(cin.flag & con0 > slack)
}
attr(JAC, "inactive.idx") <- inactive.idx
attr(JAC, "cin.idx") <- cin.idx
attr(JAC, "ceq.idx") <- ceq.idx
optim.out$con.jac <- JAC
optim.out$lambda <- rep(0, ncon)
}
}
fx <- objective_function(x) # to get "fx.group" attribute
# check convergence
warn.txt <- ""
if(converged) {
# check.gradient
if(!is.null(GRADIENT) &&
OPTIMIZER %in% c("NLMINB", "BFGS", "L-BFGS-B")) {
# compute unscaled gradient
dx <- GRADIENT(x)
# NOTE: unscaled gradient!!!
if(converged && lavoptions$check.gradient &&
any(abs(dx) > lavoptions$optim.dx.tol)) {
# ok, identify the non-zero elements
non.zero <- which(abs(dx) > lavoptions$optim.dx.tol)
# which ones are 'boundary' points, defined by lower/upper?
bound.idx <- integer(0L)
if(!is.null(lavpartable$lower)) {
bound.idx <- c(bound.idx, which(lower == x))
}
if(!is.null(lavpartable$upper)) {
bound.idx <- c(bound.idx, which(upper == x))
}
if(length(bound.idx) > 0L) {
non.zero <- non.zero[- which(non.zero %in% bound.idx) ]
}
# this has many implications ... so should be careful to
# avoid false alarm
if(length(non.zero) > 0L) {
converged <- FALSE
warn.txt <- paste("the optimizer (", OPTIMIZER, ") ",
"claimed the model converged,\n",
" but not all elements of the gradient are (near) zero;\n",
" the optimizer may not have found a local solution\n",
" use check.gradient = FALSE to skip this check.",
sep = "")
}
}
} else {
dx <- numeric(0L)
}
} else {
dx <- numeric(0L)
warn.txt <- "the optimizer warns that a solution has NOT been found!"
}
# transform back
# 3.
#if(lavoptions$optim.var.transform == "sqrt" &&
# length(lavmodel@x.free.var.idx) > 0L) {
# #x[lavmodel@x.free.var.idx] <- tan(x[lavmodel@x.free.var.idx])
# x.var <- x[lavmodel@x.free.var.idx]
# x.var.sign <- sign(x.var)
# x[lavmodel@x.free.var.idx] <- x.var.sign * (x.var * x.var) # square!
#}
# 2. unpack
if(lavmodel@eq.constraints) {
x <- as.numeric(lavmodel@eq.constraints.K %*% x) +
lavmodel@eq.constraints.k0
}
# 1. unscale
x <- x / parscale
attr(x, "converged") <- converged
attr(x, "warn.txt") <- warn.txt
attr(x, "iterations") <- iterations
attr(x, "control") <- control
attr(x, "fx") <- fx
attr(x, "dx") <- dx
attr(x, "parscale") <- parscale
if(!is.null(optim.out$con.jac)) attr(x, "con.jac") <- optim.out$con.jac
if(!is.null(optim.out$lambda)) attr(x, "con.lambda") <- optim.out$lambda
if(lavoptions$optim.partrace) {
attr(x, "partrace") <- PENV$PARTRACE
}
x
}
# backwards compatibility
# estimateModel <- lav_model_estimate
Try the lavaan package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.