Nothing
R/lav_standardize.R
In lavaan: Latent Variable Analysis
Defines functions
lav_unstandardize_ov
lav_standardize_all_nox
lav_standardize_all
lav_standardize_lv
lav_unstandardize_ov_x
lav_standardize_all_nox_x
lav_standardize_all_x
lav_standardize_lv_x
lav_standardize_lv_x <- function(x, lavobject, partable = NULL, cov.std = TRUE,
lv.var = NULL,
rotation = FALSE) {
# set new values for x
lavmodel <- lav_model_set_parameters(lavmodel = lavobject@Model, x = x)
if (rotation) {
x.unrotated <- x
lavmodel@GLIST <- lavTech(lavobject, "est.unrotated") # unrotated!
est.rot <- lav_model_efa_rotate_x(
x = x.unrotated,
lavmodel = lavmodel, # unrotated!
lavoptions = lavobject@Options,
init.rot = lavmodel@H,
type = "user",
extra = TRUE
)
GLIST <- attr(est.rot, "extra")$GLIST
attributes(est.rot) <- NULL
est <- est.rot
} else {
GLIST <- lavmodel@GLIST # if this changes, tag @TDJorgensen in commit message
est <- lav_model_get_parameters(lavmodel, type = "user")
}
x.stand.user <- lav_standardize_lv(
lavobject = lavobject,
partable = partable, est = est,
GLIST = GLIST, cov.std = cov.std,
lv.var = lv.var
)
x.stand.user
}
lav_standardize_all_x <- function(x, lavobject, partable = NULL, cov.std = TRUE,
rotation = FALSE) {
lavmodel <- lav_model_set_parameters(lavmodel = lavobject@Model, x = x)
if (rotation) {
x.unrotated <- x
lavmodel@GLIST <- lavTech(lavobject, "est.unrotated") # unrotated!
est.rot <- lav_model_efa_rotate_x(
x = x.unrotated,
lavmodel = lavmodel, # unrotated!
lavoptions = lavobject@Options,
init.rot = lavmodel@H,
type = "user",
extra = TRUE
)
GLIST <- attr(est.rot, "extra")$GLIST
attributes(est.rot) <- NULL
est <- est.rot
} else {
GLIST <- lavmodel@GLIST # if this changes, tag @TDJorgensen in commit message
est <- lav_model_get_parameters(lavmodel, type = "user")
}
x.stand.user <- lav_standardize_all(
lavobject = lavobject,
partable = partable, est = est,
est.std = NULL, GLIST = GLIST,
cov.std = cov.std
)
x.stand.user
}
lav_standardize_all_nox_x <- function(x, lavobject, partable = NULL,
cov.std = TRUE, rotation = FALSE) {
lavmodel <- lav_model_set_parameters(lavmodel = lavobject@Model, x = x)
if (rotation) {
x.unrotated <- x
lavmodel@GLIST <- lavTech(lavobject, "est.unrotated") # unrotated!
est.rot <- lav_model_efa_rotate_x(
x = x.unrotated,
lavmodel = lavmodel, # unrotated!
lavoptions = lavobject@Options,
init.rot = lavmodel@H,
type = "user",
extra = TRUE
)
GLIST <- attr(est.rot, "extra")$GLIST
attributes(est.rot) <- NULL
est <- est.rot
} else {
GLIST <- lavmodel@GLIST # if this changes, tag @TDJorgensen in commit message
est <- lav_model_get_parameters(lavmodel, type = "user")
}
x.stand.user <- lav_standardize_all_nox(
lavobject = lavobject,
partable = partable, est = est,
est.std = NULL, GLIST = GLIST,
cov.std = cov.std
)
x.stand.user
}
lav_unstandardize_ov_x <- function(x, lavobject) {
partable <- lavobject@ParTable
partable$ustart <- x
lav_unstandardize_ov(
partable = partable,
ov.var = lavobject@SampleStats@var,
cov.std = TRUE
)
}
lav_standardize_lv <- function(lavobject = NULL,
partable = NULL, est = NULL, GLIST = NULL,
cov.std = TRUE, lv.var = NULL,
lavmodel = NULL, lavpartable = NULL) {
if (is.null(lavobject)) {
stopifnot(!is.null(lavmodel))
stopifnot(!is.null(lavpartable))
if (is.null(est)) {
if (!is.null(lavpartable$est)) {
est <- lavpartable$est # if this changes, tag @TDJorgensen in commit message
} else {
lav_msg_stop(gettext("could not find `est' in lavpartable"))
}
}
} else {
lavmodel <- lavobject@Model
lavpartable <- lavobject@ParTable
if (is.null(est)) {
est <- lav_object_inspect_est(lavobject)
}
}
if (is.null(partable)) {
partable <- lavpartable
}
if (is.null(GLIST)) {
GLIST <- lavmodel@GLIST
}
out <- est
N <- length(est)
stopifnot(N == length(partable$lhs))
nmat <- lavmodel@nmat
# compute ETA
if (is.null(lv.var)) {
LV.ETA <- computeVETA(
lavmodel = lavmodel,
GLIST = GLIST
)
}
for (g in 1:lavmodel@nblocks) {
ov.names <- vnames(lavpartable, "ov", block = g) # not user,
# which may be incomplete
lv.names <- vnames(lavpartable, "lv", block = g)
# shortcut: no latents in this block, nothing to do
if (length(lv.names) == 0L) {
next
}
# which mm belong to block g?
mm.in.group <- 1:nmat[g] + cumsum(c(0, nmat))[g]
MLIST <- GLIST[mm.in.group]
if (is.null(lv.var)) {
ETA2 <- diag(LV.ETA[[g]])
} else {
ETA2 <- lv.var[[g]]
}
# change negative values to NA
ETA2[ETA2 < 0] <- as.numeric(NA)
ETA <- sqrt(ETA2)
# 1a. "=~" regular indicators
idx <- which(partable$op == "=~" & !(partable$rhs %in% lv.names) &
partable$block == g)
out[idx] <- out[idx] * ETA[match(partable$lhs[idx], lv.names)]
# 1b. "=~" regular higher-order lv indicators
idx <- which(partable$op == "=~" & !(partable$rhs %in% ov.names) &
partable$block == g)
out[idx] <- (out[idx] * ETA[match(partable$lhs[idx], lv.names)]
/ ETA[match(partable$rhs[idx], lv.names)])
# 1c. "=~" indicators that are both in ov and lv
# idx <- which(partable$op == "=~" & partable$rhs %in% ov.names
# & partable$rhs %in% lv.names &
# partable$block == g)
# 2. "~" regressions (and "<~")
idx <- which((partable$op == "~" | partable$op == "<~") &
partable$lhs %in% lv.names &
partable$block == g)
out[idx] <- out[idx] / ETA[match(partable$lhs[idx], lv.names)]
idx <- which((partable$op == "~" | partable$op == "<~") &
partable$rhs %in% lv.names &
partable$block == g)
out[idx] <- out[idx] * ETA[match(partable$rhs[idx], lv.names)]
# 3a. "~~" ov
# idx <- which(partable$op == "~~" & !(partable$lhs %in% lv.names) &
# partable$block == g)
# 3b. "~~" lv
# ATTENTION: in Mplus 4.1, the off-diagonal residual covariances
# were computed by the formula cov(i,j) / sqrt(i.var*j.var)
# were i.var and j.var where diagonal elements of ETA
#
# in Mplus 6.1 (but also AMOS and EQS), the i.var and j.var
# elements are the 'PSI' diagonal elements!!
# variances
rv.idx <- which(partable$op == "~~" & partable$rhs %in% lv.names &
partable$lhs == partable$rhs &
partable$block == g)
out[rv.idx] <- (out[rv.idx] / ETA[match(partable$lhs[rv.idx], lv.names)]
/ ETA[match(partable$rhs[rv.idx], lv.names)])
# covariances lv
# three types:
# - only lhs is LV (and fixed.x = FALSE)
# - only rhs is LV (and fixed.x = FALSE)
# - both lhs and rhs are LV (regular case)
if (cov.std) {
if (!is.complex(est[rv.idx])) {
RV <- sqrt(abs(est[rv.idx])) # abs in case of heywood cases
} else {
RV <- sqrt(est[rv.idx])
}
rv.names <- partable$lhs[rv.idx]
}
# left
idx.lhs <- which(partable$op == "~~" &
partable$lhs %in% lv.names &
partable$lhs != partable$rhs &
partable$block == g)
if (length(idx.lhs) > 0L) {
if (cov.std == FALSE) {
out[idx.lhs] <-
(out[idx.lhs] / ETA[match(partable$lhs[idx.lhs], lv.names)])
} else {
out[idx.lhs] <-
(out[idx.lhs] / RV[match(partable$lhs[idx.lhs], rv.names)])
}
}
# right
idx.rhs <- which(partable$op == "~~" &
partable$rhs %in% lv.names &
partable$lhs != partable$rhs &
partable$block == g)
if (length(idx.rhs) > 0L) {
if (cov.std == FALSE) {
out[idx.rhs] <-
(out[idx.rhs] / ETA[match(partable$rhs[idx.rhs], lv.names)])
} else {
out[idx.rhs] <-
(out[idx.rhs] / RV[match(partable$rhs[idx.rhs], rv.names)])
}
}
# 4a. "~1" ov
# idx <- which(partable$op == "~1" & !(partable$lhs %in% lv.names) &
# partable$block == g)
# 4b. "~1" lv
idx <- which(partable$op == "~1" & partable$lhs %in% lv.names &
partable$block == g)
out[idx] <- out[idx] / ETA[match(partable$lhs[idx], lv.names)]
}
# 5a ":="
idx <- which(partable$op == ":=")
if (length(idx) > 0L) {
x <- out[partable$free & !duplicated(partable$free)]
out[idx] <- lavmodel@def.function(x)
}
# 5b "=="
idx <- which(partable$op == "==")
if (length(idx) > 0L) {
x <- out[partable$free & !duplicated(partable$free)]
out[idx] <- lavmodel@ceq.function(x)
}
# 5c. "<" or ">"
idx <- which((partable$op == "<" | partable$op == ">"))
if (length(idx) > 0L) {
x <- out[partable$free & !duplicated(partable$free)]
out[idx] <- lavmodel@cin.function(x)
}
out
}
lav_standardize_all <- function(lavobject = NULL,
partable = NULL, est = NULL, est.std = NULL,
GLIST = NULL, cov.std = TRUE, ov.var = NULL,
lv.var = NULL,
lavmodel = NULL, lavpartable = NULL,
cov.x = NULL) {
if (is.null(lavobject)) {
stopifnot(!is.null(lavmodel))
stopifnot(!is.null(lavpartable))
if (is.null(est)) {
if (!is.null(lavpartable$est)) {
est <- lavpartable$est # if this changes, tag @TDJorgensen in commit message
} else {
lav_msg_stop(gettext("could not find `est' in lavpartable"))
}
}
} else {
lavmodel <- lavobject@Model
lavpartable <- lavobject@ParTable
if (is.null(est)) {
est <- lav_object_inspect_est(lavobject)
}
if (lavmodel@conditional.x) {
if (is.null(cov.x)) {
# try SampleStats slot
# if("SampleStats" %in% slotNames(lavobject)) {
# cov.x <- lavobject@SampleStats@cov.x
if (!is.null(lavobject@implied$cov.x[[1]])) {
cov.x <- lavobject@implied$cov.x # if this changes, tag @TDJorgensen in commit message
} else {
# perhaps lavaanList object
# extract it from GLIST per block
cov.x <- vector("list", length = lavmodel@nblocks)
for (b in seq_len(lavmodel@nblocks)) {
# which mm belong to block b?
mm.in.block <- (seq_len(lavmodel@nmat[b]) +
cumsum(c(0, lavmodel@nmat))[b])
MLIST <- lavmodel@GLIST[mm.in.block]
cov.x[[b]] <- MLIST[["cov.x"]]
}
}
}
}
}
if (is.null(partable)) {
partable <- lavpartable
}
if (is.null(GLIST)) {
GLIST <- lavmodel@GLIST
}
if (is.null(est.std)) {
est.std <- lav_standardize_lv(
lavobject = lavobject,
partable = partable, est = est, GLIST = GLIST,
cov.std = cov.std, lv.var = lv.var, lavmodel = lavmodel,
lavpartable = lavpartable
)
}
out <- est.std
N <- length(est.std)
stopifnot(N == length(partable$lhs))
VY <- computeVY(
lavmodel = lavmodel, GLIST = GLIST,
diagonal.only = TRUE
)
for (g in 1:lavmodel@nblocks) {
ov.names <- vnames(lavpartable, "ov", block = g) # not user
lv.names <- vnames(lavpartable, "lv", block = g)
if (is.null(ov.var)) {
OV2 <- VY[[g]]
# replace zero values by NA (but keep negative values)
zero.idx <- which(abs(OV2) < .Machine$double.eps)
if (length(zero.idx) > 0L) {
OV2[zero.idx] <- as.numeric(NA)
}
# replace negative values by NA (for sqrt)
tmp.OV2 <- OV2
neg.idx <- which(tmp.OV2 < 0)
if (length(neg.idx) > 0L) {
tmp.OV2[neg.idx] <- as.numeric(NA)
}
OV <- sqrt(tmp.OV2)
} else {
OV2 <- ov.var[[g]]
OV <- sqrt(OV2)
}
if (lavmodel@conditional.x) {
# extend OV with ov.names.x
ov.names.x <- vnames(lavpartable, "ov.x", block = g)
ov.names.nox <- vnames(lavpartable, "ov.nox", block = g)
ov.names <- c(ov.names.nox, ov.names.x)
OV2 <- c(OV2, diag(cov.x[[g]]))
OV <- c(OV, sqrt(diag(cov.x[[g]])))
}
# 1a. "=~" regular indicators
idx <- which(partable$op == "=~" & !(partable$rhs %in% lv.names) &
partable$block == g)
out[idx] <- out[idx] / OV[match(partable$rhs[idx], ov.names)]
# 1b. "=~" regular higher-order lv indicators
# 1c. "=~" indicators that are both in ov and lv
# idx <- which(partable$op == "=~" & partable$rhs %in% ov.names
# & partable$rhs %in% lv.names &
# partable$block == g)
# 2. "~" regressions (and "<~")
idx <- which((partable$op == "~" | partable$op == "<~") &
partable$lhs %in% ov.names &
partable$block == g)
out[idx] <- out[idx] / OV[match(partable$lhs[idx], ov.names)]
idx <- which((partable$op == "~" | partable$op == "<~") &
partable$rhs %in% ov.names &
partable$block == g)
out[idx] <- out[idx] * OV[match(partable$rhs[idx], ov.names)]
# 3a. "~~" ov
# ATTENTION: in Mplus 4.1, the off-diagonal residual covariances
# were computed by the formula cov(i,j) / sqrt(i.var*j.var)
# were i.var and j.var where diagonal elements of OV
#
# in Mplus 6.1 (but also AMOS and EQS), the i.var and j.var
# elements are the 'THETA' diagonal elements!!
# variances
rv.idx <- which(partable$op == "~~" & !(partable$lhs %in% lv.names) &
partable$lhs == partable$rhs &
partable$block == g)
# out[rv.idx] <- ( out[rv.idx] / OV[ match(partable$lhs[rv.idx], ov.names) ]
# / OV[ match(partable$rhs[rv.idx], ov.names) ] )
out[rv.idx] <- (out[rv.idx] /
OV2[match(partable$lhs[rv.idx], ov.names)])
# covariances ov
# three types:
# - only lhs is OV (and fixed.x = FALSE)
# - only rhs is OV (and fixed.x = FALSE)
# - both lhs and rhs are OV (regular case)
if (cov.std) {
if (!is.complex(est[rv.idx])) {
RV <- sqrt(abs(est[rv.idx]))
} else {
RV <- sqrt(est[rv.idx])
}
rv.names <- partable$lhs[rv.idx]
}
# left
idx.lhs <- which(partable$op == "~~" &
!(partable$lhs %in% lv.names) &
partable$lhs != partable$rhs &
partable$block == g)
if (length(idx.lhs) > 0L) {
if (cov.std == FALSE) {
out[idx.lhs] <-
(out[idx.lhs] / OV[match(partable$lhs[idx.lhs], ov.names)])
} else {
out[idx.lhs] <-
(out[idx.lhs] / RV[match(partable$lhs[idx.lhs], rv.names)])
}
}
# right
idx.rhs <- which(partable$op == "~~" &
!(partable$rhs %in% lv.names) &
partable$lhs != partable$rhs &
partable$block == g)
if (length(idx.rhs) > 0L) {
if (cov.std == FALSE) {
out[idx.rhs] <-
(out[idx.rhs] / OV[match(partable$rhs[idx.rhs], ov.names)])
} else {
out[idx.rhs] <-
(out[idx.rhs] / RV[match(partable$rhs[idx.rhs], rv.names)])
}
}
# 3b. "~~" lv
# idx <- which(partable$op == "~~" & partable$rhs %in% lv.names &
# partable$block == g)
# 4a. "~1" ov
idx <- which(partable$op == "~1" & !(partable$lhs %in% lv.names) &
partable$block == g)
out[idx] <- out[idx] / OV[match(partable$lhs[idx], ov.names)]
# 4b. "~1" lv
# idx <- which(partable$op == "~1" & partable$lhs %in% lv.names &
# partable$block == g)
# 4c. "|" thresholds
idx <- which(partable$op == "|" & !(partable$lhs %in% lv.names) &
partable$block == g)
out[idx] <- out[idx] / OV[match(partable$lhs[idx], ov.names)]
# 4d. "~*~" scales
idx <- which(partable$op == "~*~" & !(partable$lhs %in% lv.names) &
partable$block == g)
out[idx] <- 1.0
}
# 5a ":="
idx <- which(partable$op == ":=")
if (length(idx) > 0L) {
x <- out[partable$free & !duplicated(partable$free)]
out[idx] <- lavmodel@def.function(x)
}
# 5b "=="
idx <- which(partable$op == "==")
if (length(idx) > 0L) {
x <- out[partable$free & !duplicated(partable$free)]
out[idx] <- lavmodel@ceq.function(x)
}
# 5c. "<" or ">"
idx <- which((partable$op == "<" | partable$op == ">"))
if (length(idx) > 0L) {
x <- out[partable$free & !duplicated(partable$free)]
out[idx] <- lavmodel@cin.function(x)
}
out
}
lav_standardize_all_nox <- function(lavobject = NULL,
partable = NULL, est = NULL, est.std = NULL,
GLIST = NULL, cov.std = TRUE, ov.var = NULL,
lv.var = NULL,
lavmodel = NULL, lavpartable = NULL,
cov.x = NULL) {
if (is.null(lavobject)) {
stopifnot(!is.null(lavmodel))
stopifnot(!is.null(lavpartable))
if (is.null(est)) {
if (!is.null(lavpartable$est)) {
est <- lavpartable$est # if this changes, tag @TDJorgensen in commit message
} else {
lav_msg_stop(gettext("could not find `est' in lavpartable"))
}
}
} else {
lavmodel <- lavobject@Model
lavpartable <- lavobject@ParTable
if (is.null(est)) {
est <- lav_object_inspect_est(lavobject)
}
if (lavmodel@conditional.x) {
if (is.null(cov.x)) {
# try SampleStats slot
# if("SampleStats" %in% slotNames(lavobject)) {
# cov.x <- lavobject@SampleStats@cov.x
if (!is.null(lavobject@implied$cov.x[[1]])) {
cov.x <- lavobject@implied$cov.x # if this changes, tag @TDJorgensen in commit message
} else {
# perhaps lavaanList object
# extract it from GLIST per block
cov.x <- vector("list", length = lavmodel@nblocks)
for (b in seq_len(lavmodel@nblocks)) {
# which mm belong to block b?
mm.in.block <- (seq_len(lavmodel@nmat[b]) +
cumsum(c(0, lavmodel@nmat))[b])
MLIST <- lavmodel@GLIST[mm.in.block]
cov.x[[b]] <- MLIST[["cov.x"]]
}
}
}
}
}
if (is.null(partable)) {
partable <- lavpartable
}
if (is.null(GLIST)) {
GLIST <- lavmodel@GLIST
}
if (is.null(est.std)) {
est.std <- lav_standardize_lv(
lavobject = lavobject,
partable = partable, est = est, GLIST = GLIST,
cov.std = cov.std, lv.var = lv.var, lavmodel = lavmodel,
lavpartable = lavpartable
)
}
out <- est.std
N <- length(est.std)
stopifnot(N == length(partable$lhs))
VY <- computeVY(
lavmodel = lavmodel, GLIST = GLIST,
diagonal.only = TRUE
)
for (g in 1:lavmodel@nblocks) {
ov.names <- vnames(lavpartable, "ov", block = g)
ov.names.x <- vnames(lavpartable, "ov.x", block = g)
ov.names.nox <- vnames(lavpartable, "ov.nox", block = g)
lv.names <- vnames(lavpartable, "lv", block = g)
if (is.null(ov.var)) {
OV2 <- VY[[g]]
# replace zero values by NA (but keep negative values)
zero.idx <- which(abs(OV2) < .Machine$double.eps)
if (length(zero.idx) > 0L) {
OV2[zero.idx] <- as.numeric(NA)
}
# replace negative values by NA (for sqrt)
tmp.OV2 <- OV2
neg.idx <- which(tmp.OV2 < 0)
if (length(neg.idx) > 0L) {
tmp.OV2[neg.idx] <- as.numeric(NA)
}
OV <- sqrt(tmp.OV2)
} else {
OV2 <- ov.var[[g]]
OV <- sqrt(OV2)
}
if (lavmodel@conditional.x) {
# extend OV with ov.names.x
ov.names.x <- vnames(lavpartable, "ov.x", block = g)
ov.names <- c(ov.names.nox, ov.names.x)
OV2 <- c(OV2, diag(cov.x[[g]]))
OV <- c(OV, sqrt(diag(cov.x[[g]])))
}
# 1a. "=~" regular indicators
idx <- which(partable$op == "=~" & !(partable$rhs %in% lv.names) &
partable$block == g)
out[idx] <- out[idx] / OV[match(partable$rhs[idx], ov.names)]
# 1b. "=~" regular higher-order lv indicators
# 1c. "=~" indicators that are both in ov and lv
# idx <- which(partable$op == "=~" & partable$rhs %in% ov.names
# & partable$rhs %in% lv.names &
# partable$block == g)
# 2. "~" regressions (and "<~")
idx <- which((partable$op == "~" | partable$op == "<~") &
partable$lhs %in% ov.names &
partable$block == g)
out[idx] <- out[idx] / OV[match(partable$lhs[idx], ov.names)]
idx <- which((partable$op == "~" | partable$op == "<~") &
partable$rhs %in% ov.names.nox &
partable$block == g)
out[idx] <- out[idx] * OV[match(partable$rhs[idx], ov.names.nox)]
# 3a. "~~" ov
# ATTENTION: in Mplus 4.1, the off-diagonal residual covariances
# were computed by the formula cov(i,j) / sqrt(i.var*j.var)
# were i.var and j.var where diagonal elements of OV
#
# in Mplus 6.1 (but also AMOS and EQS), the i.var and j.var
# elements are the 'THETA' diagonal elements!!
# variances
rv.idx <- which(partable$op == "~~" & !(partable$lhs %in% lv.names) &
!(partable$lhs %in% ov.names.x) &
partable$lhs == partable$rhs &
partable$block == g)
# out[rv.idx] <- ( out[rv.idx] / OV[ match(partable$lhs[rv.idx], ov.names) ]
# / OV[ match(partable$rhs[rv.idx], ov.names) ] )
out[rv.idx] <- (out[rv.idx] /
OV2[match(partable$lhs[rv.idx], ov.names)])
# covariances ov
# three types:
# - only lhs is OV (and fixed.x = FALSE)
# - only rhs is OV (and fixed.x = FALSE)
# - both lhs and rhs are OV (regular case)
if (cov.std) {
if (!is.complex(est[rv.idx])) {
RV <- sqrt(abs(est[rv.idx]))
} else {
RV <- sqrt(est[rv.idx])
}
rv.names <- partable$lhs[rv.idx]
}
# left
idx.lhs <- which(partable$op == "~~" &
!(partable$lhs %in% lv.names) &
!(partable$lhs %in% ov.names.x) &
partable$lhs != partable$rhs &
partable$block == g)
if (length(idx.lhs) > 0L) {
if (cov.std == FALSE) {
out[idx.lhs] <-
(out[idx.lhs] / OV[match(partable$lhs[idx.lhs], ov.names)])
} else {
out[idx.lhs] <-
(out[idx.lhs] / RV[match(partable$lhs[idx.lhs], rv.names)])
}
}
# right
idx.rhs <- which(partable$op == "~~" &
!(partable$rhs %in% lv.names) &
!(partable$rhs %in% ov.names.x) &
partable$lhs != partable$rhs &
partable$block == g)
if (length(idx.rhs) > 0L) {
if (cov.std == FALSE) {
out[idx.rhs] <-
(out[idx.rhs] / OV[match(partable$rhs[idx.rhs], ov.names)])
} else {
out[idx.rhs] <-
(out[idx.rhs] / RV[match(partable$rhs[idx.rhs], rv.names)])
}
}
# 3b. "~~" lv
# idx <- which(partable$op == "~~" & partable$rhs %in% lv.names &
# partable$block == g)
# 4a. "~1" ov
idx <- which(partable$op == "~1" & !(partable$lhs %in% lv.names) &
!(partable$lhs %in% ov.names.x) &
partable$block == g)
out[idx] <- out[idx] / OV[match(partable$lhs[idx], ov.names)]
# 4b. "~1" lv
# idx <- which(partable$op == "~1" & partable$lhs %in% lv.names &
# partable$block == g)
# 4c. "|" thresholds
idx <- which(partable$op == "|" & !(partable$lhs %in% lv.names) &
partable$block == g)
out[idx] <- out[idx] / OV[match(partable$lhs[idx], ov.names)]
# 4d. "~*~" scales
idx <- which(partable$op == "~*~" & !(partable$lhs %in% lv.names) &
partable$block == g)
out[idx] <- 1.0
}
# 5a ":="
idx <- which(partable$op == ":=")
if (length(idx) > 0L) {
x <- out[partable$free & !duplicated(partable$free)]
out[idx] <- lavmodel@def.function(x)
}
# 5b "=="
idx <- which(partable$op == "==")
if (length(idx) > 0L) {
x <- out[partable$free & !duplicated(partable$free)]
out[idx] <- lavmodel@ceq.function(x)
}
# 5c. "<" or ">"
idx <- which((partable$op == "<" | partable$op == ">"))
if (length(idx) > 0L) {
x <- out[partable$free & !duplicated(partable$free)]
out[idx] <- lavmodel@cin.function(x)
}
out
}
lav_unstandardize_ov <- function(partable, ov.var = NULL, cov.std = TRUE) {
# check if ustart is missing; if so, look for est
if (is.null(partable$ustart)) {
partable$ustart <- partable$est
}
# check if block is missing
if (is.null(partable$block)) {
partable$block <- rep(1L, length(partable$ustart))
}
stopifnot(!any(is.na(partable$ustart)))
est <- out <- partable$ustart
N <- length(est)
# nblocks
nblocks <- lav_partable_nblocks(partable)
# if ov.var is NOT a list, make a list
if (!is.list(ov.var)) {
tmp <- ov.var
ov.var <- vector("list", length = nblocks)
ov.var[1:nblocks] <- list(tmp)
}
for (g in 1:nblocks) {
ov.names <- vnames(partable, "ov", block = g) # not user
lv.names <- vnames(partable, "lv", block = g)
OV <- sqrt(ov.var[[g]])
# 1a. "=~" regular indicators
idx <- which(partable$op == "=~" & !(partable$rhs %in% lv.names) &
partable$block == g)
out[idx] <- out[idx] * OV[match(partable$rhs[idx], ov.names)]
# 1b. "=~" regular higher-order lv indicators
# 1c. "=~" indicators that are both in ov and lv
# idx <- which(partable$op == "=~" & partable$rhs %in% ov.names
# & partable$rhs %in% lv.names &
# partable$block == g)
# 2. "~" regressions (and "<~")
idx <- which((partable$op == "~" | partable$op == "<~") &
partable$lhs %in% ov.names &
partable$block == g)
out[idx] <- out[idx] * OV[match(partable$lhs[idx], ov.names)]
idx <- which((partable$op == "~" | partable$op == "<~") &
partable$rhs %in% ov.names &
partable$block == g)
out[idx] <- out[idx] / OV[match(partable$rhs[idx], ov.names)]
# 3a. "~~" ov
# ATTENTION: in Mplus 4.1, the off-diagonal residual covariances
# were computed by the formula cov(i,j) / sqrt(i.var*j.var)
# were i.var and j.var where diagonal elements of OV
#
# in Mplus 6.1 (but also AMOS and EQS), the i.var and j.var
# elements are the 'THETA' diagonal elements!!
# variances
rv.idx <- which(partable$op == "~~" & !(partable$lhs %in% lv.names) &
partable$lhs == partable$rhs &
partable$block == g)
out[rv.idx] <- (out[rv.idx] * OV[match(partable$lhs[rv.idx], ov.names)]
* OV[match(partable$rhs[rv.idx], ov.names)])
# covariances
idx <- which(partable$op == "~~" & !(partable$lhs %in% lv.names) &
partable$lhs != partable$rhs &
partable$block == g)
if (length(idx) > 0L) {
if (cov.std == FALSE) {
out[idx] <- (out[idx] * OV[match(partable$lhs[idx], ov.names)]
* OV[match(partable$rhs[idx], ov.names)])
} else {
if (!is.complex(out[rv.idx])) {
RV <- sqrt(abs(out[rv.idx]))
} else {
RV <- sqrt(out[rv.idx])
}
rv.names <- partable$lhs[rv.idx]
out[idx] <- (out[idx] * RV[match(partable$lhs[idx], rv.names)]
* RV[match(partable$rhs[idx], rv.names)])
}
}
# 3b. "~~" lv
# idx <- which(partable$op == "~~" & partable$rhs %in% lv.names &
# partable$block == g)
# 4a. "~1" ov
idx <- which(partable$op == "~1" & !(partable$lhs %in% lv.names) &
partable$block == g)
out[idx] <- out[idx] * OV[match(partable$lhs[idx], ov.names)]
# 4b. "~1" lv
# idx <- which(partable$op == "~1" & partable$lhs %in% lv.names &
# partable$block == g)
}
# 5a ":="
# 5b "=="
# 5c. "<" or ">"
out
}
Try the lavaan package in your browser
Any scripts or data that you put into this service are public.
lavaan documentation built on Sept. 27, 2024, 9:07 a.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_unstandardize_ov lav_standardize_all_nox lav_standardize_all lav_standardize_lv lav_unstandardize_ov_x lav_standardize_all_nox_x lav_standardize_all_x lav_standardize_lv_x
lav_standardize_lv_x <- function(x, lavobject, partable = NULL, cov.std = TRUE,
lv.var = NULL,
rotation = FALSE) {
# set new values for x
lavmodel <- lav_model_set_parameters(lavmodel = lavobject@Model, x = x)
if (rotation) {
x.unrotated <- x
lavmodel@GLIST <- lavTech(lavobject, "est.unrotated") # unrotated!
est.rot <- lav_model_efa_rotate_x(
x = x.unrotated,
lavmodel = lavmodel, # unrotated!
lavoptions = lavobject@Options,
init.rot = lavmodel@H,
type = "user",
extra = TRUE
)
GLIST <- attr(est.rot, "extra")$GLIST
attributes(est.rot) <- NULL
est <- est.rot
} else {
GLIST <- lavmodel@GLIST # if this changes, tag @TDJorgensen in commit message
est <- lav_model_get_parameters(lavmodel, type = "user")
}
x.stand.user <- lav_standardize_lv(
lavobject = lavobject,
partable = partable, est = est,
GLIST = GLIST, cov.std = cov.std,
lv.var = lv.var
)
x.stand.user
}
lav_standardize_all_x <- function(x, lavobject, partable = NULL, cov.std = TRUE,
rotation = FALSE) {
lavmodel <- lav_model_set_parameters(lavmodel = lavobject@Model, x = x)
if (rotation) {
x.unrotated <- x
lavmodel@GLIST <- lavTech(lavobject, "est.unrotated") # unrotated!
est.rot <- lav_model_efa_rotate_x(
x = x.unrotated,
lavmodel = lavmodel, # unrotated!
lavoptions = lavobject@Options,
init.rot = lavmodel@H,
type = "user",
extra = TRUE
)
GLIST <- attr(est.rot, "extra")$GLIST
attributes(est.rot) <- NULL
est <- est.rot
} else {
GLIST <- lavmodel@GLIST # if this changes, tag @TDJorgensen in commit message
est <- lav_model_get_parameters(lavmodel, type = "user")
}
x.stand.user <- lav_standardize_all(
lavobject = lavobject,
partable = partable, est = est,
est.std = NULL, GLIST = GLIST,
cov.std = cov.std
)
x.stand.user
}
lav_standardize_all_nox_x <- function(x, lavobject, partable = NULL,
cov.std = TRUE, rotation = FALSE) {
lavmodel <- lav_model_set_parameters(lavmodel = lavobject@Model, x = x)
if (rotation) {
x.unrotated <- x
lavmodel@GLIST <- lavTech(lavobject, "est.unrotated") # unrotated!
est.rot <- lav_model_efa_rotate_x(
x = x.unrotated,
lavmodel = lavmodel, # unrotated!
lavoptions = lavobject@Options,
init.rot = lavmodel@H,
type = "user",
extra = TRUE
)
GLIST <- attr(est.rot, "extra")$GLIST
attributes(est.rot) <- NULL
est <- est.rot
} else {
GLIST <- lavmodel@GLIST # if this changes, tag @TDJorgensen in commit message
est <- lav_model_get_parameters(lavmodel, type = "user")
}
x.stand.user <- lav_standardize_all_nox(
lavobject = lavobject,
partable = partable, est = est,
est.std = NULL, GLIST = GLIST,
cov.std = cov.std
)
x.stand.user
}
lav_unstandardize_ov_x <- function(x, lavobject) {
partable <- lavobject@ParTable
partable$ustart <- x
lav_unstandardize_ov(
partable = partable,
ov.var = lavobject@SampleStats@var,
cov.std = TRUE
)
}
lav_standardize_lv <- function(lavobject = NULL,
partable = NULL, est = NULL, GLIST = NULL,
cov.std = TRUE, lv.var = NULL,
lavmodel = NULL, lavpartable = NULL) {
if (is.null(lavobject)) {
stopifnot(!is.null(lavmodel))
stopifnot(!is.null(lavpartable))
if (is.null(est)) {
if (!is.null(lavpartable$est)) {
est <- lavpartable$est # if this changes, tag @TDJorgensen in commit message
} else {
lav_msg_stop(gettext("could not find `est' in lavpartable"))
}
}
} else {
lavmodel <- lavobject@Model
lavpartable <- lavobject@ParTable
if (is.null(est)) {
est <- lav_object_inspect_est(lavobject)
}
}
if (is.null(partable)) {
partable <- lavpartable
}
if (is.null(GLIST)) {
GLIST <- lavmodel@GLIST
}
out <- est
N <- length(est)
stopifnot(N == length(partable$lhs))
nmat <- lavmodel@nmat
# compute ETA
if (is.null(lv.var)) {
LV.ETA <- computeVETA(
lavmodel = lavmodel,
GLIST = GLIST
)
}
for (g in 1:lavmodel@nblocks) {
ov.names <- vnames(lavpartable, "ov", block = g) # not user,
# which may be incomplete
lv.names <- vnames(lavpartable, "lv", block = g)
# shortcut: no latents in this block, nothing to do
if (length(lv.names) == 0L) {
next
}
# which mm belong to block g?
mm.in.group <- 1:nmat[g] + cumsum(c(0, nmat))[g]
MLIST <- GLIST[mm.in.group]
if (is.null(lv.var)) {
ETA2 <- diag(LV.ETA[[g]])
} else {
ETA2 <- lv.var[[g]]
}
# change negative values to NA
ETA2[ETA2 < 0] <- as.numeric(NA)
ETA <- sqrt(ETA2)
# 1a. "=~" regular indicators
idx <- which(partable$op == "=~" & !(partable$rhs %in% lv.names) &
partable$block == g)
out[idx] <- out[idx] * ETA[match(partable$lhs[idx], lv.names)]
# 1b. "=~" regular higher-order lv indicators
idx <- which(partable$op == "=~" & !(partable$rhs %in% ov.names) &
partable$block == g)
out[idx] <- (out[idx] * ETA[match(partable$lhs[idx], lv.names)]
/ ETA[match(partable$rhs[idx], lv.names)])
# 1c. "=~" indicators that are both in ov and lv
# idx <- which(partable$op == "=~" & partable$rhs %in% ov.names
# & partable$rhs %in% lv.names &
# partable$block == g)
# 2. "~" regressions (and "<~")
idx <- which((partable$op == "~" | partable$op == "<~") &
partable$lhs %in% lv.names &
partable$block == g)
out[idx] <- out[idx] / ETA[match(partable$lhs[idx], lv.names)]
idx <- which((partable$op == "~" | partable$op == "<~") &
partable$rhs %in% lv.names &
partable$block == g)
out[idx] <- out[idx] * ETA[match(partable$rhs[idx], lv.names)]
# 3a. "~~" ov
# idx <- which(partable$op == "~~" & !(partable$lhs %in% lv.names) &
# partable$block == g)
# 3b. "~~" lv
# ATTENTION: in Mplus 4.1, the off-diagonal residual covariances
# were computed by the formula cov(i,j) / sqrt(i.var*j.var)
# were i.var and j.var where diagonal elements of ETA
#
# in Mplus 6.1 (but also AMOS and EQS), the i.var and j.var
# elements are the 'PSI' diagonal elements!!
# variances
rv.idx <- which(partable$op == "~~" & partable$rhs %in% lv.names &
partable$lhs == partable$rhs &
partable$block == g)
out[rv.idx] <- (out[rv.idx] / ETA[match(partable$lhs[rv.idx], lv.names)]
/ ETA[match(partable$rhs[rv.idx], lv.names)])
# covariances lv
# three types:
# - only lhs is LV (and fixed.x = FALSE)
# - only rhs is LV (and fixed.x = FALSE)
# - both lhs and rhs are LV (regular case)
if (cov.std) {
if (!is.complex(est[rv.idx])) {
RV <- sqrt(abs(est[rv.idx])) # abs in case of heywood cases
} else {
RV <- sqrt(est[rv.idx])
}
rv.names <- partable$lhs[rv.idx]
}
# left
idx.lhs <- which(partable$op == "~~" &
partable$lhs %in% lv.names &
partable$lhs != partable$rhs &
partable$block == g)
if (length(idx.lhs) > 0L) {
if (cov.std == FALSE) {
out[idx.lhs] <-
(out[idx.lhs] / ETA[match(partable$lhs[idx.lhs], lv.names)])
} else {
out[idx.lhs] <-
(out[idx.lhs] / RV[match(partable$lhs[idx.lhs], rv.names)])
}
}
# right
idx.rhs <- which(partable$op == "~~" &
partable$rhs %in% lv.names &
partable$lhs != partable$rhs &
partable$block == g)
if (length(idx.rhs) > 0L) {
if (cov.std == FALSE) {
out[idx.rhs] <-
(out[idx.rhs] / ETA[match(partable$rhs[idx.rhs], lv.names)])
} else {
out[idx.rhs] <-
(out[idx.rhs] / RV[match(partable$rhs[idx.rhs], rv.names)])
}
}
# 4a. "~1" ov
# idx <- which(partable$op == "~1" & !(partable$lhs %in% lv.names) &
# partable$block == g)
# 4b. "~1" lv
idx <- which(partable$op == "~1" & partable$lhs %in% lv.names &
partable$block == g)
out[idx] <- out[idx] / ETA[match(partable$lhs[idx], lv.names)]
}
# 5a ":="
idx <- which(partable$op == ":=")
if (length(idx) > 0L) {
x <- out[partable$free & !duplicated(partable$free)]
out[idx] <- lavmodel@def.function(x)
}
# 5b "=="
idx <- which(partable$op == "==")
if (length(idx) > 0L) {
x <- out[partable$free & !duplicated(partable$free)]
out[idx] <- lavmodel@ceq.function(x)
}
# 5c. "<" or ">"
idx <- which((partable$op == "<" | partable$op == ">"))
if (length(idx) > 0L) {
x <- out[partable$free & !duplicated(partable$free)]
out[idx] <- lavmodel@cin.function(x)
}
out
}
lav_standardize_all <- function(lavobject = NULL,
partable = NULL, est = NULL, est.std = NULL,
GLIST = NULL, cov.std = TRUE, ov.var = NULL,
lv.var = NULL,
lavmodel = NULL, lavpartable = NULL,
cov.x = NULL) {
if (is.null(lavobject)) {
stopifnot(!is.null(lavmodel))
stopifnot(!is.null(lavpartable))
if (is.null(est)) {
if (!is.null(lavpartable$est)) {
est <- lavpartable$est # if this changes, tag @TDJorgensen in commit message
} else {
lav_msg_stop(gettext("could not find `est' in lavpartable"))
}
}
} else {
lavmodel <- lavobject@Model
lavpartable <- lavobject@ParTable
if (is.null(est)) {
est <- lav_object_inspect_est(lavobject)
}
if (lavmodel@conditional.x) {
if (is.null(cov.x)) {
# try SampleStats slot
# if("SampleStats" %in% slotNames(lavobject)) {
# cov.x <- lavobject@SampleStats@cov.x
if (!is.null(lavobject@implied$cov.x[[1]])) {
cov.x <- lavobject@implied$cov.x # if this changes, tag @TDJorgensen in commit message
} else {
# perhaps lavaanList object
# extract it from GLIST per block
cov.x <- vector("list", length = lavmodel@nblocks)
for (b in seq_len(lavmodel@nblocks)) {
# which mm belong to block b?
mm.in.block <- (seq_len(lavmodel@nmat[b]) +
cumsum(c(0, lavmodel@nmat))[b])
MLIST <- lavmodel@GLIST[mm.in.block]
cov.x[[b]] <- MLIST[["cov.x"]]
}
}
}
}
}
if (is.null(partable)) {
partable <- lavpartable
}
if (is.null(GLIST)) {
GLIST <- lavmodel@GLIST
}
if (is.null(est.std)) {
est.std <- lav_standardize_lv(
lavobject = lavobject,
partable = partable, est = est, GLIST = GLIST,
cov.std = cov.std, lv.var = lv.var, lavmodel = lavmodel,
lavpartable = lavpartable
)
}
out <- est.std
N <- length(est.std)
stopifnot(N == length(partable$lhs))
VY <- computeVY(
lavmodel = lavmodel, GLIST = GLIST,
diagonal.only = TRUE
)
for (g in 1:lavmodel@nblocks) {
ov.names <- vnames(lavpartable, "ov", block = g) # not user
lv.names <- vnames(lavpartable, "lv", block = g)
if (is.null(ov.var)) {
OV2 <- VY[[g]]
# replace zero values by NA (but keep negative values)
zero.idx <- which(abs(OV2) < .Machine$double.eps)
if (length(zero.idx) > 0L) {
OV2[zero.idx] <- as.numeric(NA)
}
# replace negative values by NA (for sqrt)
tmp.OV2 <- OV2
neg.idx <- which(tmp.OV2 < 0)
if (length(neg.idx) > 0L) {
tmp.OV2[neg.idx] <- as.numeric(NA)
}
OV <- sqrt(tmp.OV2)
} else {
OV2 <- ov.var[[g]]
OV <- sqrt(OV2)
}
if (lavmodel@conditional.x) {
# extend OV with ov.names.x
ov.names.x <- vnames(lavpartable, "ov.x", block = g)
ov.names.nox <- vnames(lavpartable, "ov.nox", block = g)
ov.names <- c(ov.names.nox, ov.names.x)
OV2 <- c(OV2, diag(cov.x[[g]]))
OV <- c(OV, sqrt(diag(cov.x[[g]])))
}
# 1a. "=~" regular indicators
idx <- which(partable$op == "=~" & !(partable$rhs %in% lv.names) &
partable$block == g)
out[idx] <- out[idx] / OV[match(partable$rhs[idx], ov.names)]
# 1b. "=~" regular higher-order lv indicators
# 1c. "=~" indicators that are both in ov and lv
# idx <- which(partable$op == "=~" & partable$rhs %in% ov.names
# & partable$rhs %in% lv.names &
# partable$block == g)
# 2. "~" regressions (and "<~")
idx <- which((partable$op == "~" | partable$op == "<~") &
partable$lhs %in% ov.names &
partable$block == g)
out[idx] <- out[idx] / OV[match(partable$lhs[idx], ov.names)]
idx <- which((partable$op == "~" | partable$op == "<~") &
partable$rhs %in% ov.names &
partable$block == g)
out[idx] <- out[idx] * OV[match(partable$rhs[idx], ov.names)]
# 3a. "~~" ov
# ATTENTION: in Mplus 4.1, the off-diagonal residual covariances
# were computed by the formula cov(i,j) / sqrt(i.var*j.var)
# were i.var and j.var where diagonal elements of OV
#
# in Mplus 6.1 (but also AMOS and EQS), the i.var and j.var
# elements are the 'THETA' diagonal elements!!
# variances
rv.idx <- which(partable$op == "~~" & !(partable$lhs %in% lv.names) &
partable$lhs == partable$rhs &
partable$block == g)
# out[rv.idx] <- ( out[rv.idx] / OV[ match(partable$lhs[rv.idx], ov.names) ]
# / OV[ match(partable$rhs[rv.idx], ov.names) ] )
out[rv.idx] <- (out[rv.idx] /
OV2[match(partable$lhs[rv.idx], ov.names)])
# covariances ov
# three types:
# - only lhs is OV (and fixed.x = FALSE)
# - only rhs is OV (and fixed.x = FALSE)
# - both lhs and rhs are OV (regular case)
if (cov.std) {
if (!is.complex(est[rv.idx])) {
RV <- sqrt(abs(est[rv.idx]))
} else {
RV <- sqrt(est[rv.idx])
}
rv.names <- partable$lhs[rv.idx]
}
# left
idx.lhs <- which(partable$op == "~~" &
!(partable$lhs %in% lv.names) &
partable$lhs != partable$rhs &
partable$block == g)
if (length(idx.lhs) > 0L) {
if (cov.std == FALSE) {
out[idx.lhs] <-
(out[idx.lhs] / OV[match(partable$lhs[idx.lhs], ov.names)])
} else {
out[idx.lhs] <-
(out[idx.lhs] / RV[match(partable$lhs[idx.lhs], rv.names)])
}
}
# right
idx.rhs <- which(partable$op == "~~" &
!(partable$rhs %in% lv.names) &
partable$lhs != partable$rhs &
partable$block == g)
if (length(idx.rhs) > 0L) {
if (cov.std == FALSE) {
out[idx.rhs] <-
(out[idx.rhs] / OV[match(partable$rhs[idx.rhs], ov.names)])
} else {
out[idx.rhs] <-
(out[idx.rhs] / RV[match(partable$rhs[idx.rhs], rv.names)])
}
}
# 3b. "~~" lv
# idx <- which(partable$op == "~~" & partable$rhs %in% lv.names &
# partable$block == g)
# 4a. "~1" ov
idx <- which(partable$op == "~1" & !(partable$lhs %in% lv.names) &
partable$block == g)
out[idx] <- out[idx] / OV[match(partable$lhs[idx], ov.names)]
# 4b. "~1" lv
# idx <- which(partable$op == "~1" & partable$lhs %in% lv.names &
# partable$block == g)
# 4c. "|" thresholds
idx <- which(partable$op == "|" & !(partable$lhs %in% lv.names) &
partable$block == g)
out[idx] <- out[idx] / OV[match(partable$lhs[idx], ov.names)]
# 4d. "~*~" scales
idx <- which(partable$op == "~*~" & !(partable$lhs %in% lv.names) &
partable$block == g)
out[idx] <- 1.0
}
# 5a ":="
idx <- which(partable$op == ":=")
if (length(idx) > 0L) {
x <- out[partable$free & !duplicated(partable$free)]
out[idx] <- lavmodel@def.function(x)
}
# 5b "=="
idx <- which(partable$op == "==")
if (length(idx) > 0L) {
x <- out[partable$free & !duplicated(partable$free)]
out[idx] <- lavmodel@ceq.function(x)
}
# 5c. "<" or ">"
idx <- which((partable$op == "<" | partable$op == ">"))
if (length(idx) > 0L) {
x <- out[partable$free & !duplicated(partable$free)]
out[idx] <- lavmodel@cin.function(x)
}
out
}
lav_standardize_all_nox <- function(lavobject = NULL,
partable = NULL, est = NULL, est.std = NULL,
GLIST = NULL, cov.std = TRUE, ov.var = NULL,
lv.var = NULL,
lavmodel = NULL, lavpartable = NULL,
cov.x = NULL) {
if (is.null(lavobject)) {
stopifnot(!is.null(lavmodel))
stopifnot(!is.null(lavpartable))
if (is.null(est)) {
if (!is.null(lavpartable$est)) {
est <- lavpartable$est # if this changes, tag @TDJorgensen in commit message
} else {
lav_msg_stop(gettext("could not find `est' in lavpartable"))
}
}
} else {
lavmodel <- lavobject@Model
lavpartable <- lavobject@ParTable
if (is.null(est)) {
est <- lav_object_inspect_est(lavobject)
}
if (lavmodel@conditional.x) {
if (is.null(cov.x)) {
# try SampleStats slot
# if("SampleStats" %in% slotNames(lavobject)) {
# cov.x <- lavobject@SampleStats@cov.x
if (!is.null(lavobject@implied$cov.x[[1]])) {
cov.x <- lavobject@implied$cov.x # if this changes, tag @TDJorgensen in commit message
} else {
# perhaps lavaanList object
# extract it from GLIST per block
cov.x <- vector("list", length = lavmodel@nblocks)
for (b in seq_len(lavmodel@nblocks)) {
# which mm belong to block b?
mm.in.block <- (seq_len(lavmodel@nmat[b]) +
cumsum(c(0, lavmodel@nmat))[b])
MLIST <- lavmodel@GLIST[mm.in.block]
cov.x[[b]] <- MLIST[["cov.x"]]
}
}
}
}
}
if (is.null(partable)) {
partable <- lavpartable
}
if (is.null(GLIST)) {
GLIST <- lavmodel@GLIST
}
if (is.null(est.std)) {
est.std <- lav_standardize_lv(
lavobject = lavobject,
partable = partable, est = est, GLIST = GLIST,
cov.std = cov.std, lv.var = lv.var, lavmodel = lavmodel,
lavpartable = lavpartable
)
}
out <- est.std
N <- length(est.std)
stopifnot(N == length(partable$lhs))
VY <- computeVY(
lavmodel = lavmodel, GLIST = GLIST,
diagonal.only = TRUE
)
for (g in 1:lavmodel@nblocks) {
ov.names <- vnames(lavpartable, "ov", block = g)
ov.names.x <- vnames(lavpartable, "ov.x", block = g)
ov.names.nox <- vnames(lavpartable, "ov.nox", block = g)
lv.names <- vnames(lavpartable, "lv", block = g)
if (is.null(ov.var)) {
OV2 <- VY[[g]]
# replace zero values by NA (but keep negative values)
zero.idx <- which(abs(OV2) < .Machine$double.eps)
if (length(zero.idx) > 0L) {
OV2[zero.idx] <- as.numeric(NA)
}
# replace negative values by NA (for sqrt)
tmp.OV2 <- OV2
neg.idx <- which(tmp.OV2 < 0)
if (length(neg.idx) > 0L) {
tmp.OV2[neg.idx] <- as.numeric(NA)
}
OV <- sqrt(tmp.OV2)
} else {
OV2 <- ov.var[[g]]
OV <- sqrt(OV2)
}
if (lavmodel@conditional.x) {
# extend OV with ov.names.x
ov.names.x <- vnames(lavpartable, "ov.x", block = g)
ov.names <- c(ov.names.nox, ov.names.x)
OV2 <- c(OV2, diag(cov.x[[g]]))
OV <- c(OV, sqrt(diag(cov.x[[g]])))
}
# 1a. "=~" regular indicators
idx <- which(partable$op == "=~" & !(partable$rhs %in% lv.names) &
partable$block == g)
out[idx] <- out[idx] / OV[match(partable$rhs[idx], ov.names)]
# 1b. "=~" regular higher-order lv indicators
# 1c. "=~" indicators that are both in ov and lv
# idx <- which(partable$op == "=~" & partable$rhs %in% ov.names
# & partable$rhs %in% lv.names &
# partable$block == g)
# 2. "~" regressions (and "<~")
idx <- which((partable$op == "~" | partable$op == "<~") &
partable$lhs %in% ov.names &
partable$block == g)
out[idx] <- out[idx] / OV[match(partable$lhs[idx], ov.names)]
idx <- which((partable$op == "~" | partable$op == "<~") &
partable$rhs %in% ov.names.nox &
partable$block == g)
out[idx] <- out[idx] * OV[match(partable$rhs[idx], ov.names.nox)]
# 3a. "~~" ov
# ATTENTION: in Mplus 4.1, the off-diagonal residual covariances
# were computed by the formula cov(i,j) / sqrt(i.var*j.var)
# were i.var and j.var where diagonal elements of OV
#
# in Mplus 6.1 (but also AMOS and EQS), the i.var and j.var
# elements are the 'THETA' diagonal elements!!
# variances
rv.idx <- which(partable$op == "~~" & !(partable$lhs %in% lv.names) &
!(partable$lhs %in% ov.names.x) &
partable$lhs == partable$rhs &
partable$block == g)
# out[rv.idx] <- ( out[rv.idx] / OV[ match(partable$lhs[rv.idx], ov.names) ]
# / OV[ match(partable$rhs[rv.idx], ov.names) ] )
out[rv.idx] <- (out[rv.idx] /
OV2[match(partable$lhs[rv.idx], ov.names)])
# covariances ov
# three types:
# - only lhs is OV (and fixed.x = FALSE)
# - only rhs is OV (and fixed.x = FALSE)
# - both lhs and rhs are OV (regular case)
if (cov.std) {
if (!is.complex(est[rv.idx])) {
RV <- sqrt(abs(est[rv.idx]))
} else {
RV <- sqrt(est[rv.idx])
}
rv.names <- partable$lhs[rv.idx]
}
# left
idx.lhs <- which(partable$op == "~~" &
!(partable$lhs %in% lv.names) &
!(partable$lhs %in% ov.names.x) &
partable$lhs != partable$rhs &
partable$block == g)
if (length(idx.lhs) > 0L) {
if (cov.std == FALSE) {
out[idx.lhs] <-
(out[idx.lhs] / OV[match(partable$lhs[idx.lhs], ov.names)])
} else {
out[idx.lhs] <-
(out[idx.lhs] / RV[match(partable$lhs[idx.lhs], rv.names)])
}
}
# right
idx.rhs <- which(partable$op == "~~" &
!(partable$rhs %in% lv.names) &
!(partable$rhs %in% ov.names.x) &
partable$lhs != partable$rhs &
partable$block == g)
if (length(idx.rhs) > 0L) {
if (cov.std == FALSE) {
out[idx.rhs] <-
(out[idx.rhs] / OV[match(partable$rhs[idx.rhs], ov.names)])
} else {
out[idx.rhs] <-
(out[idx.rhs] / RV[match(partable$rhs[idx.rhs], rv.names)])
}
}
# 3b. "~~" lv
# idx <- which(partable$op == "~~" & partable$rhs %in% lv.names &
# partable$block == g)
# 4a. "~1" ov
idx <- which(partable$op == "~1" & !(partable$lhs %in% lv.names) &
!(partable$lhs %in% ov.names.x) &
partable$block == g)
out[idx] <- out[idx] / OV[match(partable$lhs[idx], ov.names)]
# 4b. "~1" lv
# idx <- which(partable$op == "~1" & partable$lhs %in% lv.names &
# partable$block == g)
# 4c. "|" thresholds
idx <- which(partable$op == "|" & !(partable$lhs %in% lv.names) &
partable$block == g)
out[idx] <- out[idx] / OV[match(partable$lhs[idx], ov.names)]
# 4d. "~*~" scales
idx <- which(partable$op == "~*~" & !(partable$lhs %in% lv.names) &
partable$block == g)
out[idx] <- 1.0
}
# 5a ":="
idx <- which(partable$op == ":=")
if (length(idx) > 0L) {
x <- out[partable$free & !duplicated(partable$free)]
out[idx] <- lavmodel@def.function(x)
}
# 5b "=="
idx <- which(partable$op == "==")
if (length(idx) > 0L) {
x <- out[partable$free & !duplicated(partable$free)]
out[idx] <- lavmodel@ceq.function(x)
}
# 5c. "<" or ">"
idx <- which((partable$op == "<" | partable$op == ">"))
if (length(idx) > 0L) {
x <- out[partable$free & !duplicated(partable$free)]
out[idx] <- lavmodel@cin.function(x)
}
out
}
lav_unstandardize_ov <- function(partable, ov.var = NULL, cov.std = TRUE) {
# check if ustart is missing; if so, look for est
if (is.null(partable$ustart)) {
partable$ustart <- partable$est
}
# check if block is missing
if (is.null(partable$block)) {
partable$block <- rep(1L, length(partable$ustart))
}
stopifnot(!any(is.na(partable$ustart)))
est <- out <- partable$ustart
N <- length(est)
# nblocks
nblocks <- lav_partable_nblocks(partable)
# if ov.var is NOT a list, make a list
if (!is.list(ov.var)) {
tmp <- ov.var
ov.var <- vector("list", length = nblocks)
ov.var[1:nblocks] <- list(tmp)
}
for (g in 1:nblocks) {
ov.names <- vnames(partable, "ov", block = g) # not user
lv.names <- vnames(partable, "lv", block = g)
OV <- sqrt(ov.var[[g]])
# 1a. "=~" regular indicators
idx <- which(partable$op == "=~" & !(partable$rhs %in% lv.names) &
partable$block == g)
out[idx] <- out[idx] * OV[match(partable$rhs[idx], ov.names)]
# 1b. "=~" regular higher-order lv indicators
# 1c. "=~" indicators that are both in ov and lv
# idx <- which(partable$op == "=~" & partable$rhs %in% ov.names
# & partable$rhs %in% lv.names &
# partable$block == g)
# 2. "~" regressions (and "<~")
idx <- which((partable$op == "~" | partable$op == "<~") &
partable$lhs %in% ov.names &
partable$block == g)
out[idx] <- out[idx] * OV[match(partable$lhs[idx], ov.names)]
idx <- which((partable$op == "~" | partable$op == "<~") &
partable$rhs %in% ov.names &
partable$block == g)
out[idx] <- out[idx] / OV[match(partable$rhs[idx], ov.names)]
# 3a. "~~" ov
# ATTENTION: in Mplus 4.1, the off-diagonal residual covariances
# were computed by the formula cov(i,j) / sqrt(i.var*j.var)
# were i.var and j.var where diagonal elements of OV
#
# in Mplus 6.1 (but also AMOS and EQS), the i.var and j.var
# elements are the 'THETA' diagonal elements!!
# variances
rv.idx <- which(partable$op == "~~" & !(partable$lhs %in% lv.names) &
partable$lhs == partable$rhs &
partable$block == g)
out[rv.idx] <- (out[rv.idx] * OV[match(partable$lhs[rv.idx], ov.names)]
* OV[match(partable$rhs[rv.idx], ov.names)])
# covariances
idx <- which(partable$op == "~~" & !(partable$lhs %in% lv.names) &
partable$lhs != partable$rhs &
partable$block == g)
if (length(idx) > 0L) {
if (cov.std == FALSE) {
out[idx] <- (out[idx] * OV[match(partable$lhs[idx], ov.names)]
* OV[match(partable$rhs[idx], ov.names)])
} else {
if (!is.complex(out[rv.idx])) {
RV <- sqrt(abs(out[rv.idx]))
} else {
RV <- sqrt(out[rv.idx])
}
rv.names <- partable$lhs[rv.idx]
out[idx] <- (out[idx] * RV[match(partable$lhs[idx], rv.names)]
* RV[match(partable$rhs[idx], rv.names)])
}
}
# 3b. "~~" lv
# idx <- which(partable$op == "~~" & partable$rhs %in% lv.names &
# partable$block == g)
# 4a. "~1" ov
idx <- which(partable$op == "~1" & !(partable$lhs %in% lv.names) &
partable$block == g)
out[idx] <- out[idx] * OV[match(partable$lhs[idx], ov.names)]
# 4b. "~1" lv
# idx <- which(partable$op == "~1" & partable$lhs %in% lv.names &
# partable$block == g)
}
# 5a ":="
# 5b "=="
# 5c. "<" or ">"
out
}
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.