Skip to content

Commit 9d8398e

Browse files
yrosseelyrosseel
committed
[IV estimation] finishing vcov computating in both stages (continuous data only)
1 parent 987d2c1 commit 9d8398e

File tree

7 files changed

+444
-150
lines changed

7 files changed

+444
-150
lines changed

‎DESCRIPTION‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
Package: lavaan
22
Title: Latent Variable Analysis
3-
Version: 0.6-22.2516
3+
Version: 0.6-22.2517
44
Authors@R: c(person(given = "Yves", family = "Rosseel",
55
role = c("aut", "cre"),
66
email = "Yves.Rosseel@UGent.be",

‎R/lav_matrix.R‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1505,6 +1505,24 @@ lav_matrix_symmetric_inverse <- function(S, logdet = FALSE,
15051505
S.inv
15061506
}
15071507

1508+
# solve(A) %*% B, where A is symmetric, and B is vector of matrix
1509+
lav_matrix_symmetric_solve_spd <- function(A, B, tol = 1e-10) {
1510+
cholA <- try(chol(A), silent = TRUE)
1511+
if (!inherits(cholA, "try-error")) {
1512+
# Cholesky solve
1513+
return(backsolve(cholA, forwardsolve(t(cholA), B)))
1514+
} else {
1515+
# Eigen fallback (pseudo-inverse)
1516+
eig <- eigen(A, symmetric = TRUE)
1517+
keep <- eig$values > tol * max(eig$values)
1518+
Ainv <- eig$vectors[, keep, drop = FALSE] %*%
1519+
diag(1 / eig$values[keep], length(eig$values[keep])) %*%
1520+
t(eig$vectors[, keep, drop = FALSE])
1521+
return(Ainv %*% B)
1522+
}
1523+
}
1524+
1525+
15081526
# update inverse of A, after removing 1 or more rows (and corresponding
15091527
# colums) from A
15101528
#

‎R/lav_object_post_check.R‎

Lines changed: 12 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -6,12 +6,16 @@ lav_object_post_check <- function(object) {
66
lavdata <- object@Data
77

88
var.ov.ok <- var.lv.ok <- result.ok <- TRUE
9+
var.na <- FALSE
910

1011
# 1a. check for negative variances ov
1112
var.idx <- which(lavpartable$op == "~~" &
1213
lavpartable$lhs %in% lav_object_vnames(object, "ov") &
1314
lavpartable$lhs == lavpartable$rhs)
14-
if (length(var.idx) > 0L && any(lavpartable$est[var.idx] < 0.0)) {
15+
if (any(is.na(lavpartable$est[var.idx]))) {
16+
# perhaps estimator = "IV" + stage 1 only
17+
var.na <- TRUE
18+
} else if (length(var.idx) > 0L && any(lavpartable$est[var.idx] < 0.0)) {
1519
result.ok <- var.ov.ok <- FALSE
1620
lav_msg_warn(gettext("some estimated ov variances are negative"))
1721
}
@@ -20,14 +24,18 @@ lav_object_post_check <- function(object) {
2024
var.idx <- which(lavpartable$op == "~~" &
2125
lavpartable$lhs %in% lav_object_vnames(object, "lv") &
2226
lavpartable$lhs == lavpartable$rhs)
23-
if (length(var.idx) > 0L && any(lavpartable$est[var.idx] < 0.0)) {
27+
if (any(is.na(lavpartable$est[var.idx]))) {
28+
# perhaps estimator = "IV" + stage 1 only
29+
var.na <- TRUE
30+
} else if (length(var.idx) > 0L && any(lavpartable$est[var.idx] < 0.0)) {
2431
result.ok <- var.lv.ok <- FALSE
2532
lav_msg_warn(gettext("some estimated lv variances are negative"))
2633
}
2734

2835
# 2. is cov.lv (PSI) positive definite? (only if we did not already warn
2936
# for negative variances)
30-
if (var.lv.ok && length(lav_object_vnames(lavpartable, type = "lv.regular")) > 0L) {
37+
if (!var.na && var.lv.ok &&
38+
length(lav_object_vnames(lavpartable, type = "lv.regular")) > 0L) {
3139
ETA <- lavTech(object, "cov.lv")
3240
for (g in 1:lavdata@ngroups) {
3341
if (nrow(ETA[[g]]) == 0L) next
@@ -45,7 +53,7 @@ lav_object_post_check <- function(object) {
4553

4654
# 3. is THETA positive definite (but only for numeric variables)
4755
# and if we not already warned for negative ov variances
48-
if (var.ov.ok) {
56+
if (!var.na && var.ov.ok) {
4957
THETA <- lavTech(object, "theta")
5058
for (g in 1:lavdata@ngroups) {
5159
num.idx <- lavmodel@num.idx[[g]]

‎R/lav_options_estimator.R‎

Lines changed: 71 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -601,16 +601,25 @@ lav_options_est_iv <- function(opt) {
601601
# sample.icov not needed
602602
opt$sample.icov <- FALSE
603603

604+
# no loglik
605+
opt$loglik <- FALSE
606+
607+
# no rescale of vcov
608+
opt$sample.cov.rescale <- FALSE
609+
604610
# estimator options
605611
if (is.null(opt$estimator.args)) {
606612
# create default list
607613
opt$estimator.args <- list(iv.method = "2SLS",
608614
iv.samplestats = FALSE,
609615
iv.varcov.method = "RLS",
610-
iv.varcov.se = TRUE,
611-
iv.varcov.modelbased = TRUE,
612-
iv.varcov.jaca.numerical = FALSE,
613-
iv.varcov.jacb.numerical = FALSE)
616+
iv.sargan = TRUE,
617+
iv.vcov.stage1 <- "lm.vcov.dfres",
618+
iv.vcov.stage2 <- "delta",
619+
iv.vcov.gamma.modelbased = TRUE,
620+
iv.vcov.jack.numerical = FALSE,
621+
iv.vcov.jaca.numerical = FALSE,
622+
iv.vcov.jacb.numerical = FALSE)
614623
} else {
615624
if (is.null(opt$estimator.args$iv.method)) {
616625
opt$estimator.args$iv.method <- "2SLS"
@@ -620,30 +629,79 @@ lav_options_est_iv <- function(opt) {
620629
if (is.null(opt$estimator.args$iv.samplestats)) {
621630
opt$estimator.args$iv.samplestats <- FALSE
622631
}
632+
if (is.null(opt$estimator.args$iv.sargan)) {
633+
opt$estimator.args$iv.sargan <- TRUE
634+
}
635+
if (is.null(opt$estimator.args$iv.vcov.stage1)) {
636+
if (opt$.categorical) {
637+
opt$estimator.args$iv.vcov.stage1 <- "lm.vcov.dfres"
638+
} else {
639+
opt$estimator.args$iv.vcov.stage1 <- "lm.vcov.dfres"
640+
}
641+
if (tolower(opt$se[1]) == "none") {
642+
opt$estimator.args$iv.vcov.stage1 <- "none"
643+
}
644+
} else if (!tolower(opt$estimator.args$iv.vcov.stage1) %in%
645+
c("lm.vcov", "lm.vcov.dfres", "gamma", "none")) {
646+
lav_msg_stop(gettext("iv.vcov.stage1 should be lm.vcov, lm.vcov.dfres,
647+
gamma or none."))
648+
} else if (tolower(opt$estimator.args$iv.vcov.stage1) == "gamma" &&
649+
!opt$estimator.args$iv.samplestats) {
650+
lav_msg_stop(gettext("iv.vcov.stage1 cannot be gamma if
651+
iv.samplestats is FALSE"))
652+
}
653+
if (opt$.categorical) {
654+
opt$estimator.args$iv.samplestats <- TRUE
655+
}
656+
if (is.null(opt$estimator.args$iv.vcov.stage2)) {
657+
if (opt$.categorical ||
658+
tolower(opt$estimator.args$iv.vcov.stage1) == "gamma") {
659+
opt$estimator.args$iv.vcov.stage2 <- "h2"
660+
} else {
661+
opt$estimator.args$iv.vcov.stage2 <- "delta"
662+
}
663+
if (tolower(opt$se[1]) == "none") {
664+
opt$estimator.args$iv.vcov.stage2 <- "none"
665+
}
666+
} else if (!tolower(opt$estimator.args$iv.vcov.stage2) %in%
667+
c("h2", "delta", "none")) {
668+
lav_msg_stop(gettext("iv.vcov.stage2 should be h2, delta, or none."))
669+
} else if (tolower(opt$estimator.args$iv.vcov.stage1) %in%
670+
c("lm.vcov", "lm.vcov.dfres") &&
671+
tolower(opt$estimator.args$iv.vcov.stage2) == "h2") {
672+
lav_msg_stop(gettext("iv.vcov.stage2 should be delta (or none) if
673+
iv.vcov.stage1 is either lm.vcov or lm.vcov.dfres."))
674+
}
623675
if (opt$.categorical) {
624676
opt$estimator.args$iv.samplestats <- TRUE
625677
}
626678
if (is.null(opt$estimator.args$iv.varcov.method)) {
627679
opt$estimator.args$iv.varcov.method <- "RLS"
628680
} else if (!toupper(opt$estimator.args$iv.varcov.method) %in%
629681
c("ULS", "GLS", "2RLS", "RLS", "NONE")) {
630-
lav_msg_stop(gettext("iv.varcov.method should ULS, GLS, 2RLS, RLS
682+
lav_msg_stop(gettext("iv.varcov.method should be ULS, GLS, 2RLS, RLS
631683
or NONE."))
632684
}
633-
if (is.null(opt$estimator.args$iv.varcov.se)) {
634-
opt$estimator.args$iv.varcov.se <- TRUE
685+
if (is.null(opt$estimator.args$iv.vcov.gamma.modelbased)) {
686+
opt$estimator.args$iv.vcov.gamma.modelbased <- TRUE
635687
}
636-
if (is.null(opt$estimator.args$iv.varcov.modelbased)) {
637-
opt$estimator.args$iv.varcov.modelbased <- TRUE
688+
if (is.null(opt$estimator.args$iv.vcov.jack.numerical)) {
689+
opt$estimator.args$iv.vcov.jack.numerical <- FALSE
638690
}
639-
if (is.null(opt$estimator.args$iv.varcov.jaca.numerical)) {
640-
opt$estimator.args$iv.varcov.jaca.numerical <- FALSE
691+
if (is.null(opt$estimator.args$iv.vcov.jaca.numerical)) {
692+
opt$estimator.args$iv.vcov.jaca.numerical <- FALSE
641693
}
642-
if (is.null(opt$estimator.args$iv.varcov.jacb.numerical)) {
643-
opt$estimator.args$iv.varcov.jacb.numerical <- FALSE
694+
if (is.null(opt$estimator.args$iv.vcov.jacb.numerical)) {
695+
opt$estimator.args$iv.vcov.jacb.numerical <- FALSE
644696
}
645697
}
646698

699+
# override test if iv.varcov.method is "none"
700+
if (toupper(opt$estimator.args$iv.varcov.method) == "NONE") {
701+
opt$test <- "none"
702+
opt$standard.test <- "none"
703+
}
704+
647705
opt
648706
}
649707

‎R/lav_samplestats.R‎

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -827,7 +827,7 @@ lav_samplestats_from_data <- function(lavdata = NULL,
827827
Mplus.WLS = FALSE
828828
)
829829
}
830-
} else if (estimator %in% c("WLS", "DWLS", "ULS", "DLS", "catML", "IV")) {
830+
} else if (estimator %in% c("WLS", "DWLS", "ULS", "DLS", "catML")) {
831831
if (!categorical) {
832832
# sample size large enough?
833833
nvar <- ncol(X[[g]])
@@ -867,10 +867,13 @@ lav_samplestats_from_data <- function(lavdata = NULL,
867867
meanstructure = meanstructure
868868
)
869869
} else {
870-
if (lavoptions$se == "robust.sem.nt" || estimator == "IV") {
870+
if ("robust.sem.nt" %in% lavoptions$se ||
871+
"browne.residual.nt" %in% lavoptions$test) {
871872
NACOV[[g]] <-
872873
lav_samplestats_Gamma_NT(
873874
Y = Y,
875+
COV = cov[[g]],
876+
MEAN = mean[[g]],
874877
x.idx = x.idx[[g]],
875878
# cluster.idx = cluster.idx, # not available
876879
fixed.x = fixed.x,

0 commit comments

Comments
 (0)