R/h_modelsearch.R In psychonetrics: Structural Equation Modeling and Confirmatory Network Analysis

Documented in modelsearch

```# ggmModSelect inspired algorithm:
modelsearch <- function(x,
criterion = "bic", # Stop when criterion is no longer improved. Can also be none to ignore
matrices, # Matrices to search
prunealpha = 0.01, # Minimum p-value for edges tested to be removed
addalpha = 0.01, # Maximum p-value for edges tested to be added
verbose,
...
){
mi <- "mi"
pmi <- "pmi"

if (missing(verbose)){
verbose <- x@verbose
}

# FIXME: If number of groups > 1, stop:
if (nrow(x@sample@groups) > 1) stop("'modelsearch' is only implemented for single group models at the moment.")

# If not run, run model:
if (!x@computed){
x <- x %>% runmodel(..., verbose = verbose)
}

x <- addLog(x, "Starting modelsearch algorithm...")

# Matrices:
if (missing(matrices)){
if (x@model == "varcov"){
if (x@submodel == "ggm"){
matrices <- "omega"
} else if (x@submodel == "prec"){
matrices <- "kappa"
} else if (x@submodel == "chol"){
matrices <- "lowertri"
} else if (x@submodel == "cov"){
matrices <- "sigma"
}
} else if (x@model == "meta_varcov"){

matrices <- character(0)
if (x@types\$y == "ggm"){
matrices <- c(matrices,"omega_y")
} else if (x@types\$y == "prec"){
matrices <- c(matrices,"kappa_y")
}

# Random effects:
if (x@types\$randomEffects == "ggm"){
matrices <- c(matrices,"omega_randomEffects")
} else if (x@types\$randomEffects == "prec"){
matrices <- c(matrices,"kappa_randomEffects")
}

} else if (x@model == "lvm"){

# Only add GGM structures to search:
matrices <- character(0)
if (x@types\$latent == "ggm"){
matrices <- c(matrices,"omega_zeta")
} else if (x@types\$latent == "prec"){
matrices <- c(matrices,"kappa_zeta")
}

if (x@types\$residual == "ggm"){
matrices <- c(matrices,"omega_epsilon")
} else if (x@types\$residual == "prec"){
matrices <- c(matrices,"kappa_epsilon")
}

} else if (x@model == "lnm"){
matrices <- "omega_eta"
} else if (x@model == "rnm"){
matrices <- "omega_epsilon"
}  else if (x@model == "gvar"){
matrices <- c("beta","omega_zeta")
}  else if (x@model == "var1"){
matrices <- c("beta")
if (x@types\$zeta == "prec"){
matrices <- c(matrices,"kappa_zeta")
} else if (x@types\$zeta == "ggm"){
matrices <- c(matrices,"omega_zeta")
}

} else if (x@model == "panelvar1"){
matrices <- c("beta")
if (x@types\$contemporaneous == "prec"){
matrices <- c(matrices,"kappa_zeta")
} else if (x@types\$contemporaneous == "ggm"){
matrices <- c(matrices,"omega_zeta")
}

if (x@types\$between == "prec"){
matrices <- c(matrices,"kappa_mu")
} else if (x@types\$between == "ggm"){
matrices <- c(matrices,"omega_mu")
}

} else if (x@model %in% c("ml_lvm","dlvm1")){
if (x@model == "dlvm1") {
matrices <- c("beta")
} else {
matrices <- character(0)
}
if (x@types\$within_latent == "prec"){
matrices <- c(matrices,"kappa_zeta_within")
} else if (x@types\$within_latent == "ggm"){
matrices <- c(matrices,"omega_zeta_within")
}

if (x@types\$within_residual == "prec"){
matrices <- c(matrices,"kappa_epsilon_within")
} else if (x@types\$within_residual == "ggm"){
matrices <- c(matrices,"omega_epsilon_within")
}

if (x@types\$between_latent == "prec"){
matrices <- c(matrices,"kappa_zeta_between")
} else if (x@types\$between_latent == "ggm"){
matrices <- c(matrices,"omega_zeta_between")
}

if (x@types\$between_residual == "prec"){
matrices <- c(matrices,"kappa_epsilon_between")
} else if (x@types\$between_residual == "ggm"){
matrices <- c(matrices,"omega_epsilon_between")
}

} else if (x@model == "tsdlvm1"){
matrices <- c("beta")
if (x@types\$zeta == "prec"){
matrices <- c(matrices,"kappa_zeta")
} else if (x@types\$zeta == "ggm"){
matrices <- c(matrices,"omega_zeta")
}

if (x@types\$epsilon == "prec"){
matrices <- c(matrices,"kappa_epsilon")
} else if (x@types\$epsilon == "ggm"){
matrices <- c(matrices,"omega_epsilon")
}

}  else if (x@model == "Ising"){
matrices <- c("omega")

}  else stop("No default argument for 'matrices' for current model.")
}

# Check if MIs are added:
if (all(is.na(x@parameters[[mi]]))){
x <- x %>% addMIs(matrices = matrices)
}

# Start outer loop for evaluating ALL competing models!
repeat{
# List all parameters in matrices of interest:
ind <- which(x@parameters\$matrix %in% matrices & !x@parameters\$identified)

nona <- function(x){
x[is.na(x)] <- FALSE
x
}

allParsToConsider <- ind[
!is.na(x@parameters\$fixed[ind]) &

# Non significant edges to remove?
nona(!x@parameters\$fixed[ind]  & x@parameters\$p[ind] > prunealpha) |

]

if (length(allParsToConsider) == 0){
if (verbose){
message("No more models to search at given alpha levels.")
return(x)
}
}

# Set parameters to consider:
parsToConsider <- allParsToConsider

# Start the inner loop for subset of models:
repeat{

# Set progress bar:
if (verbose){
if (length(allParsToConsider) == length(parsToConsider)){
message("Evaluating all models at given alpha levels...")
} else {
message("Evaluating subset of models...")
}

pb <- txtProgressBar(min = 0, max = length(parsToConsider), initial = 0, style = 3)
}

# Empty list of models:
propMods <- vector("list", length(parsToConsider))

# Loop over these parameters:
for (i in seq_along(parsToConsider)){
curpar <- parsToConsider[i]

# Existing parameter?
fixed <- x@parameters\$fixed[curpar]

if (fixed){
propMods[[i]] <- x %>% freepar(matrix = x@parameters\$matrix[curpar], row = x@parameters\$row[curpar], col = x@parameters\$col[curpar],
group = x@parameters\$group_id[curpar],
start = x@parameters\$epc[curpar], verbose = FALSE)

# Run the model:
suppressWarnings(propMods[[i]] <- propMods[[i]] %>% runmodel(verbose = FALSE, addMIs = FALSE, return_improper = TRUE,...))

# Fisher information ok?
# ev <- eigen(propMods[[i]]@information)\$values
proper <- TRUE
if (!is.null(propMods[[i]]@modelmatrices[[1]]\$proper)){
propers <- sapply(propMods[[i]]@modelmatrices,"[[","proper")
proper <- all(propers)
}
# If not, try again with different starts:
# if (any(ev < -sqrt(.Machine\$double.eps))){
if (!sympd_cpp(propMods[[i]]@information) || !proper){
try <- 1
repeat{
if (try == 1){
# First try, EPC:
propMods[[i]] <- x %>% freepar(matrix = x@parameters\$matrix[curpar], row = x@parameters\$row[curpar], col = x@parameters\$col[curpar],
group = x@parameters\$group_id[curpar],
start = x@parameters\$epc[curpar], verbose = FALSE)

} else if (try == 2){
# Second try, 0.5 * EPC:
propMods[[i]] <- x %>% freepar(matrix = x@parameters\$matrix[curpar], row = x@parameters\$row[curpar], col = x@parameters\$col[curpar],
group = x@parameters\$group_id[curpar],
start = 0.5 * x@parameters\$epc[curpar], verbose = FALSE)

} else if (try == 3){
# Third try, 0.0001 * sign(EPC):
propMods[[i]] <- x %>% freepar(matrix = x@parameters\$matrix[curpar], row = x@parameters\$row[curpar], col = x@parameters\$col[curpar],
group = x@parameters\$group_id[curpar],
start = 0.0001 * sign(x@parameters\$epc[curpar]), verbose = FALSE)

} else if (try == 4){
# Final try, emergency start:
propMods[[i]] <- x %>% freepar(matrix = x@parameters\$matrix[curpar], row = x@parameters\$row[curpar], col = x@parameters\$col[curpar],
group = x@parameters\$group_id[curpar],
start = x@parameters\$epc[curpar], verbose = FALSE) %>% emergencystart

} else if (try > 4){

# Give up:
propMods[[i]] <- x
propMods[[i]]@fitmeasures[[criterion]] <- Inf
break

}

# Run the model:
suppressWarnings(propMods[[i]] <- propMods[[i]] %>% runmodel(verbose = FALSE, addMIs = FALSE,return_improper = TRUE, ...))

# Fisher information ok?
# ev <- eigen(propMods[[i]]@information)\$values

# Check for properness:
proper <- TRUE
if (!is.null(propMods[[i]]@modelmatrices[[1]]\$proper)){
propers <- sapply(propMods[[i]]@modelmatrices,"[[","proper")
proper <- all(propers)
}

# if (all(ev > -sqrt(.Machine\$double.eps))){
if (sympd_cpp(propMods[[i]]@information) && proper){
break
} else {
try <- try + 1
}

}

}

} else {

# Try to remove (prune):
propMods[[i]] <- x %>% fixpar(matrix = x@parameters\$matrix[curpar], row = x@parameters\$row[curpar], col = x@parameters\$col[curpar],
group = x@parameters\$group_id[curpar], verbose = FALSE)

# Run the model:
propMods[[i]] <- propMods[[i]] %>% runmodel(verbose = FALSE, addMIs = FALSE, return_improper = TRUE, ...)

# Fisher information ok?
# ev <- eigen(propMods[[i]]@information)\$values

# If not, try again with different starts:
# if (any(ev < -sqrt(.Machine\$double.eps))){

proper <- TRUE
if (!is.null(propMods[[i]]@modelmatrices[[1]]\$proper)){
propers <- sapply(propMods[[i]]@modelmatrices,"[[","proper")
proper <- all(propers)
}

if (!sympd_cpp(propMods[[i]]@information) || !proper){
try <- 1
repeat{
if (try == 1){
# First try, reduce all parameters in matrix a bit:
propMods[[i]] <- x %>% fixpar(matrix = x@parameters\$matrix[curpar], row = x@parameters\$row[curpar], col = x@parameters\$col[curpar],
group = x@parameters\$group_id[curpar], verbose = FALSE)

ind <- propMods[[i]]@parameters\$matrix %in% propMods[[i]]@parameters\$matrices & !propMods[[i]]@parameters\$fixed & !propMods[[i]]@parameters\$identified &
propMods[[i]]@parameters\$row != propMods[[i]]@parameters\$col

propMods[[i]]@parameters\$est[ind] <- 0.5 * propMods[[i]]@parameters\$est[ind]

} else if (try == 2){
# Second try, emergency start:
propMods[[i]] <- x %>% fixpar(matrix = x@parameters\$matrix[curpar], row = x@parameters\$row[curpar], col = x@parameters\$col[curpar],
group = x@parameters\$group_id[curpar], verbose = FALSE) %>% emergencystart

} else if (try > 2){

# Give up:
propMods[[i]] <- x
propMods[[i]]@fitmeasures[[criterion]] <- Inf
break

}

# Run the model:
propMods[[i]] <- propMods[[i]] %>% runmodel(verbose = FALSE, addMIs = FALSE, return_improper = TRUE, ...)

proper <- TRUE
if (!is.null(propMods[[i]]@modelmatrices[[1]]\$proper)){
propers <- sapply(propMods[[i]]@modelmatrices,"[[","proper")
proper <- all(propers)
}
# Fisher information ok?
# ev <- eigen(propMods[[i]]@information)\$values

# if (all(ev > -sqrt(.Machine\$double.eps))){

if (sympd_cpp(propMods[[i]]@information) && proper){
break
} else {
try <- try + 1
}

}

}

}

# Update progress:
if (verbose){
setTxtProgressBar(pb, i)
}

}

# Update progress:
if (verbose){
close(pb)
}

# Current criterion:
curCrit <- x@fitmeasures[[criterion]]

# Criterions of proposal models:
propCrits <- sapply(propMods, function(xx) xx@fitmeasures[[criterion]] )

# which improve?
whichImprove <- which(propCrits < curCrit)

# if none, break:
if (length(whichImprove) == 0){
break # FIXME
}

# Which is the best?
whichBest <- which(propCrits == min(propCrits))
if (length(whichBest) > 1){
warning("Multiple equivalent models found, selecting one at random.")
whichBest <- sample(whichBest,1)
}

# Say something:
if (verbose){
if (propMods[[whichBest]]@fitmeasures\$df < x@fitmeasures\$df){
} else {
message("Fixing one parameter...")
}
}

# Update the model:
x <- propMods[[whichBest]] %>% addMIs(verbose = FALSE)

# Update list of candidates:
parsToConsider <- parsToConsider[whichImprove[whichImprove != whichBest]]

# if none, break:
if (length(parsToConsider) == 0){
break # FIXME
}

}

# Did we test all models in the last run? If so, we are done (jaj)!
if (length(allParsToConsider) == length(parsToConsider)){
if (verbose) message("No more model found to improve fit. Returning optimal model.")
break
}

}

# Return the final model:
return(x)
}
```

Try the psychonetrics package in your browser

Any scripts or data that you put into this service are public.

psychonetrics documentation built on June 22, 2024, 10:29 a.m.