The `seqwrap` function lets the user specify which modeling algorithm to use when modeling data provided as a data frame and metadata describing each sample. Target-specific models are passed to the results together with optional user-specified summary and model-evaluation functions. The function aims not to limit the user to the preferences of specific models. This allows for comparisons between modeling frameworks and flexibility for many types of high-dimensional data. In this vignette, three popular packages are applied to the same underlying data with varying model definitions.

### Fitting count data to a negative binomial model using `glmmTBM`

We will start by simulating data from a simple experiment where count data has been gathered from 25 clusters in two conditions (e.g. treatment and control).

    library(seqwrap)
    library(glmmTMB)

    # Simulate data
    dat <- seqwrap::simcounts(n_genes = 1000,
                        n_samples = 50,
                        beta_0 = 2,
                        sigma_0 = 0.5,
                        beta_1 = 0.5,,
                        sigma_1 = 2,
                        b_0 = 1,
                        clusters = 25,
                        overdispersion_min = c(1, 10))

    # Calculate library sizes
    dat$metadata$libsize <- log(colSums(dat$data[, -1]) )

Using `seqwrap_compose` we can populate a `swcontainer` object with all needed information to fit models. The `swcontainer` objects is basically a list with data, meta data and arguments that gets passed the model algorithm. The `swcontainer` object also contain reference to the desired model function.

    swobject <- seqwrap_compose(
      data = dat$data, # A data frame with the first column as the target identifier
      metadata = dat$metadata, # Meta data bout samples, containing predictors
      samplename = "sample", # The name sample id column
      modelfun = glmmTMB::glmmTMB, # The model fitting function
      # A list of argumnets supplied to the modelfunction
      arguments = list(
        formula = y ~ x + (1 | cluster),
        family = glmmTMB::nbinom1
      ),
      summary_fun = NULL,
      eval_fun = NULL,
      additional_vars = NULL
    )

The `seqwrap` function will iterate over all targets in the data and fit the model and apply the summary and evaluation function to each model. Notice that we can use `return_models = TRUE` to save model in the output. Below we also use `subset = 1:10` to only do the first 10 targets. This is useful when we want to check if the summary and evaluation functions work as intended. The `swobject` created above is updated when we supply new arguments to `seqwrap`.

    models <- seqwrap(swobject,
                      return_models = TRUE,
                      subset = 1:10)

    ── seqwrap ─────────────────────────────────────────────────────────────────────

    Initiating clusters for parallel processing with 1 core
    Fitting 10 targets in 4 chunks of up to 3 targets
    The progress bar updates 4 times (once per 1 chunk completed), in steps of
    about 25%; use a smaller `chunk_size` for more frequent updates
    Merging and modelling data
    sep 03 2026 12:49:46: Completed model fitting and evaluation. Elapsed time was
    0.05 minutes

The resulting object can be summarized using `seqwrap_summarise`

    modelsum <- seqwrap_summarise(models)

    ── seqwrap summarise ───────────────────────────────────────────────────────────

    • A total of 50 samples and 10 targets where used in `glmmTMB::glmmTMB` with
    arguments `list(formula = y ~ x + (1 | cluster), family = glmmTMB::nbinom1)`

    • Attempting to combine results from the provided summary function.

    • Attempting to combine results from the provided evaluations function.

    • Attempting to summarise errors and warnings from the fitting process.

    ── Model summaries ─────────────────────────────────────────────────────────────

    • 10 targets have associated summaries

    ── Model evaluations ───────────────────────────────────────────────────────────

    • 10 targets have associated evaluation results

    ── Errors and warnings ─────────────────────────────────────────────────────────

    ℹ No errors or warnings were raised

    ── Combined results ────────────────────────────────────────────────────────────

    ℹ Combined results have been generated and were
    silently returned.

