Skip to content

Commit b7d7014

Browse files
yrosseelyrosseel
committed
[IV estimation] switch on iv.samplestats = TRUE if only sample statistics are provided
1 parent 4438e6d commit b7d7014

File tree

4 files changed

+161
-40
lines changed

4 files changed

+161
-40
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.2520
3+
Version: 0.6-22.2521
44
Authors@R: c(person(given = "Yves", family = "Rosseel",
55
role = c("aut", "cre"),
66
email = "Yves.Rosseel@UGent.be",

‎R/lav_lavaan_step03_data.R‎

Lines changed: 9 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -87,22 +87,26 @@ lav_lavaan_step03_data <- function(slotData = NULL, # nolint
8787
} else if (lavdata@data.type == "moment") {
8888
# check user-specified options first
8989
if (!is.null(dotdotdot$estimator)) {
90-
if (any(dotdotdot$estimator == c(
90+
if (is.list(dotdotdot$estimator)) {
91+
estimator <- toupper(dotdotdot$estimator[[1]])
92+
} else {
93+
estimator <- toupper(dotdotdot$estimator)
94+
}
95+
if (any(estimator == c(
9196
"MLM", "MLMV", "MLR", "MLR",
9297
"ULSM", "ULSMV", "ULSMVS"
9398
)) &&
9499
is.null(NACOV)) {
95100
lav_msg_stop(gettextf(
96-
"estimator %s requires full data or user-provided NACOV",
97-
dotdotdot$estimator))
98-
} else if (any(dotdotdot$estimator == c(
101+
"estimator %s requires full data or user-provided NACOV", estimator))
102+
} else if (any(estimator == c(
99103
"WLS", "WLSM", "WLSMV",
100104
"WLSMVS", "DWLS"
101105
)) &&
102106
is.null(WLS.V)) {
103107
lav_msg_stop(gettextf(
104108
"estimator %s requires full data or user-provided WLS.V and NACOV",
105-
dotdotdot$estimator))
109+
estimator))
106110
}
107111
}
108112
# catch here some options that will not work with moments

‎R/lav_matrix.R‎

Lines changed: 102 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2026,3 +2026,105 @@ lav_matrix_inverse_iminus <- function(A = NULL) {
20262026
return(IA.inv)
20272027
}
20282028
}
2029+
2030+
# let Gamma_NT = 2 * D.plus %*% (S %x% S) %*% t(D.plus), where
2031+
# D.plus is the elimination matrix: D.plus = ginv(D)
2032+
# we wish to compute R = K %*% Gamma_NT %*% t(K)
2033+
# but without computing the kronecker product explicitly
2034+
# - we allow for meanstructure TRUE/FALSE
2035+
# - we allow for fixed.x = TRUE, if x.idx is not empty
2036+
# - but NOT for conditional.x = TRUE (for now)
2037+
lav_matrix_k_gammant_kt <- function(K, S, meanstructure = FALSE,
2038+
x.idx = integer(0L)) {
2039+
2040+
# dimensions of S
2041+
nvar <- nrow(S)
2042+
pstar <- nvar * (nvar + 1L) / 2L
2043+
2044+
# dimension output
2045+
q <- nrow(K)
2046+
2047+
# sanity checks
2048+
expected_q <- if (meanstructure) nvar + pstar else pstar
2049+
if (ncol(K) != expected_q) {
2050+
lav_msg_stop(gettextf(
2051+
"ncol(K) must be %d when meanstructure = %s",
2052+
expected_q, toupper(meanstructure)
2053+
))
2054+
}
2055+
x.idx <- as.integer(x.idx)
2056+
if (length(x.idx) > 0L && any(x.idx < 1L | x.idx > nvar)) {
2057+
lav_msg_stop(gettext("x.idx contains indices outside 1:nvar."))
2058+
}
2059+
2060+
y.idx <- setdiff(seq_len(nvar), x.idx)
2061+
nvar_y <- length(y.idx)
2062+
nvar_x <- length(x.idx)
2063+
if (nvar_y == 0L) {
2064+
lav_msg_stop(gettextf("No random variables remain after removing x.idx."))
2065+
}
2066+
2067+
# Cholesky of S (no check: we assume S is pd)
2068+
L_S <- t(chol(S))
2069+
2070+
# fixed.x: conditional covariance YbarX and low-rank factor Q for R = Q Q'
2071+
if (nvar_x > 0L) {
2072+
A_yy <- S[y.idx, y.idx, drop = FALSE]
2073+
B_yx <- S[y.idx, x.idx, drop = FALSE]
2074+
C_xx <- S[x.idx, x.idx, drop = FALSE]
2075+
L_x <- t(chol(C_xx))
2076+
2077+
BL <- forwardsolve(L_x, t(B_yx))
2078+
YbarX <- A_yy - t(BL) %*% BL
2079+
2080+
F_full <- matrix(0, nvar, nvar_x)
2081+
F_full[y.idx, ] <- B_yx
2082+
F_full[x.idx, ] <- C_xx
2083+
Q <- t(forwardsolve(L_x, t(F_full)))
2084+
} else {
2085+
YbarX <- S
2086+
Q <- matrix(0, nvar, 0L)
2087+
}
2088+
L_YbarX <- t(chol(YbarX))
2089+
2090+
# build W (d^2 x p) with col j = vec(L' N_j L)
2091+
# B2 is q x pstar (rows are vech vectors); L is nvar x d
2092+
build_W <- function(B2, L) {
2093+
d <- ncol(L)
2094+
if (d == 0L) return(matrix(0, 0L, q))
2095+
2096+
# expand vech row j of B2 into nvar x nvar_q block matrix V_mat
2097+
V_mat <- matrix(0, nvar, nvar * q)
2098+
for (j in seq_len(q)) {
2099+
M <- lav_matrix_vech_reverse(B2[j, ])
2100+
diag(M) <- diag(M) * 2
2101+
V_mat[, (j - 1L) * nvar + seq_len(nvar)] <- M / 2
2102+
}
2103+
# batch compute t(L) N_j L via two (d x nvar_q) multiplies + aperm trick
2104+
LtV <- t(L) %*% V_mat
2105+
LtV3 <- array(LtV, c(d, nvar, q))
2106+
# aperm transposes each block: since M_i symmetric, t(L' M_j) = M_j L
2107+
LtVT <- matrix(aperm(LtV3, c(2L, 1L, 3L)), nvar, d * q)
2108+
matrix(t(L) %*% LtVT, d^2, q)
2109+
}
2110+
2111+
# assemble result
2112+
K2 <- if (meanstructure) {
2113+
K[, nvar + seq_len(pstar), drop = FALSE]
2114+
} else {
2115+
K
2116+
}
2117+
W_S <- build_W(K2, L_S)
2118+
W_R <- build_W(K2, Q)
2119+
out <- 2 * (crossprod(W_S) - crossprod(W_R))
2120+
2121+
if (meanstructure) {
2122+
K1 <- K[, seq_len(nvar), drop = FALSE]
2123+
Z <- K1[, y.idx, drop = FALSE] %*% L_YbarX
2124+
T1 <- tcrossprod(Z)
2125+
out <- out + T1
2126+
}
2127+
2128+
out
2129+
}
2130+

‎R/lav_sem_miiv.R‎

Lines changed: 49 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,10 @@ lav_sem_miiv_internal <- function(lavmodel = NULL, lavh1 = NULL,
1111
stopifnot(iv.method %in% "2SLS")
1212
iv.varcov.method <- toupper(lavoptions$estimator.args$iv.varcov.method)
1313
iv.samplestats <- lavoptions$estimator.args$iv.samplestats
14+
if (lavdata@data.type == "moment") {
15+
# force iv.samplestats = TRUE
16+
iv.samplestats <- TRUE
17+
}
1418
iv.vcov.stage1 <- tolower(lavoptions$estimator.args$iv.vcov.stage1)
1519
iv.sargan <- lavoptions$estimator.args$iv.sargan
1620
# just in case
@@ -814,6 +818,7 @@ lav_sem_miiv_vcov <- function(lavmodel = NULL, lavsamplestats = NULL,
814818

815819
# nblocks
816820
nblocks <- lavmodel@nblocks
821+
stopifnot(lavmodel@nblocks == 1L) # for now...
817822

818823
# eqs?
819824
if (is.null(eqs)) {
@@ -847,38 +852,36 @@ lav_sem_miiv_vcov <- function(lavmodel = NULL, lavsamplestats = NULL,
847852
}
848853

849854
# do we need gamma_big?
850-
if (iv.vcov.stage1 == "gamma" ||
851-
(iv.vcov.stage2 != "none" && length(free.undirected.idx) > 0L)) {
852-
gamma_g <- vector("list", lavmodel@ngroups)
853-
for (g in seq_len(lavmodel@ngroups)) {
854-
if (!is.null(lavsamplestats@NACOV[[g]])) {
855-
gamma_g[[g]] <- lavsamplestats@NACOV[[g]]
856-
} else {
857-
if (iv.vcov.gamma.modelbased) {
858-
mean_g <- lavimplied$mean[[g]]
859-
cov_g <- lavimplied$cov[[g]]
860-
} else {
861-
mean_g <- lavh1$implied$mean[[g]]
862-
cov_g <- lavh1$implied$cov[[g]]
863-
}
864-
# NT version (for now), model-based
865-
gamma_g[[g]] <- lav_samplestats_Gamma_NT(
866-
COV = cov_g,
867-
MEAN = mean_g,
868-
x.idx = lavsamplestats@x.idx[[g]],
869-
fixed.x = lavmodel@fixed.x,
870-
conditional.x = lavmodel@conditional.x,
871-
meanstructure = lavmodel@meanstructure,
872-
slopestructure = lavmodel@conditional.x
873-
)
874-
# weight by (group) sample size
875-
fg <- lavsamplestats@nobs[[g]] / lavsamplestats@ntotal
876-
gamma_g[[g]] <- fg * gamma_g[[g]]
855+
if (iv.vcov.stage1 == "gamma" ||
856+
(iv.vcov.stage2 != "none" && length(free.undirected.idx) > 0L)) {
857+
# gamma_g <- vector("list", lavmodel@ngroups)
858+
for (g in seq_len(lavmodel@ngroups)) {
859+
# if (!is.null(lavsamplestats@NACOV[[g]])) {
860+
# gamma_g[[g]] <- lavsamplestats@NACOV[[g]]
861+
# } else {
862+
if (iv.vcov.gamma.modelbased) {
863+
cov_g <- lavimplied$cov[[g]]
864+
} else {
865+
cov_g <- lavh1$implied$cov[[g]]
866+
}
867+
# # NT version (for now), model-based
868+
# gamma_g[[g]] <- lav_samplestats_Gamma_NT(
869+
# COV = cov_g,
870+
# MEAN = mean_g,
871+
# x.idx = lavsamplestats@x.idx[[g]],
872+
# fixed.x = lavmodel@fixed.x,
873+
# conditional.x = lavmodel@conditional.x,
874+
# meanstructure = lavmodel@meanstructure,
875+
# slopestructure = lavmodel@conditional.x
876+
# )
877+
# # weight by (group) sample size
878+
# fg <- lavsamplestats@nobs[[g]] / lavsamplestats@ntotal
879+
# gamma_g[[g]] <- fg * gamma_g[[g]]
880+
# }
877881
}
882+
# gamma_big <- lav_matrix_bdiag(gamma_g)
878883
}
879-
gamma_big <- lav_matrix_bdiag(gamma_g)
880-
}
881-
884+
882885
# stage 1: directed effects
883886
if (length(free.directed.idx) > 0L && iv.vcov.stage1 != "none") {
884887
if (iv.vcov.stage1 == "gamma") {
@@ -905,8 +908,13 @@ lav_sem_miiv_vcov <- function(lavmodel = NULL, lavsamplestats = NULL,
905908
free.directed.idx = free.directed.idx
906909
)
907910
}
911+
# K %*% Gamma_NT %*% t(K)
912+
#k_gammant_kt <- jac_k %*% gamma_big %*% t(jac_k)
913+
x.idx <- if (lavmodel@fixed.x) lavsamplestats@x.idx[[1]] else integer(0L)
914+
k_gammant_kt <- lav_matrix_k_gammant_kt(K = jac_k, S = cov_g,
915+
meanstructure = lavmodel@meanstructure, x.idx = x.idx)
908916
vcov[free.directed.idx, free.directed.idx] <-
909-
(jac_k %*% gamma_big %*% t(jac_k)) / lavsamplestats@ntotal
917+
k_gammant_kt / lavsamplestats@ntotal
910918
} else {
911919
for (b in seq_len(nblocks)) {
912920
neqs <- length(eqs[[b]])
@@ -992,8 +1000,12 @@ lav_sem_miiv_vcov <- function(lavmodel = NULL, lavsamplestats = NULL,
9921000
} else {
9931001
tmp <- h2
9941002
}
1003+
x.idx <- if (lavmodel@fixed.x) lavsamplestats@x.idx[[1]] else integer(0L)
1004+
tmp_gammant_tmpt <- lav_matrix_k_gammant_kt(K = tmp, S = cov_g,
1005+
meanstructure = lavmodel@meanstructure, x.idx = x.idx)
1006+
#tmp_gammant_tmpt_bis <- tmp %*% gamma_big %*% t(tmp)
9951007
vcov[free.undirected.idx, free.undirected.idx] <-
996-
(tmp %*% gamma_big %*% t(tmp)) / lavsamplestats@ntotal
1008+
tmp_gammant_tmpt / lavsamplestats@ntotal
9971009

9981010
# iv.vcov.stage2 == "delta"
9991011
} else {
@@ -1100,8 +1112,11 @@ lav_sem_miiv_vcov <- function(lavmodel = NULL, lavsamplestats = NULL,
11001112
iv.varcov.method = iv.varcov.method
11011113
)
11021114
}
1103-
1104-
vcov_b <- (jac_b %*% gamma_big %*% t(jac_b)) / lavsamplestats@ntotal
1115+
x.idx <- if (lavmodel@fixed.x) lavsamplestats@x.idx[[1]] else integer(0L)
1116+
jacb_gammant_jacbt <- lav_matrix_k_gammant_kt(K = jac_b, S = cov_g,
1117+
meanstructure = lavmodel@meanstructure, x.idx = x.idx)
1118+
#jacb_gammant_jacbt_bis <- jac_b %*% gamma_big %*% t(jac_b)
1119+
vcov_b <- jacb_gammant_jacbt / lavsamplestats@ntotal
11051120
vcov_ab <- vcov_a + vcov_b
11061121
vcov[free.undirected.idx, free.undirected.idx] <- vcov_ab
11071122
} # continuous

0 commit comments

Comments
 (0)