Nothing
# Latent network model creator
varcov <- function(
data, # Dataset
type = c("cov","chol","prec","ggm","cor"),
sigma = "full",
kappa = "full", # Precision
# rho = "full", # Correlations
omega = "full", # Partial correlations
lowertri = "full", # Cholesky
delta = "diag", # Used for both ggm and pcor
rho = "full", # Used for cor
SD = "full", # Used for cor
mu,
tau,
vars, # character indicating the variables Extracted if missing from data - group variable
ordered = character(0), # character indicating the variables that are ordinal
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 = "default",
optimizer,
storedata = FALSE,
WLS.W,
sampleStats, # Leave to missing
meanstructure, # Defaults to TRUE if data is used or means is used, FALSE otherwie
corinput, # Defaults to TRUE if the input is detected to consist of correlation matrix/matrices, FALSE otherwise
verbose = FALSE,
covtype = c("choose","ML","UB"),
standardize = c("none","z","quantile"),
fullFIML=FALSE,
bootstrap = FALSE,
boot_sub,
boot_resample
){
rawts = FALSE
if (rawts){
warning("'rawts' is only included for testing purposes. Please do not use!")
}
# Reset ordered if needed:
if (identical(ordered, FALSE)){
ordered <- character(0)
}
if (estimator == "default"){
if (length(ordered) > 0){
estimator <- "WLS"
} else {
estimator <- "ML"
}
}
# Check WLS for ordinal:
if (length(ordered) > 0 & !estimator %in% c("WLS","DWLS","ULS")){
stop("Ordinal data is only supported for WLS, DWLS and ULS estimators.")
}
# Check corinput FIXME: may be mixture
if (length(ordered) > 0 & missing(corinput)){
corinput <- TRUE
}
if (length(ordered) > 0 && (!missing(corinput) && !corinput)){
stop("corinput must be TRUE for ordinal data.")
}
# Disable meanstructure (FIXME: Allow for both ordinal and continuous):
if (length(ordered) > 0){
meanstructure <- FALSE
}
# Type:
type <- match.arg(type)
# Set meanstructure:
if (missing(meanstructure)){
meanstructure <- (!missing(data) || !missing(means))
}
# Check FIML:
if (!missing(data) && !meanstructure && estimator == "FIML"){
stop("meanstructure = FALSE is not yet supported for 'FIML' estimator")
}
# Obtain sample stats:
if (missing(sampleStats)){
# WLS weights:
if (missing(WLS.W)){
WLS.W <- ifelse(!estimator %in% c("WLS","ULS","DWLS"), "none",
switch(estimator,
"WLS" = "full",
"ULS" = "identity",
"DWLS" = "diag"
))
}
sampleStats <- samplestats(data = data,
vars = vars,
ordered=ordered,
groups = groups,
covs = covs,
means = means,
nobs = nobs,
missing = ifelse(estimator == "FIML","pairwise",missing),
rawts = rawts,
fimldata = estimator == "FIML",
storedata = storedata,
weightsmatrix = WLS.W,
meanstructure = meanstructure,
corinput = corinput,
covtype=covtype,
verbose=verbose,
standardize=standardize,
fullFIML=fullFIML,
bootstrap=bootstrap,
boot_sub = boot_sub,
boot_resample = boot_resample)
}
# Overwrite corinput:
corinput <- sampleStats@corinput
# Check some things:
nNode <- nrow(sampleStats@variables)
# Generate model object:
model <- generate_psychonetrics(model = "varcov",sample = sampleStats,computed = FALSE,
equal = equal,
optimizer = defaultoptimizer(), estimator = estimator, distribution = "Gaussian",
rawts = rawts, types = list(y = type),
submodel = type, meanstructure = meanstructure, verbose=verbose)
# Number of groups:
nGroup <- nrow(model@sample@groups)
# Number of means and thresholds:
nMeans <- sum(sapply(model@sample@means,function(x)sum(!is.na(x))))
if (length(ordered) > 0){
nThresh <- sum(sapply(model@sample@thresholds,function(x)sum(sapply(x,length))))
} else {
nThresh <- 0
}
# Add number of observations:
model@sample@nobs <-
nNode * (nNode-1) / 2 * nGroup + # Covariances per group
(!corinput) * nNode * nGroup + # Variances (ignored if correlation matrix is input)
meanstructure * nMeans + # Means per group
nThresh
# Model matrices:
modMatrices <- list()
# # Fix mu
# modMatrices$mu <- matrixsetup_mu(mu,nNode = nNode,nGroup = nGroup,labels = sampleStats@variables$label,equal = "mu" %in% equal,
# expmeans = model@sample@means, sampletable = sampleStats, meanstructure = meanstructure)
# Alternative, without mean structure just ignore mu:
if (meanstructure){
# Fix mu
modMatrices$mu <- matrixsetup_mu(mu,nNode = nNode,nGroup = nGroup,labels = sampleStats@variables$label,equal = "mu" %in% equal,
expmeans = model@sample@means, sampletable = sampleStats, meanstructure = meanstructure)
}
# Thresholds:
if (length(ordered) > 0){
# Ideal setup for tau is a matrix, with NA for the missing elements.
modMatrices$tau <- matrixsetup_tau(tau, nNode = nNode,nGroup = nGroup,labels = sampleStats@variables$label,
equal = "tau" %in% equal, sampleThresholds = model@sample@thresholds, sampletable = sampleStats)
}
# fixMu(mu,nGroup,nNode,"mu" %in% equal)
# Fix sigma
if (type == "cov"){
if (corinput){
stop("Correlation matrix input is not supported for type = 'cov'. Use type = 'cor' or set corinput = FALSE")
}
modMatrices$sigma <- matrixsetup_sigma(sigma,
expcov=model@sample@covs,
nNode = nNode,
nGroup = nGroup,
labels = sampleStats@variables$label,
equal = "sigma" %in% equal, sampletable = sampleStats)
} else if (type == "chol"){
if (corinput){
stop("Correlation matrix input is not supported for type = 'chol'.")
}
modMatrices$lowertri <- matrixsetup_lowertri(lowertri,
expcov=model@sample@covs,
nNode = nNode,
nGroup = nGroup,
labels = sampleStats@variables$label,
equal = "lowertri" %in% equal, sampletable = sampleStats)
} else if (type == "ggm"){
# Add omega matrix:
modMatrices$omega <- matrixsetup_omega(omega,
expcov=model@sample@covs,
nNode = nNode,
nGroup = nGroup,
labels = sampleStats@variables$label,
equal = "omega" %in% equal, sampletable = sampleStats)
if (!corinput){
# Add delta matrix (ingored if corinput == TRUE):
modMatrices$delta <- matrixsetup_delta(delta,
expcov=model@sample@covs,
nNode = nNode,
nGroup = nGroup,
labels = sampleStats@variables$label,
equal = "delta" %in% equal, sampletable = sampleStats,
omegaStart = modMatrices$omega$start)
}
} else if (type == "prec"){
if (corinput){
stop("Correlation matrix input is not supported for type = 'prec'. Use type = 'ggm' or set corinput = FALSE")
}
# Add omega matrix:
modMatrices$kappa <- matrixsetup_kappa(kappa,
expcov=model@sample@covs,
nNode = nNode,
nGroup = nGroup,
labels = sampleStats@variables$label,
equal = "kappa" %in% equal, sampletable = sampleStats)
} else if (type == "cor"){
# Add rho matrix:
modMatrices$rho <- matrixsetup_rho(rho,
expcov=model@sample@covs,
nNode = nNode,
nGroup = nGroup,
labels = sampleStats@variables$label,
equal = "rho" %in% equal, sampletable = sampleStats)
if (!corinput){
# Add SD matrix (ignored if corinput == TRUE):
modMatrices$SD <- matrixsetup_SD(SD,
expcov=model@sample@covs,
nNode = nNode,
nGroup = nGroup,
labels = sampleStats@variables$label,
equal = "SD" %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
# if (type == "cov" || type == "prec"){
# model@extramatrices <- list(
# D = psychonetrics::duplicationMatrix(nNode), # non-strict duplciation matrix
# L = psychonetrics::eliminationMatrix(nNode), # Elinimation matrix
# In = as(diag(nNode),"dMatrix") # Identity of dim n
# )
# } else if (type == "chol"){
# model@extramatrices <- list(
# D = psychonetrics::duplicationMatrix(nNode), # non-strict duplciation matrix
# L = psychonetrics::eliminationMatrix(nNode), # Elinimation matrix
# In = as(diag(nNode),"dMatrix"), # Identity of dim n
# # C = as(lavaan::lav_matrix_commutation(nNode,nNode),"sparseMatrix")
# C = as(lavaan::lav_matrix_commutation(nNode,nNode),"dMatrix")
# )
# } else if (type == "ggm" || type == "cor"){
# model@extramatrices <- list(
# D = psychonetrics::duplicationMatrix(nNode), # non-strict duplciation matrix
# L = psychonetrics::eliminationMatrix(nNode), # Elinimation matrix
# Dstar = psychonetrics::duplicationMatrix(nNode,diag = FALSE), # Strict duplicaton matrix
# In = as(diag(nNode),"dMatrix"), # Identity of dim n
# A = psychonetrics::diagonalizationMatrix(nNode)
# )
# }
model@extramatrices <- list(
D = psychonetrics::duplicationMatrix(nNode), # non-strict duplciation matrix
L = psychonetrics::eliminationMatrix(nNode), # Elinimation matrix
In = as(diag(nNode),"dMatrix"), # Identity of dim n
C = as(lavaan::lav_matrix_commutation(nNode,nNode),"dMatrix"),
Dstar = psychonetrics::duplicationMatrix(nNode,diag = FALSE), # Strict duplicaton matrix
A = psychonetrics::diagonalizationMatrix(nNode)
)
# Form the model matrices
model@modelmatrices <- formModelMatrices(model)
### Baseline model ###
if (baseline_saturated){
# Form baseline model:
if ("omega" %in% equal | "sigma" %in% equal | "kappa" %in% equal){
equal <- c(equal,"lowertri")
}
if (!corinput){
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,
meanstructure=meanstructure,
corinput = corinput,
ordered = ordered,
baseline_saturated = FALSE,sampleStats=sampleStats)
} else {
model@baseline_saturated$baseline <- varcov(data,
type = "cor",
rho = "zero",
vars = vars,
groups = groups,
covs = covs,
means = means,
nobs = nobs,
missing = missing,
equal = equal,
estimator = estimator,
meanstructure=meanstructure,
corinput = corinput,
ordered = ordered,
baseline_saturated = FALSE,sampleStats=sampleStats)
}
# Add model:
# model@baseline_saturated$baseline@fitfunctions$extramatrices$M <- Mmatrix(model@baseline_saturated$baseline@parameters)
### Saturated model ###
if (!corinput){
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,
meanstructure=meanstructure,
corinput = corinput,
ordered = ordered,
baseline_saturated = FALSE,sampleStats=sampleStats)
} else {
model@baseline_saturated$saturated <- varcov(data,
type = "cor",
lowertri = "full",
vars = vars,
groups = groups,
covs = covs,
means = means,
nobs = nobs,
missing = missing,
equal = equal,
estimator = estimator,
meanstructure=meanstructure,
corinput = corinput,
ordered = ordered,
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)
}
}
if (missing(optimizer)){
model <- setoptimizer(model, "default")
} else {
model <- setoptimizer(model, optimizer)
}
# Return model:
return(model)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.