We have left the arguments `summary_fun` and `eval_fun` as NULL. This means that `seqwrap` will combine model parameter estimates using `broom.mixed::`/`broom::tidy`, and report the convergence diagnostics that the fitting algorithm itself supplies. The evaluation is algorithm-sensitive: `lme4` reports an optimiser code together with singularity messages, `nlme` reports whether an approximate variance-covariance matrix could be formed, and `glmmTMB` reports both an optimiser code and whether the Hessian was positive definite. These are reduced to a common set of columns so that they can be combined across targets. These evaluations are thought of as minimal summaries of the model and might not completely satisfy the needs of the user. From our example above, the fitted parameters from the models are returned in `modelsum$summaries` and model evaluations are returned in `modelsum$evaluations`.

The functions are exported from the seqwrap package and we can inspect them:

    seqwrap::generic_summary

    function(x) {
      out <- broom.mixed::tidy(x)
      return(out)
    }
    <bytecode: 0x00000298ca53a708>
    <environment: namespace:seqwrap>

    seqwrap::generic_evaluation

    function(x, ...) {
      UseMethod("generic_evaluation")
    }
    <bytecode: 0x00000298c747e6c8>
    <environment: namespace:seqwrap>

We encourage users to define their own custom functions.

### Fitting transformed count data to a Gaussian model using `lme4`

Instead of counts, let’s say we are interested in analyzing counts per million, expressed as $\text{CPM} = \frac{\text{Count}}{\text{Library size}} \times 10^6$. We will convert our simulated count data to CPM.

    lib_sizes <- colSums(dat$data[,-1]) / 1e6

    cpm <- sweep(dat$data[,-1], MARGIN = 2,
                 STATS = lib_sizes,
                 FUN = "/")

    dat$data <- data.frame(cbind(dat$data[,1], cpm))

We will create a new `swcontainer` using `lme4::lmer` as our fitting algorithm.

    swobject <- seqwrap_compose(
      data = dat$data, # A data frame with the first column as the target identifier
      metadata = dat$metadata, # Meta data bout samples, containing predictors
      samplename = "sample", # The name sample id column
      modelfun = lme4::lmer, # The model fitting function
      # A list of argumnets supplied to the modelfunction
      arguments = list(
        formula = y ~ x + (1 | cluster)
      ),
      summary_fun = NULL,
      eval_fun = NULL,
      additional_vars = NULL
    )

Next we can run `seqwrap` to get estimates for our genes, to save some time we will start by using a subset of the simulated genes.

    models <- seqwrap(swobject,
                      return_models = TRUE,
                      subset = 1:10)

    ── seqwrap ─────────────────────────────────────────────────────────────────────

    Initiating clusters for parallel processing with 1 core
    Fitting 10 targets in 4 chunks of up to 3 targets
    The progress bar updates 4 times (once per 1 chunk completed), in steps of
    about 25%; use a smaller `chunk_size` for more frequent updates
    Merging and modelling data
    sep 03 2026 12:49:48: Completed model fitting and evaluation. Elapsed time was
    0.03 minutes

