Learn R Programming

psychonetrics (version 0.17.8)

ml_lvm: Multi-level latent variable model family

Description

This family is the two-level random intercept variant of the lvm model family. It is mostly a special case of the dlvm1 family, with the addition of structural effects rather than temporal effects in the beta matrix.

Usage

ml_lnm(...)
ml_rnm(...)
ml_lrnm(...)
ml_lvm(data, lambda, clusters, within_latent = c("cov",
                   "chol", "prec", "ggm", "cor"), within_residual = c("cov",
                   "chol", "prec", "ggm", "cor"), between_latent = c("cov",
                   "chol", "prec", "ggm", "cor"), between_residual = c("cov",
                   "chol", "prec", "ggm", "cor"), beta_within = "zero",
                   beta_between = "zero", omega_zeta_within = "full",
                   delta_zeta_within = "full", kappa_zeta_within =
                   "full", sigma_zeta_within = "full",
                   lowertri_zeta_within = "full", omega_epsilon_within =
                   "zero", delta_epsilon_within = "diag",
                   kappa_epsilon_within = "diag", sigma_epsilon_within =
                   "diag", lowertri_epsilon_within = "diag",
                   omega_zeta_between = "full", delta_zeta_between =
                   "full", kappa_zeta_between = "full",
                   sigma_zeta_between = "full", lowertri_zeta_between =
                   "full", omega_epsilon_between = "zero",
                   delta_epsilon_between = "diag", kappa_epsilon_between
                   = "diag", sigma_epsilon_between = "diag",
                   lowertri_epsilon_between = "diag", nu, nu_eta,
                   identify = TRUE, identification = c("loadings",
                   "variance"), vars, latents, groups, equal = "none",
                   baseline_saturated = TRUE, estimator = c("default",
                   "FIML", "ML"),
                   optimizer, storedata = FALSE, verbose =
                   FALSE, standardize = c("none", "z", "quantile"),
                   sampleStats, bootstrap = FALSE, boot_sub,
                   boot_resample,
                   rho_zeta_within = "full", SD_zeta_within = "full",
                   rho_epsilon_within = "zero", SD_epsilon_within = "diag",
                   rho_zeta_between = "full", SD_zeta_between = "full",
                   rho_epsilon_between = "zero", SD_epsilon_between = "diag")

Value

An object of the class psychonetrics (psychonetrics-class)

Arguments

data

A data frame encoding the data used in the analysis. Must be a raw dataset.

lambda

A model matrix encoding the factor loading structure. Each row indicates an indicator and each column a latent. A 0 encodes a fixed to zero element, a 1 encoding a free to estimate element, and higher integers encoding equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix. Could also be the result of simplestructure.

clusters

A string indicating the variable in the dataset that describes group membership.

within_latent

The type of within-person latent contemporaneous model to be used.

within_residual

The type of within-person residual model to be used.

between_latent

The type of between-person latent model to be used.

between_residual

The type of between-person residual model to be used.

beta_within

A model matrix encoding the within-cluster structural. A 0 encodes a fixed to zero element, a 1 encoding a free to estimate element, and higher integers encoding equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix. Defaults to "zero".

beta_between

A model matrix encoding the between-cluster structural. A 0 encodes a fixed to zero element, a 1 encoding a free to estimate element, and higher integers encoding equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix. Defaults to "zero".

omega_zeta_within

Only used when within_latent = "ggm". Can be "full", "zero", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.

delta_zeta_within

Only used when within_latent = "ggm". Can be "diag", "zero" (not recommended), or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.

kappa_zeta_within

Only used when within_latent = "prec". Can be "full", "diag", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.

sigma_zeta_within

Only used when within_latent = "cov". Can be "full", "diag", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.

lowertri_zeta_within

Only used when within_latent = "chol". Can be "full", "diag", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.

omega_epsilon_within

Only used when within_residual = "ggm". Can be "full", "zero", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.

delta_epsilon_within

Only used when within_residual = "ggm". Can be "diag", "zero" (not recommended), or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.

kappa_epsilon_within

Only used when within_residual = "prec". Can be "full", "diag", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.

sigma_epsilon_within

Only used when within_residual = "cov". Can be "full", "diag", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.

lowertri_epsilon_within

Only used when within_residual = "chol". Can be "full", "diag", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.

omega_zeta_between

Only used when between_latent = "ggm". Can be "full", "zero", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.

