# Multiple analysis functions

#### Phil Chalmers

#### 2025-03-28

This vignette demonstrates one of the newer features in the `SimDesign` package pertaining to multiple analysis function definitions that can be selected for each `Design` condition or whenever they are applicable. The purpose of providing multiple analysis functions is to

1. Remove the less readable if-then-else combinations that can appear when writing simulation code, where specific analysis functions are not intended to be used on a given generated dataset,
2. Provide better automatic naming of the analysis results across independent subroutines,
3. Construct more readable code when the `Analyse()` function contains too much code to easily track, and to
4. Create a more modular approach to isolating the analysis functions for the purpose of redistribution or reusing in related projects

Functionality speaking, this type of organization does not change how `SimDesign` generally operates. For that reason, the coding style presented in this vignette can be considered optional. However, if any of the above points resonate well with you then following the details of this coding organization style may prove useful.

# 1 Description of structure

The usual work-flow with `SimDesign` requires first calling `SimFunctions()` to generate a working template, such as the following.

    SimDesign::SimFunctions()

    ## #-------------------------------------------------------------------
    ##
    ## library(SimDesign)
    ##
    ## Design <- createDesign(factor1 = NA,
    ##                        factor2 = NA)
    ##
    ## #-------------------------------------------------------------------
    ##
    ## Generate <- function(condition, fixed_objects) {
    ##     dat <- data.frame()
    ##     dat
    ## }
    ##
    ## Analyse <- function(condition, dat, fixed_objects) {
    ##     ret <- nc(stat1 = NaN, stat2 = NaN)
    ##     ret
    ## }
    ##
    ## Summarise <- function(condition, results, fixed_objects) {
    ##     ret <- c(bias = NaN, RMSE = NaN)
    ##     ret
    ## }
    ##
    ## #-------------------------------------------------------------------
    ##
    ## res <- runSimulation(design=Design, replications=2, generate=Generate,
    ##                      analyse=Analyse, summarise=Summarise)
    ## res

which uses the default `nAnalyses=1` to generate only a single `Analyse()` function. In the context of multiple analysis functions, however, users may be more interested in passing the number of analysis functions they believe they will need in their simulation (e.g., if analyzing a \(t\)-test setup to compare the Welch versus independent samples t-test, then two analysis functions should be used). Passing `nAnalyses=2` to `SimFunctions()` creates the following template:

    SimDesign::SimFunctions(nAnalyses = 2)

    ## #-------------------------------------------------------------------
    ##
    ## library(SimDesign)
    ##
    ## Design <- createDesign(factor1 = NA,
    ##                        factor2 = NA)
    ##
    ## #-------------------------------------------------------------------
    ##
    ## Generate <- function(condition, fixed_objects) {
    ##     dat <- data.frame()
    ##     dat
    ## }
    ##
    ## Analyse.A1 <- function(condition, dat, fixed_objects) {
    ##     ret <- nc(stat1 = NaN, stat2 = NaN)
    ##     ret
    ## }
    ##
    ## Analyse.A2 <- function(condition, dat, fixed_objects) {
    ##     ret <- nc(stat1 = NaN, stat2 = NaN)
    ##     ret
    ## }
    ##
    ## #-------------------------------------------------------------------
    ##
    ## Summarise <- function(condition, results, fixed_objects) {
    ##     ret <- c(bias = NaN, RMSE = NaN)
    ##     ret
    ## }
    ##
    ## #-------------------------------------------------------------------
    ##
    ## res <- runSimulation(design=Design, replications=2,generate=Generate,
    ##                      analyse=list(A1=Analyse.A1, A2=Analyse.A2),
    ##                      summarise=Summarise)
    ## res

Notice in this case that there are two `Analyse.#()` definitions constructed, and when passed to `runSimulation()` these are organized as a named `list`. The names of the list will ultimately be attached to the names of the analysis objects so that there is no ambiguity in the outputted information. However, the inputs to the `Analyse()` functions will always be the same, as the `dat` object formed by the `Generate()` call will be passed to each of these Analyse definitions (hence, the generate data is held constant across the respective analyses).