The generic summary function will provide fitted parameter statistics.

    modelsum <- seqwrap_summarise(models)

    ── seqwrap summarise ───────────────────────────────────────────────────────────

    • A total of 50 samples and 10 targets where used in `lme4::lmer` with
    arguments `list(formula = y ~ x + (1 | cluster))`

    • Attempting to combine results from the provided summary function.

    • Attempting to combine results from the provided evaluations function.

    • Attempting to summarise errors and warnings from the fitting process.

    ── Model summaries ─────────────────────────────────────────────────────────────

    • 10 targets have associated summaries

    ── Model evaluations ───────────────────────────────────────────────────────────

    • 10 targets have associated evaluation results

    ── Errors and warnings ─────────────────────────────────────────────────────────

    ℹ No errors or warnings were raised

    ── Combined results ────────────────────────────────────────────────────────────

    ℹ Combined results have been generated and were
    silently returned.

    modelsum$summaries

       target   effect    group            term      estimate std.error   statistic
    1   gene1    fixed     <NA>     (Intercept)  5.399693e+02  62.68847   8.6135345
    2   gene1    fixed     <NA>               x -5.343370e+02  88.65489  -6.0271576
    3   gene1 ran_pars  cluster sd__(Intercept)  2.205648e-05        NA          NA
    4   gene1 ran_pars Residual sd__Observation  3.134424e+02        NA          NA
    5  gene10    fixed     <NA>     (Intercept)  5.701182e+02  50.01939  11.3979438
    6  gene10    fixed     <NA>               x -5.362535e+02  70.73810  -7.5808289
    7  gene10 ran_pars  cluster sd__(Intercept)  0.000000e+00        NA          NA
    8  gene10 ran_pars Residual sd__Observation  2.500970e+02        NA          NA
    9   gene2    fixed     <NA>     (Intercept)  8.464301e+02 119.98004   7.0547574
    10  gene2    fixed     <NA>               x -8.322984e+02 168.48314  -4.9399509
    11  gene2 ran_pars  cluster sd__(Intercept)  7.105034e+01        NA          NA
    12  gene2 ran_pars Residual sd__Observation  5.956779e+02        NA          NA
    13  gene3    fixed     <NA>     (Intercept)  1.912564e+03 113.71921  16.8183000
    14  gene3    fixed     <NA>               x -1.723231e+03 160.82325 -10.7150633
    15  gene3 ran_pars  cluster sd__(Intercept)  3.976711e-07        NA          NA
    16  gene3 ran_pars Residual sd__Observation  5.685961e+02        NA          NA
    17  gene4    fixed     <NA>     (Intercept)  8.382400e+02  55.13350  15.2038242
    18  gene4    fixed     <NA>               x -7.413738e+02  77.78148  -9.5314947
    19  gene4 ran_pars  cluster sd__(Intercept)  1.918541e+01        NA          NA
    20  gene4 ran_pars Residual sd__Observation  2.749991e+02        NA          NA
    21  gene5    fixed     <NA>     (Intercept)  1.172927e+03 150.55430   7.7907245
    22  gene5    fixed     <NA>               x -1.170599e+03 212.91592  -5.4979381
    23  gene5 ran_pars  cluster sd__(Intercept)  2.199980e-01        NA          NA
    24  gene5 ran_pars Residual sd__Observation  7.527714e+02        NA          NA
    25  gene6    fixed     <NA>     (Intercept)  1.869216e+03 188.33816   9.9247854
    26  gene6    fixed     <NA>               x  1.428747e+02 266.35038   0.5364164
    27  gene6 ran_pars  cluster sd__(Intercept)  0.000000e+00        NA          NA
    28  gene6 ran_pars Residual sd__Observation  9.416908e+02        NA          NA
    29  gene7    fixed     <NA>     (Intercept)  1.311955e+03 122.79115  10.6844401
    30  gene7    fixed     <NA>               x -1.096413e+03 173.65292  -6.3138163
    31  gene7 ran_pars  cluster sd__(Intercept)  1.772548e-07        NA          NA
    32  gene7 ran_pars Residual sd__Observation  6.139558e+02        NA          NA
    33  gene8    fixed     <NA>     (Intercept)  5.167672e+02 395.56412   1.3064056
    34  gene8    fixed     <NA>               x  5.926088e+03 545.06122  10.8723352
    35  gene8 ran_pars  cluster sd__(Intercept)  4.451153e+02        NA          NA
    36  gene8 ran_pars Residual sd__Observation  1.927082e+03        NA          NA
    37  gene9    fixed     <NA>     (Intercept)  7.447068e+02  90.65968   8.2143108
    38  gene9    fixed     <NA>               x -4.356680e+02 121.43194  -3.5877546
    39  gene9 ran_pars  cluster sd__(Intercept)  1.454578e+02        NA          NA
    40  gene9 ran_pars Residual sd__Observation  4.293267e+02        NA          NA

### Using `nlme::lme` and `nlme::gls`, need for `additional_vars`