delta_zeta_between

Only used when between_latent = "ggm". Can be "diag", "zero" (not recommended), or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.

kappa_zeta_between

Only used when between_latent = "prec". Can be "full", "diag", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.

sigma_zeta_between

Only used when between_latent = "cov". Can be "full", "diag", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.

lowertri_zeta_between

Only used when between_latent = "chol". Can be "full", "diag", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.

omega_epsilon_between

Only used when between_residual = "ggm". Can be "full", "zero", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.

delta_epsilon_between

Only used when between_residual = "ggm". Can be "diag", "zero" (not recommended), or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.

kappa_epsilon_between

Only used when between_residual = "prec". Can be "full", "diag", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.

sigma_epsilon_between

Only used when between_residual = "cov". Can be "full", "diag", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.

lowertri_epsilon_between

Only used when between_residual = "chol". Can be "full", "diag", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.

nu

Optional vector encoding the intercepts of the observed variables. Set elements to 0 to indicate fixed to zero constrains, 1 to indicate free intercepts, and higher integers to indicate equality constraints. For multiple groups, this argument can be a list or array with each element/column encoding such a vector.

nu_eta

Optional vector encoding the intercepts of the latent variables. Set elements to 0 to indicate fixed to zero constrains, 1 to indicate free intercepts, and higher integers to indicate equality constraints. For multiple groups, this argument can be a list or array with each element/column encoding such a vector.

identify

Logical, should the model be automatically identified?

identification

Type of identification used. "loadings" to fix the first factor loadings to 1, and "variance" to fix the diagonal of the latent variable model matrix (sigma_zeta, lowertri_zeta, delta_zeta or kappa_zeta) to 1.

vars

An optional character vector with names of the variables used.

latents

An optional character vector with names of the latent variables.

groups

An optional string indicating the name of the group variable in data.

equal

A character vector indicating which matrices should be constrained equal across groups.

baseline_saturated

A logical indicating if the baseline and saturated model should be included. Mostly used internally and NOT Recommended to be used manually.

estimator

The estimator to be used. "ML" and "FIML" maximize the same (exact) two-level maximum-likelihood objective and converge to the same solution; they differ only in how that likelihood is computed, exactly as for the single-level model families in this package, where "ML" fits the model to summary statistics and "FIML" fits it row by row to the raw data. "ML" (experimental) evaluates the likelihood from two-level sufficient statistics (pooled within-cluster covariance and per-cluster-size between moments; McDonald & Goldstein, 1989), which for complete data is much faster and does not slow down as clusters grow. "FIML" (wide-format full-information maximum likelihood) stacks each cluster's members into one wide observation; it was the only estimator before version 0.15.31. Both handle within-cluster missing data ("ML" through a per-pattern, per-cluster likelihood; see Details). This is NOT the approximate "MUML" / pseudo-balanced estimator: the sufficient-statistics evaluation is exact for unbalanced designs, and the two estimators' fit-function values coincide to numerical precision (verified to ~1e-13). "default" selects "ML" when the data are complete and the largest cluster contains more than 5 units, and "FIML" otherwise; with missing data the default stays on the long-standing "FIML" path, so the (R-only) missing-data "ML" estimator is used only when estimator = "ML" is requested explicitly.

optimizer

The optimizer to be used. Can be one of "nlminb" (the default R nlminb function), "ucminf" (from the optimr package), "nloptr_TNEWTON" (preconditioned truncated Newton via nloptr), and "LBFGS++" (pure C++ L-BFGS-B). Defaults to "nlminb".

storedata

Logical, should the raw data be stored? Needed for bootstrapping (see bootstrap).

verbose

Logical, should progress be printed to the console?

standardize

Which standardization method should be used? "none" (default) for no standardization, "z" for z-scores, and "quantile" for a non-parametric transformation to the quantiles of the marginal standard normal distribution.

sampleStats

An optional sample statistics object. Mostly used internally.

bootstrap

Should the data be bootstrapped? If TRUE the data are resampled and a bootstrap sample is created. These must be aggregated using aggregate_bootstraps! Can be TRUE or FALSE. Can also be "nonparametric" (which sets boot_sub = 1 and boot_resample = TRUE) or "case" (which sets boot_sub = 0.75 and boot_resample = FALSE).

boot_sub

Proportion of cases to be subsampled (round(boot_sub * N)).

boot_resample