The above template should of course be modified to replace the less useful names of the `Analyse.#()` components. By default users will want to change these to something like `Analyse.some_statistic`, `Analyse.some_other_statistic`, …, `Analyse.some_other_other_statistic`, and so on, where the number of `Analyse.#` function definitions will ultimately end up in the `runSimulation(..., Analyse=list())` input. Supplying better names to the named `list` component is also recommended as these will be used to name the associated results in the `Summarise()` step.

Finally, note that all the rules about objects and object naming from the typical single Analyse function still apply and are properly checked internally for suitable names and consistency. The independently defined Analyse functions are also interchangable and removable/replaceable, which makes the structure of the Generate-Analyse-Summarise setup more modular with respect to the analysis components.

## 1.1 An example

The following code is adopted from the Wiki, and so details about the simulation should be obtained from that source.

    library(SimDesign)
    # SimFunctions(nAnalyses = 2)

    sample_sizes <- c(250, 500, 1000)
    nitems <- c(10, 20)
    Design <- createDesign(sample_size = sample_sizes,
                           nitems = nitems)

    # create list of additional parameters which are fixed across conditions
    set.seed(1)
    pars_10 <- rbind(a = round(rlnorm(10, .3, .5)/1.702, 2),
                     d = round(rnorm(10, 0, .5)/1.702, 2))
    pars_20 <- rbind(a = round(rlnorm(20, .3, .5)/1.702, 2),
                     d = round(rnorm(20, 0, .5)/1.702, 2))
    pars <- list(ten=pars_10, twenty=pars_20)

    P_logit <- function(a, d, Theta) exp(a * Theta + d) / (1 + exp(a * Theta + d))
    P_ogive <- function(a, d, Theta) pnorm(a * Theta + d)

    Generate <- function(condition, fixed_objects) {

        N <- condition$sample_size
        nitems <- condition$nitems
        nitems_name <- ifelse(nitems == 10, 'ten', 'twenty')

        #extract objects from fixed_objects
        a <- fixed_objects[[nitems_name]]['a', ]
        d <- fixed_objects[[nitems_name]]['d', ]

        dat <- matrix(NA, N, nitems)
        colnames(dat) <- paste0('item_', 1:nitems)
        Theta <- rnorm(N)
        for(j in 1:nitems){
            p <- P_ogive(a[j], d[j], Theta)
            for(i in 1:N)
                dat[i,j] <- sample(c(1,0), 1, prob = c(p[i], 1 - p[i]))
        }
        as.data.frame(dat) #data.frame works nicer with lavaan
    }

    Analyse.FIML <- function(condition, dat, fixed_objects) {
        mod <- mirt(dat, 1L, verbose=FALSE)
        if(!extract.mirt(mod, 'converged')) stop('mirt did not converge')
        cfs <- mirt::coef(mod, simplify = TRUE, digits = Inf)
        FIML_as <- cfs$items[,1L] / 1.702

        ret <- c(as=unname(FIML_as))
        ret
    }

    Analyse.DWLS <- function(condition, dat, fixed_objects) {
        nitems <- condition$nitems
        lavmod <- paste0('F =~ ', paste0('NA*', colnames(dat)[1L], ' + '),
                         paste0(colnames(dat)[-1L], collapse = ' + '),
                         '\nF ~~ 1*F')
        lmod <- sem(lavmod, dat, ordered = colnames(dat))
        if(!lavInspect(lmod, 'converged')) stop('lavaan did not converge')
        cfs2 <- lavaan::coef(lmod)
        DWLS_alpha <- cfs2[1L:nitems]
        const <- sqrt(1 - DWLS_alpha^2)
        DWLS_as <- DWLS_alpha / const

        ret <- c(as=unname(DWLS_as))
        ret
    }

    Summarise <- function(condition, results, fixed_objects) {
        nitems <- condition$nitems
        nitems_name <- ifelse(nitems == 10, 'ten', 'twenty')

        #extract objects from fixed_objects
        a <- fixed_objects[[nitems_name]]['a', ]
        pop <- c(a, a)

        obt_bias <- bias(results, pop)
        obt_RMSE <- RMSE(results, pop)
        ret <- c(bias=obt_bias, RMSE=obt_RMSE)
        ret
    }

    res <- runSimulation(Design, replications=100, verbose=FALSE, parallel=TRUE,
                         generate=Generate,
                         analyse=list(FIML=Analyse.FIML, DWLS=Analyse.DWLS),
                         summarise=Summarise, filename = 'mirt_lavaan',
                         packages=c('mirt', 'lavaan'), fixed_objects=pars)

    res

    ## # A tibble: 6 × 86
    ##   sample_size nitems bias.FIML.as1 bias.FIML.as2 bias.FIML.as3 bias.FIML.as4
    ##         <dbl>  <dbl>         <dbl>         <dbl>         <dbl>         <dbl>
    ## 1         250     10     0.0099887     0.013457     -0.020161      0.18616
    ## 2         500     10    -0.021136      0.018720     -0.018173      0.10424
    ## 3        1000     10    -0.0039473    -0.0014778    -0.013393      0.077890
    ## 4         250     20     0.060647      0.011787      0.064688     -0.026314
    ## 5         500     20     0.019202      0.045376      0.021837     -0.0063587
    ## 6        1000     20     0.019976      0.018608      0.0079817    -0.011498
    ## # ℹ 80 more variables: bias.FIML.as5 <dbl>, bias.FIML.as6 <dbl>,
    ## #   bias.FIML.as7 <dbl>, bias.FIML.as8 <dbl>, bias.FIML.as9 <dbl>,
    ## #   bias.FIML.as10 <dbl>, bias.DWLS.as1 <dbl>, bias.DWLS.as2 <dbl>,
    ## #   bias.DWLS.as3 <dbl>, bias.DWLS.as4 <dbl>, bias.DWLS.as5 <dbl>,
    ## #   bias.DWLS.as6 <dbl>, bias.DWLS.as7 <dbl>, bias.DWLS.as8 <dbl>,
    ## #   bias.DWLS.as9 <dbl>, bias.DWLS.as10 <dbl>, RMSE.FIML.as1 <dbl>,
    ## #   RMSE.FIML.as2 <dbl>, RMSE.FIML.as3 <dbl>, RMSE.FIML.as4 <dbl>, …

