diff --git a/.github/workflows/test_package_ubuntu.yml b/.github/workflows/test_package_ubuntu.yml index 464c259..0450eca 100644 --- a/.github/workflows/test_package_ubuntu.yml +++ b/.github/workflows/test_package_ubuntu.yml @@ -30,11 +30,11 @@ jobs: # install two helper packages we need sudo apt install --no-install-recommends software-properties-common dirmngr # add the signing key (by Michael Rutter) for these repos - # To verify key, run gpg --show-keys /etc/apt/trusted.gpg.d/cran_ubuntu_key.asc + # To verify key, run gpg --show-keys /etc/apt/trusted.gpg.d/cran_ubuntu_key.asc # Fingerprint: E298A3A825C0D65DFD57CBB651716619E084DAB9 wget -qO- https://cloud.r-project.org/bin/linux/ubuntu/marutter_pubkey.asc | sudo tee -a /etc/apt/trusted.gpg.d/cran_ubuntu_key.asc # add the R 4.0 repo from CRAN -- adjust 'focal' to 'groovy' or 'bionic' as needed - sudo add-apt-repository "deb https://cloud.r-project.org/bin/linux/ubuntu $(lsb_release -cs)-cran40/" + sudo add-apt-repository "deb https://cloud.r-project.org/bin/linux/ubuntu $(lsb_release -cs)-cran40/" - name: install compiled packages run: | sudo add-apt-repository ppa:c2d4u.team/c2d4u4.0+ @@ -46,6 +46,7 @@ jobs: sudo apt install --no-install-recommends r-cran-rcpp sudo apt install --no-install-recommends r-cran-testthat sudo apt install --no-install-recommends r-cran-rcmdcheck + sudo apt install --no-install-recommends r-cran-dplyr - name: Check run: rcmdcheck::rcmdcheck(path = ".", args = "--no-manual", error_on = "warning") shell: Rscript {0} diff --git a/DESCRIPTION b/DESCRIPTION index 651c1e9..f0063ed 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: mxsem Type: Package Title: Specify 'OpenMx' Models with a 'lavaan'-Style Syntax -Version: 0.0.5 +Version: 0.0.6 Date: 2023-07-12 Authors@R: c(person(given = "Jannik H.", family = "Orzek", role = c("aut", "cre", "cph"), @@ -17,7 +17,9 @@ Depends: OpenMx Imports: Rcpp (>= 1.0.10), stats, - methods + methods, + dplyr, + utils LinkingTo: Rcpp RoxygenNote: 7.2.3 Encoding: UTF-8 diff --git a/NAMESPACE b/NAMESPACE index 64f5bbc..fbe2b62 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,10 +1,15 @@ # Generated by roxygen2: do not edit by hand +S3method(print,multi_group_parameters) +export(get_groups) +export(get_individual_algebra_results) export(mxsem) +export(mxsem_group_by) export(parameters) export(set_starting_values) export(simulate_latent_growth_curve) export(simulate_moderated_nonlinear_factor_analysis) +export(summarize_multi_group_model) export(unicode_directed) export(unicode_undirected) import(OpenMx) @@ -13,4 +18,6 @@ importFrom(methods,is) importFrom(stats,cov) importFrom(stats,rnorm) importFrom(stats,runif) +importFrom(utils,setTxtProgressBar) +importFrom(utils,txtProgressBar) useDynLib(mxsem, .registration = TRUE) diff --git a/R/get_individual_algebra_results.R b/R/get_individual_algebra_results.R new file mode 100644 index 0000000..dfa946b --- /dev/null +++ b/R/get_individual_algebra_results.R @@ -0,0 +1,99 @@ +#' get_individual_algebra_results +#' +#' evaluates algebras for each subject in the data set. This function is +#' useful if you have algebras with definition variables (e.g., in mnlfa). +#' @param mxModel mxModel with algebras +#' @param algebra_names optional: Only compute individual algebras for a subset +#' of the parameters +#' @param progress_bar should a progress bar be shown? +#' @returns a list of data frames. The list contains data frames for each of the algebras. +#' The data frames contain the individual specific algebra results as well as all +#' definition variables used to predict said algebra +#' @export +#' @importFrom utils txtProgressBar +#' @importFrom utils setTxtProgressBar +#' @examples +#' library(mxsem) +#' +#' set.seed(123) +#' dataset <- simulate_moderated_nonlinear_factor_analysis(N = 50) +#' +#' model <- " +#' xi =~ x1 + x2 + x3 +#' eta =~ y1 + y2 + y3 +#' eta ~ {a := a0 + data.k*a1}*xi +#' " +#' fit <- mxsem(model = model, +#' data = dataset) |> +#' mxTryHard() +#' +#' algebra_results <- get_individual_algebra_results(mxModel = fit, +#' progress_bar = FALSE) +#' +#' # the following plot will only show two data points because there is only +#' # two values for the definition variable k (0 or 1). +#' +#' plot(x = algebra_results[["a"]]$k, +#' y = algebra_results[["a"]]$algebra_result) +get_individual_algebra_results <- function(mxModel, + algebra_names = NULL, + progress_bar = TRUE){ + n_subjects <- mxModel$data$numObs + if(is.null(algebra_names)){ + algebra_names <- names(mxModel$algebras) + }else{ + if(any(! algebra_names %in% names(mxModel$algebras))) + stop("Could not find the following algebra(s) in your model: ", + paste0(algebra_names[! algebra_names %in% names(mxModel$algebras)], collapse = ", "), + ". The algebras are: ", + paste0(names(mxModel$algebras), collapse = ", "), ".") + } + + if(is.null(algebra_names) | (length(algebra_names) == 0)) + stop("Could not find any algebras in your OpenMx model.") + + + algebra_results <- vector("list", length(algebra_names)) + names(algebra_results) <- algebra_names + + if(progress_bar) + pb <- utils::txtProgressBar(min = 0, + max = length(algebra_names) * n_subjects, + initial = 0, + style = 3) + + it <- 0 + + for(algebra_name in algebra_names){ + + # find definition variables used in this algebra + algebra_elements <- extract_algebra_elements(mxAlgebra_formula = mxModel$algebras[[algebra_name]]$formula) + definition_variables <- algebra_elements[grepl("^data\\.", x = algebra_elements)] |> + gsub(pattern = "data\\.", replacement = "", x = _) + + algebra_result <- data.frame(person = 1:n_subjects, + mxModel$data$observed[,definition_variables, drop = FALSE], + algebra_result = NA) + + for(i in 1:n_subjects){ + it <- it + 1 + if(progress_bar) + utils::setTxtProgressBar(pb = pb, + value = it) + + algebra_result_i <- mxEvalByName(name = algebra_name, + model = mxModel, + compute = TRUE, + defvar.row = i) + if((nrow(algebra_result_i) != 1) | + (ncol(algebra_result_i) != 1)) + stop("This function cannot handle algebras with non-scalar outcomes.") + + algebra_result$algebra_result[i] <- algebra_result_i[1,1] + } + + algebra_results[[algebra_name]] <- algebra_result + } + + return(algebra_results) +} diff --git a/R/mxsem.R b/R/mxsem.R index b1d75c4..c6ee73f 100644 --- a/R/mxsem.R +++ b/R/mxsem.R @@ -314,7 +314,7 @@ mxsem <- function(model, mx_data <- data }else{ if(!add_intercepts){ - mx_data <- OpenMx::mxData(observed = stats::cov(data), + mx_data <- OpenMx::mxData(observed = stats::cov(as.matrix(data[,parameter_table$variables$manifests])), type = "cov", numObs = nrow(data)) }else{ diff --git a/R/mxsem_group_by.R b/R/mxsem_group_by.R new file mode 100644 index 0000000..41b0b59 --- /dev/null +++ b/R/mxsem_group_by.R @@ -0,0 +1,252 @@ +#' mxsem_group_by +#' +#' creates a multi-group model from an OpenMx model. +#' +#' mxsem_group_by creates a multi-group model by splitting the data found +#' in an mxModel object using dplyr's group_by function. The general idea is +#' as follows: +#' +#' 1. The function extracts the data from mxModel +#' 2. The data is split using the group_by function of dplyr with the variables +#' in grouping_variables +#' 3. a separate model is set up for each group. All parameters that match +#' those specified in the parameters argument are group specific +#' +#' **Warning**: The multi-group model may differ from **lavaan**! For instance, +#' **lavaan** will automatically set the latent variances for all but the first +#' group free if the loadings are fixed to equality. Such automatic procedures +#' are not yet implemented in **mxsem**. +#' +#' +#' @param mxModel mxModel with the entire data +#' @param grouping_variables Variables used to split the data in groups +#' @param parameters the parameters that should be group specific. By default +#' all parameters are group specific. +#' @param use_grepl if set to TRUE, grepl is used to check which parameters are +#' group-specific. For instance, if parameters = "a" and use_grepl = TRUE, all parameters +#' whose label contains the letter "a" will be group specific. If use_grep = FALSE +#' only the parameter that has the label "a" is group specific. +#' @return mxModel with multiple groups. Use get_groups to extract the groups +#' @export +#' @examples +#' # THE FOLLOWING EXAMPLE IS ADAPTED FROM +#' # https://openmx.ssri.psu.edu/docs/OpenMx/latest/_static/Rdoc/mxModel.html +#' library(mxsem) +#' +#' model <- 'spatial =~ visual + cubes + paper +#' verbal =~ general + paragrap + sentence +#' math =~ numeric + series + arithmet' +#' +#' mg_model <- mxsem(model = model, +#' data = OpenMx::HS.ability.data) |> +#' # we want separate models for all combinations of grades and schools: +#' mxsem_group_by(grouping_variables = "school") |> +#' mxTryHard() +#' +#' # let's summarize the results: +#' summarize_multi_group_model(mg_model) +mxsem_group_by <- function(mxModel, + grouping_variables, + parameters = c(".*"), + use_grepl = TRUE){ + + if(!is(mxModel, "MxRAMModel")) + stop("mxModel must be of class MxRAMModel.") + if(mxModel$data$type != "raw") + stop("The data of mxModel must be of type 'raw'") + + splitted_data <- mxModel$data$observed |> + dplyr::group_by(dplyr::across(dplyr::all_of(grouping_variables))) |> + dplyr::group_split() + + n_groups <- length(splitted_data) + + names(splitted_data) <- paste0("group_", seq_len(n_groups)) + + group_specific_parameters <- c() + parameter_labels <- names(OpenMx::omxGetParameters(mxModel)) + for(p in parameters){ + if(use_grepl){ + group_specific_parameters <- c(group_specific_parameters, + parameter_labels[grepl(p, parameter_labels)]) + }else{ + group_specific_parameters <- c(group_specific_parameters, + parameter_labels[parameter_labels == p] + ) + } + } + common_parameters <- parameter_labels[!parameter_labels %in% group_specific_parameters] + + message(paste0("The following parameters will be the same across groups: ", + paste0(common_parameters, collapse = ", "))) + message(paste0("The following parameters will be group specific: ", + paste0(group_specific_parameters, collapse = ", "))) + + mg_model <- OpenMx::mxModel(name = mxModel$name) + + fitfunction <- c() + + for(gr in seq_len(n_groups)){ + + fitfunction <- c(fitfunction, + paste0(mxModel$name, "_group_", gr, ".fitfunction")) + + current_model <- OpenMx::mxModel( + name = paste0(mxModel$name, "_group_", gr), + mxModel, + OpenMx::mxData(splitted_data[[gr]], type = "raw") + ) + + current_parameter_labels <- names(OpenMx::omxGetParameters(mxModel)) + + for(p in parameters){ + if(use_grepl){ + current_parameter_labels[grepl(p, current_parameter_labels)] <- paste0( + current_parameter_labels[grepl(p, current_parameter_labels)], + "_group_", gr + ) + }else{ + current_parameter_labels[current_parameter_labels == p] <- paste0( + current_parameter_labels[current_parameter_labels == p], + "_group_", gr + ) + } + } + + current_model <- OpenMx::omxSetParameters(model = current_model, + labels = names(OpenMx::omxGetParameters(mxModel)), + values = OpenMx::omxGetParameters(mxModel), + newlabels = current_parameter_labels) + mg_model <- OpenMx::mxModel( + mg_model, + current_model + ) + } + + mg_model <- OpenMx::mxModel( + mg_model, + # define fitfunction as the sum of the fitfunctions of the group models + OpenMx::mxAlgebraFromString(paste0(fitfunction, collapse = " + "), + name = "mg_fitfunction"), + OpenMx::mxFitFunctionAlgebra(algebra = "mg_fitfunction") + ) + + attr(mg_model, which = "groups") <- splitted_data + attr(mg_model, which = "grouping_variables") <- grouping_variables + attr(mg_model, which = "common_parameters") <- common_parameters + attr(mg_model, which = "group_specific_parameters") <- group_specific_parameters + + return(mg_model) +} + +#' get_groups +#' +#' returns a list of groups for a multi group model +#' @param multi_group_model multi group model created with mxsem_group_by +#' @return list with data for each group +#' @export +#' @examples +#' # THE FOLLOWING EXAMPLE IS ADAPTED FROM +#' # https://openmx.ssri.psu.edu/docs/OpenMx/latest/_static/Rdoc/mxModel.html +#' library(mxsem) +#' +#' model <- 'spatial =~ visual + cubes + paper +#' verbal =~ general + paragrap + sentence +#' math =~ numeric + series + arithmet' +#' +#' mg_model <- mxsem(model = model, +#' data = OpenMx::HS.ability.data) |> +#' # we want separate models for all combinations of grades and schools: +#' mxsem_group_by(grouping_variables = "school") |> +#' mxTryHard() +#' +#' # let's summarize the results: +#' summarize_multi_group_model(mg_model) +#' +#' # let's get the groups: +#' get_groups(mg_model) +get_groups <- function(multi_group_model){ + return(attr(multi_group_model, "groups")) +} + + +#' summarize_multi_group_model +#' +#' summarize the results of a multi group model created with mxsem_group_by +#' @param multi_group_model multi group model created with mxsem_group_by +#' @return list with goup specific parameters and common parameters +#' @export +#' @examples +#' # THE FOLLOWING EXAMPLE IS ADAPTED FROM +#' # https://openmx.ssri.psu.edu/docs/OpenMx/latest/_static/Rdoc/mxModel.html +#' library(mxsem) +#' +#' model <- 'spatial =~ visual + cubes + paper +#' verbal =~ general + paragrap + sentence +#' math =~ numeric + series + arithmet' +#' +#' mg_model <- mxsem(model = model, +#' data = OpenMx::HS.ability.data) |> +#' # we want separate models for all combinations of grades and schools: +#' mxsem_group_by(grouping_variables = "school") |> +#' mxTryHard() +#' +#' # let's summarize the results: +#' summarize_multi_group_model(mg_model) +summarize_multi_group_model <- function(multi_group_model){ + + groups <- attr(multi_group_model, which = "groups") + grouping_variables <- attr(multi_group_model, which = "grouping_variables") + common_parameters <- attr(multi_group_model, "common_parameters") + parameters <- OpenMx::omxGetParameters(model = multi_group_model) + summmarized <- summary(multi_group_model) + + parameter_table <- vector("list", length(groups) + 1) + names(parameter_table) <- c("common parameters", + names(groups)) + + parameter_table[["common parameters"]] <- summmarized$parameters[summmarized$parameters$name %in% common_parameters,] + + for(gr in seq_len(length(groups))){ + + group_name <- names(groups)[gr] + parameter_table[[group_name]] <- list( + parameters = summmarized$parameters[grepl(pattern = paste0(group_name, "$"), + x = names(parameters)),], + group_setting = groups[[gr]][,grouping_variables] |> + unique() + ) + } + + class(parameter_table) <- "multi_group_parameters" + return(parameter_table) +} + + +#' print the multi_group_parameters +#' @param x object from summarize_multi_group_model +#' @param ... not used +#' @method print multi_group_parameters +#' @export +print.multi_group_parameters <- function(x, ...){ + console_width <- getOption("width") + cat("\n") + cat("\n") + cat(paste0(rep("-", console_width), collapse = "")) + cat("\n") + cat("Common Parameters:") + cat("\n") + print(x$`common parameters`) + + for(gr in 2:length(x)){ + cat(paste0(rep("-", console_width), collapse = "")) + cat("\n") + cat(paste0(names(x[gr])[1], ":")) + cat("\n") + print(x[[gr]][["group_setting"]]) + cat("\n") + print(x[[gr]][["parameters"]]) + } + cat("\n") + cat(paste0(rep("-", console_width), collapse = "")) +} diff --git a/man/get_groups.Rd b/man/get_groups.Rd new file mode 100644 index 0000000..f45a984 --- /dev/null +++ b/man/get_groups.Rd @@ -0,0 +1,38 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/mxsem_group_by.R +\name{get_groups} +\alias{get_groups} +\title{get_groups} +\usage{ +get_groups(multi_group_model) +} +\arguments{ +\item{multi_group_model}{multi group model created with mxsem_group_by} +} +\value{ +list with data for each group +} +\description{ +returns a list of groups for a multi group model +} +\examples{ +# THE FOLLOWING EXAMPLE IS ADAPTED FROM +# https://openmx.ssri.psu.edu/docs/OpenMx/latest/_static/Rdoc/mxModel.html +library(mxsem) + +model <- 'spatial =~ visual + cubes + paper + verbal =~ general + paragrap + sentence + math =~ numeric + series + arithmet' + +mg_model <- mxsem(model = model, + data = OpenMx::HS.ability.data) |> + # we want separate models for all combinations of grades and schools: + mxsem_group_by(grouping_variables = "school") |> + mxTryHard() + +# let's summarize the results: +summarize_multi_group_model(mg_model) + +# let's get the groups: +get_groups(mg_model) +} diff --git a/man/get_individual_algebra_results.Rd b/man/get_individual_algebra_results.Rd new file mode 100644 index 0000000..a570e64 --- /dev/null +++ b/man/get_individual_algebra_results.Rd @@ -0,0 +1,53 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/get_individual_algebra_results.R +\name{get_individual_algebra_results} +\alias{get_individual_algebra_results} +\title{get_individual_algebra_results} +\usage{ +get_individual_algebra_results( + mxModel, + algebra_names = NULL, + progress_bar = TRUE +) +} +\arguments{ +\item{mxModel}{mxModel with algebras} + +\item{algebra_names}{optional: Only compute individual algebras for a subset +of the parameters} + +\item{progress_bar}{should a progress bar be shown?} +} +\value{ +a list of data frames. The list contains data frames for each of the algebras. +The data frames contain the individual specific algebra results as well as all +definition variables used to predict said algebra +} +\description{ +evaluates algebras for each subject in the data set. This function is +useful if you have algebras with definition variables (e.g., in mnlfa). +} +\examples{ +library(mxsem) + +set.seed(123) +dataset <- simulate_moderated_nonlinear_factor_analysis(N = 50) + +model <- " + xi =~ x1 + x2 + x3 + eta =~ y1 + y2 + y3 + eta ~ {a := a0 + data.k*a1}*xi + " +fit <- mxsem(model = model, + data = dataset) |> + mxTryHard() + +algebra_results <- get_individual_algebra_results(mxModel = fit, + progress_bar = FALSE) + +# the following plot will only show two data points because there is only +# two values for the definition variable k (0 or 1). + +plot(x = algebra_results[["a"]]$k, + y = algebra_results[["a"]]$algebra_result) +} diff --git a/man/mxsem_group_by.Rd b/man/mxsem_group_by.Rd new file mode 100644 index 0000000..32a8e46 --- /dev/null +++ b/man/mxsem_group_by.Rd @@ -0,0 +1,66 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/mxsem_group_by.R +\name{mxsem_group_by} +\alias{mxsem_group_by} +\title{mxsem_group_by} +\usage{ +mxsem_group_by( + mxModel, + grouping_variables, + parameters = c(".*"), + use_grepl = TRUE +) +} +\arguments{ +\item{mxModel}{mxModel with the entire data} + +\item{grouping_variables}{Variables used to split the data in groups} + +\item{parameters}{the parameters that should be group specific. By default +all parameters are group specific.} + +\item{use_grepl}{if set to TRUE, grepl is used to check which parameters are +group-specific. For instance, if parameters = "a" and use_grepl = TRUE, all parameters +whose label contains the letter "a" will be group specific. If use_grep = FALSE +only the parameter that has the label "a" is group specific.} +} +\value{ +mxModel with multiple groups. Use get_groups to extract the groups +} +\description{ +creates a multi-group model from an OpenMx model. +} +\details{ +mxsem_group_by creates a multi-group model by splitting the data found +in an mxModel object using dplyr's group_by function. The general idea is +as follows: + +1. The function extracts the data from mxModel +2. The data is split using the group_by function of dplyr with the variables +in grouping_variables +3. a separate model is set up for each group. All parameters that match +those specified in the parameters argument are group specific + +**Warning**: The multi-group model may differ from **lavaan**! For instance, +**lavaan** will automatically set the latent variances for all but the first +group free if the loadings are fixed to equality. Such automatic procedures +are not yet implemented in **mxsem**. +} +\examples{ +# THE FOLLOWING EXAMPLE IS ADAPTED FROM +# https://openmx.ssri.psu.edu/docs/OpenMx/latest/_static/Rdoc/mxModel.html +library(mxsem) + +model <- 'spatial =~ visual + cubes + paper + verbal =~ general + paragrap + sentence + math =~ numeric + series + arithmet' + +mg_model <- mxsem(model = model, + data = OpenMx::HS.ability.data) |> + # we want separate models for all combinations of grades and schools: + mxsem_group_by(grouping_variables = "school") |> + mxTryHard() + +# let's summarize the results: +summarize_multi_group_model(mg_model) +} diff --git a/man/print.multi_group_parameters.Rd b/man/print.multi_group_parameters.Rd new file mode 100644 index 0000000..f5655f8 --- /dev/null +++ b/man/print.multi_group_parameters.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/mxsem_group_by.R +\name{print.multi_group_parameters} +\alias{print.multi_group_parameters} +\title{print the multi_group_parameters} +\usage{ +\method{print}{multi_group_parameters}(x, ...) +} +\arguments{ +\item{x}{object from summarize_multi_group_model} + +\item{...}{not used} +} +\description{ +print the multi_group_parameters +} diff --git a/man/summarize_multi_group_model.Rd b/man/summarize_multi_group_model.Rd new file mode 100644 index 0000000..23a93e6 --- /dev/null +++ b/man/summarize_multi_group_model.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/mxsem_group_by.R +\name{summarize_multi_group_model} +\alias{summarize_multi_group_model} +\title{summarize_multi_group_model} +\usage{ +summarize_multi_group_model(multi_group_model) +} +\arguments{ +\item{multi_group_model}{multi group model created with mxsem_group_by} +} +\value{ +list with goup specific parameters and common parameters +} +\description{ +summarize the results of a multi group model created with mxsem_group_by +} +\examples{ +# THE FOLLOWING EXAMPLE IS ADAPTED FROM +# https://openmx.ssri.psu.edu/docs/OpenMx/latest/_static/Rdoc/mxModel.html +library(mxsem) + +model <- 'spatial =~ visual + cubes + paper + verbal =~ general + paragrap + sentence + math =~ numeric + series + arithmet' + +mg_model <- mxsem(model = model, + data = OpenMx::HS.ability.data) |> + # we want separate models for all combinations of grades and schools: + mxsem_group_by(grouping_variables = "school") |> + mxTryHard() + +# let's summarize the results: +summarize_multi_group_model(mg_model) +} diff --git a/tests/testthat/test-Political-Democracy.R b/tests/testthat/test-Political-Democracy.R index db57f5a..71467a3 100644 --- a/tests/testthat/test-Political-Democracy.R +++ b/tests/testthat/test-Political-Democracy.R @@ -1,6 +1,7 @@ test_that("Political Democracy works", { set.seed(123) library(lavaan) + library(mxsem) model <- ' # latent variable definitions diff --git a/tests/testthat/test-get_individual_algebras.R b/tests/testthat/test-get_individual_algebras.R new file mode 100644 index 0000000..fdee313 --- /dev/null +++ b/tests/testthat/test-get_individual_algebras.R @@ -0,0 +1,57 @@ +test_that("testing individual algebras", { + set.seed(123) + dataset <- simulate_moderated_nonlinear_factor_analysis(N = 50) + + model <- " + xi =~ x1 + x2 + x3 + eta =~ y1 + y2 + y3 + eta ~ a*xi + + !a0 + !a1 + a := a0 + data.k*a1 + " + + mod <- mxsem(model = model, + data = dataset) |> + mxTryHard() + + ind_alg <- get_individual_algebra_results(mxModel = mod) + + testthat::expect_true(length(ind_alg) == 1) + testthat::expect_true(names(ind_alg) == c("a")) + testthat::expect_true(all(colnames(ind_alg$a) == c("person", "k", "algebra_result"))) + + testthat::expect_error(get_individual_algebra_results(mxModel = mod, + algebra_names = "b")) + + model <- " + xi =~ x1 + x2 + x3 + eta =~ y1 + y2 + y3 + eta ~ {a := a0 + data.k*a1}*xi + " + + mod2 <- mxsem(model = model, + data = dataset) |> + mxTryHard() + + ind_alg <- get_individual_algebra_results(mxModel = mod2) + + testthat::expect_true(length(ind_alg) == 1) + testthat::expect_true(names(ind_alg) == c("a")) + testthat::expect_true(all(colnames(ind_alg$a) == c("person", "k", "algebra_result"))) + + testthat::expect_error(get_individual_algebra_results(mxModel = mod2, + algebra_names = "b")) + + model <- " + xi =~ x1 + x2 + x3 + eta =~ y1 + y2 + y3 + eta ~ xi + " + mod3 <- mxsem(model = model, + data = dataset) |> + mxTryHard() + + testthat::expect_error(get_individual_algebra_results(mxModel = mod3)) +}) diff --git a/tests/testthat/test-mnlfa.R b/tests/testthat/test-mnlfa.R index 651332c..cbe060a 100644 --- a/tests/testthat/test-mnlfa.R +++ b/tests/testthat/test-mnlfa.R @@ -22,4 +22,22 @@ test_that("mnlfa works", { testthat::expect_true(abs(omxGetParameters(mod)["a0"] - .7) < .1) testthat::expect_true(abs(omxGetParameters(mod)["a1"] - -.2) < .1) + model <- " + xi =~ x1 + x2 + x3 + eta =~ y1 + y2 + y3 + eta ~ {a0 + data.k*a1}*xi + " + + mod2 <- mxsem(model = model, + data = dataset) |> + mxTryHard() + + omxGetParameters(mod2) + + testthat::expect_true(abs(omxGetParameters(mod2)["a0"] - .7) < .1) + testthat::expect_true(abs(omxGetParameters(mod2)["a1"] - -.2) < .1) + + testthat::expect_true(abs(omxGetParameters(mod)["a0"] - omxGetParameters(mod2)["a0"]) < .001) + testthat::expect_true(abs(omxGetParameters(mod)["a1"] - omxGetParameters(mod2)["a1"]) < .001) + }) diff --git a/tests/testthat/test-multi-group.R b/tests/testthat/test-multi-group.R new file mode 100644 index 0000000..4bbc9ef --- /dev/null +++ b/tests/testthat/test-multi-group.R @@ -0,0 +1,75 @@ +test_that("multiplication works", { + library(lavaan) + library(mxsem) + set.seed(123) + HS.model <- 'visual =~ x1 + x2 + x3 + textual =~ x4 + x5 + x6 + speed =~ x7 + x8 + x9' + + fit <- cfa(HS.model, + data = HolzingerSwineford1939, + group = "school") + + mg_model <- mxsem(model = HS.model, + data = HolzingerSwineford1939) |> + # we want separate models for all combinations of grades and schools: + mxsem_group_by(grouping_variables = c("school")) |> + mxRun() + + testthat::expect_true(abs( + mg_model$fitfunction$result[[1]] - + -2*logLik(fit) + ) < .1) + + HS.model <- 'visual =~ x1 + x2 + x3 + textual =~ x4 + x5 + x6 + speed =~ x7 + x8 + x9 + # mxsem differs from lavaan in that it will not + # automatically set intercepts freed for all but one group + # once other parameters are set to equality + visual ~ 0*1 + textual ~ 0*1 + speed ~ 0*1 + ' + + fit <- cfa(HS.model, + data = HolzingerSwineford1939, + group = "school", + group.equal = c("intercepts", "loadings")) + + mg_model <- mxsem(model = HS.model, + data = HolzingerSwineford1939) |> + # we want separate models for all combinations of grades and schools: + mxsem_group_by(grouping_variables = c("school"), + parameters = mxsem::unicode_undirected()) |> + mxTryHard() + + testthat::expect_true(abs( + mg_model$fitfunction$result[[1]] - + -2*logLik(fit) + ) < .1) + + + HS.model <- 'visual =~ x1 + x2 + x3 + textual =~ x4 + x5 + x6 + speed =~ x7 + x8 + x9 + ' + + fit <- cfa(HS.model, + data = HolzingerSwineford1939, + group = "school", + group.equal = c("loadings")) + + mg_model <- mxsem(model = HS.model, + data = HolzingerSwineford1939) |> + # we want separate models for all combinations of grades and schools: + mxsem_group_by(grouping_variables = c("school"), + parameters = c(mxsem::unicode_undirected(), + "one")) |> + mxTryHard() + + testthat::expect_true(abs( + mg_model$fitfunction$result[[1]] - + -2*logLik(fit) + ) < .1) +}) diff --git a/vignettes/Moderated-Nonlinear-Factor-Analysis.Rmd b/vignettes/Moderated-Nonlinear-Factor-Analysis.Rmd index b357812..5593ad6 100644 --- a/vignettes/Moderated-Nonlinear-Factor-Analysis.Rmd +++ b/vignettes/Moderated-Nonlinear-Factor-Analysis.Rmd @@ -41,9 +41,6 @@ are taken directly from that script: # at https://osf.io/527zr library(mokken) -#> Loading required package: poLCA -#> Loading required package: scatterplot3d -#> Loading required package: MASS ## Load data data("DS14", package="mokken") @@ -161,7 +158,7 @@ With **mxsem**, this can be implemented as follows: ```r -mlnfa_syntax <- " +mnlfa_syntax <- " ==== MNLFA ==== SI =~ {lSi_1 := lSi0_1 + lSi1_1*data.Age + lSi2_1*data.Male }*Si1 + @@ -232,7 +229,7 @@ Finally, we pass the syntax to the `mxsem()`-function to create an **OpenMx** mo ```r library(mxsem) -mlnfa_model <- mxsem(model = mlnfa_syntax, +mnlfa_model <- mxsem(model = mnlfa_syntax, data = DS14, # we scaled the latent variables manually, # so we will set all automatic scalings to FALSE: @@ -248,8 +245,8 @@ The model can now be fitted using `mxRun()` or `mxTryHard()`. ```r -mlnfa_model <- mxRun(mlnfa_model) -summary(mlnfa_model) +mnlfa_model <- mxRun(mnlfa_model) +summary(mnlfa_model) ```
@@ -406,8 +403,8 @@ summary(mlnfa_model) #> RMSEA: 0 [95% CI (NA, NA)] #> Prob(RMSEA <= 0.05): NA #> To get additional fit indices, see help(mxRefModels) -#> timestamp: 2023-08-02 17:35:53 -#> Wall clock time: 226.1492 secs +#> timestamp: 2023-08-25 14:36:18 +#> Wall clock time: 292.0618 secs #> optimizer: SLSQP #> OpenMx version number: 2.21.8 #> Need help? See help(mxSummary) @@ -418,4 +415,48 @@ Checking the regression coefficients `lSi1_1`, `lSi1_2`, ... will tell us if there is a linear change across age or if individuals with `Male = 0` differ from individuals with `Male = 1`. +## Plotting Individual Parameters + +MNLFA predicts individual parameter values (e.g., `lSi_1`) using definition +variables. To get a better picture of the individual parameters, **mxsem** +provides the `get_individual_algebra_results` function. This function will compute +for each algebra the individual parameters. Depending on the sample size and +the number of algebras in the model, this may take some time. Therefore, we will +only extract the individual parameter values for `lSi_1` as an example. + + +```r +lSi_1 <- get_individual_algebra_results(mxModel = mnlfa_model, + algebra_names = "lSi_1", + progress_bar = FALSE) +head(lSi_1$lSi_1) +#> person Age Male algebra_result +#> 1 1 0.03136864 1 0.8699212 +#> 2 2 -0.44266587 1 0.9246368 +#> 3 3 -0.82189349 1 0.9684094 +#> 4 4 -0.63227968 1 0.9465231 +#> 5 5 0.03136864 1 0.8699212 +#> 6 6 -0.06343826 0 0.7408356 +``` + +The function will return a list with data frames for all requested algebras. Here, +the list has only one element: `lSi_1`. The data frame will have fields for the +person, the definition variables used in the algebra (`Age` and `Male` in this case) +and the person specific parameter (`algebra_result`). We can plot these results +as follows: + + +```r +library(ggplot2) +ggplot(data = lSi_1$lSi_1, + aes(x = Age, + y = algebra_result, + color = factor(Male))) + + ylab("Individual Parameter Value for lSi_1") + + geom_point() +``` + + +![](figures/mnlfa_lS_1.png) + ## References diff --git a/vignettes/Moderated-Nonlinear-Factor-Analysis.mxsemmd b/vignettes/Moderated-Nonlinear-Factor-Analysis.mxsemmd index 1b512e0..88b24e5 100644 --- a/vignettes/Moderated-Nonlinear-Factor-Analysis.mxsemmd +++ b/vignettes/Moderated-Nonlinear-Factor-Analysis.mxsemmd @@ -18,6 +18,7 @@ knitr::opts_chunk$set( ) library(mxsem) library(lavaan) +library(ggplot2) set.seed(123) # to knit this vignette, use # knitr::knit(input = "Moderated-Nonlinear-Factor-Analysis.mxsemmd", @@ -222,7 +223,7 @@ We want to predict each and every parameter in this model using `Age` and `Male` With **mxsem**, this can be implemented as follows: ```{r} -mlnfa_syntax <- " +mnlfa_syntax <- " ==== MNLFA ==== SI =~ {lSi_1 := lSi0_1 + lSi1_1*data.Age + lSi2_1*data.Male }*Si1 + @@ -292,7 +293,7 @@ Finally, we pass the syntax to the `mxsem()`-function to create an **OpenMx** mo ```{r} library(mxsem) -mlnfa_model <- mxsem(model = mlnfa_syntax, +mnlfa_model <- mxsem(model = mnlfa_syntax, data = DS14, # we scaled the latent variables manually, # so we will set all automatic scalings to FALSE: @@ -306,17 +307,17 @@ mlnfa_model <- mxsem(model = mlnfa_syntax, The model can now be fitted using `mxRun()` or `mxTryHard()`. ```{r, include = FALSE} -mlnfa_model <- mxRun(mlnfa_model) +mnlfa_model <- mxRun(mnlfa_model) ``` ```{r, eval = FALSE} -mlnfa_model <- mxRun(mlnfa_model) -summary(mlnfa_model) +mnlfa_model <- mxRun(mnlfa_model) +summary(mnlfa_model) ```
Show summary ```{r, echo = FALSE} -summary(mlnfa_model) +summary(mnlfa_model) ```
@@ -324,4 +325,49 @@ Checking the regression coefficients `lSi1_1`, `lSi1_2`, ... will tell us if there is a linear change across age or if individuals with `Male = 0` differ from individuals with `Male = 1`. +## Plotting Individual Parameters + +MNLFA predicts individual parameter values (e.g., `lSi_1`) using definition +variables. To get a better picture of the individual parameters, **mxsem** +provides the `get_individual_algebra_results` function. This function will compute +for each algebra the individual parameters. Depending on the sample size and +the number of algebras in the model, this may take some time. Therefore, we will +only extract the individual parameter values for `lSi_1` as an example. + +```{r} +lSi_1 <- get_individual_algebra_results(mxModel = mnlfa_model, + algebra_names = "lSi_1", + progress_bar = FALSE) +head(lSi_1$lSi_1) +``` + +The function will return a list with data frames for all requested algebras. Here, +the list has only one element: `lSi_1`. The data frame will have fields for the +person, the definition variables used in the algebra (`Age` and `Male` in this case) +and the person specific parameter (`algebra_result`). We can plot these results +as follows: + +```{r, eval = FALSE} +library(ggplot2) +ggplot(data = lSi_1$lSi_1, + aes(x = Age, + y = algebra_result, + color = factor(Male))) + + ylab("Individual Parameter Value for lSi_1") + + geom_point() +``` + +```{r, include = FALSE} +library(ggplot2) +png(filename = "figures/mnlfa_lS_1.png", width = 10, height = 10, units = "cm", res = 200) +ggplot(data = lSi_1$lSi_1, + aes(x = Age, + y = algebra_result, + color = factor(Male))) + + ylab("Individual Parameter Value for lSi_1") + + geom_point() +dev.off() +``` +![](figures/mnlfa_lS_1.png) + ## References diff --git a/vignettes/figures/mnlfa_lS_1.png b/vignettes/figures/mnlfa_lS_1.png new file mode 100644 index 0000000..21293a5 Binary files /dev/null and b/vignettes/figures/mnlfa_lS_1.png differ