# Latent variable model (lvm)
lvm <- function(
data, # Dataset
lambda, # only required non-missing matrix
latent = c("cov","chol","prec","ggm"), # Maybe add cor at some point, but not now
residual = c("cov","chol","prec","ggm"),
# Latent matrices:
sigma_zeta = "full",
kappa_zeta = "full", # Precision
omega_zeta = "full", # Partial correlations
lowertri_zeta = "full", # Cholesky
delta_zeta = "full", # Used for both ggm and pcor
# Residual matrices:
sigma_epsilon = "diag", #
kappa_epsilon = "diag", # Precision
omega_epsilon = "zero", # Partial correlations
lowertri_epsilon = "diag", # Cholesky
delta_epsilon = "diag", # Used for both ggm and pcor
# Beta:
beta = "zero",
# Mean structure:
nu,
nu_eta,
# Identification:
identify = TRUE,
identification = c("loadings","variance"),
# Rest:
vars, # character indicating the variables Extracted if missing from data - group variable
latents, # Name of latent varianles
groups, # ignored if missing. Can be character indicating groupvar, or vector with names of groups
covs, # alternative covs (array nvar * nvar * ngroup)
means, # alternative means (matrix nvar * ngroup)
nobs, # Alternative if data is missing (length ngroup)
missing = "listwise",
equal = "none", # Can also be any of the matrices
baseline_saturated = TRUE, # Leave to TRUE! Only used to stop recursive calls
estimator = "ML",
optimizer,
storedata = FALSE,
WLS.W,
covtype = c("choose","ML","UB"),
standardize = c("none","z","quantile"),
sampleStats,
verbose=FALSE,
simplelambdastart = FALSE
){
rawts = FALSE
if (rawts){
warning("'rawts' is only included for testing purposes. Please do not use!")
}
# Type:
latent <- match.arg(latent)
residual <- match.arg(residual)
# Identification:
identification <- match.arg(identification)
# WLS weights:
if (missing(WLS.W)){
WLS.W <- ifelse(!estimator %in% c("WLS","ULS","DWLS"), "none",
switch(estimator,
"WLS" = "full",
"ULS" = "identity",
"DWLS" = "diag"
))
}
# Obtain sample stats:
if (missing(sampleStats)){
sampleStats <- samplestats(data = data,
vars = vars,
groups = groups,
covs = covs,
means = means,
nobs = nobs,
missing = ifelse(estimator == "FIML","pairwise",missing),
rawts = rawts,
fimldata = estimator == "FIML",
storedata = storedata,
covtype=covtype,
weightsmatrix = WLS.W,
verbose=verbose,
standardize=standardize)
}
# Check some things:
nNode <- nrow(sampleStats@variables)
# Generate model object:
model <- generate_psychonetrics(model = "lvm",sample = sampleStats,computed = FALSE,
equal = equal,identification=identification,
optimizer = defaultoptimizer(), estimator = estimator, distribution = "Gaussian",
rawts = rawts, types = list(latent = latent, residual = residual),
verbose=verbose)
# Submodel:
latentCov <- latent %in% c("cov","chol")
residCov <- residual %in% c("cov", "chol")
if (latentCov & residCov){
model@submodel <- "sem"
} else if (latentCov & !residCov){
model@submodel <- "rnm"
} else if (!latentCov & residCov){
model@submodel <- "lnm"
} else if (!latentCov & !residCov){
model@submodel <- "lrnm"
}
# Number of groups:
nGroup <- nrow(model@sample@groups)
# Add number of observations:
model@sample@nobs <-
nNode * (nNode+1) / 2 * nGroup + # Covariances per group
nNode * nGroup # Means per group
# Stop if lambda is missing or a character:
if (missing(lambda)){
stop("'lambda' may not be missing")
}
if (is.character(lambda)){
stop("'lambda' may not be a string")
}
# Number of latents:
nLatent <- ncol(lambda)
# If latents is not provided, make it:
if (missing(latents)){
latents <- paste0("Eta_",seq_len(nLatent))
}
if (length(latents) != nLatent){
stop("Length of 'latents' is not equal to number of latent variables in model.")
}
# Model matrices:
modMatrices <- list()
# Fix nu
modMatrices$nu <- matrixsetup_mu(nu,nNode = nNode,nGroup = nGroup,labels = sampleStats@variables$label,equal = "nu" %in% equal,
expmeans = model@sample@means, sampletable = sampleStats, name = "nu")
# Fix nu_eta
modMatrices$nu_eta <- matrixsetup_mu(nu_eta,nNode = nLatent,nGroup = nGroup,labels = latents,equal = "nu_eta" %in% equal,
expmeans = lapply(seq_len(nGroup),function(x)rep(0,nLatent)), sampletable = sampleStats, name = "nu_eta")
# Setup lambda:
modMatrices$lambda <- matrixsetup_lambda(lambda, expcov=model@sample@covs, nGroup = nGroup,
observednames = sampleStats@variables$label, latentnames = latents,
sampletable = sampleStats, identification = identification, simple = simplelambdastart)
# Setup beta:
modMatrices$beta <- matrixsetup_beta(beta, nNode = nLatent, nGroup = nGroup, labels = latents, sampletable = sampleStats)
# Compute the expected latent and residual cov matrices:
expLatSigma <- lapply(1:nGroup,function(x)matrix(0,nLatent,nLatent))
expResidSigma <- lapply(1:nGroup,function(x)matrix(0,nNode,nNode))
# For each group:
for (g in 1:nGroup){
expResidSigma[[g]] <- modMatrices$lambda$sigma_epsilon_start[,,g]
expLatSigma[[g]] <- modMatrices$lambda$sigma_zeta_start[,,g]
# # Current cov estimate:
# curcov <- as.matrix(sampleStats@covs[[g]])
#
#
# # Cur loadings:
# curLambda <- modMatrices$lambda$start[,,g]
# if (!is.matrix(curLambda)){
# curLambda <- as.matrix(curLambda)
# }
#
# ### Run factor start ###:
# # facStart <- factorstart(curcov, modMatrices$lambda[[1]][,,1])
#
#
# ###
#
#
# # Residual variances, let's start by putting the vars on 1/4 times the observed variances:
# Theta <- diag(diag(curcov)/4)
# # Theta <- modMatrices$lambda$thetaStart[,,g]
# fa <- factanal(covmat = curcov, factors = nLatent, covar = TRUE)
#
# fa <- fa(r = curcov, nfactors = nLatent, rotate = "promax", covar = TRUE)
#
# Theta <- diag(fa$uniquenesses)
#
# # Check if this is positive definite:
# ev <- eigen(curcov - Theta)$values
#
# # Shrink until it is positive definite:
# loop <- 0
# repeat{
# ev <- eigen(curcov - Theta)$values
# if (loop == 100){
# # give up...
#
# Theta <- diag(nrow(Theta))
# break
# }
# if (all(ev>0)){
# break
# }
# Theta <- Theta * 0.9
# loop <- loop + 1
# }
#
# # Expected residual sigma:
# expResidSigma[[g]] <- Theta
#
# # This means that the factor-part is expected to be:
# factorPart <- curcov - Theta
#
# # Let's take a pseudoinverse:
# inv <- corpcor::pseudoinverse(kronecker(curLambda,curLambda))
#
# # And obtain psi estimate:
# expLatSigma[[g]] <- matrix(inv %*% as.vector(factorPart),nLatent,nLatent)
}
# Latent varcov:
if (latent == "cov"){
modMatrices$sigma_zeta <- matrixsetup_sigma(sigma_zeta,
name = "sigma_zeta",
expcov=expLatSigma,
nNode = nLatent,
nGroup = nGroup,
labels = latents,
equal = "sigma_zeta" %in% equal, sampletable = sampleStats,
beta = modMatrices$beta[[1]])
} else if (latent == "chol"){
modMatrices$lowertri_zeta <- matrixsetup_lowertri(lowertri_zeta,
name = "lowertri_zeta",
expcov=expLatSigma,
nNode = nLatent,
nGroup = nGroup,
labels = latents,
equal = "lowertri_zeta" %in% equal, sampletable = sampleStats,
beta = modMatrices$beta[[1]])
} else if (latent == "ggm"){
# Add omega matrix:
modMatrices$omega_zeta <- matrixsetup_omega(omega_zeta,
name = "omega_zeta",
expcov=expLatSigma,
nNode = nLatent,
nGroup = nGroup,
labels = latents,
equal = "lowertri_zeta" %in% equal, sampletable = sampleStats,
beta = modMatrices$beta[[1]])
# Add delta matrix:
modMatrices$delta_zeta <- matrixsetup_delta(delta_zeta,
name = "delta_zeta",
expcov=expLatSigma,
nNode = nLatent,
nGroup = nGroup,
labels = latents,
equal = "delta_zeta" %in% equal, sampletable = sampleStats,
omegaStart = modMatrices$omega_zeta$start)
} else if (latent == "prec"){
# Add omega matrix:
modMatrices$kappa_zeta <- matrixsetup_kappa(kappa_zeta,
name = "kappa_zeta",
expcov=expLatSigma,
nNode = nLatent,
nGroup = nGroup,
labels = latents,
equal = "kappa_zeta" %in% equal, sampletable = sampleStats,
beta = modMatrices$beta[[1]])
}
### Residual varcov ###
if (residual == "cov"){
modMatrices$sigma_epsilon <- matrixsetup_sigma(sigma_epsilon,
name = "sigma_epsilon",
expcov=expResidSigma,
nNode = nNode,
nGroup = nGroup,
labels = sampleStats@variables$label,
equal = "sigma_epsilon" %in% equal, sampletable = sampleStats)
} else if (residual == "chol"){
modMatrices$lowertri_epsilon <- matrixsetup_lowertri(lowertri_epsilon,
name = "lowertri_epsilon",
expcov=expResidSigma,
nNode = nNode,
nGroup = nGroup,
labels = sampleStats@variables$label,
equal = "lowertri_epsilon" %in% equal, sampletable = sampleStats)
} else if (residual == "ggm"){
# Add omega matrix:
modMatrices$omega_epsilon <- matrixsetup_omega(omega_epsilon,
name = "omega_epsilon",
expcov=expResidSigma,
nNode = nNode,
nGroup = nGroup,
labels = sampleStats@variables$label,
equal = "omega_epsilon" %in% equal, sampletable = sampleStats)
# Add delta matrix:
modMatrices$delta_epsilon <- matrixsetup_delta(delta_epsilon,
name = "delta_epsilon",
expcov=expResidSigma,
nNode = nNode,
nGroup = nGroup,
labels = sampleStats@variables$label,
equal = "delta_epsilon" %in% equal, sampletable = sampleStats,
omegaStart = modMatrices$omega_epsilon$start)
} else if (residual == "prec"){
# Add omega matrix:
modMatrices$kappa_epsilon <- matrixsetup_kappa(kappa_epsilon,
name = "kappa_epsilon",
expcov=expResidSigma,
nNode = nNode,
nGroup = nGroup,
labels = sampleStats@variables$label,
equal = "kappa_epsilon" %in% equal, sampletable = sampleStats)
}
# Generate the full parameter table:
pars <- do.call(generateAllParameterTables, modMatrices)
# Store in model:
model@parameters <- pars$partable
model@matrices <- pars$mattable
# Extra matrices:
model@extramatrices <- list(
D = psychonetrics::duplicationMatrix(nNode), # non-strict duplciation matrix
Deta = psychonetrics::duplicationMatrix(nLatent), # non-strict duplciation matrix
L = psychonetrics::eliminationMatrix(nNode), # Elinimation matrix
L_eta = psychonetrics::eliminationMatrix(nLatent), # Elinimation matrix
Dstar = psychonetrics::duplicationMatrix(nNode,diag = FALSE), # Strict duplicaton matrix
Dstar_eta = psychonetrics::duplicationMatrix(nLatent,diag = FALSE), # Strict duplicaton matrix
In = as(diag(nNode),"dMatrix"), # Identity of dim n
Inlatent = as(diag(nLatent),"dMatrix"),
C = as(lavaan::lav_matrix_commutation(nNode, nLatent),"dMatrix"),
Cbeta = as(lavaan::lav_matrix_commutation(nLatent, nLatent),"dMatrix"),
C_chol = as(lavaan::lav_matrix_commutation(nNode, nNode),"dMatrix"),
A = psychonetrics::diagonalizationMatrix(nNode),
Aeta = psychonetrics::diagonalizationMatrix(nLatent)
)
# Form the model matrices
model@modelmatrices <- formModelMatrices(model)
### Baseline model ###
if (baseline_saturated){
# Form baseline model:
model@baseline_saturated$baseline <- varcov(data,
type = "chol",
lowertri = "diag",
vars = vars,
groups = groups,
covs = covs,
means = means,
nobs = nobs,
missing = missing,
equal = equal,
estimator = estimator,
baseline_saturated = FALSE, sampleStats = sampleStats)
# Add model:
# model@baseline_saturated$baseline@fitfunctions$extramatrices$M <- Mmatrix(model@baseline_saturated$baseline@parameters)
### Saturated model ###
model@baseline_saturated$saturated <- varcov(data,
type = "chol",
lowertri = "full",
vars = vars,
groups = groups,
covs = covs,
means = means,
nobs = nobs,
missing = missing,
equal = equal,
estimator = estimator,
baseline_saturated = FALSE, sampleStats = sampleStats)
# if not FIML, Treat as computed:
if (estimator != "FIML"){
model@baseline_saturated$saturated@computed <- TRUE
# FIXME: TODO
model@baseline_saturated$saturated@objective <- psychonetrics_fitfunction(parVector(model@baseline_saturated$saturated),model@baseline_saturated$saturated)
}
}
# Identify model:
if (identify){
model <- identify(model)
}
if (missing(optimizer)){
model <- setoptimizer(model, "default")
} else {
model <- setoptimizer(model, optimizer)
}
# Return model:
return(model)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.