Logical, should the bootstrap be with replacement (TRUE) or without replacement (FALSE)

rho_zeta_within

Only used when within_latent = "cor". Correlation matrix specification ("full", "zero", or a pattern matrix; see omega_zeta_within).

SD_zeta_within

Only used when within_latent = "cor". Diagonal standard-deviation matrix specification ("full", "diag", or a pattern matrix; see delta_zeta_within).

rho_epsilon_within

Only used when within_residual = "cor". Correlation matrix specification ("zero", "full", or a pattern matrix).

SD_epsilon_within

Only used when within_residual = "cor". Diagonal standard-deviation matrix specification ("diag", "full", or a pattern matrix).

rho_zeta_between

Only used when between_latent = "cor". Correlation matrix specification ("full", "zero", or a pattern matrix).

SD_zeta_between

Only used when between_latent = "cor". Diagonal standard-deviation matrix specification ("full", "diag", or a pattern matrix).

rho_epsilon_between

Only used when between_residual = "cor". Correlation matrix specification ("zero", "full", or a pattern matrix).

SD_epsilon_between

Only used when between_residual = "cor". Diagonal standard-deviation matrix specification ("diag", "full", or a pattern matrix).

...

Arguments sent to 'ml_lvm'

Author

Sacha Epskamp <mail@sachaepskamp.com>

Details

The "ML" and "FIML" estimators are two computational routes to the same exact two-level maximum-likelihood estimate, not two different statistical methods and not an exact-versus-approximate pair: they maximize the same log-likelihood and (up to optimizer tolerance) return identical estimates, standard errors and fit. The naming follows the convention used throughout psychonetrics -- "ML" fits to summary statistics while "FIML" fits row by row to the raw data -- and does not by itself indicate complete versus missing data, since both routes can accommodate within-cluster missingness. Choose "ML" for speed on (near-)complete data with non-trivial cluster sizes and "FIML" (or "default") otherwise.

With estimator = "FIML", every cluster is treated as one wide-format observation (members stacked side by side) and the model is estimated with full-information maximum likelihood over the missingness patterns. This handles missing data, but its cost grows quickly with the cluster size.

With estimator = "ML" (experimental), the same two-level random-intercept likelihood is instead evaluated from sufficient statistics (McDonald & Goldstein, 1989; Muthen, 1990) when the data are complete. Writing \(\Sigma_W\) and \(\Sigma_B\) for the implied within- and between-cluster covariance matrices and \(\mu\) for the implied mean vector, minus twice the log-likelihood equals $$(N - J)\left[\ln|\Sigma_W| + \mathrm{tr}(\Sigma_W^{-1} S_{PW})\right] + \sum_s m_s \left[\ln|\Sigma_s| + n_s \mathrm{tr}(\Sigma_s^{-1} A_s)\right] + N p \ln(2\pi),$$ in which \(N\) is the number of units, \(J\) the number of clusters, \(S_{PW}\) the pooled within-cluster covariance matrix, and, per distinct cluster size \(n_s\) with \(m_s\) clusters, \(\Sigma_s = \Sigma_W + n_s \Sigma_B\) and \(A_s\) the second-moment matrix of the corresponding cluster means around \(\mu\). The two estimators optimize numerically identical objectives (the fit function values coincide at identical parameter values), so they converge to the same solution; the same likelihood is also used by lavaan for two-level models (Rosseel, 2012). The sufficient statistics are computed once, making "ML" much faster than "FIML" when clusters are large, and its cost does not depend on the cluster sizes. For complete data, standard errors and modification indices for "ML" use the analytic EXPECTED Fisher information. Note that lavaan uses the OBSERVED information by default for multilevel models, so standard errors reported by lavaan typically differ by a few percent; with lavaan::sem(..., information = "expected") the standard errors agree.

