functions to plot individual parameters when using definition variabl…

…es in algebras (e.g., mnlfa)

DESCRIPTION

@@ -18,7 +18,8 @@ Imports:
   Rcpp (>= 1.0.10), 
   stats,
   methods,
-  dplyr
+  dplyr,
+  utils
 LinkingTo: Rcpp
 RoxygenNote: 7.2.3
 Encoding: UTF-8

NAMESPACE

@@ -2,6 +2,7 @@
 
 S3method(print,multi_group_parameters)
 export(get_groups)
+export(get_individual_algebra_results)
 export(mxsem)
 export(mxsem_group_by)
 export(parameters)
@@ -17,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)

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)
+}

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))
+})

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)
+
 })

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)
 ```
 
 <details>
@@ -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:20:14 
+#> Wall clock time: 357.9646 secs 
 #> optimizer:  SLSQP 
 #> OpenMx version number: 2.21.8 
 #> Need help?  See help(mxSummary)
@@ -418,4 +415,47 @@ 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()
+```
+
+![plot of chunk unnamed-chunk-25](figure/unnamed-chunk-25-1.png)
+
 ## References