Nothing
R/lav_model_gradient.R
In lavaan: Latent Variable Analysis
Defines functions
lav_model_gradient_DD
computeOmega
computeDeltaDx
computeDelta
lav_model_gradient
# model gradient
lav_model_gradient <- function(lavmodel = NULL,
GLIST = NULL,
lavsamplestats = NULL,
lavdata = NULL,
lavcache = NULL,
type = "free",
verbose = FALSE,
group.weight = TRUE,
Delta = NULL,
m.el.idx = NULL,
x.el.idx = NULL,
ceq.simple = FALSE) {
nmat <- lavmodel@nmat
estimator <- lavmodel@estimator
representation <- lavmodel@representation
meanstructure <- lavmodel@meanstructure
categorical <- lavmodel@categorical
group.w.free <- lavmodel@group.w.free
fixed.x <- lavmodel@fixed.x
conditional.x <- lavmodel@conditional.x
num.idx <- lavmodel@num.idx
th.idx <- lavmodel@th.idx
nx.free <- lavmodel@nx.free
if(.hasSlot(lavmodel, "estimator.args")) {
estimator.args <- lavmodel@estimator.args
} else {
estimator.args <- list()
}
# state or final?
if(is.null(GLIST)) GLIST <- lavmodel@GLIST
if(estimator == "REML") warning("analytical gradient not implement; use numerical approximation")
# group.weight
# FIXME --> block.weight
if(group.weight) {
if(estimator %in% c("ML","PML","FML","MML","REML","NTRLS","catML")) {
group.w <- (unlist(lavsamplestats@nobs)/lavsamplestats@ntotal)
} else if(estimator == "DLS") {
if(estimator.args$dls.FtimesNminus1) {
group.w <- ((unlist(lavsamplestats@nobs)-1)/lavsamplestats@ntotal)
} else {
group.w <- (unlist(lavsamplestats@nobs)/lavsamplestats@ntotal)
}
} else {
# FIXME: double check!
group.w <- ((unlist(lavsamplestats@nobs)-1)/lavsamplestats@ntotal)
}
} else {
group.w <- rep(1.0, lavmodel@nblocks)
}
# do we need WLS.est?
if(estimator %in% c("WLS", "DWLS", "ULS", "GLS", "NTRLS", "DLS")) {
# always compute WLS.est
WLS.est <- lav_model_wls_est(lavmodel = lavmodel, GLIST = GLIST) #,
# cov.x = lavsamplestats@cov.x)
}
if(estimator %in% c("ML", "PML", "FML", "REML", "NTRLS", "catML")) {
# compute moments for all groups
#if(conditional.x) {
# Sigma.hat <- computeSigmaHatJoint(lavmodel = lavmodel,
# GLIST = GLIST,
# extra = (estimator %in% c("ML", "REML","NTRLS")))
#} else {
Sigma.hat <- computeSigmaHat(lavmodel = lavmodel, GLIST = GLIST,
extra = (estimator %in% c("ML", "REML",
"NTRLS", "catML")))
#}
if(meanstructure) {
#if(conditional.x) {
# Mu.hat <- computeMuHat(lavmodel = lavmodel, GLIST = GLIST)
#} else {
Mu.hat <- computeMuHat(lavmodel = lavmodel, GLIST = GLIST)
#}
}
if(categorical) {
TH <- computeTH(lavmodel = lavmodel, GLIST = GLIST)
}
if(conditional.x) {
PI <- computePI(lavmodel = lavmodel, GLIST = GLIST)
} else if(estimator == "PML") {
PI <- vector("list", length = lavmodel@nblocks)
}
if(group.w.free) {
GW <- computeGW(lavmodel = lavmodel, GLIST = GLIST)
}
} else if(estimator == "DLS" && estimator.args$dls.GammaNT == "model") {
Sigma.hat <- computeSigmaHat(lavmodel = lavmodel, GLIST = GLIST,
extra = FALSE)
Mu.hat <- computeMuHat(lavmodel = lavmodel, GLIST = GLIST)
} else if(estimator == "MML") {
TH <- computeTH( lavmodel = lavmodel, GLIST = GLIST)
THETA <- computeTHETA(lavmodel = lavmodel, GLIST = GLIST)
GW <- computeGW( lavmodel = lavmodel, GLIST = GLIST)
}
# four approaches (FIXME!!!! merge this!)
# - ML approach: using Omega (and Omega.mu)
# Omega = 'POST' = Sigma.inv %*% (S - Sigma) %*% t(Sigma.inv)
# (still 2x faster than Delta method)
# - WLS/DWLS/GLS: using Delta + WLS.V; support for fixed.x, conditional.x
# - (ML)/NTRLS: using Delta, no support for fixed.x, conditional.x
# - PML/FML/MML: custom
# 1. ML approach
if( (estimator == "ML" || estimator == "REML" || estimator == "catML") &&
lavdata@nlevels == 1L &&
!lavmodel@conditional.x ) {
if(meanstructure) {
Omega <- computeOmega(Sigma.hat=Sigma.hat, Mu.hat=Mu.hat,
lavsamplestats=lavsamplestats,
estimator=estimator,
meanstructure=TRUE,
conditional.x = conditional.x)
Omega.mu <- attr(Omega, "mu")
} else {
Omega <- computeOmega(Sigma.hat=Sigma.hat, Mu.hat=NULL,
lavsamplestats=lavsamplestats,
estimator=estimator,
meanstructure=FALSE,
conditional.x = conditional.x)
Omega.mu <- vector("list", length = lavmodel@nblocks)
}
# compute DX (for all elements in every model matrix)
DX <- vector("list", length=length(GLIST))
for(g in 1:lavmodel@nblocks) {
# which mm belong to group g?
mm.in.group <- 1:nmat[g] + cumsum(c(0,nmat))[g]
mm.names <- names( GLIST[mm.in.group] )
if(representation == "LISREL") {
DX.group <- derivative.F.LISREL(GLIST[mm.in.group],
Omega[[g]],
Omega.mu[[g]])
# FIXME!!!
# add empty gamma
if(lavmodel@conditional.x) {
DX.group$gamma <- lavmodel@GLIST$gamma
}
# only save what we need
DX[mm.in.group] <- DX.group[ mm.names ]
} else if(representation == "RAM") {
DX.group <- lav_ram_df(GLIST[mm.in.group],
Omega[[g]],
Omega.mu[[g]])
# only save what we need
DX[mm.in.group] <- DX.group[ mm.names ]
} else {
stop("only LISREL and RAM representation has been implemented for now")
}
# weight by group
if(lavmodel@nblocks > 1L) {
for(mm in mm.in.group) {
DX[[mm]] <- group.w[g] * DX[[mm]]
}
}
}
# extract free parameters
if(type == "free") {
if(.hasSlot(lavmodel, "ceq.simple.only") &&
lavmodel@ceq.simple.only) { # new in 0.6-11
dx <- numeric( lavmodel@nx.unco )
for(g in 1:lavmodel@nblocks) {
mm.in.group <- 1:nmat[g] + cumsum(c(0,nmat))[g]
for(mm in mm.in.group) {
m.free.idx <- lavmodel@m.free.idx[[mm]]
x.unco.idx <- lavmodel@x.unco.idx[[mm]]
dx[x.unco.idx] <- DX[[mm]][m.free.idx]
}
}
if(ceq.simple) {
dx <- drop( crossprod(lavmodel@ceq.simple.K, dx) )
}
} else {
dx <- numeric( nx.free )
for(g in 1:lavmodel@nblocks) {
mm.in.group <- 1:nmat[g] + cumsum(c(0,nmat))[g]
for(mm in mm.in.group) {
m.free.idx <- lavmodel@m.free.idx[[mm]]
x.free.idx <- lavmodel@x.free.idx[[mm]]
dx[x.free.idx] <- DX[[mm]][m.free.idx]
}
}
}
} else {
dx <- DX
# handle equality constraints
### FIXME!!!! TODO!!!!
}
} else # ML
# 2. using Delta - *LS family
if(estimator %in% c("WLS", "DWLS", "ULS", "GLS", "NTGLS", "DLS")) {
if(type != "free") {
if(is.null(Delta))
stop("FIXME: Delta should be given if type != free")
#stop("FIXME: WLS gradient with type != free needs fixing!")
} else {
Delta <- computeDelta(lavmodel = lavmodel, GLIST. = GLIST,
ceq.simple = ceq.simple)
}
for(g in 1:lavmodel@nblocks) {
#diff <- as.matrix(lavsamplestats@WLS.obs[[g]] - WLS.est[[g]])
#group.dx <- -1 * ( t(Delta[[g]]) %*% lavsamplestats@WLS.V[[g]] %*% diff)
# 0.5-17: use crossprod twice; treat DWLS/ULS special
if(estimator == "WLS" ||
estimator == "GLS" ||
estimator == "DLS" ||
estimator == "NTRLS") {
# full weight matrix
diff <- lavsamplestats@WLS.obs[[g]] - WLS.est[[g]]
# full weight matrix
if(estimator == "GLS" || estimator == "WLS") {
WLS.V <- lavsamplestats@WLS.V[[g]]
group.dx <- -1 * crossprod(Delta[[g]],
crossprod(WLS.V, diff))
} else if(estimator == "DLS") {
if(estimator.args$dls.GammaNT == "sample") {
WLS.V <- lavsamplestats@WLS.V[[g]] # for now
} else {
dls.a <- estimator.args$dls.a
GammaNT <- lav_samplestats_Gamma_NT(
COV = Sigma.hat[[g]],
MEAN = Mu.hat[[g]],
rescale = FALSE,
x.idx = lavsamplestats@x.idx[[g]],
fixed.x = lavmodel@fixed.x,
conditional.x = lavmodel@conditional.x,
meanstructure = lavmodel@meanstructure,
slopestructure = lavmodel@conditional.x)
W.DLS <- (1 - dls.a)*lavsamplestats@NACOV[[g]] + dls.a*GammaNT
WLS.V <- lav_matrix_symmetric_inverse(W.DLS)
}
group.dx <- -1 * crossprod(Delta[[g]],
crossprod(WLS.V, diff))
} else if(estimator == "NTRLS") {
stopifnot(!conditional.x)
#WLS.V <- lav_samplestats_Gamma_inverse_NT(
# ICOV = attr(Sigma.hat[[g]],"inv")[,,drop=FALSE],
# COV = Sigma.hat[[g]][,,drop=FALSE],
# MEAN = Mu.hat[[g]],
# x.idx = lavsamplestats@x.idx[[g]],
# fixed.x = fixed.x,
# conditional.x = conditional.x,
# meanstructure = meanstructure,
# slopestructure = conditional.x)
S <- lavsamplestats@cov[[g]]
Sigma <- Sigma.hat[[g]]
Sigma.inv <- attr(Sigma, "inv")
nvar <- NROW(Sigma)
if(meanstructure) {
MEAN <- lavsamplestats@mean[[g]]; Mu <- Mu.hat[[g]]
POST.Sigma <- lav_matrix_duplication_pre(
matrix((Sigma.inv %*% (S - Sigma) %*% t(Sigma.inv)) %*%
(diag(nvar) + (S - Sigma) %*% Sigma.inv) +
(Sigma.inv %*% tcrossprod(MEAN - Mu) %*% Sigma.inv),
ncol = 1) )
POST.Mu <- as.numeric(2 * Sigma.inv %*% (MEAN - Mu))
POST <- c(POST.Mu, POST.Sigma)
} else {
POST <- lav_matrix_duplication_pre(
matrix((Sigma.inv %*% (S - Sigma) %*% t(Sigma.inv)) %*%
(diag(nvar) + (S - Sigma) %*% Sigma.inv), ncol = 1))
}
group.dx <- as.numeric( -1 * crossprod(Delta[[g]], POST) )
}
} else
if(estimator == "DWLS" || estimator == "ULS") {
# diagonal weight matrix
diff <- lavsamplestats@WLS.obs[[g]] - WLS.est[[g]]
group.dx <- -1 * crossprod(Delta[[g]],
lavsamplestats@WLS.VD[[g]] * diff)
}
group.dx <- group.w[g] * group.dx
if(g == 1) {
dx <- group.dx
} else {
dx <- dx + group.dx
}
} # g
if(type == "free") {
# nothing to do
} else {
# make a GLIST
dx <- lav_model_x2GLIST(lavmodel = lavmodel, x = dx,
type = "custom", setDelta = FALSE,
m.el.idx = m.el.idx,
x.el.idx = x.el.idx)
}
} # WLS
# ML + conditional.x
else if(estimator %in% c("ML", "catML") && lavmodel@conditional.x
&& lavdata@nlevels == 1L) {
if(type != "free") {
if(is.null(Delta))
stop("FIXME: Delta should be given if type != free")
#stop("FIXME: WLS gradient with type != free needs fixing!")
} else {
Delta <- computeDelta(lavmodel = lavmodel, GLIST. = GLIST,
ceq.simple = ceq.simple)
}
for(g in 1:lavmodel@nblocks) {
# augmented mean.x + cov.x matrix
mean.x <- lavsamplestats@mean.x[[g]]
cov.x <- lavsamplestats@cov.x[[g]]
C3 <- rbind(c(1,mean.x),
cbind(mean.x, cov.x + tcrossprod(mean.x)))
Sigma <- Sigma.hat[[g]]
Mu.g <- Mu.hat[[g]]
PI.g <- PI[[g]]
Sigma.inv <- attr(Sigma, "inv")
nvar <- NROW(Sigma)
S <- lavsamplestats@res.cov[[g]]
# beta
OBS <- t( cbind(lavsamplestats@res.int[[g]],
lavsamplestats@res.slopes[[g]]) )
EST <- t( cbind(Mu.g, PI.g) )
#obs.beta <- c(lavsamplestats@res.int[[g]],
# lav_matrix_vec(lavsamplestats@res.slopes[[g]]))
#est.beta <- c(Mu.g, lav_matrix_vec(PI.g))
#beta.COV <- C3 %x% Sigma.inv
#a <- t(obs.beta - est.beta)
#b <- as.matrix(obs.beta - est.beta)
#K <- lav_matrix_commutation(m = nvar, n = nvar)
#AB <- (K %x% diag(NROW(C3)*NROW(C3))) %*%
# (diag(nvar) %x% lav_matrix_vec(C3) %x% diag(nvar))
#K <- lav_matrix_commutation(m = nvar, n = NROW(C3))
#AB <- ( diag(NROW(C3)) %x% K %x% diag(nvar) ) %*%
# (lav_matrix_vec(C3) %x% diag( nvar * nvar) )
#POST.beta <- 2 * beta.COV %*% (obs.beta - est.beta)
d.BETA <- C3 %*% (OBS - EST) %*% Sigma.inv
# NOTE: the vecr here, unlike lav_mvreg_dlogl_beta
# this is because DELTA has used vec(t(BETA)),
# instead of vec(BETA)
#POST.beta <- 2 * lav_matrix_vecr(d.BETA)
# NOT any longer, since 0.6-1!!!
POST.beta <- 2 * lav_matrix_vec(d.BETA)
#POST.sigma1 <- lav_matrix_duplication_pre(
# (Sigma.inv %x% Sigma.inv) %*% t(AB) %*% (t(a) %x% b) )
# Sigma
#POST.sigma2 <- lav_matrix_duplication_pre(
# matrix( lav_matrix_vec(
# Sigma.inv %*% (S - Sigma) %*% t(Sigma.inv)), ncol = 1L))
W.tilde <- S + t(OBS - EST) %*% C3 %*% (OBS - EST)
d.SIGMA <- (Sigma.inv - Sigma.inv %*% W.tilde %*% Sigma.inv)
d.vechSigma <- as.numeric( lav_matrix_duplication_pre(
as.matrix(lav_matrix_vec(d.SIGMA)) ) )
POST.sigma <- -1 * d.vechSigma
#POST <- c(POST.beta, POST.sigma1 + POST.sigma2)
POST <- c(POST.beta, POST.sigma)
group.dx <- as.numeric( -1 * crossprod(Delta[[g]], POST) )
# because we still use obj/2, we need to divide by 2!
group.dx <- group.dx / 2 # fixed in 0.6-1
group.dx <- group.w[g] * group.dx
if(g == 1) {
dx <- group.dx
} else {
dx <- dx + group.dx
}
} # g
if(type == "free") {
# nothing to do
} else {
# make a GLIST
dx <- lav_model_x2GLIST(lavmodel = lavmodel, x = dx,
type = "custom", setDelta = FALSE,
m.el.idx = m.el.idx,
x.el.idx = x.el.idx)
}
} # ML + conditional.x
else if(estimator == "ML" && lavdata@nlevels > 1L) {
if(type != "free") {
stop("FIXME: type != free in lav_model_gradient for estimator ML for nlevels > 1")
} else {
Delta <- computeDelta(lavmodel = lavmodel, GLIST. = GLIST,
ceq.simple = ceq.simple)
}
# for each upper-level group....
for(g in 1:lavmodel@ngroups) {
if(!lavsamplestats@missing.flag) { # complete data
if(lavmodel@conditional.x) {
DX <- lav_mvreg_cluster_dlogl_2l_samplestats(
YLp = lavsamplestats@YLp[[g]],
Lp = lavdata@Lp[[g]],
Res.Sigma.W = Sigma.hat[[(g-1)*2 + 1]],
Res.Int.W = Mu.hat[[(g-1)*2 + 1]],
Res.Pi.W = PI[[(g-1)*2 + 1]],
Res.Sigma.B = Sigma.hat[[(g-1)*2 + 2]],
Res.Int.B = Mu.hat[[(g-1)*2 + 2]],
Res.Pi.B = PI[[(g-1)*2 + 2]],
Sinv.method = "eigen")
} else {
DX <- lav_mvnorm_cluster_dlogl_2l_samplestats(
YLp = lavsamplestats@YLp[[g]],
Lp = lavdata@Lp[[g]],
Mu.W = Mu.hat[[(g-1)*2 + 1]],
Sigma.W = Sigma.hat[[(g-1)*2 + 1]],
Mu.B = Mu.hat[[(g-1)*2 + 2]],
Sigma.B = Sigma.hat[[(g-1)*2 + 2]],
Sinv.method = "eigen")
}
} else {
# missing data
if(lavmodel@conditional.x) {
stop("lavaan ERROR: gradient for twolevel + conditional.x + fiml is not ready; use optim.gradient = \"numerical\"")
} else {
DX <- lav_mvnorm_cluster_missing_dlogl_2l_samplestats(
Y1 = lavdata@X[[g]],
Y2 = lavsamplestats@YLp[[g]][[2]]$Y2,
Lp = lavdata@Lp[[g]],
Mp = lavdata@Mp[[g]],
Mu.W = Mu.hat[[(g-1)*2 + 1]],
Sigma.W = Sigma.hat[[(g-1)*2 + 1]],
Mu.B = Mu.hat[[(g-1)*2 + 2]],
Sigma.B = Sigma.hat[[(g-1)*2 + 2]],
Sinv.method = "eigen")
}
}
group.dx <- as.numeric( DX %*% Delta[[g]] )
# group weights (if any)
group.dx <- group.w[g] * group.dx
if(g == 1) {
dx <- group.dx
} else {
dx <- dx + group.dx
}
} # g
# divide by 2 * N
dx <- dx / (2 * lavsamplestats@ntotal)
#cat("dx1 (numerical) = \n"); print( zapsmall(dx1) )
#cat("dx (analytic) = \n"); print( zapsmall(dx ) )
} # ML + two-level
else if(estimator == "PML" || estimator == "FML" ||
estimator == "MML") {
if(type != "free") {
stop("FIXME: type != free in lav_model_gradient for estimator PML")
} else {
Delta <- computeDelta(lavmodel = lavmodel, GLIST. = GLIST,
ceq.simple = ceq.simple)
}
for(g in 1:lavmodel@nblocks) {
#print(GLIST)
#print(lav_model_get_parameters(lavmodel = lavmodel, GLIST = GLIST))
#print(Sigma.hat[[g]])
#print(TH[[g]])
#cat("*****\n")
# compute partial derivative of logLik with respect to
# thresholds/means, slopes, variances, correlations
if(estimator == "PML") {
if(lavdata@nlevels > 1L) {
stop("lavaan ERROR: PL gradient + multilevel not implemented; try optim.gradient = \"numerical\"")
} else if(conditional.x) {
d1 <- pml_deriv1(Sigma.hat = Sigma.hat[[g]],
Mu.hat = Mu.hat[[g]],
TH = TH[[g]],
th.idx = th.idx[[g]],
num.idx = num.idx[[g]],
X = lavdata@X[[g]],
lavcache = lavcache[[g]],
eXo = lavdata@eXo[[g]],
PI = PI[[g]],
missing = lavdata@missing)
} else {
d1 <- pml_deriv1(Sigma.hat = Sigma.hat[[g]],
Mu.hat = Mu.hat[[g]],
TH = TH[[g]],
th.idx = th.idx[[g]],
num.idx = num.idx[[g]],
X = lavdata@X[[g]],
lavcache = lavcache[[g]],
eXo = NULL,
PI = NULL,
missing = lavdata@missing)
} # not conditional.x
# chain rule (fmin)
group.dx <-
as.numeric(t(d1) %*% Delta[[g]])
} # PML
else if(estimator == "FML") {
d1 <- fml_deriv1(Sigma.hat = Sigma.hat[[g]],
TH = TH[[g]],
th.idx = th.idx[[g]],
num.idx = num.idx[[g]],
X = lavdata@X[[g]],
lavcache = lavcache[[g]])
# chain rule (fmin)
group.dx <-
as.numeric(t(d1) %*% Delta[[g]])/lavsamplestats@nobs[[g]]
} else if(estimator == "MML") {
group.dx <-
lav_model_gradient_mml(lavmodel = lavmodel,
GLIST = GLIST,
THETA = THETA[[g]],
TH = TH[[g]],
group = g,
lavdata = lavdata,
sample.mean = lavsamplestats@mean[[g]],
sample.mean.x = lavsamplestats@mean.x[[g]],
lavcache = lavcache)
}
# group weights (if any)
group.dx <- group.w[g] * group.dx
if(g == 1) {
dx <- group.dx
} else {
dx <- dx + group.dx
}
} # g
} else {
stop("lavaan ERROR: no analytical gradient available for estimator ",
estimator)
}
# group.w.free for ML
if(lavmodel@group.w.free &&
estimator %in% c("ML","MML","FML","PML","REML","catML")) {
#est.prop <- unlist( computeGW(lavmodel = lavmodel, GLIST = GLIST) )
#obs.prop <- unlist(lavsamplestats@group.w)
# FIXME: G2 based -- ML and friends only!!
#dx.GW <- - (obs.prop - est.prop)
# poisson version
est.freq <- exp(unlist(computeGW(lavmodel = lavmodel, GLIST = GLIST)))
obs.freq <- unlist(lavsamplestats@group.w) * lavsamplestats@ntotal
dx.GW <- - (obs.freq - est.freq)
# divide by N (to be consistent with the rest of lavaan)
dx.GW <- dx.GW / lavsamplestats@ntotal
# remove last element (fixed LAST group to zero)
# dx.GW <- dx.GW[-length(dx.GW)]
# fill in in dx
gw.mat.idx <- which(names(lavmodel@GLIST) == "gw")
gw.x.idx <- unlist( lavmodel@x.free.idx[gw.mat.idx] )
dx[gw.x.idx] <- dx.GW
}
# dx is 1xnpar matrix of LIST (type != "free")
if(is.matrix(dx)) {
dx <- as.numeric(dx)
}
dx
}
# for testing purposes only
# computeDeltaNumerical <- function(lavmodel = NULL, GLIST = NULL, g = 1L) {
#
# # state or final?
# if(is.null(GLIST)) GLIST <- lavmodel@GLIST
#
# compute.moments <- function(x) {
# GLIST <- lav_model_x2GLIST(lavmodel = lavmodel, x=x, type="free")
# Sigma.hat <- computeSigmaHat(lavmodel = lavmodel, GLIST = GLIST)
# S.vec <- lav_matrix_vech(Sigma.hat[[g]])
# if(lavmodel@meanstructure) {
# Mu.hat <- computeMuHat(lavmodel = lavmodel, GLIST=GLIST)
# out <- c(Mu.hat[[g]], S.vec)
# } else {
# out <- S.vec
# }
# out
# }
#
# x <- lav_model_get_parameters(lavmodel = lavmodel, GLIST=GLIST, type="free")
# Delta <- lav_func_jacobian_complex(func=compute.moments, x = x)
#
# Delta
# }
### FIXME: should we here also:
### - weight for groups? (no, for now)
### - handle equality constraints? (yes, for now)
computeDelta <- function(lavmodel = NULL, GLIST. = NULL,
m.el.idx. = NULL, x.el.idx. = NULL,
ceq.simple = FALSE,
force.conditional.x.false = FALSE) {
representation <- lavmodel@representation
categorical <- lavmodel@categorical
if(.hasSlot(lavmodel, "correlation")) {
correlation <- lavmodel@correlation
} else {
correlation <- FALSE
}
conditional.x <- lavmodel@conditional.x
group.w.free <- lavmodel@group.w.free
nmat <- lavmodel@nmat
nblocks <- lavmodel@nblocks
nvar <- lavmodel@nvar
num.idx <- lavmodel@num.idx
th.idx <- lavmodel@th.idx
nexo <- lavmodel@nexo
parameterization <- lavmodel@parameterization
# number of thresholds per group (if any)
nth <- sapply(th.idx, function(x) sum(x > 0L))
# state or final?
if(is.null(GLIST.))
GLIST <- lavmodel@GLIST
else
GLIST <- GLIST.
# type = "free" or something else?
type <- "nonfree"
m.el.idx <- m.el.idx.; x.el.idx <- x.el.idx.
if(is.null(m.el.idx) && is.null(x.el.idx))
type <- "free"
# number of rows in DELTA.group
pstar <- integer(nblocks)
for(g in 1:nblocks) {
pstar[g] <- as.integer(nvar[g] * (nvar[g] + 1) / 2)
if(lavmodel@meanstructure) {
pstar[g] <- nvar[g] + pstar[g] # first the means, then sigma
}
if(categorical) {
pstar[g] <- pstar[g] - nvar[g] # remove variances
pstar[g] <- pstar[g] - nvar[g] # remove means
pstar[g] <- pstar[g] + nth[g] # add thresholds
pstar[g] <- pstar[g] + length(num.idx[[g]]) # add num means
pstar[g] <- pstar[g] + length(num.idx[[g]]) # add num vars
} else if(correlation) {
pstar[g] <- pstar[g] - nvar[g] # remove variances
}
if(conditional.x && nexo[g] > 0L) {
pstar[g] <- pstar[g] + (nvar[g] * nexo[g]) # add slopes
}
if(group.w.free) {
pstar[g] <- pstar[g] + 1L # add group weight
}
}
# number of columns in DELTA + m.el.idx/x.el.idx
if(type == "free") {
if(.hasSlot(lavmodel, "ceq.simple.only") && lavmodel@ceq.simple.only) {
NCOL <- lavmodel@nx.unco
} else {
NCOL <- lavmodel@nx.free
}
m.el.idx <- x.el.idx <- vector("list", length=length(GLIST))
for(mm in 1:length(GLIST)) {
m.el.idx[[mm]] <- lavmodel@m.free.idx[[mm]]
if(.hasSlot(lavmodel, "ceq.simple.only") &&
lavmodel@ceq.simple.only) {
x.el.idx[[mm]] <- lavmodel@x.unco.idx[[mm]]
} else {
x.el.idx[[mm]] <- lavmodel@x.free.idx[[mm]]
}
# handle symmetric matrices
if(lavmodel@isSymmetric[mm]) {
# since we use 'x.free.idx', only symmetric elements
# are duplicated (not the equal ones, only in x.free.free)
dix <- duplicated(x.el.idx[[mm]])
if(any(dix)) {
m.el.idx[[mm]] <- m.el.idx[[mm]][!dix]
x.el.idx[[mm]] <- x.el.idx[[mm]][!dix]
}
}
}
} else {
## FIXME: this does *not* take into account symmetric
## matrices; hence NCOL will be too large, and empty
## columns will be added
## this is ugly, but it doesn't hurt
## alternative could be:
## NCOL <- sum(unlist(lapply(x.el.idx, function(x) length(unique(x)))))
#NCOL <- sum(unlist(lapply(m.el.idx, length)))
NCOL <- sum(unlist(lapply(x.el.idx, function(x) length(unique(x)))))
# sanity check
#nx <- sum(unlist(lapply(x.el.idx, length)))
#stopifnot(NCOL == nx)
}
# compute Delta
Delta <- vector("list", length=nblocks)
for(g in 1:nblocks) {
Delta.group <- matrix(0, nrow=pstar[g], ncol=NCOL)
# which mm belong to group g?
mm.in.group <- 1:nmat[g] + cumsum(c(0,nmat))[g]
# label rows of Delta.group --- FIXME!!!
#if(categorical) {
# # 1. th (means interleaved?)
# # 2. pi
# # 3. var num + cor
#} else {
# if(meanstructure) {
# }
#}
#if(group.w.free) {
#}
# if theta, do some preparation
if(representation == "LISREL" && parameterization == "theta") {
sigma.hat <- computeSigmaHat.LISREL(MLIST=GLIST[mm.in.group],
delta=FALSE)
dsigma <- diag(sigma.hat)
# dcor/dcov for sigma
R <- lav_deriv_cov2cor(sigma.hat, num.idx = lavmodel@num.idx[[g]])
theta.var.idx <- lav_matrix_diagh_idx(nvar[g])
}
for(mm in mm.in.group) {
mname <- names(lavmodel@GLIST)[mm]
# skip empty ones
if(!length(m.el.idx[[mm]])) next
# get Delta columns for this model matrix
if(representation == "LISREL") {
# Sigma
DELTA <- dxSigma <-
derivative.sigma.LISREL(m = mname,
idx = m.el.idx[[mm]],
MLIST = GLIST[ mm.in.group ],
delta = parameterization == "delta")
if(categorical && parameterization == "theta") {
DELTA <- R %*% DELTA
}
if(categorical) {
# reorder: first variances (of numeric), then covariances
cov.idx <- lav_matrix_vech_idx(nvar[g])
covd.idx <- lav_matrix_vech_idx(nvar[g], diagonal = FALSE)
var.idx <- which(is.na(match(cov.idx,
covd.idx)))[num.idx[[g]]]
cor.idx <- match(covd.idx, cov.idx)
DELTA <- rbind(DELTA[var.idx,,drop=FALSE],
DELTA[cor.idx,,drop=FALSE])
}
# correlation structure?
if(!categorical && correlation) {
rm.idx <- lav_matrix_diagh_idx(nvar[g])
DELTA <- DELTA[-rm.idx, , drop = FALSE]
}
if(!categorical) {
if(conditional.x) {
# means/intercepts
DELTA.mu <- derivative.mu.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
# slopes
if(lavmodel@nexo[g] > 0L) {
DELTA.pi <- derivative.pi.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
if(lavmodel@multilevel) {
DELTA <- rbind(DELTA.mu, DELTA.pi, DELTA)
} else {
# ATTENTION: we need to change the order here
# lav_mvreg_scores_* uses 'Beta' where the
# the intercepts are just the first row
# using the col-major approach, we need to
# interweave the intercepts with the slopes!
nEls <- NROW(DELTA.mu) + NROW(DELTA.pi)
# = (nexo + 1 int) * nvar
# intercepts on top
tmp <- rbind(DELTA.mu, DELTA.pi)
# change row index
row.idx <- lav_matrix_vec(matrix(seq.int(nEls),
nrow = lavmodel@nexo[g] + 1L,
ncol = lavmodel@nvar[g], byrow = TRUE))
DELTA.beta <- tmp[row.idx,,drop = FALSE]
DELTA <- rbind(DELTA.beta, DELTA)
}
} else {
DELTA <- rbind(DELTA.mu, DELTA)
}
} else if(!conditional.x && lavmodel@meanstructure) {
DELTA.mu <- derivative.mu.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
DELTA <- rbind(DELTA.mu, DELTA)
}
}
else if(categorical) {
DELTA.th <- derivative.th.LISREL(m=mname,
idx=m.el.idx[[mm]],
th.idx=th.idx[[g]],
MLIST=GLIST[ mm.in.group ],
delta = TRUE)
if(parameterization == "theta") {
# dy/ddsigma = -0.5/(ddsigma*sqrt(ddsigma))
dDelta.dx <-
( dxSigma[theta.var.idx,,drop=FALSE] *
-0.5 / (dsigma*sqrt(dsigma)) )
dth.dDelta <-
derivative.th.LISREL(m = "delta",
idx = 1:nvar[g],
MLIST = GLIST[ mm.in.group ],
th.idx = th.idx[[g]])
# add dth.dDelta %*% dDelta.dx
no.num.idx <- which(th.idx[[g]] > 0)
DELTA.th[no.num.idx,] <-
DELTA.th[no.num.idx,,drop=FALSE] +
(dth.dDelta %*% dDelta.dx)[no.num.idx,,drop=FALSE]
}
if(conditional.x && lavmodel@nexo[g] > 0L) {
DELTA.pi <-
derivative.pi.LISREL(m=mname,
idx=m.el.idx[[mm]],
MLIST=GLIST[ mm.in.group ])
if(parameterization == "theta") {
dpi.dDelta <-
derivative.pi.LISREL(m = "delta",
idx = 1:nvar[g],
MLIST = GLIST[ mm.in.group ])
# add dpi.dDelta %*% dDelta.dx
no.num.idx <-
which(!seq.int(1L,nvar[g]) %in% num.idx[[g]])
no.num.idx <- rep(seq.int(0,nexo[g]-1) * nvar[g],
each=length(no.num.idx)) + no.num.idx
DELTA.pi[no.num.idx,] <-
DELTA.pi[no.num.idx,,drop=FALSE] +
(dpi.dDelta %*% dDelta.dx)[no.num.idx,,drop=FALSE]
}
DELTA <- rbind(DELTA.th, DELTA.pi, DELTA)
} else {
DELTA <- rbind(DELTA.th, DELTA)
}
}
if(group.w.free) {
DELTA.gw <- derivative.gw.LISREL(m=mname,
idx=m.el.idx[[mm]],
MLIST=GLIST[ mm.in.group ])
DELTA <- rbind(DELTA.gw, DELTA)
}
} else if(representation == "RAM") {
DELTA <- dxSigma <-
lav_ram_dsigma(m = mname,
idx = m.el.idx[[mm]],
MLIST = GLIST[ mm.in.group ])
if(lavmodel@meanstructure) {
DELTA.mu <- lav_ram_dmu(m = mname,
idx = m.el.idx[[mm]],
MLIST = GLIST[ mm.in.group ])
DELTA <- rbind(DELTA.mu, DELTA)
}
} else {
stop("representation ", representation, " not implemented yet")
}
Delta.group[ ,x.el.idx[[mm]]] <- DELTA
} # mm
# if type == "free" take care of equality constraints
if(type == "free" && ceq.simple &&
.hasSlot(lavmodel, "ceq.simple.only") && lavmodel@ceq.simple.only) {
Delta.group <- Delta.group %*% lavmodel@ceq.simple.K
}
Delta[[g]] <- Delta.group
} # g
# if multilevel, rbind levels within group
if(.hasSlot(lavmodel, "multilevel") && lavmodel@multilevel) {
DELTA <- vector("list", length = lavmodel@ngroups)
for(g in 1:lavmodel@ngroups) {
DELTA[[g]] <- rbind( Delta[[(g-1)*2 + 1]],
Delta[[(g-1)*2 + 2]] )
}
Delta <- DELTA
}
Delta
}
computeDeltaDx <- function(lavmodel = NULL, GLIST = NULL, target = "lambda",
ceq.simple = FALSE) {
# state or final?
if(is.null(GLIST)) GLIST <- lavmodel@GLIST
representation <- lavmodel@representation
nmat <- lavmodel@nmat
nblocks <- lavmodel@nblocks
num.idx <- lavmodel@num.idx
th.idx <- lavmodel@th.idx
# number of columns in DELTA + m.el.idx/x.el.idx
type <- "free"
#if(type == "free") {
if(.hasSlot(lavmodel, "ceq.simple.only") && lavmodel@ceq.simple.only) {
NCOL <- lavmodel@nx.unco
} else {
NCOL <- lavmodel@nx.free
}
m.el.idx <- x.el.idx <- vector("list", length=length(GLIST))
for(mm in 1:length(GLIST)) {
m.el.idx[[mm]] <- lavmodel@m.free.idx[[mm]]
if(.hasSlot(lavmodel, "ceq.simple.only") &&
lavmodel@ceq.simple.only) {
x.el.idx[[mm]] <- lavmodel@x.unco.idx[[mm]]
} else {
x.el.idx[[mm]] <- lavmodel@x.free.idx[[mm]]
}
# handle symmetric matrices
if(lavmodel@isSymmetric[mm]) {
# since we use 'x.free.idx', only symmetric elements
# are duplicated (not the equal ones, only in x.free.free)
dix <- duplicated(x.el.idx[[mm]])
if(any(dix)) {
m.el.idx[[mm]] <- m.el.idx[[mm]][!dix]
x.el.idx[[mm]] <- x.el.idx[[mm]][!dix]
}
}
}
#} else {
# NCOL <- sum(unlist(lapply(x.el.idx, function(x) length(unique(x)))))
#}
# compute Delta per group
Delta <- vector("list", length=nblocks)
for(g in 1:nblocks) {
mm.in.group <- 1:nmat[g] + cumsum(c(0,nmat))[g]
Delta.group <- NULL
for(mm in mm.in.group) {
mname <- names(lavmodel@GLIST)[mm]
# skip empty ones
if(!length(m.el.idx[[mm]])) next
# get Delta columns for this model matrix
if(representation == "LISREL") {
if(target == "lambda") {
DELTA <- derivative.lambda.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
} else if(target == "th") {
DELTA <- derivative.th.LISREL(m=mname, th.idx = th.idx[[g]],
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ],
delta=TRUE)
} else if(target == "mu") {
DELTA <- derivative.mu.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
} else if(target == "nu") {
DELTA <- derivative.nu.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
} else if(target == "tau") {
DELTA <- derivative.tau.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
} else if(target == "theta") {
DELTA <- derivative.theta.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
} else if(target == "gamma") {
DELTA <- derivative.gamma.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
} else if(target == "beta") {
DELTA <- derivative.beta.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
} else if(target == "alpha") {
DELTA <- derivative.alpha.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
} else if(target == "psi") {
DELTA <- derivative.psi.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
} else if(target == "sigma") {
DELTA <- derivative.sigma.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ],
delta=TRUE)
} else {
stop("lavaan ERROR: target ", target, " not implemented yet")
}
# initialize?
if(is.null(Delta.group)) {
Delta.group <- matrix(0, nrow=nrow(DELTA), ncol=NCOL)
}
Delta.group[ ,x.el.idx[[mm]]] <- DELTA
}
} # mm
if(type == "free" && ceq.simple &&
.hasSlot(lavmodel, "ceq.simple.only") && lavmodel@ceq.simple.only) {
Delta.group <- Delta.group %*% lavmodel@ceq.simple.K
}
Delta[[g]] <- Delta.group
} # g
Delta
}
computeOmega <- function(Sigma.hat=NULL, Mu.hat=NULL,
lavsamplestats=NULL, estimator="ML",
meanstructure=FALSE, conditional.x = FALSE) {
# nblocks
nblocks <- length(Sigma.hat)
Omega <- vector("list", length = nblocks)
Omega.mu <- vector("list", length = nblocks)
for(g in 1:nblocks) {
# ML
if(estimator %in% c("ML", "REML", "catML")) {
if(attr(Sigma.hat[[g]], "po") == FALSE) {
# FIXME: WHAT IS THE BEST THING TO DO HERE??
# CURRENTLY: stop
warning("lav_model_gradient: Sigma.hat is not positive definite\n")
Sigma.hat.inv <- MASS::ginv(Sigma.hat[[g]])
Sigma.hat.log.det <- log(.Machine$double.eps)
} else {
Sigma.hat.inv <- attr(Sigma.hat[[g]], "inv")
Sigma.hat.log.det <- attr(Sigma.hat[[g]], "log.det")
}
if(!lavsamplestats@missing.flag) { # complete data
if(meanstructure) {
if(conditional.x) {
diff <- lavsamplestats@res.int[[g]] - Mu.hat[[g]]
W.tilde <- lavsamplestats@res.cov[[g]] + tcrossprod(diff)
} else {
diff <- lavsamplestats@mean[[g]] - Mu.hat[[g]]
W.tilde <- lavsamplestats@cov[[g]] + tcrossprod(diff)
}
# Browne 1995 eq 4.55
Omega.mu[[g]] <- t(t(diff) %*% Sigma.hat.inv)
Omega[[g]] <-
( Sigma.hat.inv %*% (W.tilde - Sigma.hat[[g]]) %*%
Sigma.hat.inv )
} else {
if(conditional.x) {
W.tilde <- lavsamplestats@res.cov[[g]]
} else {
W.tilde <- lavsamplestats@cov[[g]]
}
Omega[[g]] <-
( Sigma.hat.inv %*% (W.tilde - Sigma.hat[[g]]) %*%
Sigma.hat.inv )
}
} else { # missing data
M <- lavsamplestats@missing[[g]]
nvar <- ncol(lavsamplestats@cov[[g]])
OMEGA <- matrix(0, nvar, nvar)
OMEGA.MU <- matrix(0, nvar, 1)
for(p in 1:length(M)) {
SX <- M[[p]][["SY"]]
MX <- M[[p]][["MY"]]
nobs <- M[[p]][["freq"]]
var.idx <- M[[p]][["var.idx"]]
Sigma.inv <- inv.chol(Sigma.hat[[g]][var.idx, var.idx],
logdet=FALSE)
Mu <- Mu.hat[[g]][var.idx]
W.tilde <- SX + tcrossprod(MX - Mu)
OMEGA.MU[var.idx, 1] <-
( OMEGA.MU[var.idx, 1] + nobs/lavsamplestats@ntotal *
t(t(MX - Mu) %*% Sigma.inv) )
OMEGA[var.idx, var.idx] <-
( OMEGA[var.idx, var.idx] + nobs/lavsamplestats@ntotal *
(Sigma.inv %*%
(W.tilde - Sigma.hat[[g]][var.idx,var.idx]) %*%
Sigma.inv ) )
}
Omega.mu[[g]] <- OMEGA.MU
Omega[[g]] <- OMEGA
} # missing
# GLS
} else if(estimator == "GLS") {
W.inv <- lavsamplestats@icov[[g]]
W <- lavsamplestats@cov[[g]]
Omega[[g]] <- (lavsamplestats@nobs[[g]]-1)/lavsamplestats@nobs[[g]] *
(W.inv %*% (W - Sigma.hat[[g]]) %*% W.inv)
if(meanstructure) {
diff <- as.matrix(lavsamplestats@mean[[g]] - Mu.hat[[g]])
Omega.mu[[g]] <- t( t(diff) %*% W.inv )
}
}
} # g
if(meanstructure) attr(Omega, "mu") <- Omega.mu
Omega
}
lav_model_gradient_DD <- function(lavmodel, GLIST = NULL, group = 1L) {
if(is.null(GLIST)) GLIST <- lavmodel@GLIST
#### FIX th + mu!!!!!
Delta.lambda <- computeDeltaDx(lavmodel, GLIST=GLIST, target="lambda")[[group]]
Delta.tau <- computeDeltaDx(lavmodel, GLIST=GLIST, target="tau" )[[group]]
Delta.nu <- computeDeltaDx(lavmodel, GLIST=GLIST, target="nu" )[[group]]
Delta.theta <- computeDeltaDx(lavmodel, GLIST=GLIST, target="theta" )[[group]]
Delta.beta <- computeDeltaDx(lavmodel, GLIST=GLIST, target="beta" )[[group]]
Delta.psi <- computeDeltaDx(lavmodel, GLIST=GLIST, target="psi" )[[group]]
Delta.alpha <- computeDeltaDx(lavmodel, GLIST=GLIST, target="alpha" )[[group]]
Delta.gamma <- computeDeltaDx(lavmodel, GLIST=GLIST, target="gamma" )[[group]]
ov.y.dummy.ov.idx <- lavmodel@ov.y.dummy.ov.idx[[group]]
ov.x.dummy.ov.idx <- lavmodel@ov.x.dummy.ov.idx[[group]]
ov.y.dummy.lv.idx <- lavmodel@ov.y.dummy.lv.idx[[group]]
ov.x.dummy.lv.idx <- lavmodel@ov.x.dummy.lv.idx[[group]]
ov.dummy.idx <- c(ov.y.dummy.ov.idx, ov.x.dummy.ov.idx)
lv.dummy.idx <- c(ov.y.dummy.lv.idx, ov.x.dummy.lv.idx)
th.idx <- lavmodel@th.idx[[group]]
num.idx <- lavmodel@num.idx[[group]]
ord.idx <- unique( th.idx[th.idx > 0L] )
# fix Delta's...
mm.in.group <- 1:lavmodel@nmat[group] + cumsum(c(0,lavmodel@nmat))[group]
MLIST <- GLIST[ mm.in.group ]
DD <- list()
nvar <- lavmodel@nvar
nfac <- ncol(MLIST$lambda) - length(lv.dummy.idx)
# DD$theta
theta.idx <- lav_matrix_diagh_idx(nvar)
DD$theta <- Delta.theta[theta.idx,,drop=FALSE]
if(length(ov.dummy.idx) > 0L) {
psi.idx <- lav_matrix_diagh_idx( ncol(MLIST$psi) )[lv.dummy.idx]
DD$theta[ov.dummy.idx,] <- Delta.psi[psi.idx,,drop=FALSE]
}
# num only? FIXME or just all of them?
DD$theta <- DD$theta[num.idx,,drop=FALSE]
# DD$nu
DD$nu <- Delta.nu
if(length(ov.dummy.idx) > 0L) {
DD$nu[ov.dummy.idx,] <- Delta.alpha[lv.dummy.idx,]
}
DD$nu <- DD$nu[num.idx,,drop=FALSE] # needed?
# DD$lambda
nr <- nvar; nc <- nfac
lambda.idx <- nr*((1:nc) - 1L) + rep(1:nvar, each=nc)
DD$lambda <- Delta.lambda[lambda.idx,,drop=FALSE]
if(length(ov.dummy.idx) > 0L) {
nr <- nrow(MLIST$beta); nc <- nfac # only the first 1:nfac columns
# beta.idx <- rep(nr*((1:nc) - 1L), each=length(lv.dummy.idx)) + rep(lv.dummy.idx, times=nc) ## FIXME
beta.idx <- rep(nr*((1:nc) - 1L), times=length(lv.dummy.idx)) + rep(lv.dummy.idx, each=nc)
#l.idx <- inr*((1:nc) - 1L) + rep(ov.dummy.idx, each=nc) ## FIXME
# l.idx <- rep(nr*((1:nc) - 1L), each=length(ov.dummy.idx)) + rep(ov.dummy.idx, times=nc)
l.idx <- rep(nr*((1:nc) - 1L), times=length(ov.dummy.idx)) + rep(ov.dummy.idx, each=nc)
DD$lambda[match(l.idx, lambda.idx),] <- Delta.beta[beta.idx,,drop=FALSE]
}
# DD$KAPPA
DD$kappa <- Delta.gamma
if(length(ov.dummy.idx) > 0L) {
nr <- nrow(MLIST$gamma); nc <- ncol(MLIST$gamma)
kappa.idx <- nr*((1:nc) - 1L) + rep(lv.dummy.idx, each=nc)
DD$kappa <- DD$kappa[kappa.idx,,drop=FALSE]
}
# DD$GAMMA
if(!is.null(MLIST$gamma)) {
nr <- nrow(MLIST$gamma); nc <- ncol(MLIST$gamma)
lv.idx <- 1:nfac
# MUST BE ROWWISE!
gamma.idx <- rep(nr*((1:nc) - 1L), times=length(lv.idx)) + rep(lv.idx, each=nc)
DD$gamma <- Delta.gamma[gamma.idx,,drop=FALSE]
}
# DD$BETA
if(!is.null(MLIST$beta)) {
nr <- nc <- nrow(MLIST$beta)
lv.idx <- 1:nfac
# MUST BE ROWWISE!
beta.idx <- rep(nr*((1:nfac) - 1L), times=nfac) + rep(lv.idx, each=nfac)
DD$beta <- Delta.beta[beta.idx,,drop=FALSE]
}
## DD$psi
DD$psi <- Delta.psi
if(length(lv.dummy.idx) > 0L) {
nr <- nc <- nrow(MLIST$psi)
lv.idx <- 1:nfac
# MUST BE ROWWISE!
psi.idx <- rep(nr*((1:nfac) - 1L), times=nfac) + rep(lv.idx, each=nfac)
DD$psi <- DD$psi[psi.idx,,drop=FALSE]
}
## DD$tau
if(!is.null(MLIST$tau)) {
DD$tau <- Delta.tau
}
DD
}
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_gradient_DD computeOmega computeDeltaDx computeDelta lav_model_gradient
# model gradient
lav_model_gradient <- function(lavmodel = NULL,
GLIST = NULL,
lavsamplestats = NULL,
lavdata = NULL,
lavcache = NULL,
type = "free",
verbose = FALSE,
group.weight = TRUE,
Delta = NULL,
m.el.idx = NULL,
x.el.idx = NULL,
ceq.simple = FALSE) {
nmat <- lavmodel@nmat
estimator <- lavmodel@estimator
representation <- lavmodel@representation
meanstructure <- lavmodel@meanstructure
categorical <- lavmodel@categorical
group.w.free <- lavmodel@group.w.free
fixed.x <- lavmodel@fixed.x
conditional.x <- lavmodel@conditional.x
num.idx <- lavmodel@num.idx
th.idx <- lavmodel@th.idx
nx.free <- lavmodel@nx.free
if(.hasSlot(lavmodel, "estimator.args")) {
estimator.args <- lavmodel@estimator.args
} else {
estimator.args <- list()
}
# state or final?
if(is.null(GLIST)) GLIST <- lavmodel@GLIST
if(estimator == "REML") warning("analytical gradient not implement; use numerical approximation")
# group.weight
# FIXME --> block.weight
if(group.weight) {
if(estimator %in% c("ML","PML","FML","MML","REML","NTRLS","catML")) {
group.w <- (unlist(lavsamplestats@nobs)/lavsamplestats@ntotal)
} else if(estimator == "DLS") {
if(estimator.args$dls.FtimesNminus1) {
group.w <- ((unlist(lavsamplestats@nobs)-1)/lavsamplestats@ntotal)
} else {
group.w <- (unlist(lavsamplestats@nobs)/lavsamplestats@ntotal)
}
} else {
# FIXME: double check!
group.w <- ((unlist(lavsamplestats@nobs)-1)/lavsamplestats@ntotal)
}
} else {
group.w <- rep(1.0, lavmodel@nblocks)
}
# do we need WLS.est?
if(estimator %in% c("WLS", "DWLS", "ULS", "GLS", "NTRLS", "DLS")) {
# always compute WLS.est
WLS.est <- lav_model_wls_est(lavmodel = lavmodel, GLIST = GLIST) #,
# cov.x = lavsamplestats@cov.x)
}
if(estimator %in% c("ML", "PML", "FML", "REML", "NTRLS", "catML")) {
# compute moments for all groups
#if(conditional.x) {
# Sigma.hat <- computeSigmaHatJoint(lavmodel = lavmodel,
# GLIST = GLIST,
# extra = (estimator %in% c("ML", "REML","NTRLS")))
#} else {
Sigma.hat <- computeSigmaHat(lavmodel = lavmodel, GLIST = GLIST,
extra = (estimator %in% c("ML", "REML",
"NTRLS", "catML")))
#}
if(meanstructure) {
#if(conditional.x) {
# Mu.hat <- computeMuHat(lavmodel = lavmodel, GLIST = GLIST)
#} else {
Mu.hat <- computeMuHat(lavmodel = lavmodel, GLIST = GLIST)
#}
}
if(categorical) {
TH <- computeTH(lavmodel = lavmodel, GLIST = GLIST)
}
if(conditional.x) {
PI <- computePI(lavmodel = lavmodel, GLIST = GLIST)
} else if(estimator == "PML") {
PI <- vector("list", length = lavmodel@nblocks)
}
if(group.w.free) {
GW <- computeGW(lavmodel = lavmodel, GLIST = GLIST)
}
} else if(estimator == "DLS" && estimator.args$dls.GammaNT == "model") {
Sigma.hat <- computeSigmaHat(lavmodel = lavmodel, GLIST = GLIST,
extra = FALSE)
Mu.hat <- computeMuHat(lavmodel = lavmodel, GLIST = GLIST)
} else if(estimator == "MML") {
TH <- computeTH( lavmodel = lavmodel, GLIST = GLIST)
THETA <- computeTHETA(lavmodel = lavmodel, GLIST = GLIST)
GW <- computeGW( lavmodel = lavmodel, GLIST = GLIST)
}
# four approaches (FIXME!!!! merge this!)
# - ML approach: using Omega (and Omega.mu)
# Omega = 'POST' = Sigma.inv %*% (S - Sigma) %*% t(Sigma.inv)
# (still 2x faster than Delta method)
# - WLS/DWLS/GLS: using Delta + WLS.V; support for fixed.x, conditional.x
# - (ML)/NTRLS: using Delta, no support for fixed.x, conditional.x
# - PML/FML/MML: custom
# 1. ML approach
if( (estimator == "ML" || estimator == "REML" || estimator == "catML") &&
lavdata@nlevels == 1L &&
!lavmodel@conditional.x ) {
if(meanstructure) {
Omega <- computeOmega(Sigma.hat=Sigma.hat, Mu.hat=Mu.hat,
lavsamplestats=lavsamplestats,
estimator=estimator,
meanstructure=TRUE,
conditional.x = conditional.x)
Omega.mu <- attr(Omega, "mu")
} else {
Omega <- computeOmega(Sigma.hat=Sigma.hat, Mu.hat=NULL,
lavsamplestats=lavsamplestats,
estimator=estimator,
meanstructure=FALSE,
conditional.x = conditional.x)
Omega.mu <- vector("list", length = lavmodel@nblocks)
}
# compute DX (for all elements in every model matrix)
DX <- vector("list", length=length(GLIST))
for(g in 1:lavmodel@nblocks) {
# which mm belong to group g?
mm.in.group <- 1:nmat[g] + cumsum(c(0,nmat))[g]
mm.names <- names( GLIST[mm.in.group] )
if(representation == "LISREL") {
DX.group <- derivative.F.LISREL(GLIST[mm.in.group],
Omega[[g]],
Omega.mu[[g]])
# FIXME!!!
# add empty gamma
if(lavmodel@conditional.x) {
DX.group$gamma <- lavmodel@GLIST$gamma
}
# only save what we need
DX[mm.in.group] <- DX.group[ mm.names ]
} else if(representation == "RAM") {
DX.group <- lav_ram_df(GLIST[mm.in.group],
Omega[[g]],
Omega.mu[[g]])
# only save what we need
DX[mm.in.group] <- DX.group[ mm.names ]
} else {
stop("only LISREL and RAM representation has been implemented for now")
}
# weight by group
if(lavmodel@nblocks > 1L) {
for(mm in mm.in.group) {
DX[[mm]] <- group.w[g] * DX[[mm]]
}
}
}
# extract free parameters
if(type == "free") {
if(.hasSlot(lavmodel, "ceq.simple.only") &&
lavmodel@ceq.simple.only) { # new in 0.6-11
dx <- numeric( lavmodel@nx.unco )
for(g in 1:lavmodel@nblocks) {
mm.in.group <- 1:nmat[g] + cumsum(c(0,nmat))[g]
for(mm in mm.in.group) {
m.free.idx <- lavmodel@m.free.idx[[mm]]
x.unco.idx <- lavmodel@x.unco.idx[[mm]]
dx[x.unco.idx] <- DX[[mm]][m.free.idx]
}
}
if(ceq.simple) {
dx <- drop( crossprod(lavmodel@ceq.simple.K, dx) )
}
} else {
dx <- numeric( nx.free )
for(g in 1:lavmodel@nblocks) {
mm.in.group <- 1:nmat[g] + cumsum(c(0,nmat))[g]
for(mm in mm.in.group) {
m.free.idx <- lavmodel@m.free.idx[[mm]]
x.free.idx <- lavmodel@x.free.idx[[mm]]
dx[x.free.idx] <- DX[[mm]][m.free.idx]
}
}
}
} else {
dx <- DX
# handle equality constraints
### FIXME!!!! TODO!!!!
}
} else # ML
# 2. using Delta - *LS family
if(estimator %in% c("WLS", "DWLS", "ULS", "GLS", "NTGLS", "DLS")) {
if(type != "free") {
if(is.null(Delta))
stop("FIXME: Delta should be given if type != free")
#stop("FIXME: WLS gradient with type != free needs fixing!")
} else {
Delta <- computeDelta(lavmodel = lavmodel, GLIST. = GLIST,
ceq.simple = ceq.simple)
}
for(g in 1:lavmodel@nblocks) {
#diff <- as.matrix(lavsamplestats@WLS.obs[[g]] - WLS.est[[g]])
#group.dx <- -1 * ( t(Delta[[g]]) %*% lavsamplestats@WLS.V[[g]] %*% diff)
# 0.5-17: use crossprod twice; treat DWLS/ULS special
if(estimator == "WLS" ||
estimator == "GLS" ||
estimator == "DLS" ||
estimator == "NTRLS") {
# full weight matrix
diff <- lavsamplestats@WLS.obs[[g]] - WLS.est[[g]]
# full weight matrix
if(estimator == "GLS" || estimator == "WLS") {
WLS.V <- lavsamplestats@WLS.V[[g]]
group.dx <- -1 * crossprod(Delta[[g]],
crossprod(WLS.V, diff))
} else if(estimator == "DLS") {
if(estimator.args$dls.GammaNT == "sample") {
WLS.V <- lavsamplestats@WLS.V[[g]] # for now
} else {
dls.a <- estimator.args$dls.a
GammaNT <- lav_samplestats_Gamma_NT(
COV = Sigma.hat[[g]],
MEAN = Mu.hat[[g]],
rescale = FALSE,
x.idx = lavsamplestats@x.idx[[g]],
fixed.x = lavmodel@fixed.x,
conditional.x = lavmodel@conditional.x,
meanstructure = lavmodel@meanstructure,
slopestructure = lavmodel@conditional.x)
W.DLS <- (1 - dls.a)*lavsamplestats@NACOV[[g]] + dls.a*GammaNT
WLS.V <- lav_matrix_symmetric_inverse(W.DLS)
}
group.dx <- -1 * crossprod(Delta[[g]],
crossprod(WLS.V, diff))
} else if(estimator == "NTRLS") {
stopifnot(!conditional.x)
#WLS.V <- lav_samplestats_Gamma_inverse_NT(
# ICOV = attr(Sigma.hat[[g]],"inv")[,,drop=FALSE],
# COV = Sigma.hat[[g]][,,drop=FALSE],
# MEAN = Mu.hat[[g]],
# x.idx = lavsamplestats@x.idx[[g]],
# fixed.x = fixed.x,
# conditional.x = conditional.x,
# meanstructure = meanstructure,
# slopestructure = conditional.x)
S <- lavsamplestats@cov[[g]]
Sigma <- Sigma.hat[[g]]
Sigma.inv <- attr(Sigma, "inv")
nvar <- NROW(Sigma)
if(meanstructure) {
MEAN <- lavsamplestats@mean[[g]]; Mu <- Mu.hat[[g]]
POST.Sigma <- lav_matrix_duplication_pre(
matrix((Sigma.inv %*% (S - Sigma) %*% t(Sigma.inv)) %*%
(diag(nvar) + (S - Sigma) %*% Sigma.inv) +
(Sigma.inv %*% tcrossprod(MEAN - Mu) %*% Sigma.inv),
ncol = 1) )
POST.Mu <- as.numeric(2 * Sigma.inv %*% (MEAN - Mu))
POST <- c(POST.Mu, POST.Sigma)
} else {
POST <- lav_matrix_duplication_pre(
matrix((Sigma.inv %*% (S - Sigma) %*% t(Sigma.inv)) %*%
(diag(nvar) + (S - Sigma) %*% Sigma.inv), ncol = 1))
}
group.dx <- as.numeric( -1 * crossprod(Delta[[g]], POST) )
}
} else
if(estimator == "DWLS" || estimator == "ULS") {
# diagonal weight matrix
diff <- lavsamplestats@WLS.obs[[g]] - WLS.est[[g]]
group.dx <- -1 * crossprod(Delta[[g]],
lavsamplestats@WLS.VD[[g]] * diff)
}
group.dx <- group.w[g] * group.dx
if(g == 1) {
dx <- group.dx
} else {
dx <- dx + group.dx
}
} # g
if(type == "free") {
# nothing to do
} else {
# make a GLIST
dx <- lav_model_x2GLIST(lavmodel = lavmodel, x = dx,
type = "custom", setDelta = FALSE,
m.el.idx = m.el.idx,
x.el.idx = x.el.idx)
}
} # WLS
# ML + conditional.x
else if(estimator %in% c("ML", "catML") && lavmodel@conditional.x
&& lavdata@nlevels == 1L) {
if(type != "free") {
if(is.null(Delta))
stop("FIXME: Delta should be given if type != free")
#stop("FIXME: WLS gradient with type != free needs fixing!")
} else {
Delta <- computeDelta(lavmodel = lavmodel, GLIST. = GLIST,
ceq.simple = ceq.simple)
}
for(g in 1:lavmodel@nblocks) {
# augmented mean.x + cov.x matrix
mean.x <- lavsamplestats@mean.x[[g]]
cov.x <- lavsamplestats@cov.x[[g]]
C3 <- rbind(c(1,mean.x),
cbind(mean.x, cov.x + tcrossprod(mean.x)))
Sigma <- Sigma.hat[[g]]
Mu.g <- Mu.hat[[g]]
PI.g <- PI[[g]]
Sigma.inv <- attr(Sigma, "inv")
nvar <- NROW(Sigma)
S <- lavsamplestats@res.cov[[g]]
# beta
OBS <- t( cbind(lavsamplestats@res.int[[g]],
lavsamplestats@res.slopes[[g]]) )
EST <- t( cbind(Mu.g, PI.g) )
#obs.beta <- c(lavsamplestats@res.int[[g]],
# lav_matrix_vec(lavsamplestats@res.slopes[[g]]))
#est.beta <- c(Mu.g, lav_matrix_vec(PI.g))
#beta.COV <- C3 %x% Sigma.inv
#a <- t(obs.beta - est.beta)
#b <- as.matrix(obs.beta - est.beta)
#K <- lav_matrix_commutation(m = nvar, n = nvar)
#AB <- (K %x% diag(NROW(C3)*NROW(C3))) %*%
# (diag(nvar) %x% lav_matrix_vec(C3) %x% diag(nvar))
#K <- lav_matrix_commutation(m = nvar, n = NROW(C3))
#AB <- ( diag(NROW(C3)) %x% K %x% diag(nvar) ) %*%
# (lav_matrix_vec(C3) %x% diag( nvar * nvar) )
#POST.beta <- 2 * beta.COV %*% (obs.beta - est.beta)
d.BETA <- C3 %*% (OBS - EST) %*% Sigma.inv
# NOTE: the vecr here, unlike lav_mvreg_dlogl_beta
# this is because DELTA has used vec(t(BETA)),
# instead of vec(BETA)
#POST.beta <- 2 * lav_matrix_vecr(d.BETA)
# NOT any longer, since 0.6-1!!!
POST.beta <- 2 * lav_matrix_vec(d.BETA)
#POST.sigma1 <- lav_matrix_duplication_pre(
# (Sigma.inv %x% Sigma.inv) %*% t(AB) %*% (t(a) %x% b) )
# Sigma
#POST.sigma2 <- lav_matrix_duplication_pre(
# matrix( lav_matrix_vec(
# Sigma.inv %*% (S - Sigma) %*% t(Sigma.inv)), ncol = 1L))
W.tilde <- S + t(OBS - EST) %*% C3 %*% (OBS - EST)
d.SIGMA <- (Sigma.inv - Sigma.inv %*% W.tilde %*% Sigma.inv)
d.vechSigma <- as.numeric( lav_matrix_duplication_pre(
as.matrix(lav_matrix_vec(d.SIGMA)) ) )
POST.sigma <- -1 * d.vechSigma
#POST <- c(POST.beta, POST.sigma1 + POST.sigma2)
POST <- c(POST.beta, POST.sigma)
group.dx <- as.numeric( -1 * crossprod(Delta[[g]], POST) )
# because we still use obj/2, we need to divide by 2!
group.dx <- group.dx / 2 # fixed in 0.6-1
group.dx <- group.w[g] * group.dx
if(g == 1) {
dx <- group.dx
} else {
dx <- dx + group.dx
}
} # g
if(type == "free") {
# nothing to do
} else {
# make a GLIST
dx <- lav_model_x2GLIST(lavmodel = lavmodel, x = dx,
type = "custom", setDelta = FALSE,
m.el.idx = m.el.idx,
x.el.idx = x.el.idx)
}
} # ML + conditional.x
else if(estimator == "ML" && lavdata@nlevels > 1L) {
if(type != "free") {
stop("FIXME: type != free in lav_model_gradient for estimator ML for nlevels > 1")
} else {
Delta <- computeDelta(lavmodel = lavmodel, GLIST. = GLIST,
ceq.simple = ceq.simple)
}
# for each upper-level group....
for(g in 1:lavmodel@ngroups) {
if(!lavsamplestats@missing.flag) { # complete data
if(lavmodel@conditional.x) {
DX <- lav_mvreg_cluster_dlogl_2l_samplestats(
YLp = lavsamplestats@YLp[[g]],
Lp = lavdata@Lp[[g]],
Res.Sigma.W = Sigma.hat[[(g-1)*2 + 1]],
Res.Int.W = Mu.hat[[(g-1)*2 + 1]],
Res.Pi.W = PI[[(g-1)*2 + 1]],
Res.Sigma.B = Sigma.hat[[(g-1)*2 + 2]],
Res.Int.B = Mu.hat[[(g-1)*2 + 2]],
Res.Pi.B = PI[[(g-1)*2 + 2]],
Sinv.method = "eigen")
} else {
DX <- lav_mvnorm_cluster_dlogl_2l_samplestats(
YLp = lavsamplestats@YLp[[g]],
Lp = lavdata@Lp[[g]],
Mu.W = Mu.hat[[(g-1)*2 + 1]],
Sigma.W = Sigma.hat[[(g-1)*2 + 1]],
Mu.B = Mu.hat[[(g-1)*2 + 2]],
Sigma.B = Sigma.hat[[(g-1)*2 + 2]],
Sinv.method = "eigen")
}
} else {
# missing data
if(lavmodel@conditional.x) {
stop("lavaan ERROR: gradient for twolevel + conditional.x + fiml is not ready; use optim.gradient = \"numerical\"")
} else {
DX <- lav_mvnorm_cluster_missing_dlogl_2l_samplestats(
Y1 = lavdata@X[[g]],
Y2 = lavsamplestats@YLp[[g]][[2]]$Y2,
Lp = lavdata@Lp[[g]],
Mp = lavdata@Mp[[g]],
Mu.W = Mu.hat[[(g-1)*2 + 1]],
Sigma.W = Sigma.hat[[(g-1)*2 + 1]],
Mu.B = Mu.hat[[(g-1)*2 + 2]],
Sigma.B = Sigma.hat[[(g-1)*2 + 2]],
Sinv.method = "eigen")
}
}
group.dx <- as.numeric( DX %*% Delta[[g]] )
# group weights (if any)
group.dx <- group.w[g] * group.dx
if(g == 1) {
dx <- group.dx
} else {
dx <- dx + group.dx
}
} # g
# divide by 2 * N
dx <- dx / (2 * lavsamplestats@ntotal)
#cat("dx1 (numerical) = \n"); print( zapsmall(dx1) )
#cat("dx (analytic) = \n"); print( zapsmall(dx ) )
} # ML + two-level
else if(estimator == "PML" || estimator == "FML" ||
estimator == "MML") {
if(type != "free") {
stop("FIXME: type != free in lav_model_gradient for estimator PML")
} else {
Delta <- computeDelta(lavmodel = lavmodel, GLIST. = GLIST,
ceq.simple = ceq.simple)
}
for(g in 1:lavmodel@nblocks) {
#print(GLIST)
#print(lav_model_get_parameters(lavmodel = lavmodel, GLIST = GLIST))
#print(Sigma.hat[[g]])
#print(TH[[g]])
#cat("*****\n")
# compute partial derivative of logLik with respect to
# thresholds/means, slopes, variances, correlations
if(estimator == "PML") {
if(lavdata@nlevels > 1L) {
stop("lavaan ERROR: PL gradient + multilevel not implemented; try optim.gradient = \"numerical\"")
} else if(conditional.x) {
d1 <- pml_deriv1(Sigma.hat = Sigma.hat[[g]],
Mu.hat = Mu.hat[[g]],
TH = TH[[g]],
th.idx = th.idx[[g]],
num.idx = num.idx[[g]],
X = lavdata@X[[g]],
lavcache = lavcache[[g]],
eXo = lavdata@eXo[[g]],
PI = PI[[g]],
missing = lavdata@missing)
} else {
d1 <- pml_deriv1(Sigma.hat = Sigma.hat[[g]],
Mu.hat = Mu.hat[[g]],
TH = TH[[g]],
th.idx = th.idx[[g]],
num.idx = num.idx[[g]],
X = lavdata@X[[g]],
lavcache = lavcache[[g]],
eXo = NULL,
PI = NULL,
missing = lavdata@missing)
} # not conditional.x
# chain rule (fmin)
group.dx <-
as.numeric(t(d1) %*% Delta[[g]])
} # PML
else if(estimator == "FML") {
d1 <- fml_deriv1(Sigma.hat = Sigma.hat[[g]],
TH = TH[[g]],
th.idx = th.idx[[g]],
num.idx = num.idx[[g]],
X = lavdata@X[[g]],
lavcache = lavcache[[g]])
# chain rule (fmin)
group.dx <-
as.numeric(t(d1) %*% Delta[[g]])/lavsamplestats@nobs[[g]]
} else if(estimator == "MML") {
group.dx <-
lav_model_gradient_mml(lavmodel = lavmodel,
GLIST = GLIST,
THETA = THETA[[g]],
TH = TH[[g]],
group = g,
lavdata = lavdata,
sample.mean = lavsamplestats@mean[[g]],
sample.mean.x = lavsamplestats@mean.x[[g]],
lavcache = lavcache)
}
# group weights (if any)
group.dx <- group.w[g] * group.dx
if(g == 1) {
dx <- group.dx
} else {
dx <- dx + group.dx
}
} # g
} else {
stop("lavaan ERROR: no analytical gradient available for estimator ",
estimator)
}
# group.w.free for ML
if(lavmodel@group.w.free &&
estimator %in% c("ML","MML","FML","PML","REML","catML")) {
#est.prop <- unlist( computeGW(lavmodel = lavmodel, GLIST = GLIST) )
#obs.prop <- unlist(lavsamplestats@group.w)
# FIXME: G2 based -- ML and friends only!!
#dx.GW <- - (obs.prop - est.prop)
# poisson version
est.freq <- exp(unlist(computeGW(lavmodel = lavmodel, GLIST = GLIST)))
obs.freq <- unlist(lavsamplestats@group.w) * lavsamplestats@ntotal
dx.GW <- - (obs.freq - est.freq)
# divide by N (to be consistent with the rest of lavaan)
dx.GW <- dx.GW / lavsamplestats@ntotal
# remove last element (fixed LAST group to zero)
# dx.GW <- dx.GW[-length(dx.GW)]
# fill in in dx
gw.mat.idx <- which(names(lavmodel@GLIST) == "gw")
gw.x.idx <- unlist( lavmodel@x.free.idx[gw.mat.idx] )
dx[gw.x.idx] <- dx.GW
}
# dx is 1xnpar matrix of LIST (type != "free")
if(is.matrix(dx)) {
dx <- as.numeric(dx)
}
dx
}
# for testing purposes only
# computeDeltaNumerical <- function(lavmodel = NULL, GLIST = NULL, g = 1L) {
#
# # state or final?
# if(is.null(GLIST)) GLIST <- lavmodel@GLIST
#
# compute.moments <- function(x) {
# GLIST <- lav_model_x2GLIST(lavmodel = lavmodel, x=x, type="free")
# Sigma.hat <- computeSigmaHat(lavmodel = lavmodel, GLIST = GLIST)
# S.vec <- lav_matrix_vech(Sigma.hat[[g]])
# if(lavmodel@meanstructure) {
# Mu.hat <- computeMuHat(lavmodel = lavmodel, GLIST=GLIST)
# out <- c(Mu.hat[[g]], S.vec)
# } else {
# out <- S.vec
# }
# out
# }
#
# x <- lav_model_get_parameters(lavmodel = lavmodel, GLIST=GLIST, type="free")
# Delta <- lav_func_jacobian_complex(func=compute.moments, x = x)
#
# Delta
# }
### FIXME: should we here also:
### - weight for groups? (no, for now)
### - handle equality constraints? (yes, for now)
computeDelta <- function(lavmodel = NULL, GLIST. = NULL,
m.el.idx. = NULL, x.el.idx. = NULL,
ceq.simple = FALSE,
force.conditional.x.false = FALSE) {
representation <- lavmodel@representation
categorical <- lavmodel@categorical
if(.hasSlot(lavmodel, "correlation")) {
correlation <- lavmodel@correlation
} else {
correlation <- FALSE
}
conditional.x <- lavmodel@conditional.x
group.w.free <- lavmodel@group.w.free
nmat <- lavmodel@nmat
nblocks <- lavmodel@nblocks
nvar <- lavmodel@nvar
num.idx <- lavmodel@num.idx
th.idx <- lavmodel@th.idx
nexo <- lavmodel@nexo
parameterization <- lavmodel@parameterization
# number of thresholds per group (if any)
nth <- sapply(th.idx, function(x) sum(x > 0L))
# state or final?
if(is.null(GLIST.))
GLIST <- lavmodel@GLIST
else
GLIST <- GLIST.
# type = "free" or something else?
type <- "nonfree"
m.el.idx <- m.el.idx.; x.el.idx <- x.el.idx.
if(is.null(m.el.idx) && is.null(x.el.idx))
type <- "free"
# number of rows in DELTA.group
pstar <- integer(nblocks)
for(g in 1:nblocks) {
pstar[g] <- as.integer(nvar[g] * (nvar[g] + 1) / 2)
if(lavmodel@meanstructure) {
pstar[g] <- nvar[g] + pstar[g] # first the means, then sigma
}
if(categorical) {
pstar[g] <- pstar[g] - nvar[g] # remove variances
pstar[g] <- pstar[g] - nvar[g] # remove means
pstar[g] <- pstar[g] + nth[g] # add thresholds
pstar[g] <- pstar[g] + length(num.idx[[g]]) # add num means
pstar[g] <- pstar[g] + length(num.idx[[g]]) # add num vars
} else if(correlation) {
pstar[g] <- pstar[g] - nvar[g] # remove variances
}
if(conditional.x && nexo[g] > 0L) {
pstar[g] <- pstar[g] + (nvar[g] * nexo[g]) # add slopes
}
if(group.w.free) {
pstar[g] <- pstar[g] + 1L # add group weight
}
}
# number of columns in DELTA + m.el.idx/x.el.idx
if(type == "free") {
if(.hasSlot(lavmodel, "ceq.simple.only") && lavmodel@ceq.simple.only) {
NCOL <- lavmodel@nx.unco
} else {
NCOL <- lavmodel@nx.free
}
m.el.idx <- x.el.idx <- vector("list", length=length(GLIST))
for(mm in 1:length(GLIST)) {
m.el.idx[[mm]] <- lavmodel@m.free.idx[[mm]]
if(.hasSlot(lavmodel, "ceq.simple.only") &&
lavmodel@ceq.simple.only) {
x.el.idx[[mm]] <- lavmodel@x.unco.idx[[mm]]
} else {
x.el.idx[[mm]] <- lavmodel@x.free.idx[[mm]]
}
# handle symmetric matrices
if(lavmodel@isSymmetric[mm]) {
# since we use 'x.free.idx', only symmetric elements
# are duplicated (not the equal ones, only in x.free.free)
dix <- duplicated(x.el.idx[[mm]])
if(any(dix)) {
m.el.idx[[mm]] <- m.el.idx[[mm]][!dix]
x.el.idx[[mm]] <- x.el.idx[[mm]][!dix]
}
}
}
} else {
## FIXME: this does *not* take into account symmetric
## matrices; hence NCOL will be too large, and empty
## columns will be added
## this is ugly, but it doesn't hurt
## alternative could be:
## NCOL <- sum(unlist(lapply(x.el.idx, function(x) length(unique(x)))))
#NCOL <- sum(unlist(lapply(m.el.idx, length)))
NCOL <- sum(unlist(lapply(x.el.idx, function(x) length(unique(x)))))
# sanity check
#nx <- sum(unlist(lapply(x.el.idx, length)))
#stopifnot(NCOL == nx)
}
# compute Delta
Delta <- vector("list", length=nblocks)
for(g in 1:nblocks) {
Delta.group <- matrix(0, nrow=pstar[g], ncol=NCOL)
# which mm belong to group g?
mm.in.group <- 1:nmat[g] + cumsum(c(0,nmat))[g]
# label rows of Delta.group --- FIXME!!!
#if(categorical) {
# # 1. th (means interleaved?)
# # 2. pi
# # 3. var num + cor
#} else {
# if(meanstructure) {
# }
#}
#if(group.w.free) {
#}
# if theta, do some preparation
if(representation == "LISREL" && parameterization == "theta") {
sigma.hat <- computeSigmaHat.LISREL(MLIST=GLIST[mm.in.group],
delta=FALSE)
dsigma <- diag(sigma.hat)
# dcor/dcov for sigma
R <- lav_deriv_cov2cor(sigma.hat, num.idx = lavmodel@num.idx[[g]])
theta.var.idx <- lav_matrix_diagh_idx(nvar[g])
}
for(mm in mm.in.group) {
mname <- names(lavmodel@GLIST)[mm]
# skip empty ones
if(!length(m.el.idx[[mm]])) next
# get Delta columns for this model matrix
if(representation == "LISREL") {
# Sigma
DELTA <- dxSigma <-
derivative.sigma.LISREL(m = mname,
idx = m.el.idx[[mm]],
MLIST = GLIST[ mm.in.group ],
delta = parameterization == "delta")
if(categorical && parameterization == "theta") {
DELTA <- R %*% DELTA
}
if(categorical) {
# reorder: first variances (of numeric), then covariances
cov.idx <- lav_matrix_vech_idx(nvar[g])
covd.idx <- lav_matrix_vech_idx(nvar[g], diagonal = FALSE)
var.idx <- which(is.na(match(cov.idx,
covd.idx)))[num.idx[[g]]]
cor.idx <- match(covd.idx, cov.idx)
DELTA <- rbind(DELTA[var.idx,,drop=FALSE],
DELTA[cor.idx,,drop=FALSE])
}
# correlation structure?
if(!categorical && correlation) {
rm.idx <- lav_matrix_diagh_idx(nvar[g])
DELTA <- DELTA[-rm.idx, , drop = FALSE]
}
if(!categorical) {
if(conditional.x) {
# means/intercepts
DELTA.mu <- derivative.mu.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
# slopes
if(lavmodel@nexo[g] > 0L) {
DELTA.pi <- derivative.pi.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
if(lavmodel@multilevel) {
DELTA <- rbind(DELTA.mu, DELTA.pi, DELTA)
} else {
# ATTENTION: we need to change the order here
# lav_mvreg_scores_* uses 'Beta' where the
# the intercepts are just the first row
# using the col-major approach, we need to
# interweave the intercepts with the slopes!
nEls <- NROW(DELTA.mu) + NROW(DELTA.pi)
# = (nexo + 1 int) * nvar
# intercepts on top
tmp <- rbind(DELTA.mu, DELTA.pi)
# change row index
row.idx <- lav_matrix_vec(matrix(seq.int(nEls),
nrow = lavmodel@nexo[g] + 1L,
ncol = lavmodel@nvar[g], byrow = TRUE))
DELTA.beta <- tmp[row.idx,,drop = FALSE]
DELTA <- rbind(DELTA.beta, DELTA)
}
} else {
DELTA <- rbind(DELTA.mu, DELTA)
}
} else if(!conditional.x && lavmodel@meanstructure) {
DELTA.mu <- derivative.mu.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
DELTA <- rbind(DELTA.mu, DELTA)
}
}
else if(categorical) {
DELTA.th <- derivative.th.LISREL(m=mname,
idx=m.el.idx[[mm]],
th.idx=th.idx[[g]],
MLIST=GLIST[ mm.in.group ],
delta = TRUE)
if(parameterization == "theta") {
# dy/ddsigma = -0.5/(ddsigma*sqrt(ddsigma))
dDelta.dx <-
( dxSigma[theta.var.idx,,drop=FALSE] *
-0.5 / (dsigma*sqrt(dsigma)) )
dth.dDelta <-
derivative.th.LISREL(m = "delta",
idx = 1:nvar[g],
MLIST = GLIST[ mm.in.group ],
th.idx = th.idx[[g]])
# add dth.dDelta %*% dDelta.dx
no.num.idx <- which(th.idx[[g]] > 0)
DELTA.th[no.num.idx,] <-
DELTA.th[no.num.idx,,drop=FALSE] +
(dth.dDelta %*% dDelta.dx)[no.num.idx,,drop=FALSE]
}
if(conditional.x && lavmodel@nexo[g] > 0L) {
DELTA.pi <-
derivative.pi.LISREL(m=mname,
idx=m.el.idx[[mm]],
MLIST=GLIST[ mm.in.group ])
if(parameterization == "theta") {
dpi.dDelta <-
derivative.pi.LISREL(m = "delta",
idx = 1:nvar[g],
MLIST = GLIST[ mm.in.group ])
# add dpi.dDelta %*% dDelta.dx
no.num.idx <-
which(!seq.int(1L,nvar[g]) %in% num.idx[[g]])
no.num.idx <- rep(seq.int(0,nexo[g]-1) * nvar[g],
each=length(no.num.idx)) + no.num.idx
DELTA.pi[no.num.idx,] <-
DELTA.pi[no.num.idx,,drop=FALSE] +
(dpi.dDelta %*% dDelta.dx)[no.num.idx,,drop=FALSE]
}
DELTA <- rbind(DELTA.th, DELTA.pi, DELTA)
} else {
DELTA <- rbind(DELTA.th, DELTA)
}
}
if(group.w.free) {
DELTA.gw <- derivative.gw.LISREL(m=mname,
idx=m.el.idx[[mm]],
MLIST=GLIST[ mm.in.group ])
DELTA <- rbind(DELTA.gw, DELTA)
}
} else if(representation == "RAM") {
DELTA <- dxSigma <-
lav_ram_dsigma(m = mname,
idx = m.el.idx[[mm]],
MLIST = GLIST[ mm.in.group ])
if(lavmodel@meanstructure) {
DELTA.mu <- lav_ram_dmu(m = mname,
idx = m.el.idx[[mm]],
MLIST = GLIST[ mm.in.group ])
DELTA <- rbind(DELTA.mu, DELTA)
}
} else {
stop("representation ", representation, " not implemented yet")
}
Delta.group[ ,x.el.idx[[mm]]] <- DELTA
} # mm
# if type == "free" take care of equality constraints
if(type == "free" && ceq.simple &&
.hasSlot(lavmodel, "ceq.simple.only") && lavmodel@ceq.simple.only) {
Delta.group <- Delta.group %*% lavmodel@ceq.simple.K
}
Delta[[g]] <- Delta.group
} # g
# if multilevel, rbind levels within group
if(.hasSlot(lavmodel, "multilevel") && lavmodel@multilevel) {
DELTA <- vector("list", length = lavmodel@ngroups)
for(g in 1:lavmodel@ngroups) {
DELTA[[g]] <- rbind( Delta[[(g-1)*2 + 1]],
Delta[[(g-1)*2 + 2]] )
}
Delta <- DELTA
}
Delta
}
computeDeltaDx <- function(lavmodel = NULL, GLIST = NULL, target = "lambda",
ceq.simple = FALSE) {
# state or final?
if(is.null(GLIST)) GLIST <- lavmodel@GLIST
representation <- lavmodel@representation
nmat <- lavmodel@nmat
nblocks <- lavmodel@nblocks
num.idx <- lavmodel@num.idx
th.idx <- lavmodel@th.idx
# number of columns in DELTA + m.el.idx/x.el.idx
type <- "free"
#if(type == "free") {
if(.hasSlot(lavmodel, "ceq.simple.only") && lavmodel@ceq.simple.only) {
NCOL <- lavmodel@nx.unco
} else {
NCOL <- lavmodel@nx.free
}
m.el.idx <- x.el.idx <- vector("list", length=length(GLIST))
for(mm in 1:length(GLIST)) {
m.el.idx[[mm]] <- lavmodel@m.free.idx[[mm]]
if(.hasSlot(lavmodel, "ceq.simple.only") &&
lavmodel@ceq.simple.only) {
x.el.idx[[mm]] <- lavmodel@x.unco.idx[[mm]]
} else {
x.el.idx[[mm]] <- lavmodel@x.free.idx[[mm]]
}
# handle symmetric matrices
if(lavmodel@isSymmetric[mm]) {
# since we use 'x.free.idx', only symmetric elements
# are duplicated (not the equal ones, only in x.free.free)
dix <- duplicated(x.el.idx[[mm]])
if(any(dix)) {
m.el.idx[[mm]] <- m.el.idx[[mm]][!dix]
x.el.idx[[mm]] <- x.el.idx[[mm]][!dix]
}
}
}
#} else {
# NCOL <- sum(unlist(lapply(x.el.idx, function(x) length(unique(x)))))
#}
# compute Delta per group
Delta <- vector("list", length=nblocks)
for(g in 1:nblocks) {
mm.in.group <- 1:nmat[g] + cumsum(c(0,nmat))[g]
Delta.group <- NULL
for(mm in mm.in.group) {
mname <- names(lavmodel@GLIST)[mm]
# skip empty ones
if(!length(m.el.idx[[mm]])) next
# get Delta columns for this model matrix
if(representation == "LISREL") {
if(target == "lambda") {
DELTA <- derivative.lambda.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
} else if(target == "th") {
DELTA <- derivative.th.LISREL(m=mname, th.idx = th.idx[[g]],
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ],
delta=TRUE)
} else if(target == "mu") {
DELTA <- derivative.mu.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
} else if(target == "nu") {
DELTA <- derivative.nu.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
} else if(target == "tau") {
DELTA <- derivative.tau.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
} else if(target == "theta") {
DELTA <- derivative.theta.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
} else if(target == "gamma") {
DELTA <- derivative.gamma.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
} else if(target == "beta") {
DELTA <- derivative.beta.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
} else if(target == "alpha") {
DELTA <- derivative.alpha.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
} else if(target == "psi") {
DELTA <- derivative.psi.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ])
} else if(target == "sigma") {
DELTA <- derivative.sigma.LISREL(m=mname,
idx=m.el.idx[[mm]], MLIST=GLIST[ mm.in.group ],
delta=TRUE)
} else {
stop("lavaan ERROR: target ", target, " not implemented yet")
}
# initialize?
if(is.null(Delta.group)) {
Delta.group <- matrix(0, nrow=nrow(DELTA), ncol=NCOL)
}
Delta.group[ ,x.el.idx[[mm]]] <- DELTA
}
} # mm
if(type == "free" && ceq.simple &&
.hasSlot(lavmodel, "ceq.simple.only") && lavmodel@ceq.simple.only) {
Delta.group <- Delta.group %*% lavmodel@ceq.simple.K
}
Delta[[g]] <- Delta.group
} # g
Delta
}
computeOmega <- function(Sigma.hat=NULL, Mu.hat=NULL,
lavsamplestats=NULL, estimator="ML",
meanstructure=FALSE, conditional.x = FALSE) {
# nblocks
nblocks <- length(Sigma.hat)
Omega <- vector("list", length = nblocks)
Omega.mu <- vector("list", length = nblocks)
for(g in 1:nblocks) {
# ML
if(estimator %in% c("ML", "REML", "catML")) {
if(attr(Sigma.hat[[g]], "po") == FALSE) {
# FIXME: WHAT IS THE BEST THING TO DO HERE??
# CURRENTLY: stop
warning("lav_model_gradient: Sigma.hat is not positive definite\n")
Sigma.hat.inv <- MASS::ginv(Sigma.hat[[g]])
Sigma.hat.log.det <- log(.Machine$double.eps)
} else {
Sigma.hat.inv <- attr(Sigma.hat[[g]], "inv")
Sigma.hat.log.det <- attr(Sigma.hat[[g]], "log.det")
}
if(!lavsamplestats@missing.flag) { # complete data
if(meanstructure) {
if(conditional.x) {
diff <- lavsamplestats@res.int[[g]] - Mu.hat[[g]]
W.tilde <- lavsamplestats@res.cov[[g]] + tcrossprod(diff)
} else {
diff <- lavsamplestats@mean[[g]] - Mu.hat[[g]]
W.tilde <- lavsamplestats@cov[[g]] + tcrossprod(diff)
}
# Browne 1995 eq 4.55
Omega.mu[[g]] <- t(t(diff) %*% Sigma.hat.inv)
Omega[[g]] <-
( Sigma.hat.inv %*% (W.tilde - Sigma.hat[[g]]) %*%
Sigma.hat.inv )
} else {
if(conditional.x) {
W.tilde <- lavsamplestats@res.cov[[g]]
} else {
W.tilde <- lavsamplestats@cov[[g]]
}
Omega[[g]] <-
( Sigma.hat.inv %*% (W.tilde - Sigma.hat[[g]]) %*%
Sigma.hat.inv )
}
} else { # missing data
M <- lavsamplestats@missing[[g]]
nvar <- ncol(lavsamplestats@cov[[g]])
OMEGA <- matrix(0, nvar, nvar)
OMEGA.MU <- matrix(0, nvar, 1)
for(p in 1:length(M)) {
SX <- M[[p]][["SY"]]
MX <- M[[p]][["MY"]]
nobs <- M[[p]][["freq"]]
var.idx <- M[[p]][["var.idx"]]
Sigma.inv <- inv.chol(Sigma.hat[[g]][var.idx, var.idx],
logdet=FALSE)
Mu <- Mu.hat[[g]][var.idx]
W.tilde <- SX + tcrossprod(MX - Mu)
OMEGA.MU[var.idx, 1] <-
( OMEGA.MU[var.idx, 1] + nobs/lavsamplestats@ntotal *
t(t(MX - Mu) %*% Sigma.inv) )
OMEGA[var.idx, var.idx] <-
( OMEGA[var.idx, var.idx] + nobs/lavsamplestats@ntotal *
(Sigma.inv %*%
(W.tilde - Sigma.hat[[g]][var.idx,var.idx]) %*%
Sigma.inv ) )
}
Omega.mu[[g]] <- OMEGA.MU
Omega[[g]] <- OMEGA
} # missing
# GLS
} else if(estimator == "GLS") {
W.inv <- lavsamplestats@icov[[g]]
W <- lavsamplestats@cov[[g]]
Omega[[g]] <- (lavsamplestats@nobs[[g]]-1)/lavsamplestats@nobs[[g]] *
(W.inv %*% (W - Sigma.hat[[g]]) %*% W.inv)
if(meanstructure) {
diff <- as.matrix(lavsamplestats@mean[[g]] - Mu.hat[[g]])
Omega.mu[[g]] <- t( t(diff) %*% W.inv )
}
}
} # g
if(meanstructure) attr(Omega, "mu") <- Omega.mu
Omega
}
lav_model_gradient_DD <- function(lavmodel, GLIST = NULL, group = 1L) {
if(is.null(GLIST)) GLIST <- lavmodel@GLIST
#### FIX th + mu!!!!!
Delta.lambda <- computeDeltaDx(lavmodel, GLIST=GLIST, target="lambda")[[group]]
Delta.tau <- computeDeltaDx(lavmodel, GLIST=GLIST, target="tau" )[[group]]
Delta.nu <- computeDeltaDx(lavmodel, GLIST=GLIST, target="nu" )[[group]]
Delta.theta <- computeDeltaDx(lavmodel, GLIST=GLIST, target="theta" )[[group]]
Delta.beta <- computeDeltaDx(lavmodel, GLIST=GLIST, target="beta" )[[group]]
Delta.psi <- computeDeltaDx(lavmodel, GLIST=GLIST, target="psi" )[[group]]
Delta.alpha <- computeDeltaDx(lavmodel, GLIST=GLIST, target="alpha" )[[group]]
Delta.gamma <- computeDeltaDx(lavmodel, GLIST=GLIST, target="gamma" )[[group]]
ov.y.dummy.ov.idx <- lavmodel@ov.y.dummy.ov.idx[[group]]
ov.x.dummy.ov.idx <- lavmodel@ov.x.dummy.ov.idx[[group]]
ov.y.dummy.lv.idx <- lavmodel@ov.y.dummy.lv.idx[[group]]
ov.x.dummy.lv.idx <- lavmodel@ov.x.dummy.lv.idx[[group]]
ov.dummy.idx <- c(ov.y.dummy.ov.idx, ov.x.dummy.ov.idx)
lv.dummy.idx <- c(ov.y.dummy.lv.idx, ov.x.dummy.lv.idx)
th.idx <- lavmodel@th.idx[[group]]
num.idx <- lavmodel@num.idx[[group]]
ord.idx <- unique( th.idx[th.idx > 0L] )
# fix Delta's...
mm.in.group <- 1:lavmodel@nmat[group] + cumsum(c(0,lavmodel@nmat))[group]
MLIST <- GLIST[ mm.in.group ]
DD <- list()
nvar <- lavmodel@nvar
nfac <- ncol(MLIST$lambda) - length(lv.dummy.idx)
# DD$theta
theta.idx <- lav_matrix_diagh_idx(nvar)
DD$theta <- Delta.theta[theta.idx,,drop=FALSE]
if(length(ov.dummy.idx) > 0L) {
psi.idx <- lav_matrix_diagh_idx( ncol(MLIST$psi) )[lv.dummy.idx]
DD$theta[ov.dummy.idx,] <- Delta.psi[psi.idx,,drop=FALSE]
}
# num only? FIXME or just all of them?
DD$theta <- DD$theta[num.idx,,drop=FALSE]
# DD$nu
DD$nu <- Delta.nu
if(length(ov.dummy.idx) > 0L) {
DD$nu[ov.dummy.idx,] <- Delta.alpha[lv.dummy.idx,]
}
DD$nu <- DD$nu[num.idx,,drop=FALSE] # needed?
# DD$lambda
nr <- nvar; nc <- nfac
lambda.idx <- nr*((1:nc) - 1L) + rep(1:nvar, each=nc)
DD$lambda <- Delta.lambda[lambda.idx,,drop=FALSE]
if(length(ov.dummy.idx) > 0L) {
nr <- nrow(MLIST$beta); nc <- nfac # only the first 1:nfac columns
# beta.idx <- rep(nr*((1:nc) - 1L), each=length(lv.dummy.idx)) + rep(lv.dummy.idx, times=nc) ## FIXME
beta.idx <- rep(nr*((1:nc) - 1L), times=length(lv.dummy.idx)) + rep(lv.dummy.idx, each=nc)
#l.idx <- inr*((1:nc) - 1L) + rep(ov.dummy.idx, each=nc) ## FIXME
# l.idx <- rep(nr*((1:nc) - 1L), each=length(ov.dummy.idx)) + rep(ov.dummy.idx, times=nc)
l.idx <- rep(nr*((1:nc) - 1L), times=length(ov.dummy.idx)) + rep(ov.dummy.idx, each=nc)
DD$lambda[match(l.idx, lambda.idx),] <- Delta.beta[beta.idx,,drop=FALSE]
}
# DD$KAPPA
DD$kappa <- Delta.gamma
if(length(ov.dummy.idx) > 0L) {
nr <- nrow(MLIST$gamma); nc <- ncol(MLIST$gamma)
kappa.idx <- nr*((1:nc) - 1L) + rep(lv.dummy.idx, each=nc)
DD$kappa <- DD$kappa[kappa.idx,,drop=FALSE]
}
# DD$GAMMA
if(!is.null(MLIST$gamma)) {
nr <- nrow(MLIST$gamma); nc <- ncol(MLIST$gamma)
lv.idx <- 1:nfac
# MUST BE ROWWISE!
gamma.idx <- rep(nr*((1:nc) - 1L), times=length(lv.idx)) + rep(lv.idx, each=nc)
DD$gamma <- Delta.gamma[gamma.idx,,drop=FALSE]
}
# DD$BETA
if(!is.null(MLIST$beta)) {
nr <- nc <- nrow(MLIST$beta)
lv.idx <- 1:nfac
# MUST BE ROWWISE!
beta.idx <- rep(nr*((1:nfac) - 1L), times=nfac) + rep(lv.idx, each=nfac)
DD$beta <- Delta.beta[beta.idx,,drop=FALSE]
}
## DD$psi
DD$psi <- Delta.psi
if(length(lv.dummy.idx) > 0L) {
nr <- nc <- nrow(MLIST$psi)
lv.idx <- 1:nfac
# MUST BE ROWWISE!
psi.idx <- rep(nr*((1:nfac) - 1L), times=nfac) + rep(lv.idx, each=nfac)
DD$psi <- DD$psi[psi.idx,,drop=FALSE]
}
## DD$tau
if(!is.null(MLIST$tau)) {
DD$tau <- Delta.tau
}
DD
}
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.