When the data contain within-cluster missing values (assumed missing at random), estimator = "ML" switches to a per-pattern, per-cluster evaluation of the same two-level likelihood: for each missingness pattern the within-cluster precision and log-determinant are obtained by a symmetric (Schur) inverse update of \(\Sigma_W^{-1}\), and these are accumulated per cluster into the quantities entering minus twice the log-likelihood, $$\sum_P f_P \ln|\Sigma_W^{(P)}| + \sum_j \ln|I + \Sigma_B A_j| + \sum_i (y_i - \mu)'\Sigma_W^{(P_i)-1}(y_i - \mu) - \sum_j p_j'(I + \Sigma_B A_j)^{-1}\Sigma_B\, p_j + n_{\mathrm{obs}} \ln(2\pi),$$ with \(A_j = \sum_{i \in j}\Sigma_W^{(P_i)-1}\), \(p_j = \sum_{i \in j}\Sigma_W^{(P_i)-1}(y_i - \mu)\) (observed coordinates only), \(f_P\) the number of units in pattern \(P\), and \(n_{\mathrm{obs}}\) the number of observed values. This is the same two-level missing-data likelihood that lavaan::sem(..., cluster = , missing = "ml") maximizes (Rosseel, 2012), and the two agree to optimizer tolerance. The missing-data path is evaluated in R (no C++ acceleration) and reduces exactly to the sufficient-statistics likelihood above when no values are missing; its standard errors use the numeric (observed-information-based) Fisher information.

The saturated (unrestricted) and baseline (independence) reference models used by runmodel for the fit measures are themselves estimated as two-level models: the saturated model freely estimates the means and the full within- and between-cluster covariance matrices, and the baseline model only estimates means and the variances at both levels. These are the same reference models lavaan uses for two-level SEM, so the saturated and baseline log-likelihoods, the model and baseline chi-square, the degrees of freedom, the CFI and the TLI agree with lavaan (up to optimization tolerances), provided the lavaan model uses the same parameterization (ml_lvm shares the factor loadings across the within and between levels, so the matching lavaan model must constrain its level-1 and level-2 loadings to be equal). One convention difference remains: psychonetrics evaluates sample-size-dependent fit measures, such as the RMSEA (with its confidence interval) and the BIC, using the number of independent observations, which in these models is the number of clusters \(J\), whereas lavaan uses the total number of level-1 units \(N\). Because the RMSEA equals \(\sqrt{\max((\chi^2/n)/\mathit{df} - 1/n,\, 0)}\) with \(n\) the number of independent observations, and the chi-square and degrees of freedom otherwise agree, the RMSEA reported by psychonetrics is larger than the value reported by lavaan by exactly a factor \(\sqrt{N/J}\). This applies to both estimators of the ml_lvm family.

References

McDonald, R. P., & Goldstein, H. (1989). Balanced versus unbalanced designs for linear structural relations in two-level data. British Journal of Mathematical and Statistical Psychology, 42(2), 215-232.

Muthen, B. O. (1990). Mean and covariance structure analysis of hierarchical data. UCLA Statistics Series #62.

Rosseel, Y. (2012). lavaan: An R Package for Structural Equation Modeling. Journal of Statistical Software, 48(2), 1-36. tools:::Rd_expr_doi("10.18637/jss.v048.i02")

Examples

Run this code
# \donttest{
# Simulate two-level data with 2 latent variables on both levels:
set.seed(1)
J <- 50 # clusters
nj <- 10 # units per cluster
p <- 6
Lambda <- matrix(0, p, 2)
Lambda[1:3,1] <- c(1, 0.8, 0.7)
Lambda[4:6,2] <- c(1, 0.9, 0.6)
PsiW <- matrix(c(1,0.3,0.3,1), 2)
PsiB <- matrix(c(0.5,0.2,0.2,0.4), 2)
data <- do.call(rbind, lapply(seq_len(J), function(j){
  bj <- MASS::mvrnorm(1, rep(0,2), PsiB)
  ej <- rnorm(p, sd = 0.4) # between-level residuals
  t(replicate(nj, as.vector(Lambda %*% (bj + MASS::mvrnorm(1, rep(0,2), PsiW))) +
      ej + rnorm(p, sd = 0.7)))
}))
data <- as.data.frame(data)
names(data) <- paste0("y", 1:p)
data$cluster <- rep(seq_len(J), each = nj)

# Lambda model matrix (shared across levels):
lambda <- matrix(0, p, 2)
lambda[1:3,1] <- lambda[4:6,2] <- 1

# Two-level ML estimator (fast, complete data):
mod <- ml_lvm(data, lambda = lambda, clusters = "cluster", estimator = "ML")
mod <- runmodel(mod)
mod

# The FIML estimator (slower; supports missing data) converges to the
# same solution:
mod_fiml <- ml_lvm(data, lambda = lambda, clusters = "cluster",
                   estimator = "FIML")
mod_fiml <- runmodel(mod_fiml)
# }

Run the code above in your browser using DataLab