We will now try to fit the same model as previously fitted in `lme4::lmer` using `nlme::lme`. We can use the same `swcontainer` object and update the relevant lines. Notice that `nlme::lme` requires a different syntax for defining the model. Unsurprisingly, the results are the same Figure 1.

    models_lme <- seqwrap(swobject,
                      modelfun = nlme::lme,

                      arguments = alist(fixed = y ~ x,
                                        random = ~ 1|cluster),
                      return_models = TRUE,
                      eval_fun = NULL,
                      subset = 1:10)

    ── seqwrap ─────────────────────────────────────────────────────────────────────

    Initiating clusters for parallel processing with 1 core
    Fitting 10 targets in 4 chunks of up to 3 targets
    The progress bar updates 4 times (once per 1 chunk completed), in steps of
    about 25%; use a smaller `chunk_size` for more frequent updates
    Merging and modelling data
    sep 03 2026 12:49:49: Completed model fitting and evaluation. Elapsed time was
    0.01 minutes

    modelsum_lme <- seqwrap_summarise(models_lme)

    ── seqwrap summarise ───────────────────────────────────────────────────────────
    • A total of 50 samples and 10 targets where used in `The model function has
    been updated` with arguments `list(formula = y ~ x + (1 | cluster))`
    • Attempting to combine results from the provided summary function.
    • Attempting to combine results from the provided evaluations function.
    • Attempting to summarise errors and warnings from the fitting process.

    ── Model summaries ─────────────────────────────────────────────────────────────
    • 10 targets have associated summaries

    ── Model evaluations ───────────────────────────────────────────────────────────
    • 10 targets have associated evaluation results

    ── Errors and warnings ─────────────────────────────────────────────────────────
    ℹ No errors or warnings were raised

    ── Combined results ────────────────────────────────────────────────────────────
    ℹ Combined results have been generated and were
    silently returned.

    lme_est <- data.frame(modelsum_lme$summaries[, c(1,4,5, 6,9)])
    lmer_est <- data.frame(modelsum$summaries[, c(1,4,5, 6)])

Figure 1: Similar estimates between `lme4::lmer` and `nlme::lme` models

`seqwrap` will attempt to extract relevant metadata for use in modeling. This is done by extracting variable names from model formulations. In some cases, this can be difficult because we use additional functionality in some packages. Models constructed in `nlme::gls` are an example.

We can construct a similar model to `nlme::lme` using `nlme::CorCompSymm`, resulting in a model where the correlations between observations within clusters are uniform (see `?nlme::CorCompSymm`). To help `seqwrap` understand that we want to include `cluster` from the metadata in all fitting objects, we need to include it in `additional_vars`.

    models_gls <- seqwrap(swobject,
                      modelfun = nlme::gls,
                      arguments = alist(
                      model = y ~ x,
                      correlation = nlme::corCompSymm(form = ~ 1|cluster) ),
                      additional_vars = "cluster",
                      return_models = TRUE,
                      eval_fun = NULL,
                      subset = 1:10)

    ── seqwrap ─────────────────────────────────────────────────────────────────────

    Initiating clusters for parallel processing with 1 core
    Fitting 10 targets in 4 chunks of up to 3 targets
    The progress bar updates 4 times (once per 1 chunk completed), in steps of
    about 25%; use a smaller `chunk_size` for more frequent updates
    Merging and modelling data
    sep 03 2026 12:49:50: Completed model fitting and evaluation. Elapsed time was
    0.01 minutes

    modelsum_gls <- seqwrap_summarise(models_gls)

    ── seqwrap summarise ───────────────────────────────────────────────────────────
    • A total of 50 samples and 10 targets where used in `The model function has
    been updated` with arguments `list(formula = y ~ x + (1 | cluster))`
    • Attempting to combine results from the provided summary function.
    • Attempting to combine results from the provided evaluations function.
    • Attempting to summarise errors and warnings from the fitting process.

    ── Model summaries ─────────────────────────────────────────────────────────────
    • 10 targets have associated summaries

    ── Model evaluations ───────────────────────────────────────────────────────────
    • 10 targets have associated evaluation results

    ── Errors and warnings ─────────────────────────────────────────────────────────
    ℹ No errors or warnings were raised

    ── Combined results ────────────────────────────────────────────────────────────
    ℹ Combined results have been generated and were
    silently returned.

As evident from Figure 2 we can see that the two models will produce slightly different p-values.

    gls_est <- data.frame(modelsum_gls$summaries[, c(1,2,3,4,6)])

Figure 2: Similar estimates between `lme4::lmer` and `nlme::lme` models