Nothing
R/lav_partable_bounds.R
In lavaan: Latent Variable Analysis
Defines functions
lav_partable_add_bounds
lav_partable_add_bounds <- function(partable = NULL,
lavpta = NULL,
lavh1 = NULL,
lavdata = NULL,
lavsamplestats = NULL,
lavoptions = NULL) {
# no support (yet) for multilevel
if(lav_partable_nlevels(partable) > 1L) {
return(partable)
}
# check optim.bounds
if(is.null(lavoptions$optim.bounds)) {
# <0.6-6 version
return(partable)
} else {
if(!is.null(lavoptions$bounds) && lavoptions$bounds == "none") {
# no bounds needed
return(partable)
}
# no support from effect.coding (for now)
if(!is.null(lavoptions$effect.coding) &&
nchar(lavoptions$effect.coding[1L]) > 0L) {
warning("lavaan WARNING: automatic bounds not available (yet) if effect.coding is used")
return(partable)
}
optim.bounds <- lavoptions$optim.bounds
# check the elements
if(is.null(optim.bounds$lower)) {
optim.bounds$lower <- character(0L)
} else {
optim.bounds$lower <- as.character(optim.bounds$lower)
}
if(is.null(optim.bounds$upper)) {
optim.bounds$upper <- character(0L)
} else {
optim.bounds$upper <- as.character(optim.bounds$upper)
}
if(is.null(optim.bounds$min.reliability.marker)) {
optim.bounds$min.reliability.marker <- 0.0
} else {
if(optim.bounds$min.reliability.marker < 0 ||
optim.bounds$min.reliability.marker > 1.0) {
stop("lavaan ERROR: optim.bounds$min.reliability.marker ",
"is out of range: ", optim.bounds$min.reliability.marker)
}
}
if(is.null(optim.bounds$min.var.ov)) {
optim.bounds$min.var.ov <- -Inf
}
if(is.null(optim.bounds$min.var.lv.exo)) {
optim.bounds$min.var.lv.exo <- 0.0
}
if(is.null(optim.bounds$min.var.lv.endo)) {
optim.bounds$min.var.lv.endo <- 0.0
}
if(is.null(optim.bounds$max.r2.lv.endo)) {
optim.bounds$max.r2.lv.endo <- 1.0
}
if(is.null(optim.bounds$lower.factor)) {
optim.bounds$lower.factor <- rep(1.0, length(optim.bounds$lower))
} else {
if(length(optim.bounds$lower.factor) == 1L &&
is.numeric(optim.bounds$lower.factor)) {
optim.bounds$lower.factor <- rep(optim.bounds$lower.factor,
length(optim.bounds$lower))
} else if(length(optim.bounds$lower.factor) !=
length(optim.bounds$lower)) {
stop("lavaan ERROR: length(optim.bounds$lower.factor) is not ",
"equal to length(optim.bounds$lower)")
}
}
lower.factor <- optim.bounds$lower.factor
if(is.null(optim.bounds$upper.factor)) {
optim.bounds$upper.factor <- rep(1.0, length(optim.bounds$upper))
} else {
if(length(optim.bounds$upper.factor) == 1L &&
is.numeric(optim.bounds$upper.factor)) {
optim.bounds$upper.factor <- rep(optim.bounds$upper.factor,
length(optim.bounds$upper))
} else if(length(optim.bounds$upper.factor) !=
length(optim.bounds$upper)) {
stop("lavaan ERROR: length(optim.bounds$upper.factor) is not ",
"equal to length(optim.bounds$upper)")
}
}
upper.factor <- optim.bounds$upper.factor
}
# shortcut
REL <- optim.bounds$min.reliability.marker
# nothing to do
if(length(optim.bounds$lower) == 0L &&
length(optim.bounds$upper) == 0L) {
return(partable)
} else {
# we compute ALL bounds, then we select what we need
# (otherwise, we can not use the 'factor')
if(!is.null(partable$lower)) {
lower.user <- partable$lower
} else {
partable$lower <- lower.user <- rep(-Inf, length(partable$lhs))
}
if(!is.null(partable$upper)) {
upper.user <- partable$upper
} else {
partable$upper <- upper.user <- rep(+Inf, length(partable$lhs))
}
# the 'automatic' bounds
lower.auto <- rep(-Inf, length(partable$lhs))
upper.auto <- rep(+Inf, length(partable$lhs))
}
# make sure we have lavpta
if(is.null(lavpta)) {
lavpta <- lav_partable_attributes(partable)
}
# check blocks
if(is.null(partable$block)) {
partable$block <- rep(1L, length(partable$lhs))
}
block.values <- lav_partable_block_values(partable)
# check groups
if(is.null(partable$group)) {
partable$group <- rep(1L, length(partable$lhs))
}
group.values <- lav_partable_group_values(partable)
ngroups <- length(group.values)
# compute bounds per group ### TODO: add levels/classes/...
b <- 0L
for(g in seq_len(ngroups)) {
# next block
b <- b + 1L
# for this block
ov.names <- lavpta$vnames$ov[[b]]
lv.names <- lavpta$vnames$lv[[b]]
lv.names.x <- lavpta$vnames$lv.x[[b]]
if(length(lv.names.x) > 0L) {
lv.names.endo <- lv.names[! lv.names %in% lv.names.x ]
} else {
lv.names.endo <- lv.names
}
lv.marker <- lavpta$vnames$lv.marker[[b]]
# OV.VAR for this group
if(lavsamplestats@missing.flag && lavdata@nlevels == 1L) {
OV.VAR <- diag(lavsamplestats@missing.h1[[g]]$sigma)
} else {
if(lavoptions$conditional.x) {
OV.VAR <- diag(lavsamplestats@res.cov[[g]])
} else {
OV.VAR <- diag(lavsamplestats@cov[[g]])
}
}
# we 'process' the parameters per 'type', so we can choose
# to apply (or not) upper/lower bounds for each type separately
################################
## 1. (residual) ov variances ##
################################
par.idx <- which(partable$group == group.values[g] &
partable$op == "~~" &
partable$lhs %in% ov.names &
partable$lhs == partable$rhs)
if(length(par.idx) > 0L) {
# lower == 0
lower.auto[par.idx] <- 0
# upper == var(ov)
var.idx <- match(partable$lhs[par.idx], ov.names)
upper.auto[par.idx] <- OV.VAR[var.idx]
# if reliability > 0, adapt marker indicators only
if(REL > 0) {
marker.idx <- which(partable$group == group.values[g] &
partable$op == "~~" &
partable$lhs %in% lv.marker &
partable$lhs == partable$rhs)
marker.var.idx <- match(partable$lhs[marker.idx], ov.names)
# upper = (1-REL)*OVAR
upper.auto[marker.idx] <- (1 - REL) * OV.VAR[marker.var.idx]
}
# range
bound.range <- upper.auto[par.idx] - pmax(lower.auto[par.idx], 0)
# enlarge lower?
if("ov.var" %in% optim.bounds$lower) {
factor <- lower.factor[ which(optim.bounds$lower == "ov.var") ]
if( is.finite(factor) && factor != 1.0 ) {
new.range <- bound.range * factor
diff <- abs(new.range - bound.range)
lower.auto[par.idx] <- lower.auto[par.idx] - diff
}
}
# enlarge upper?
if("ov.var" %in% optim.bounds$upper) {
factor <- upper.factor[ which(optim.bounds$upper == "ov.var") ]
if( is.finite(factor) && factor != 1.0 ) {
new.range <- bound.range * factor
diff <- abs(new.range - bound.range)
upper.auto[par.idx] <- upper.auto[par.idx] + diff
}
}
# min.var.ov?
min.idx <- which(lower.auto[par.idx] < optim.bounds$min.var.ov)
if(length(min.idx) > 0L) {
lower.auto[par.idx[min.idx]] <- optim.bounds$min.var.ov
}
# requested?
if("ov.var" %in% optim.bounds$lower) {
partable$lower[par.idx] <- lower.auto[par.idx]
}
if("ov.var" %in% optim.bounds$upper) {
partable$upper[par.idx] <- upper.auto[par.idx]
}
} # (res) ov variances
################################
## 2. (residual) lv variances ##
################################
# first collect lower/upper bounds for TOTAL variances in lv.names
LV.VAR.LB <- numeric( length(lv.names) )
LV.VAR.UB <- numeric( length(lv.names) )
if(lavoptions$std.lv) {
LV.VAR.LB <- rep(1.0, length(lv.names))
LV.VAR.UB <- rep(1.0, length(lv.names))
} else {
for(i in seq_len(length(lv.names))) {
this.lv.name <- lv.names[i]
this.lv.marker <- lv.marker[i]
if(nchar(this.lv.marker) > 0L && this.lv.marker %in% ov.names) {
marker.var <- OV.VAR[ match(this.lv.marker, ov.names) ]
LOWER <- marker.var - (1 - REL)*marker.var
LV.VAR.LB[i] <- max(LOWER, optim.bounds$min.var.lv.exo)
#LV.VAR.UB[i] <- marker.var - REL*marker.var
LV.VAR.UB[i] <- marker.var
} else {
LV.VAR.LB[i] <- optim.bounds$min.var.lv.exo
LV.VAR.UB[i] <- max(OV.VAR)
}
}
}
# use these bounds for the free parameters
par.idx <- which(partable$group == group.values[g] &
partable$op == "~~" &
partable$lhs %in% lv.names &
partable$lhs == partable$rhs)
if(length(par.idx) > 0L) {
# adjust for endogenenous lv
LV.VAR.LB2 <- LV.VAR.LB
endo.idx <- which(lv.names %in% lv.names.endo)
if(length(endo.idx) > 0L) {
LV.VAR.LB2[endo.idx] <- optim.bounds$min.var.lv.endo
if(optim.bounds$max.r2.lv.endo != 1) {
LV.VAR.LB2[endo.idx] <- (1 - optim.bounds$max.r2.lv.endo) * LV.VAR.UB[endo.idx]
}
}
exo.idx <- which(!lv.names %in% lv.names.endo)
if(length(exo.idx) > 0L && optim.bounds$min.var.lv.exo != 0) {
LV.VAR.LB2[exo.idx] <- optim.bounds$min.var.lv.exo
}
lower.auto[par.idx] <- LV.VAR.LB2[ match(partable$lhs[par.idx],
lv.names) ]
upper.auto[par.idx] <- LV.VAR.UB[ match(partable$lhs[par.idx],
lv.names) ]
# range
bound.range <- upper.auto[par.idx] - pmax(lower.auto[par.idx], 0)
# enlarge lower?
if("lv.var" %in% optim.bounds$lower) {
factor <- lower.factor[ which(optim.bounds$lower == "lv.var") ]
if( is.finite(factor) && factor != 1.0 ) {
new.range <- bound.range * factor
diff <- abs(new.range - bound.range)
lower.auto[par.idx] <- lower.auto[par.idx] - diff
}
}
# enlarge upper?
if("lv.var" %in% optim.bounds$upper) {
factor <- upper.factor[ which(optim.bounds$upper == "lv.var") ]
if( is.finite(factor) && factor != 1.0 ) {
new.range <- bound.range * factor
diff <- abs(new.range - bound.range)
upper.auto[par.idx] <- upper.auto[par.idx] + diff
}
}
# requested?
if("lv.var" %in% optim.bounds$lower) {
partable$lower[par.idx] <- lower.auto[par.idx]
}
if("lv.var" %in% optim.bounds$upper) {
partable$upper[par.idx] <- upper.auto[par.idx]
}
} # lv variances
#############################################
## 3. factor loadings (ov indicators only) ##
#############################################
# lambda_p^(u) = sqrt( upper(res.var.indicators_p) /
# lower(var.factor) )
ov.ind.names <- lavpta$vnames$ov.ind[[b]]
par.idx <- which(partable$group == group.values[g] &
partable$op == "=~" &
partable$lhs %in% lv.names &
partable$rhs %in% ov.ind.names)
if(length(par.idx) > 0L) {
# if negative LV variances are allowed (due to factor > 1)
# make them equal to zero
LV.VAR.LB[ LV.VAR.LB < 0 ] <- 0.0
var.all <- OV.VAR[ match(partable$rhs[par.idx], ov.names) ]
tmp <- LV.VAR.LB[ match(partable$lhs[par.idx], lv.names) ]
tmp[is.na(tmp)] <- 0 # just in case...
lower.auto[par.idx] <- -1 * sqrt(var.all/tmp) # -Inf if tmp==0
upper.auto[par.idx] <- +1 * sqrt(var.all/tmp) # +Inf if tmp==0
# if std.lv = TRUE, force 'first' loading to be positive?
#if(lavoptions$std.lv) {
# # get index 'first' indicators
# first.idx <- which(!duplicated(partable$lhs[par.idx]))
# lower.auto[par.idx][first.idx] <- 0
#}
# range
bound.range <- upper.auto[par.idx] - lower.auto[par.idx]
# enlarge lower?
if("loadings" %in% optim.bounds$lower) {
factor <- lower.factor[ which(optim.bounds$lower=="loadings") ]
if( is.finite(factor) && factor != 1.0 ) {
new.range <- bound.range * factor
ok.idx <- is.finite(new.range)
if(length(ok.idx) > 0L) {
diff <- abs(new.range[ok.idx] - bound.range[ok.idx])
lower.auto[par.idx][ok.idx] <-
lower.auto[par.idx][ok.idx] - diff
}
}
}
# enlarge upper?
if("loadings" %in% optim.bounds$upper) {
factor <- upper.factor[ which(optim.bounds$upper=="loadings") ]
if( is.finite(factor) && factor != 1.0 ) {
new.range <- bound.range * factor
ok.idx <- is.finite(new.range)
if(length(ok.idx) > 0L) {
diff <- abs(new.range[ok.idx] - bound.range[ok.idx])
upper.auto[par.idx][ok.idx] <-
upper.auto[par.idx][ok.idx] + diff
}
}
}
# requested?
if("loadings" %in% optim.bounds$lower) {
partable$lower[par.idx] <- lower.auto[par.idx]
}
if("loadings" %in% optim.bounds$upper) {
partable$upper[par.idx] <- upper.auto[par.idx]
}
} # lambda
####################
## 4. covariances ##
####################
# | sqrt(var(x)) sqrt(var(y)) | <= cov(x,y)
par.idx <- which(partable$group == group.values[g] &
partable$op == "~~" &
partable$lhs != partable$rhs)
if(length(par.idx) > 0L) {
for(i in seq_len(length(par.idx))) {
# this lhs/rhs
this.lhs <- partable$lhs[ par.idx[i] ]
this.rhs <- partable$rhs[ par.idx[i] ]
# 2 possibilities:
# - variances are free parameters
# - variances are fixed (eg std.lv = TRUE)
# var idx
lhs.var.idx <- which(partable$group == group.values[g] &
partable$op == "~~" &
partable$lhs == this.lhs &
partable$lhs == partable$rhs)
rhs.var.idx <- which(partable$group == group.values[g] &
partable$op == "~~" &
partable$lhs == this.rhs &
partable$lhs == partable$rhs)
# upper bounds
lhs.upper <- upper.auto[lhs.var.idx]
rhs.upper <- upper.auto[rhs.var.idx]
# compute upper bounds for this cov (assuming >0 vars)
if(is.finite(lhs.upper) && is.finite(rhs.upper)) {
upper.cov <- sqrt(lhs.upper) * sqrt(rhs.upper)
upper.auto[par.idx[i]] <- +1 * upper.cov
lower.auto[par.idx[i]] <- -1 * upper.cov
}
}
# range
bound.range <- upper.auto[par.idx] - lower.auto[par.idx]
# enlarge lower?
if("covariances" %in% optim.bounds$lower) {
factor <-
lower.factor[ which(optim.bounds$lower=="covariances") ]
if( is.finite(factor) && factor != 1.0 ) {
new.range <- bound.range * factor
ok.idx <- is.finite(new.range)
if(length(ok.idx) > 0L) {
diff <- new.range[ok.idx] - bound.range[ok.idx]
lower.auto[par.idx][ok.idx] <-
lower.auto[par.idx][ok.idx] - diff
}
}
}
# enlarge upper?
if("covariances" %in% optim.bounds$upper) {
factor <-
upper.factor[ which(optim.bounds$upper=="covariances") ]
if( is.finite(factor) && factor != 1.0 ) {
new.range <- bound.range * factor
ok.idx <- is.finite(new.range)
if(length(ok.idx) > 0L) {
diff <- new.range[ok.idx] - bound.range[ok.idx]
upper.auto[par.idx][ok.idx] <-
upper.auto[par.idx][ok.idx] + diff
}
}
}
# requested?
if("covariances" %in% optim.bounds$lower) {
partable$lower[par.idx] <- lower.auto[par.idx]
}
if("covariances" %in% optim.bounds$upper) {
partable$upper[par.idx] <- upper.auto[par.idx]
}
} # covariances
} # g
# overwrite with lower.user (except -Inf)
not.inf.idx <- which(lower.user > -Inf)
if(length(not.inf.idx) > 0L) {
partable$lower[not.inf.idx] <- lower.user[not.inf.idx]
}
# overwrite with upper.user (except +Inf)
not.inf.idx <- which(upper.user < +Inf)
if(length(not.inf.idx) > 0L) {
partable$upper[not.inf.idx] <- upper.user[not.inf.idx]
}
# non-free
non.free.idx <- which(partable$free == 0L)
if(length(non.free.idx) > 0L && !is.null(partable$ustart)) {
partable$lower[non.free.idx] <- partable$ustart[non.free.idx]
partable$upper[non.free.idx] <- partable$ustart[non.free.idx]
}
partable
}
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_partable_add_bounds
lav_partable_add_bounds <- function(partable = NULL,
lavpta = NULL,
lavh1 = NULL,
lavdata = NULL,
lavsamplestats = NULL,
lavoptions = NULL) {
# no support (yet) for multilevel
if(lav_partable_nlevels(partable) > 1L) {
return(partable)
}
# check optim.bounds
if(is.null(lavoptions$optim.bounds)) {
# <0.6-6 version
return(partable)
} else {
if(!is.null(lavoptions$bounds) && lavoptions$bounds == "none") {
# no bounds needed
return(partable)
}
# no support from effect.coding (for now)
if(!is.null(lavoptions$effect.coding) &&
nchar(lavoptions$effect.coding[1L]) > 0L) {
warning("lavaan WARNING: automatic bounds not available (yet) if effect.coding is used")
return(partable)
}
optim.bounds <- lavoptions$optim.bounds
# check the elements
if(is.null(optim.bounds$lower)) {
optim.bounds$lower <- character(0L)
} else {
optim.bounds$lower <- as.character(optim.bounds$lower)
}
if(is.null(optim.bounds$upper)) {
optim.bounds$upper <- character(0L)
} else {
optim.bounds$upper <- as.character(optim.bounds$upper)
}
if(is.null(optim.bounds$min.reliability.marker)) {
optim.bounds$min.reliability.marker <- 0.0
} else {
if(optim.bounds$min.reliability.marker < 0 ||
optim.bounds$min.reliability.marker > 1.0) {
stop("lavaan ERROR: optim.bounds$min.reliability.marker ",
"is out of range: ", optim.bounds$min.reliability.marker)
}
}
if(is.null(optim.bounds$min.var.ov)) {
optim.bounds$min.var.ov <- -Inf
}
if(is.null(optim.bounds$min.var.lv.exo)) {
optim.bounds$min.var.lv.exo <- 0.0
}
if(is.null(optim.bounds$min.var.lv.endo)) {
optim.bounds$min.var.lv.endo <- 0.0
}
if(is.null(optim.bounds$max.r2.lv.endo)) {
optim.bounds$max.r2.lv.endo <- 1.0
}
if(is.null(optim.bounds$lower.factor)) {
optim.bounds$lower.factor <- rep(1.0, length(optim.bounds$lower))
} else {
if(length(optim.bounds$lower.factor) == 1L &&
is.numeric(optim.bounds$lower.factor)) {
optim.bounds$lower.factor <- rep(optim.bounds$lower.factor,
length(optim.bounds$lower))
} else if(length(optim.bounds$lower.factor) !=
length(optim.bounds$lower)) {
stop("lavaan ERROR: length(optim.bounds$lower.factor) is not ",
"equal to length(optim.bounds$lower)")
}
}
lower.factor <- optim.bounds$lower.factor
if(is.null(optim.bounds$upper.factor)) {
optim.bounds$upper.factor <- rep(1.0, length(optim.bounds$upper))
} else {
if(length(optim.bounds$upper.factor) == 1L &&
is.numeric(optim.bounds$upper.factor)) {
optim.bounds$upper.factor <- rep(optim.bounds$upper.factor,
length(optim.bounds$upper))
} else if(length(optim.bounds$upper.factor) !=
length(optim.bounds$upper)) {
stop("lavaan ERROR: length(optim.bounds$upper.factor) is not ",
"equal to length(optim.bounds$upper)")
}
}
upper.factor <- optim.bounds$upper.factor
}
# shortcut
REL <- optim.bounds$min.reliability.marker
# nothing to do
if(length(optim.bounds$lower) == 0L &&
length(optim.bounds$upper) == 0L) {
return(partable)
} else {
# we compute ALL bounds, then we select what we need
# (otherwise, we can not use the 'factor')
if(!is.null(partable$lower)) {
lower.user <- partable$lower
} else {
partable$lower <- lower.user <- rep(-Inf, length(partable$lhs))
}
if(!is.null(partable$upper)) {
upper.user <- partable$upper
} else {
partable$upper <- upper.user <- rep(+Inf, length(partable$lhs))
}
# the 'automatic' bounds
lower.auto <- rep(-Inf, length(partable$lhs))
upper.auto <- rep(+Inf, length(partable$lhs))
}
# make sure we have lavpta
if(is.null(lavpta)) {
lavpta <- lav_partable_attributes(partable)
}
# check blocks
if(is.null(partable$block)) {
partable$block <- rep(1L, length(partable$lhs))
}
block.values <- lav_partable_block_values(partable)
# check groups
if(is.null(partable$group)) {
partable$group <- rep(1L, length(partable$lhs))
}
group.values <- lav_partable_group_values(partable)
ngroups <- length(group.values)
# compute bounds per group ### TODO: add levels/classes/...
b <- 0L
for(g in seq_len(ngroups)) {
# next block
b <- b + 1L
# for this block
ov.names <- lavpta$vnames$ov[[b]]
lv.names <- lavpta$vnames$lv[[b]]
lv.names.x <- lavpta$vnames$lv.x[[b]]
if(length(lv.names.x) > 0L) {
lv.names.endo <- lv.names[! lv.names %in% lv.names.x ]
} else {
lv.names.endo <- lv.names
}
lv.marker <- lavpta$vnames$lv.marker[[b]]
# OV.VAR for this group
if(lavsamplestats@missing.flag && lavdata@nlevels == 1L) {
OV.VAR <- diag(lavsamplestats@missing.h1[[g]]$sigma)
} else {
if(lavoptions$conditional.x) {
OV.VAR <- diag(lavsamplestats@res.cov[[g]])
} else {
OV.VAR <- diag(lavsamplestats@cov[[g]])
}
}
# we 'process' the parameters per 'type', so we can choose
# to apply (or not) upper/lower bounds for each type separately
################################
## 1. (residual) ov variances ##
################################
par.idx <- which(partable$group == group.values[g] &
partable$op == "~~" &
partable$lhs %in% ov.names &
partable$lhs == partable$rhs)
if(length(par.idx) > 0L) {
# lower == 0
lower.auto[par.idx] <- 0
# upper == var(ov)
var.idx <- match(partable$lhs[par.idx], ov.names)
upper.auto[par.idx] <- OV.VAR[var.idx]
# if reliability > 0, adapt marker indicators only
if(REL > 0) {
marker.idx <- which(partable$group == group.values[g] &
partable$op == "~~" &
partable$lhs %in% lv.marker &
partable$lhs == partable$rhs)
marker.var.idx <- match(partable$lhs[marker.idx], ov.names)
# upper = (1-REL)*OVAR
upper.auto[marker.idx] <- (1 - REL) * OV.VAR[marker.var.idx]
}
# range
bound.range <- upper.auto[par.idx] - pmax(lower.auto[par.idx], 0)
# enlarge lower?
if("ov.var" %in% optim.bounds$lower) {
factor <- lower.factor[ which(optim.bounds$lower == "ov.var") ]
if( is.finite(factor) && factor != 1.0 ) {
new.range <- bound.range * factor
diff <- abs(new.range - bound.range)
lower.auto[par.idx] <- lower.auto[par.idx] - diff
}
}
# enlarge upper?
if("ov.var" %in% optim.bounds$upper) {
factor <- upper.factor[ which(optim.bounds$upper == "ov.var") ]
if( is.finite(factor) && factor != 1.0 ) {
new.range <- bound.range * factor
diff <- abs(new.range - bound.range)
upper.auto[par.idx] <- upper.auto[par.idx] + diff
}
}
# min.var.ov?
min.idx <- which(lower.auto[par.idx] < optim.bounds$min.var.ov)
if(length(min.idx) > 0L) {
lower.auto[par.idx[min.idx]] <- optim.bounds$min.var.ov
}
# requested?
if("ov.var" %in% optim.bounds$lower) {
partable$lower[par.idx] <- lower.auto[par.idx]
}
if("ov.var" %in% optim.bounds$upper) {
partable$upper[par.idx] <- upper.auto[par.idx]
}
} # (res) ov variances
################################
## 2. (residual) lv variances ##
################################
# first collect lower/upper bounds for TOTAL variances in lv.names
LV.VAR.LB <- numeric( length(lv.names) )
LV.VAR.UB <- numeric( length(lv.names) )
if(lavoptions$std.lv) {
LV.VAR.LB <- rep(1.0, length(lv.names))
LV.VAR.UB <- rep(1.0, length(lv.names))
} else {
for(i in seq_len(length(lv.names))) {
this.lv.name <- lv.names[i]
this.lv.marker <- lv.marker[i]
if(nchar(this.lv.marker) > 0L && this.lv.marker %in% ov.names) {
marker.var <- OV.VAR[ match(this.lv.marker, ov.names) ]
LOWER <- marker.var - (1 - REL)*marker.var
LV.VAR.LB[i] <- max(LOWER, optim.bounds$min.var.lv.exo)
#LV.VAR.UB[i] <- marker.var - REL*marker.var
LV.VAR.UB[i] <- marker.var
} else {
LV.VAR.LB[i] <- optim.bounds$min.var.lv.exo
LV.VAR.UB[i] <- max(OV.VAR)
}
}
}
# use these bounds for the free parameters
par.idx <- which(partable$group == group.values[g] &
partable$op == "~~" &
partable$lhs %in% lv.names &
partable$lhs == partable$rhs)
if(length(par.idx) > 0L) {
# adjust for endogenenous lv
LV.VAR.LB2 <- LV.VAR.LB
endo.idx <- which(lv.names %in% lv.names.endo)
if(length(endo.idx) > 0L) {
LV.VAR.LB2[endo.idx] <- optim.bounds$min.var.lv.endo
if(optim.bounds$max.r2.lv.endo != 1) {
LV.VAR.LB2[endo.idx] <- (1 - optim.bounds$max.r2.lv.endo) * LV.VAR.UB[endo.idx]
}
}
exo.idx <- which(!lv.names %in% lv.names.endo)
if(length(exo.idx) > 0L && optim.bounds$min.var.lv.exo != 0) {
LV.VAR.LB2[exo.idx] <- optim.bounds$min.var.lv.exo
}
lower.auto[par.idx] <- LV.VAR.LB2[ match(partable$lhs[par.idx],
lv.names) ]
upper.auto[par.idx] <- LV.VAR.UB[ match(partable$lhs[par.idx],
lv.names) ]
# range
bound.range <- upper.auto[par.idx] - pmax(lower.auto[par.idx], 0)
# enlarge lower?
if("lv.var" %in% optim.bounds$lower) {
factor <- lower.factor[ which(optim.bounds$lower == "lv.var") ]
if( is.finite(factor) && factor != 1.0 ) {
new.range <- bound.range * factor
diff <- abs(new.range - bound.range)
lower.auto[par.idx] <- lower.auto[par.idx] - diff
}
}
# enlarge upper?
if("lv.var" %in% optim.bounds$upper) {
factor <- upper.factor[ which(optim.bounds$upper == "lv.var") ]
if( is.finite(factor) && factor != 1.0 ) {
new.range <- bound.range * factor
diff <- abs(new.range - bound.range)
upper.auto[par.idx] <- upper.auto[par.idx] + diff
}
}
# requested?
if("lv.var" %in% optim.bounds$lower) {
partable$lower[par.idx] <- lower.auto[par.idx]
}
if("lv.var" %in% optim.bounds$upper) {
partable$upper[par.idx] <- upper.auto[par.idx]
}
} # lv variances
#############################################
## 3. factor loadings (ov indicators only) ##
#############################################
# lambda_p^(u) = sqrt( upper(res.var.indicators_p) /
# lower(var.factor) )
ov.ind.names <- lavpta$vnames$ov.ind[[b]]
par.idx <- which(partable$group == group.values[g] &
partable$op == "=~" &
partable$lhs %in% lv.names &
partable$rhs %in% ov.ind.names)
if(length(par.idx) > 0L) {
# if negative LV variances are allowed (due to factor > 1)
# make them equal to zero
LV.VAR.LB[ LV.VAR.LB < 0 ] <- 0.0
var.all <- OV.VAR[ match(partable$rhs[par.idx], ov.names) ]
tmp <- LV.VAR.LB[ match(partable$lhs[par.idx], lv.names) ]
tmp[is.na(tmp)] <- 0 # just in case...
lower.auto[par.idx] <- -1 * sqrt(var.all/tmp) # -Inf if tmp==0
upper.auto[par.idx] <- +1 * sqrt(var.all/tmp) # +Inf if tmp==0
# if std.lv = TRUE, force 'first' loading to be positive?
#if(lavoptions$std.lv) {
# # get index 'first' indicators
# first.idx <- which(!duplicated(partable$lhs[par.idx]))
# lower.auto[par.idx][first.idx] <- 0
#}
# range
bound.range <- upper.auto[par.idx] - lower.auto[par.idx]
# enlarge lower?
if("loadings" %in% optim.bounds$lower) {
factor <- lower.factor[ which(optim.bounds$lower=="loadings") ]
if( is.finite(factor) && factor != 1.0 ) {
new.range <- bound.range * factor
ok.idx <- is.finite(new.range)
if(length(ok.idx) > 0L) {
diff <- abs(new.range[ok.idx] - bound.range[ok.idx])
lower.auto[par.idx][ok.idx] <-
lower.auto[par.idx][ok.idx] - diff
}
}
}
# enlarge upper?
if("loadings" %in% optim.bounds$upper) {
factor <- upper.factor[ which(optim.bounds$upper=="loadings") ]
if( is.finite(factor) && factor != 1.0 ) {
new.range <- bound.range * factor
ok.idx <- is.finite(new.range)
if(length(ok.idx) > 0L) {
diff <- abs(new.range[ok.idx] - bound.range[ok.idx])
upper.auto[par.idx][ok.idx] <-
upper.auto[par.idx][ok.idx] + diff
}
}
}
# requested?
if("loadings" %in% optim.bounds$lower) {
partable$lower[par.idx] <- lower.auto[par.idx]
}
if("loadings" %in% optim.bounds$upper) {
partable$upper[par.idx] <- upper.auto[par.idx]
}
} # lambda
####################
## 4. covariances ##
####################
# | sqrt(var(x)) sqrt(var(y)) | <= cov(x,y)
par.idx <- which(partable$group == group.values[g] &
partable$op == "~~" &
partable$lhs != partable$rhs)
if(length(par.idx) > 0L) {
for(i in seq_len(length(par.idx))) {
# this lhs/rhs
this.lhs <- partable$lhs[ par.idx[i] ]
this.rhs <- partable$rhs[ par.idx[i] ]
# 2 possibilities:
# - variances are free parameters
# - variances are fixed (eg std.lv = TRUE)
# var idx
lhs.var.idx <- which(partable$group == group.values[g] &
partable$op == "~~" &
partable$lhs == this.lhs &
partable$lhs == partable$rhs)
rhs.var.idx <- which(partable$group == group.values[g] &
partable$op == "~~" &
partable$lhs == this.rhs &
partable$lhs == partable$rhs)
# upper bounds
lhs.upper <- upper.auto[lhs.var.idx]
rhs.upper <- upper.auto[rhs.var.idx]
# compute upper bounds for this cov (assuming >0 vars)
if(is.finite(lhs.upper) && is.finite(rhs.upper)) {
upper.cov <- sqrt(lhs.upper) * sqrt(rhs.upper)
upper.auto[par.idx[i]] <- +1 * upper.cov
lower.auto[par.idx[i]] <- -1 * upper.cov
}
}
# range
bound.range <- upper.auto[par.idx] - lower.auto[par.idx]
# enlarge lower?
if("covariances" %in% optim.bounds$lower) {
factor <-
lower.factor[ which(optim.bounds$lower=="covariances") ]
if( is.finite(factor) && factor != 1.0 ) {
new.range <- bound.range * factor
ok.idx <- is.finite(new.range)
if(length(ok.idx) > 0L) {
diff <- new.range[ok.idx] - bound.range[ok.idx]
lower.auto[par.idx][ok.idx] <-
lower.auto[par.idx][ok.idx] - diff
}
}
}
# enlarge upper?
if("covariances" %in% optim.bounds$upper) {
factor <-
upper.factor[ which(optim.bounds$upper=="covariances") ]
if( is.finite(factor) && factor != 1.0 ) {
new.range <- bound.range * factor
ok.idx <- is.finite(new.range)
if(length(ok.idx) > 0L) {
diff <- new.range[ok.idx] - bound.range[ok.idx]
upper.auto[par.idx][ok.idx] <-
upper.auto[par.idx][ok.idx] + diff
}
}
}
# requested?
if("covariances" %in% optim.bounds$lower) {
partable$lower[par.idx] <- lower.auto[par.idx]
}
if("covariances" %in% optim.bounds$upper) {
partable$upper[par.idx] <- upper.auto[par.idx]
}
} # covariances
} # g
# overwrite with lower.user (except -Inf)
not.inf.idx <- which(lower.user > -Inf)
if(length(not.inf.idx) > 0L) {
partable$lower[not.inf.idx] <- lower.user[not.inf.idx]
}
# overwrite with upper.user (except +Inf)
not.inf.idx <- which(upper.user < +Inf)
if(length(not.inf.idx) > 0L) {
partable$upper[not.inf.idx] <- upper.user[not.inf.idx]
}
# non-free
non.free.idx <- which(partable$free == 0L)
if(length(non.free.idx) > 0L && !is.null(partable$ustart)) {
partable$lower[non.free.idx] <- partable$ustart[non.free.idx]
partable$upper[non.free.idx] <- partable$ustart[non.free.idx]
}
partable
}
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.