Skip to content

Commit

Permalink
better handling of missing-by-design settings when using multilevel S…
Browse files Browse the repository at this point in the history
…EM + FIML
  • Loading branch information
yrosseel committed May 24, 2024
1 parent f5ff179 commit 808a0be
Show file tree
Hide file tree
Showing 7 changed files with 1,011 additions and 352 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Package: lavaan
Title: Latent Variable Analysis
Version: 0.6-18.2112
Version: 0.6-18.2113
Authors@R: c(person(given = "Yves", family = "Rosseel",
role = c("aut", "cre"),
email = "Yves.Rosseel@UGent.be",
Expand Down
18 changes: 10 additions & 8 deletions R/lav_data.R
Original file line number Diff line number Diff line change
Expand Up @@ -1153,14 +1153,16 @@ lav_data_full <- function(data = NULL, # data.frame
paste(empty.case.idx, collapse = " ")))
}
}
if (warn && any(Mp[[g]]$coverage == 0)) {
lav_msg_warn(gettext(
"due to missing values, some pairwise combinations have 0% coverage;
use lavInspect(fit, \"coverage\") to investigate."))
} else if (warn && any(Mp[[g]]$coverage < 0.1)) {
lav_msg_warn(gettext(
"due to missing values, some pairwise combinations have less than
10% coverage; use lavInspect(fit, \"coverage\") to investigate."))
if (warn && any(Mp[[g]]$coverage < 0.1)) {
coverage.vech <- lav_matrix_vech(Mp[[g]]$coverage, diagonal = FALSE)
small.idx <- which(coverage.vech < 0.1)
if (all(coverage.vech[small.idx] == 0)) {
# 0.6-18: no warning --> this could be due to missing by design
} else {
lav_msg_warn(gettext(
"due to missing values, some pairwise combinations have less than
10% coverage; use lavInspect(fit, \"coverage\") to investigate."))
}
}
# in case we had observations with only missings
nobs[[g]] <- NROW(X[[g]]) - length(Mp[[g]]$empty.idx)
Expand Down
75 changes: 54 additions & 21 deletions R/lav_model_h1_information.R
Original file line number Diff line number Diff line change
Expand Up @@ -572,16 +572,29 @@ lav_model_h1_information_observed <- function(lavobject = NULL,
SIGMA.W <- implied$cov[[(g - 1) * lavdata@nlevels + 1L]]
SIGMA.B <- implied$cov[[(g - 1) * lavdata@nlevels + 2L]]

# clustered data
A1[[g]] <- lav_mvnorm_cluster_information_observed(
Lp = lavdata@Lp[[g]],
YLp = lavsamplestats@YLp[[g]],
Mu.W = MU.W,
Sigma.W = SIGMA.W,
Mu.B = MU.B,
Sigma.B = SIGMA.B,
x.idx = lavsamplestats@x.idx[[g]]
)
if (lavdata@missing == "ml") {
A1[[g]] <- lav_mvnorm_cluster_missing_information_observed(
Y1 = lavdata@X[[g]],
Y2 = lavsamplestats@YLp[[g]][[2]]$Y2,
Lp = lavdata@Lp[[g]],
Mp = lavdata@Mp[[g]],
Mu.W = MU.W,
Sigma.W = SIGMA.W,
Mu.B = MU.B,
Sigma.B = SIGMA.B,
x.idx = lavsamplestats@x.idx[[g]]
)
} else {
A1[[g]] <- lav_mvnorm_cluster_information_observed(
Lp = lavdata@Lp[[g]],
YLp = lavsamplestats@YLp[[g]],
Mu.W = MU.W,
Sigma.W = SIGMA.W,
Mu.B = MU.B,
Sigma.B = SIGMA.B,
x.idx = lavsamplestats@x.idx[[g]]
)
}
} # g
} # ML + multilevel

Expand Down Expand Up @@ -761,6 +774,10 @@ lav_model_h1_information_firstorder <- function(lavobject = NULL,
# if not-structured, we use lavh1, and that is always
# 'unconditional' (for now)
if (lavmodel@conditional.x && structured) {
if (lavdata@missing == "ml") {
lav_msg_stop(gettext("firstorder information matrix not available ",
"(yet) if conditional.x + fiml"))
}
Res.Sigma.W <- implied$res.cov[[(g - 1) * lavdata@nlevels + 1L]]
Res.Int.W <- implied$res.int[[(g - 1) * lavdata@nlevels + 1L]]
Res.Pi.W <- implied$res.slopes[[(g - 1) * lavdata@nlevels + 1L]]
Expand All @@ -785,17 +802,33 @@ lav_model_h1_information_firstorder <- function(lavobject = NULL,
MU.B <- implied$mean[[(g - 1) * lavdata@nlevels + 2L]]
SIGMA.W <- implied$cov[[(g - 1) * lavdata@nlevels + 1L]]
SIGMA.B <- implied$cov[[(g - 1) * lavdata@nlevels + 2L]]
B1[[g]] <- lav_mvnorm_cluster_information_firstorder(
Y1 = lavdata@X[[g]],
YLp = lavsamplestats@YLp[[g]],
Lp = lavdata@Lp[[g]],
Mu.W = MU.W,
Sigma.W = SIGMA.W,
Mu.B = MU.B,
Sigma.B = SIGMA.B,
x.idx = lavsamplestats@x.idx[[g]],
divide.by.two = TRUE
)
if (lavdata@missing == "ml") {
B1[[g]] <- lav_mvnorm_cluster_missing_information_firstorder(
Y1 = lavdata@X[[g]],
Y2 = lavsamplestats@YLp[[g]][[2]]$Y2,
Lp = lavdata@Lp[[g]],
Mp = lavdata@Mp[[g]],
Mu.W = MU.W,
Sigma.W = SIGMA.W,
Mu.B = MU.B,
Sigma.B = SIGMA.B,
x.idx = lavsamplestats@x.idx[[g]],
divide.by.two = TRUE
)
} else {
# no missing values
B1[[g]] <- lav_mvnorm_cluster_information_firstorder(
Y1 = lavdata@X[[g]],
YLp = lavsamplestats@YLp[[g]],
Lp = lavdata@Lp[[g]],
Mu.W = MU.W,
Sigma.W = SIGMA.W,
Mu.B = MU.B,
Sigma.B = SIGMA.B,
x.idx = lavsamplestats@x.idx[[g]],
divide.by.two = TRUE
)
}
}
} else if (estimator == "ML" && lavdata@nlevels == 1L) {
if (length(lavdata@cluster) > 0L) {
Expand Down
5 changes: 4 additions & 1 deletion R/lav_mvnorm_cluster.R
Original file line number Diff line number Diff line change
Expand Up @@ -56,13 +56,15 @@ lav_mvnorm_cluster_implied22l <- function(Lp = NULL,
Mu.B.tilde <- numeric(p.tilde)
Mu.B.tilde[ov.idx[[2]]] <- Mu.B


# add Mu.W[within.idx] to Mu.B
Mu.WB.tilde <- numeric(p.tilde)
Mu.WB.tilde[within.idx] <- Mu.W.tilde[within.idx]
Mu.WB.tilde[both.idx] <- (Mu.B.tilde[both.idx] +
Mu.W.tilde[both.idx])

# set Mu.W[both.idx] to zero (after we added tot WB)
Mu.W.tilde[both.idx] <- 0

# map to matrices needed for loglik
if (length(within.idx) > 0L) {
Mu.B.tilde[within.idx] <- 0
Expand Down Expand Up @@ -106,6 +108,7 @@ lav_mvnorm_cluster_2l2implied <- function(Lp,
# between.idx
between.idx <- Lp$between.idx[[2]]
within.idx <- Lp$within.idx[[2]]
both.idx <- Lp$both.idx[[2]]

# ov.idx per level
ov.idx <- Lp$ov.idx
Expand Down
Loading

0 comments on commit 808a0be

Please sign in to comment.