In this particular formulation the `mirt` and `lavaan` package analyses have been completely isolated into their own respective functions, and in principle could therefore be analyzed independently in future simulation studies. This adds a nicer layer of potential modularity to the Analyse portion of the `SimDesign` framework, where re-using or modifying previous `SimDesign` code should be less painful.

# 2 `AnalyseIf()`

In situations where analysis functions defined in the `analyse` list should only be applied in certain design conditions, users can include an `AnalyseIf()` definition at the beginning of their respective functions to ensure the analyses are only executed when the provided logical is `TRUE`. This logical ensures the data generation conditions are suitable for the analysis function to be investigated; otherwise, it is skipped over in the generate-analyse-summarise work-flow.

As a continuation from above, say that an investigator was also interested in recovering the slope parameters of a factor analysis model where the observed indicator variable were continuous as well as discrete. The `Design` definition may therefore look like the following.

    Design <- createDesign(sample_size = sample_sizes,
                           nitems = nitems,
                           indicators = c('discrete', 'continuous'))
    Design

    ## # A tibble: 12 × 3
    ##    sample_size nitems indicators
    ##          <dbl>  <dbl> <chr>
    ##  1         250     10 discrete
    ##  2         500     10 discrete
    ##  3        1000     10 discrete
    ##  4         250     20 discrete
    ##  5         500     20 discrete
    ##  6        1000     20 discrete
    ##  7         250     10 continuous
    ##  8         500     10 continuous
    ##  9        1000     10 continuous
    ## 10         250     20 continuous
    ## 11         500     20 continuous
    ## 12        1000     20 continuous

