New Upstream Release - r-cran-lavaan
Ready changes
Summary
Merged new upstream version: 0.6.15 (was: 0.6.14).
Diff
diff --git a/DESCRIPTION b/DESCRIPTION
index 540db9f..e57443b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: lavaan
Title: Latent Variable Analysis
-Version: 0.6-14
+Version: 0.6-15
Authors@R: c(person(given = "Yves", family = "Rosseel",
role = c("aut", "cre"),
email = "Yves.Rosseel@UGent.be",
@@ -57,7 +57,7 @@ LazyData: yes
ByteCompile: true
URL: https://lavaan.ugent.be
NeedsCompilation: no
-Packaged: 2023-02-09 12:46:29 UTC; yves
+Packaged: 2023-03-14 15:12:27 UTC; yves
Author: Yves Rosseel [aut, cre] (<https://orcid.org/0000-0002-4129-4477>),
Terrence D. Jorgensen [aut] (<https://orcid.org/0000-0001-5111-6773>),
Nicholas Rockwood [aut] (<https://orcid.org/0000-0001-5931-183X>),
@@ -74,4 +74,4 @@ Author: Yves Rosseel [aut, cre] (<https://orcid.org/0000-0002-4129-4477>),
Han Du [ctb]
Maintainer: Yves Rosseel <Yves.Rosseel@UGent.be>
Repository: CRAN
-Date/Publication: 2023-02-09 17:50:02 UTC
+Date/Publication: 2023-03-14 18:50:02 UTC
diff --git a/MD5 b/MD5
index 278e22b..f890966 100644
--- a/MD5
+++ b/MD5
@@ -1,5 +1,5 @@
-21712371c94d55d7f3769f66d1f487e9 *DESCRIPTION
-7bdc3539cab05c7d06b1ffdcf8e94145 *NAMESPACE
+ceccc6c1de4f09cd871a5dfe003eefe5 *DESCRIPTION
+0e50b3573d5997effa05c508cbd7f923 *NAMESPACE
8c07055a2f6d2c4fed52c259cdeef334 *R/00class.R
d3ff28a992eb9ad9f3302dab2d1db831 *R/00generic.R
38752e01bebf37f3e6a30e121565a942 *R/ctr_estfun.R
@@ -13,7 +13,7 @@ f8449cb947eaeb0f474f7e6a71cf90db *R/ctr_pml_plrt.R
7037519cf3c4973a139602463cda59f6 *R/ctr_pml_plrt2.R
1f6abdeeee0f828057b0550be0a27575 *R/ctr_pml_plrt_nested.R
6635c66e092e2bb8609674f3a1d43924 *R/ctr_pml_utils.R
-14e65418e7e7176e4cf707946632b2fd *R/lav_bootstrap.R
+1d53c8c7f7f87a09d2d10b8ae8167a2b *R/lav_bootstrap.R
802062ef84ed7cb27563939c9cfc4caf *R/lav_bootstrap_lrt.R
786955839de8160c7e81d23c636ee109 *R/lav_bvmix.R
4db8c88bc456d03f7e5eb93ca13f2bad *R/lav_bvord.R
@@ -22,9 +22,10 @@ de4d729b1c11b91d149d37468efd7c92 *R/lav_cfa_1fac.R
6ab321fbdc2199123a68d06d32507ea1 *R/lav_cfa_bentler1982.R
0b4b0e7c54eb33ea29bb597b29a6efb7 *R/lav_cfa_fabin.R
8c8196806e4831bb4973d7613130a6ea *R/lav_cfa_guttman1952.R
-9d094e4f1360b839869c079fb0e2b5d8 *R/lav_cfa_utils.R
+459371fe3ee8767385dc3b50f485abc8 *R/lav_cfa_jamesstein.R
+89362108a9aa09bd8a9bb31e8d7b01f8 *R/lav_cfa_utils.R
7e34cae25007369f6eefb62e15fc2956 *R/lav_constraints.R
-5f1943d5b50b829ff175d2e4060afa3e *R/lav_cor.R
+c54ee15382ab153c7e84ed15ec78fe01 *R/lav_cor.R
b6e7214ae9c97cac7d0aa998682ea112 *R/lav_data.R
12fc53360224efb5cb05d7d1362913b8 *R/lav_data_patterns.R
510d1f1e5455693f07d76f501885497c *R/lav_data_print.R
@@ -42,7 +43,7 @@ e6633a2398d5837013b0841cceb4a714 *R/lav_efa_print.R
5b1408371dd3b6303e579d6bee6d6ed5 *R/lav_fit_aic.R
a7c4280094e9cba742ed82edccc78db4 *R/lav_fit_cfi.R
debace4e95c32d5c76190c5541b4040b *R/lav_fit_gfi.R
-84e8ef21698150669b2e893a16e2fdb3 *R/lav_fit_measures.R
+e57c23e463d14413998b62a815d2763f *R/lav_fit_measures.R
373616ebbeb1182310e96fc77a84d748 *R/lav_fit_other.R
239ab22eccc8c6d08d387378020cfe30 *R/lav_fit_rmsea.R
0cbb5b20e43624b28d68a2ef3046c72f *R/lav_fit_srmr.R
@@ -58,7 +59,7 @@ da1f299fdd19b249f16557bffe447c6c *R/lav_lavaanList_methods.R
04e57c141fc9ecb958491d8a187dfc98 *R/lav_lavaanList_multipleGroups.R
a4603d9dd7d900db3e0ac06aa3912d57 *R/lav_lavaanList_multipleImputation.R
b0a9f7eff4f795b955e1bc8d7e44387b *R/lav_lavaanList_simulate.R
-3d3397a2e2084d4a76c9f0c8fb94c0e3 *R/lav_matrix.R
+86d7db55ecd08909714f9e1915bd33ae *R/lav_matrix.R
69d1d8697b94d2882c1852e5d4323e0d *R/lav_matrix_rotate.R
a8db43f516ad7cee63261fd4a640f321 *R/lav_matrix_rotate_methods.R
341752928a4366761365b66e681ead51 *R/lav_matrix_rotate_utils.R
@@ -101,8 +102,8 @@ bdb77c93e5ef0c514220940b9b858c98 *R/lav_object_inspect.R
0eb2e28617f98430eb94cdf9bc70bae9 *R/lav_objective.R
7458ef10c42164365d6cdd606ea21968 *R/lav_optim_gn.R
c1eea693e94eaa1d21a33db58870c970 *R/lav_optim_nlminb_constr.R
-97dcb35a552d219c584cb048165f60a8 *R/lav_optim_noniter.R
-cdb56a180cba732fa3cd20e5044b62e7 *R/lav_options.R
+5b705348c5131ab60e472465eea36134 *R/lav_optim_noniter.R
+9a0bf7b810dc94b64dd76505facc4b84 *R/lav_options.R
3198011aa92ba0079caa41d183b9d47a *R/lav_partable.R
00b1df081c75af43613d181e7121e83e *R/lav_partable_attributes.R
fea6713eaefb45aeeaafe4025cb11dad *R/lav_partable_bounds.R
@@ -119,11 +120,11 @@ ca3ea92f0ede0741192173e18ab62efc *R/lav_partable_full.R
990d96f3e550f1808807707dd314b96d *R/lav_partable_subset.R
b75fd0632684d788dc731e54de35f8fe *R/lav_partable_unrestricted.R
c8f24afbd6e6c8afef3285e4f2abcf4e *R/lav_partable_utils.R
-d48f3e258822c1f24f233b13aa821455 *R/lav_partable_vnames.R
+73d46bd54614166b41b58865f169c897 *R/lav_partable_vnames.R
e94c70b6ff32491c3c7a3d547de866e6 *R/lav_predict.R
a32a85b7962efe002db3c6f983d65625 *R/lav_predict_y.R
ef4aab312133fcecd6b33c10e6ee03e3 *R/lav_prelis.R
-a3261cb1f08618c1bc7765f73db1a9a9 *R/lav_print.R
+5f20fc8fc5d52ec550a34d78cb180d52 *R/lav_print.R
704a67f366f62b1597c514c2742d6369 *R/lav_representation.R
5ef4a5b05de392cf249738eb931c9d07 *R/lav_representation_lisrel.R
9cb82137123d946b609a213dec35b2ba *R/lav_representation_ram.R
@@ -186,7 +187,7 @@ b88073fe6a281aa96c0fc126e3fecf8c *man/bootstrap.Rd
e10ce2f596f1e9105b12d1e9860aacfd *man/getCov.Rd
c09580c5994692adf594e66741cfce54 *man/growth.Rd
73e825fa761127c81205c806be420ffa *man/inspectSampleCov.Rd
-a9db35f788a7250e1661e72710bc908b *man/lavCor.Rd
+bdeed4c0e6f40ea31a4c7caab104accd *man/lavCor.Rd
225af6a64d1104380844614ca4191caa *man/lavExport.Rd
a08b3a571a4a07694853667a7a6a7372 *man/lavInspect.Rd
26119d475ad895f46006afe1e0587060 *man/lavListInspect.Rd
diff --git a/NAMESPACE b/NAMESPACE
index f1bf1ab..4884fa5 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -18,7 +18,7 @@ importFrom("stats",
"dnorm", "lm.fit", "na.omit", "nlminb", "optim", "pchisq",
"plogis", "pnorm", "qchisq", "qnorm", "quantile", "rnorm", "runif",
"sd", "terms", "uniroot", "var", "weighted.mean", "aggregate",
- "dlogis", "qlogis", "optimize",
+ "dlogis", "qlogis", "optimize", "lm",
# generics
"coef", "residuals", "resid", "fitted.values", "fitted",
"predict", "update", "anova", "vcov")
diff --git a/R/lav_bootstrap.R b/R/lav_bootstrap.R
index 47e4865..5bee90c 100644
--- a/R/lav_bootstrap.R
+++ b/R/lav_bootstrap.R
@@ -496,7 +496,7 @@ lav_bootstrap_internal <- function(object = NULL,
}
# fill in container
- t.star <- do.call("rbind", res)
+ t.star[] <- do.call("rbind", res)
# handle errors
error.idx <- which(sapply(res, function(x) is.na(x[1L])))
diff --git a/R/lav_cfa_jamesstein.R b/R/lav_cfa_jamesstein.R
new file mode 100644
index 0000000..b7a2988
--- /dev/null
+++ b/R/lav_cfa_jamesstein.R
@@ -0,0 +1,288 @@
+# James-Stein estimator
+#
+# Burghgraeve, E., De Neve, J., & Rosseel, Y. (2021). Estimating structural
+# equation models using James-Stein type shrinkage estimators. Psychometrika,
+# 86(1), 96-130.
+#
+# YR 08 Feb 2023: - first version in lavaan, cfa only (for now)
+
+lav_cfa_jamesstein <- function(S,
+ Y = NULL, # raw data
+ marker.idx = NULL,
+ lambda.nonzero.idx = NULL,
+ theta = NULL, # vector!
+ theta.bounds = TRUE,
+ aggregated = FALSE) { # aggregated?
+ # dimensions
+ nvar <- ncol(S); nfac <- length(marker.idx)
+ stopifnot(length(theta) == nvar)
+ N <- nrow(Y)
+ stopifnot(ncol(Y) == nvar)
+
+ # overview of lambda structure
+ B <- LAMBDA <- B.nomarker <- matrix(0, nvar, nfac)
+ lambda.marker.idx <- (seq_len(nfac) - 1L)*nvar + marker.idx
+ B[lambda.marker.idx ] <- LAMBDA[lambda.marker.idx ] <- 1L
+ B[lambda.nonzero.idx] <- B.nomarker[lambda.nonzero.idx] <- 1L
+
+ # Nu
+ NU <- numeric(nvar)
+
+ # do we first 'clip' the theta values so they are within standard bounds?
+ # (Question: do we need the 0.01 and 0.99 multipliers?)
+ diagS <- diag(S)
+ if(theta.bounds) {
+ # lower bound
+ lower.bound <- diagS * 0 # * 0.01
+ too.small.idx <- which(theta < lower.bound)
+ if(length(too.small.idx) > 0L) {
+ theta[ too.small.idx ] <- lower.bound[ too.small.idx ]
+ }
+
+ # upper bound
+ upper.bound <- diagS * 1 # * 0.99
+ too.large.idx <- which(theta > upper.bound)
+ if(length(too.large.idx) > 0L) {
+ theta[ too.large.idx ] <- upper.bound[ too.large.idx ]
+ }
+ }
+
+ # compute conditional expectation conditional on the scaling indicator
+ E.JS1 <- lav_cfa_jamesstein_ce(Y = Y, marker.idx = marker.idx,
+ resvars.markers = theta[marker.idx])
+
+ # compute LAMBDA
+ for(f in seq_len(nfac)) {
+ nomarker.idx <- which(B.nomarker[,f] == 1)
+ Y.nomarker.f <- Y[, nomarker.idx, drop = FALSE]
+
+ # regress no.marker.idx data on E(\eta|Y)
+ fit <- lm(Y.nomarker.f ~ E.JS1[ ,f, drop = FALSE])
+
+ # extract 'lambda' values
+ LAMBDA[nomarker.idx, f] <- drop(coef(fit)[-1,])
+
+ # (optional) extract means
+ # NU[nomarker.idx] <- drop(coef(fit)[1,])
+
+ if(aggregated) {
+
+ # local copy of 'scaling' LAMBDA
+ LAMBDA.scaling <- LAMBDA
+
+ J <- length(nomarker.idx)
+ for(j in seq_len(J)) {
+ # data without this indicator
+ j.idx <- nomarker.idx[j]
+ no.j.idx <- c(marker.idx[f],nomarker.idx[-j])
+ Y.agg <- Y[, no.j.idx, drop = FALSE]
+ Y.j <- Y[, j.idx, drop = FALSE]
+
+ # retrieve estimated values scaling JS
+ lambda.JS.scaling <- LAMBDA.scaling[no.j.idx, f, drop = FALSE]
+
+ # optimize the weights
+ starting.weights <- rep(1/J, times = J)
+ w <- optim(par = starting.weights,
+ fn = lav_cfa_jamesstein_rel,
+ data = Y.agg,
+ resvars = theta[no.j.idx])$par
+
+ # make sure the weights sum up to 1
+ w.optim <- w/sum(w)
+
+ # compute aggregated indicator using the optimal weights
+ y_agg <- t(t(w.optim) %*% t(Y.agg))
+
+ # compute error variance of the aggregated indicator
+ var_eps_agg <- drop(t(w.optim) %*%
+ diag(theta[no.j.idx], nrow = length(no.j.idx)) %*% w.optim)
+
+ # compute conditional expectation using aggregated indicator
+ tmp <- lav_cfa_jamesstein_ce(Y = y_agg, marker.idx = 1L,
+ resvars.markers = var_eps_agg)
+ CE_agg <- tmp / drop(w.optim %*% lambda.JS.scaling)
+
+ # compute factor loading
+ fit <- lm(Y.j ~ CE_agg)
+ LAMBDA[j.idx, f] <- drop(coef(fit)[-1])
+
+ # (optional) extract means
+ # NU[j.idx] <- drop(coef(fit)[1,])
+ } # j
+ } # aggregate
+
+ } # f
+
+ out <- list(lambda = LAMBDA, nu = NU)
+}
+
+# internal function to be used inside lav_optim_noniter
+# return 'x', the estimated vector of free parameters
+lav_cfa_jamesstein_internal <- function(lavobject = NULL, # convenience
+ # internal slot
+ lavmodel = NULL,
+ lavsamplestats = NULL,
+ lavpartable = NULL,
+ lavpta = NULL,
+ lavdata = NULL,
+ lavoptions = NULL,
+ theta.bounds = TRUE) {
+
+ if(!is.null(lavobject)) {
+ stopifnot(inherits(lavobject, "lavaan"))
+
+ # extract slots
+ lavmodel <- lavobject@Model
+ lavsamplestats <- lavobject@SampleStats
+ lavpartable <- lavobject@ParTable
+ lavpta <- lavobject@pta
+ lavdata <- lavobject@Data
+ lavoptions <- lavobject@Options
+ }
+
+ # no structural part!
+ if(any(lavpartable$op == "~")) {
+ stop("lavaan ERROR: JS(A) estimator only available for CFA models (for now)")
+ }
+ # no BETA matrix! (i.e., no higher-order factors)
+ if(!is.null(lavmodel@GLIST$beta)) {
+ stop("lavaan ERROR: JS(A) estimator not available for models the require a BETA matrix")
+ }
+ # no std.lv = TRUE for now
+ if(lavoptions$std.lv) {
+ stop("lavaan ERROR: JS(A) estimator not available if std.lv = TRUE")
+ }
+
+ nblocks <- lav_partable_nblocks(lavpartable)
+ stopifnot(nblocks == 1L) # for now
+ b <- 1L
+ sample.cov <- lavsamplestats@cov[[b]]
+ nvar <- nrow(sample.cov)
+ lv.names <- lavpta$vnames$lv.regular[[b]]
+ nfac <- length(lv.names)
+ marker.idx <- lavpta$vidx$lv.marker[[b]]
+ lambda.idx <- which(names(lavmodel@GLIST) == "lambda")
+ lambda.nonzero.idx <- lavmodel@m.free.idx[[lambda.idx]]
+ # only diagonal THETA for now...
+ theta.idx <- which(names(lavmodel@GLIST) == "theta") # usually '2'
+ m.theta <- lavmodel@m.free.idx[[theta.idx]]
+ nondiag.idx <- m.theta[!m.theta %in% lav_matrix_diag_idx(nvar)]
+ if(length(nondiag.idx) > 0L) {
+ warning("lavaan WARNING: this implementation of JS/JSA does not handle correlated residuals yet!")
+ }
+
+
+ # 1. obtain estimate for (diagonal elements of) THETA
+ # for now we use Spearman per factor
+ B <- matrix(0, nvar, nfac)
+ lambda.marker.idx <- (seq_len(nfac) - 1L)*nvar + marker.idx
+ B[lambda.marker.idx ] <- 1L
+ B[lambda.nonzero.idx] <- 1L
+ theta <- numeric(nvar)
+ for(f in seq_len(nfac)) {
+ ov.idx <- which(B[,f] == 1L)
+ S.fac <- sample.cov[ov.idx, ov.idx, drop = FALSE]
+ theta[ov.idx] <- lav_cfa_theta_spearman(S.fac, bounds = "wide")
+ }
+ THETA <- diag(theta, nrow = nvar)
+
+ # 2. run James-Stein algorithm
+ Y <- lavdata@X[[1]] # raw data
+ aggregated <- FALSE
+ if(lavoptions$estimator == "JSA") {
+ aggregated <- TRUE
+ }
+ out <- lav_cfa_jamesstein(S = sample.cov, Y = Y, marker.idx = marker.idx,
+ lambda.nonzero.idx = lambda.nonzero.idx,
+ theta = theta,
+ # experimental
+ theta.bounds = theta.bounds,
+ #
+ aggregated = aggregated)
+ LAMBDA <- out$lambda
+
+ # 3. PSI
+ PSI <- lav_cfa_lambdatheta2psi(lambda = LAMBDA, theta = theta,
+ S = sample.cov, mapping = "ML")
+
+ # store matrices in lavmodel@GLIST
+ lavmodel@GLIST$lambda <- LAMBDA
+ lavmodel@GLIST$theta <- THETA
+ lavmodel@GLIST$psi <- PSI
+
+ # extract free parameters only
+ x <- lav_model_get_parameters(lavmodel)
+
+ # apply bounds (if any)
+ if(!is.null(lavpartable$lower)) {
+ lower.x <- lavpartable$lower[lavpartable$free > 0]
+ too.small.idx <- which(x < lower.x)
+ if(length(too.small.idx) > 0L) {
+ x[ too.small.idx ] <- lower.x[ too.small.idx ]
+ }
+ }
+ if(!is.null(lavpartable$upper)) {
+ upper.x <- lavpartable$upper[lavpartable$free > 0]
+ too.large.idx <- which(x > upper.x)
+ if(length(too.large.idx) > 0L) {
+ x[ too.large.idx ] <- upper.x[ too.large.idx ]
+ }
+ }
+
+ x
+}
+
+
+# Conditional expectation (Section 2.1, eq. 10)
+lav_cfa_jamesstein_ce <- function(Y = NULL,
+ marker.idx = NULL,
+ resvars.markers = NULL) {
+
+ Y <- as.matrix(Y)
+
+ # sample size
+ N <- nrow(Y)
+ N1 <- N - 1
+ N3 <- N - 3
+
+ # markers only
+ Y.marker <- Y[ ,marker.idx, drop = FALSE]
+
+ # means and variances
+ MEAN <- colMeans(Y.marker, na.rm = TRUE)
+ VAR <- apply(Y.marker, 2, var, na.rm = TRUE)
+
+ # 1 - R per maker
+ oneminR <- N3 * resvars.markers / (N1 * VAR)
+
+ # R per marker
+ R <- 1 - oneminR
+
+ # create E(\eta | Y)
+ E.eta.cond.Y <- t( t(Y.marker) * R + oneminR * MEAN )
+
+ E.eta.cond.Y
+}
+
+# Reliability function used to obtain the weights (Section 4, Aggregation)
+lav_cfa_jamesstein_rel <- function(w = NULL, data = NULL, resvars = NULL) {
+
+ # construct weight vector
+ w <- matrix(w, ncol = 1)
+
+ # construct aggregated indicator: y_agg = t(w) %*% y_i
+ y_agg <- t(t(w) %*% t(data))
+
+ # calculate variance of aggregated indicator
+ var_y_agg <- var(y_agg)
+
+ # calculate error variance of the aggregated indicator
+ var_eps_agg <- t(w) %*% diag(resvars) %*% w
+
+ # reliability function to be maximized
+ rel <- (var_y_agg - var_eps_agg) %*% solve(var_y_agg)
+
+ # return value
+ return(-rel)
+}
diff --git a/R/lav_cfa_utils.R b/R/lav_cfa_utils.R
index ab7c6c4..958fb9d 100644
--- a/R/lav_cfa_utils.R
+++ b/R/lav_cfa_utils.R
@@ -101,6 +101,80 @@ lav_cfa_lambda2thetapsi <- function(lambda = NULL, S = NULL, S.inv = NULL,
list(lambda = LAMBDA, theta = theta.nobounds, psi = PSI)
}
+# compute PSI, given lambda and theta using either ULS, GLS, ML
+# this function assumes:
+# - THETA is diagonal
+# - PSI is unrestricted
+#
+# YR 08 Mar 2023: - first version
+lav_cfa_lambdatheta2psi <- function(lambda = NULL, theta = NULL, # vector!
+ S = NULL, S.inv = NULL,
+ mapping = "ML", nobs = 20L) {
+ LAMBDA <- as.matrix(lambda)
+ nvar <- nrow(LAMBDA); nfac <- ncol(LAMBDA)
+
+ theta.nobounds <- theta
+
+ # ALWAYS check bounds for theta to compute PSI
+ diagS <- diag(S)
+ # lower bound
+ lower.bound <- diagS * 0 # * 0.01
+ too.small.idx <- which(theta < lower.bound)
+ if(length(too.small.idx) > 0L) {
+ theta[ too.small.idx ] <- lower.bound[ too.small.idx ]
+ }
+
+ # upper bound
+ upper.bound <- diagS * 1 # * 0.99
+ too.large.idx <- which(theta > upper.bound)
+ if(length(too.large.idx) > 0L) {
+ theta[ too.large.idx ] <- upper.bound[ too.large.idx ]
+ }
+
+ # psi
+ diag.theta <- diag(theta, nvar)
+ lambda <- try(lav_matrix_symmetric_diff_smallest_root(S, diag.theta),
+ silent = TRUE)
+ if(inherits(lambda, "try-error")) {
+ warning("lavaan WARNING: failed to compute lambda")
+ SminTheta <- S - diag.theta # and hope for the best
+ } else {
+ cutoff <- 1 + 1/(nobs - 1)
+ if(lambda < cutoff) {
+ lambda.star <- lambda - 1/(nobs - 1)
+ SminTheta <- S - lambda.star * diag.theta
+ } else {
+ SminTheta <- S - diag.theta
+ }
+ }
+
+ # mapping matrix
+ if(mapping == "ML") {
+ Ti <- 1/theta
+ zero.theta.idx <- which(abs(theta) < 0.01) # be conservative
+ if(length(zero.theta.idx) > 0L) {
+ Ti[zero.theta.idx] <- 1
+ }
+ M <- solve(t(LAMBDA) %*% diag(Ti, nvar) %*% LAMBDA) %*% t(LAMBDA) %*% diag(Ti, nvar)
+ } else if(mapping == "GLS") {
+ if(is.null(S.inv)) {
+ S.inv <- try(solve(S), silent = TRUE)
+ }
+ if(inherits(S.inv, "try-error")) {
+ M <- tcrossprod(solve(crossprod(LAMBDA)), LAMBDA)
+ } else {
+ M <- solve(t(LAMBDA) %*% S.inv %*% LAMBDA) %*% t(LAMBDA) %*% S.inv
+ }
+ } else if(mapping == "ULS") {
+ M <- tcrossprod(solve(crossprod(LAMBDA)), LAMBDA)
+ }
+
+ # compute PSI
+ PSI <- M %*% SminTheta %*% t(M)
+
+ PSI
+}
+
# compute theta elements for a 1-factor model
lav_cfa_theta_spearman <- function(S, bounds = "wide") {
p <- ncol(S); out <- numeric(p)
diff --git a/R/lav_cor.R b/R/lav_cor.R
index a7bdfaa..d1707e4 100644
--- a/R/lav_cor.R
+++ b/R/lav_cor.R
@@ -12,6 +12,7 @@ lavCor <- function(object,
group = NULL,
missing = "listwise",
ov.names.x = NULL,
+ sampling.weights = NULL,
# lavaan options
se = "none",
test = "none",
@@ -44,6 +45,14 @@ lavCor <- function(object,
}
}
+ # extract sampling.weights.normalization from dots (for lavData() call)
+ dots <- list(...)
+ sampling.weights.normalization <- "total"
+ if (!is.null(dots$sampling.weights.normalization)) {
+ sampling.weights.normalization <- dots$sampling.weights.normalization
+ }
+
+
# check object class
if(inherits(object, "lavData")) {
lav.data <- object
@@ -54,6 +63,9 @@ lavCor <- function(object,
if(!is.null(group)) {
NAMES <- NAMES[- match(group, NAMES)]
}
+ if(!is.null(sampling.weights)) {
+ NAMES <- NAMES[- match(sampling.weights, NAMES)]
+ }
if(is.logical(ordered)) {
ordered.flag <- ordered
if(ordered.flag) {
@@ -81,8 +93,10 @@ lavCor <- function(object,
}
lav.data <- lavData(data = object, group = group,
ov.names = NAMES, ordered = ordered,
+ sampling.weights = sampling.weights,
ov.names.x = ov.names.x,
- lavoptions = list(missing = missing))
+ lavoptions = list(missing = missing,
+ sampling.weights.normalization = sampling.weights.normalization))
} else {
stop("lavaan ERROR: lavCor can not handle objects of class ",
paste(class(object), collapse= " "))
@@ -98,8 +112,7 @@ lavCor <- function(object,
}
}
- # extract partable options from dots
- dots <- list(...)
+ # extract more partable options from dots
meanstructure <- FALSE
fixed.x <- FALSE
mimic <- "lavaan"
@@ -120,7 +133,6 @@ lavCor <- function(object,
conditional.x <- dots$conditional.x
}
-
# override, only for backwards compatibility (eg moments() in JWileymisc)
#if(missing %in% c("ml", "fiml")) {
# meanstructure = TRUE
@@ -133,6 +145,7 @@ lavCor <- function(object,
lavoptions = list(meanstructure = meanstructure,
fixed.x = fixed.x,
conditional.x = conditional.x,
+ # sampling.weights.normalization = sampling.weights.normalization,
group.w.free = FALSE,
missing = missing,
correlation = categorical,
diff --git a/R/lav_fit_measures.R b/R/lav_fit_measures.R
index b537705..04c52bb 100644
--- a/R/lav_fit_measures.R
+++ b/R/lav_fit_measures.R
@@ -54,7 +54,7 @@ fitMeasures.efaList <- fitmeasures.efaList <- function(object,
rmsea.close.h0 = 0.05,
rmsea.notclose.h0 = 0.08,
cat.check.pd = TRUE),
- vector = "list", ...) {
+ output = "list", ...) {
# kill object$loadings if present
object[["loadings"]] <- NULL
diff --git a/R/lav_matrix.R b/R/lav_matrix.R
index aa42299..6072157 100644
--- a/R/lav_matrix.R
+++ b/R/lav_matrix.R
@@ -741,8 +741,10 @@ lav_matrix_commutation <- .com1
# = permuting the rows of A
lav_matrix_commutation_pre <- function(A = matrix(0,0,0)) {
+ A <- as.matrix(A)
+
# number of rows of A
- n2 <- NROW(A)
+ n2 <- nrow(A)
# K_nn only (n2 = m * n)
stopifnot(sqrt(n2) == round(sqrt(n2)))
@@ -758,6 +760,56 @@ lav_matrix_commutation_pre <- function(A = matrix(0,0,0)) {
OUT
}
+# compute A %*% K_n without explicitly computing K
+# K_n = K_nn, so sqrt(ncol(A)) must be an integer!
+# = permuting the columns of A
+lav_matrix_commutation_post <- function(A = matrix(0,0,0)) {
+
+ A <- as.matrix(A)
+
+ # number of columns of A
+ n2 <- ncol(A)
+
+ # K_nn only (n2 = m * n)
+ stopifnot(sqrt(n2) == round(sqrt(n2)))
+
+ # dimension
+ n <- sqrt(n2)
+
+ # compute col indices
+ #row.idx <- as.integer(t(matrix(1:n2, n, n)))
+ col.idx <- rep(1:n, each = n) + (0:(n-1L))*n
+
+ OUT <- A[, col.idx, drop = FALSE]
+ OUT
+}
+
+# compute K_n %*% A %*% K_n without explicitly computing K
+# K_n = K_nn, so sqrt(ncol(A)) must be an integer!
+# = permuting both the rows AND columns of A
+lav_matrix_commutation_pre_post <- function(A = matrix(0,0,0)) {
+
+ A <- as.matrix(A)
+
+ # number of columns of A
+ n2 <- NCOL(A)
+
+ # K_nn only (n2 = m * n)
+ stopifnot(sqrt(n2) == round(sqrt(n2)))
+
+ # dimension
+ n <- sqrt(n2)
+
+ # compute col indices
+ row.idx <- rep(1:n, each = n) + (0:(n-1L))*n
+ col.idx <- row.idx
+
+ OUT <- A[row.idx, col.idx, drop = FALSE]
+ OUT
+}
+
+
+
# compute K_mn %*% A without explicitly computing K
# = permuting the rows of A
lav_matrix_commutation_mn_pre <- function(A, m = 1L, n = 1L) {
diff --git a/R/lav_optim_noniter.R b/R/lav_optim_noniter.R
index 6346dd0..910cb70 100644
--- a/R/lav_optim_noniter.R
+++ b/R/lav_optim_noniter.R
@@ -50,7 +50,7 @@ lav_optim_noniter <- function(lavmodel = NULL, lavsamplestats = NULL,
x <- try(lav_cfa_fabin_internal(lavmodel = lavmodel,
lavsamplestats = lavsamplestats, lavpartable = lavpartable,
lavpta = lavpta, lavoptions = lavoptions), silent = TRUE)
- } else if(lavoptions$estimator == "GUTTMAN1952") {
+ } else if(lavoptions$estimator == "MGM") {
x <- try(lav_cfa_guttman1952_internal(lavmodel = lavmodel,
lavsamplestats = lavsamplestats, lavpartable = lavpartable,
lavpta = lavpta, lavoptions = lavoptions), silent = TRUE)
@@ -58,6 +58,15 @@ lav_optim_noniter <- function(lavmodel = NULL, lavsamplestats = NULL,
x <- try(lav_cfa_bentler1982_internal(lavmodel = lavmodel,
lavsamplestats = lavsamplestats, lavpartable = lavpartable,
lavpta = lavpta, lavoptions = lavoptions), silent = TRUE)
+ } else if(lavoptions$estimator %in% c("JS", "JSA")) {
+ x <- try(lav_cfa_jamesstein_internal(lavmodel = lavmodel,
+ lavsamplestats = lavsamplestats, lavpartable = lavpartable,
+ lavdata = lavdata,
+ lavpta = lavpta, lavoptions = lavoptions), silent = TRUE)
+ } else if(lavoptions$estimator == "BENTLER1982") {
+ x <- try(lav_cfa_bentler1982_internal(lavmodel = lavmodel,
+ lavsamplestats = lavsamplestats, lavpartable = lavpartable,
+ lavpta = lavpta, lavoptions = lavoptions), silent = TRUE)
} else {
warning("lavaan WARNING: unknown (noniterative) estimator: ",
lavoptions$estimator, " (returning starting values)")
diff --git a/R/lav_options.R b/R/lav_options.R
index 9ce61bd..15865c2 100644
--- a/R/lav_options.R
+++ b/R/lav_options.R
@@ -1203,11 +1203,14 @@ lav_options_set <- function(opt = NULL) {
##################################################################
- # FABIN, GUTTMAN1952, BENTLER, ... #
+ # FABIN, MULTIPLE-GROUP-METHOD (MGM( BENTLER, ... #
##################################################################
} else if(opt$estimator %in% c("fabin", "fabin2", "fabin3",
- "guttman", "gutman", "guttman1952",
- "bentler", "bentler1982")) {
+ "mgm", "guttman", "gutman", "guttman1952",
+ "js", "jsa", "james-stein", "james.stein",
+ "james-stein-aggregated",
+ "james.stein.aggregated",
+ "bentler", "bentler1982")) {
# experimental, for cfa or sam only
# sample.cov.rescale
@@ -1222,10 +1225,16 @@ lav_options_set <- function(opt = NULL) {
# estimator
if(opt$estimator == "fabin") {
opt$estimator <- "FABIN2"
- } else if(opt$estimator %in% c("guttman", "gutman", "guttmann")) {
- opt$estimator <- "GUTTMAN1952"
+ } else if(opt$estimator %in% c("mgm", "guttman", "gutman",
+ "guttmann", "guttman1952")) {
+ opt$estimator <- "MGM"
} else if(opt$estimator %in% c("bentler", "bentler1982")) {
opt$estimator <- "BENTLER1982"
+ } else if(opt$estimator %in% c("js", "james-stein", "james.stein")) {
+ opt$estimator <- "JS"
+ } else if(opt$estimator %in% c("jsa", "james-stein-aggregated",
+ "james.stein.aggregated")) {
+ opt$estimator <- "JSA"
} else {
opt$estimator <- toupper(opt$estimator)
}
@@ -1284,8 +1293,8 @@ lav_options_set <- function(opt = NULL) {
}
}
- # options for guttman1952
- if(opt$estimator == "GUTTMAN1952") {
+ # options for guttman1952 multiple group method
+ if(opt$estimator == "MGM") {
if(is.null(opt$estimator.args)) {
opt$estimator.args <- list(psi.mapping = FALSE)
} else {
diff --git a/R/lav_partable_vnames.R b/R/lav_partable_vnames.R
index 0178762..5892bb1 100644
--- a/R/lav_partable_vnames.R
+++ b/R/lav_partable_vnames.R
@@ -682,7 +682,8 @@ lav_partable_vnames <- function(partable, type = NULL, ...,
target.idx <- which(x[[b]] %in% ov.names.data)
if(length(target.idx) > 0L) {
new.ov <-
- ov.names.data[match(x[[b]], ov.names.data)]
+ ov.names.data[sort(match(x[[b]],
+ ov.names.data))]
# rm NA's (eg lv's in eqs.y)
na.idx <- which(is.na(new.ov))
if(length(na.idx) > 0L) {
diff --git a/R/lav_print.R b/R/lav_print.R
index 805567e..3bc0257 100644
--- a/R/lav_print.R
+++ b/R/lav_print.R
@@ -52,6 +52,7 @@ print.lavaan.matrix.symmetric <- function(x, ..., nd = 3L, shift = 0L,
# in package nlme
x <- as.matrix(x) # just in case
y <- x; y <- unclass(y)
+ attributes(y)[c("header","footer")] <- NULL
ll <- lower.tri(x, diag = TRUE)
y[ ll] <- format(round(x[ll], digits = nd))
y[!ll] <- ""
@@ -74,7 +75,17 @@ print.lavaan.matrix.symmetric <- function(x, ..., nd = 3L, shift = 0L,
rownames(y) <- empty.string
}
}
+
+ if(!is.null(attr(x, "header"))) {
+ cat("\n", attr(x, "header"), "\n\n", sep = "")
+ }
+
print(y, ..., quote = FALSE, right = TRUE)
+
+ if(!is.null(attr(x, "footer"))) {
+ cat("\n", attr(x, "footer"), "\n\n", sep = "")
+ }
+
invisible(x)
}
@@ -82,6 +93,7 @@ print.lavaan.matrix.symmetric <- function(x, ..., nd = 3L, shift = 0L,
print.lavaan.matrix <- function(x, ..., nd = 3L, shift = 0L) {
x <- as.matrix(x) # just in case
y <- unclass(x)
+ attributes(y)[c("header","footer")] <- NULL
if(!is.null(colnames(x))) {
colnames(y) <- abbreviate(colnames(x), minlength = nd + 3L)
}
@@ -93,15 +105,29 @@ print.lavaan.matrix <- function(x, ..., nd = 3L, shift = 0L) {
rownames(y) <- empty.string
}
}
+ if(!is.null(attr(x, "header"))) {
+ cat("\n", attr(x, "header"), "\n\n", sep = "")
+ }
+
print( round(y, nd), right = TRUE, ... )
+
+ if(!is.null(attr(x, "footer"))) {
+ cat("\n", attr(x, "footer"), "\n\n", sep = "")
+ }
+
invisible(x)
}
print.lavaan.vector <- function(x, ..., nd = 3L, shift = 0L) {
y <- unclass(x)
+ attributes(y)[c("header","footer")] <- NULL
#if(!is.null(names(x))) {
# names(y) <- abbreviate(names(x), minlength = nd + 3)
#}
+ if(!is.null(attr(x, "header"))) {
+ cat("\n", attr(x, "header"), "\n\n", sep = "")
+ }
+
if(shift > 0L) {
empty.string <- strrep(x = " ", times = shift)
tmp <- format(y, digits = nd, width = 2L + nd)
@@ -110,6 +136,11 @@ print.lavaan.vector <- function(x, ..., nd = 3L, shift = 0L) {
} else {
print( round(y, nd), right = TRUE, ... )
}
+
+ if(!is.null(attr(x, "footer"))) {
+ cat("\n", attr(x, "footer"), "\n\n", sep = "")
+ }
+
invisible(x)
}
diff --git a/debian/changelog b/debian/changelog
index bc3d87d..ad632f3 100644
--- a/debian/changelog
+++ b/debian/changelog
@@ -1,3 +1,9 @@
+r-cran-lavaan (0.6.15-1) UNRELEASED; urgency=low
+
+ * New upstream release.
+
+ -- Debian Janitor <janitor@jelmer.uk> Tue, 20 Jun 2023 03:49:21 -0000
+
r-cran-lavaan (0.6.14-1) unstable; urgency=medium
* Team upload.
diff --git a/man/lavCor.Rd b/man/lavCor.Rd
index 6386bb3..cbfeddf 100644
--- a/man/lavCor.Rd
+++ b/man/lavCor.Rd
@@ -6,7 +6,8 @@ Fit an unrestricted model to compute polychoric, polyserial and/or Pearson
correlations.}
\usage{
lavCor(object, ordered = NULL, group = NULL, missing = "listwise",
- ov.names.x = NULL, se = "none", test = "none",
+ ov.names.x = NULL, sampling.weights = NULL,
+ se = "none", test = "none",
estimator = "two.step", baseline = FALSE, ...,
cor.smooth = FALSE, cor.smooth.tol = 1e-04, output = "cor")
}
@@ -28,6 +29,8 @@ lavCor(object, ordered = NULL, group = NULL, missing = "listwise",
(and mean vector). If \code{"pairwise"}, pairwise deletion is used. If
\code{"default"}, the value is set depending on the estimator and the
mimic option.}
+\item{sampling.weights}{Only used if \code{object} is a \code{data.frame}.
+ Specify a variable containing sampling weights.}
\item{ov.names.x}{Only used if \code{object} is a \code{data.frame}. Specify
variables that need to be treated as exogenous. Only used if at least
one variable is declared as ordered.}