Provided that the `Generate()` step utilized this `indicators` character vector, this would imply that the `dat` object returned from `Generate()` could consist of discrete or continuous data. In the case of continuous indicator variables, `lavaan` could be used as it supports such indicator types; however, `mirt` cannot. So, to ensure that only the analysis function pertaining to `lavaan` is used one could include the following replacement definition that used `mirt`, but now includes an `AnalyseIf()` logical given the `indicators` variable’s state.

    Analyse.FIML <- function(condition, dat, fixed_objects) {
        AnalyseIf(condition$indicators == 'discrete')
        # equivalently:
        #   AnalyseIf(indicators == 'discrete', condition)
        #   with(condition, AnalyseIf(indicators == 'discrete'))
        mod <- mirt(dat, 1L, verbose=FALSE)
        if(!extract.mirt(mod, 'converged')) stop('mirt did not converge')
        cfs <- coef(mod, simplify = TRUE, digits = Inf)
        FIML_as <- cfs$items[,1L] / 1.702

        ret <- c(as=unname(FIML_as))
        ret
    }

Using this definition the final object returned by `runSimulation()` will provide suitable `NA` placeholders (where appropriate). For continuous indicators the results will be presented as though `mirt` was never used for the continuous indicator conditions controlled by the `Design` object.

Note that a similar logic can be made for multiple `Generate()` functions, templated via `SimDesign::SimFunctions(nGenerate = 2)`, and requiring specific `GenerateIf()` tests to indicate which generation function should be used. For those interested in this type of structure, the following could be used to separate the discrete/continuous data generation per `Design` row:

    Generate.G1 <- function(condition, fixed_objects) {
        GenerateIf(condition$indicators == 'discrete')
        ...
        dat
    }

    Generate.G2 <- function(condition, fixed_objects) {
        GenerateIf(condition$indicators == 'continuous')
        ...
        dat
    }

    res <- runSimulation(design=Design, replications=1000,
                         generate=list(G1=Generate.G1, G2=Generate.G2),
                         analyse=list(DWLS=Analyse.DWLS, FIML=Analyse.FIML),
                         summarise=Summarise)

## 2.1 Applying one analyse function per-condition

Interestingly, `AnalyseIf()` could also be used to select only one analysis function at a time given the components in the `Design` object. For instance, if the `Design` definition were constructed using

    Design <- createDesign(sample_size = sample_sizes,
                           nitems = nitems,
                           method = c('FIML', 'DWLS'))
    Design

    ## # A tibble: 12 × 3
    ##    sample_size nitems method
    ##          <dbl>  <dbl> <chr>
    ##  1         250     10 FIML
    ##  2         500     10 FIML
    ##  3        1000     10 FIML
    ##  4         250     20 FIML
    ##  5         500     20 FIML
    ##  6        1000     20 FIML
    ##  7         250     10 DWLS
    ##  8         500     10 DWLS
    ##  9        1000     10 DWLS
    ## 10         250     20 DWLS
    ## 11         500     20 DWLS
    ## 12        1000     20 DWLS

and the analysis functions above were supplied defined as

    Analyse.FIML <- function(condition, dat, fixed_objects) {
        AnalyseIf(method == 'FIML', condition)
        #...
    }

    Analyse.DWLS <- function(condition, dat, fixed_objects) {
        AnalyseIf(method == 'DWLS', condition)
        # ...
    }

    # ...
    res <- runSimulation(Design, replications=100, verbose=FALSE, parallel=TRUE,
                         generate=Generate,
                         analyse=list(Analyse.FIML, Analyse.DWLS),
                         summarise=Summarise, filename = 'mirt_lavaan',
                         packages=c('mirt', 'lavaan'), fixed_objects=pars)

then only one analysis function will be applied at a time in the simulation experiment. Note that in this case there is no need to append ‘MML’ or ‘DWLS’ to the `results` objects as this becomes redundant with the `method` column in the `Design` object, and so the `analyse` list input is specified as an unnamed list (cf. earlier when the input was named, which appended `MML.` and `DWLS.` to the `results` output in `Summarise()`).

Users may find this a more natural setup than having to merge all analysis information into a single `Analyse()` definition. The downside, however, is that the analysis function will be applied to different generated datasets, which while theoretically unbiased could have ramifications should the analysis functions throw errors at different rates (even when explicitly supplying a `seed` vector input to `runSimulation()`).