R/sigma_model.R

In biospi/deamass: Statistical modelling for mass spectrometry omics

#' @import data.table
#' @include generics.R
setMethod("model", "sigma_block", function(object, input, chain = 1) {
  ctrl <- control(object)
  cat(paste0("[", Sys.time(), "]  SIGMA-", toupper(input), " block=", name(object), " chain=", chain, "/", ctrl@nchain, "\n"))

  # load metadata
  DT.index <- fst::read.fst(file.path(filepath(object), input, "input.index.fst"), as.data.table = T)

  nitt <- ctrl@nwarmup + (ctrl@nsample * ctrl@thin) / ctrl@nchain
  DT.priors <- priors(object, input = input, as.data.table = T)
  if (!is.null(DT.priors)) {
    DT.priors[, Assay := as.integer(Assay)]

    if (ctrl@eb.model == "") {
      DT.priors <- NULL
    } else if (ctrl@eb.model == "fit") {
      DT.priors[Effect == "Measurements"]$s <- DT.priors[Effect == "Measurements", s0]
      DT.priors[Effect == "Measurements"]$df <- DT.priors[Effect == "Measurements", df0]
      DT.priors[Effect == "Components"]$s <- DT.priors[Effect == "Components", s0]
      DT.priors[Effect == "Components"]$df <- DT.priors[Effect == "Components", df0]
    }
  }

  # create subdirs
  dir.create(file.path(filepath(object), input, "timings"), showWarnings = F)
  if ("summaries" %in% ctrl@keep) dir.create(file.path(filepath(object), input, "summaries"), showWarnings = F)

  unlink(file.path(filepath(object), input, "dist_*.measurement.stdevs.fst"))
  unlink(file.path(filepath(object), input, "dist_*.component.stdevs.fst"))
  unlink(file.path(filepath(object), input, "dist_*.assay.deviations.fst"))
  dir.create(file.path(filepath(object), input, "measurement.stdevs"), showWarnings = F)
  dir.create(file.path(filepath(object), input, "component.stdevs"), showWarnings = F)
  if (input == "model0") dir.create(file.path(filepath(object), input, "assay.deviations"), showWarnings = F)

  if (input != "model0" || "model0" %in% ctrl@keep) {
    unlink(file.path(filepath(object), input, "dist_*.group.quants.fst"))
    unlink(file.path(filepath(object), input, "dist_*.group.means.fst"))
    unlink(file.path(filepath(object), input, "dist_*.measurement.means.fst"))
    unlink(file.path(filepath(object), input, "dist_*.component.means.fst"))
    unlink(file.path(filepath(object), input, "dist_*.component.deviations.fst"))
    dir.create(file.path(filepath(object), input, "group.quants"), showWarnings = F)
    dir.create(file.path(filepath(object), input, "group.means"), showWarnings = F)
    dir.create(file.path(filepath(object), input, "measurement.means"), showWarnings = F)
    dir.create(file.path(filepath(object), input, "component.means"), showWarnings = F)
    dir.create(file.path(filepath(object), input, "component.deviations"), showWarnings = F)
  }

  cat(paste0("[", Sys.time(), "]   modelling ngroup=", nrow(DT.index), " nitt=", nitt, "...\n"))

  # run model
  DT.index <- merge(DT.index, groups(object, as.data.table = T)[, .(Group, pred.time)], by = "Group", sort = F)
  DT.index[, Group := as.integer(Group)]
  outputs <- rbindlists(parallel_lapply(split(DT.index, by = "Group", drop = T), function(item, object, input, chain, DT.priors) {
    ctrl <- control(object)
    nitt <- ctrl@nwarmup + (ctrl@nsample * ctrl@thin) / ctrl@nchain

    # load data
    group <- item[1, Group]
    DT <- fst::read.fst(file.path(filepath(object), input, item[1, file]), from = item[1, from], to = item[1, to], as.data.table = T)[, !"Group"]

    # prepare DT for MCMCglmm
    DT[, Component := factor(Component)]
    DT[, Measurement := interaction(Component, factor(Measurement), drop = T, lex.order = T)]
    DT[, Assay := factor(Assay)]
    nC <- nlevels(DT$Component)
    nM <- nlevels(DT$Measurement)
    nA <- nlevels(DT$Assay)

    # fixed effects
    fixed <- as.formula(paste(ifelse("Count1" %in% colnames(DT), "c(Count0, Count1)", "Count0"), "~ ", ifelse(nM == 1, "", "Measurement-1 +"), ifelse(nA == 1, " 1", " Assay")))

    # measurement rcov
    if (nM == 1 || ctrl@measurement.model == "single") {
      rcov <- as.formula("~units")
      if (is.null(DT.priors)) {
        prior.rcov <- list(V = 1, nu = 0.02)
      } else {
        prior.rcov <- list(V = (log(2) * DT.priors[Effect == "Measurements", s])^2, nu = DT.priors[Effect == "Measurements", df])
      }
    } else {
      rcov <- as.formula("~idh(Measurement):units")
      if (is.null(DT.priors)) {
        prior.rcov <- list(V = diag(nM), nu = 0.02)
      } else {
        prior.rcov <- list(V = (log(2) * DT.priors[Effect == "Measurements", s])^2 * diag(nM), nu = DT.priors[Effect == "Measurements", df])
      }
    }

    # component random effect
    if (ctrl@component.model == "") {
      random.component <- NULL
    } else if (ctrl@component.model == "single" || nC == 1) {
      random.component <- "Component:Assay"
      if (is.null(DT.priors)) {
        prior.component <- list(Component = list(V = 1, nu = 1, alpha.mu = 0, alpha.V = 25^2))
      } else {
        prior.component <- list(Component = list(V = (log(2) * DT.priors[Effect == "Components", s])^2, nu = DT.priors[Effect == "Components", df]))
      }
    } else {
      random.component <- "idh(Component):Assay"
      if (is.null(DT.priors)) {
        prior.component <- list(Component = list(V = diag(nC), nu = nC, alpha.mu = rep(0, nC), alpha.V = diag(25^2, nC)))
      } else {
        prior.component <- list(Component = list(V = (log(2) * DT.priors[Effect == "Components", s])^2 * diag(nC), nu = DT.priors[Effect == "Components", df]))
      }
    }

    # assay random effect
    if (ctrl@assay.model == "") {
      random.assay <- NULL
    } else if (ctrl@assay.model == "measurement") {
      for (l in levels(DT$Assay)) DT[, paste0("Assay", l) := ifelse(Assay == l, 1, 0)]
      random.assay <- paste(paste0("idh(Assay", levels(DT$Assay), "):Component:Measurement"), collapse = "+")
      if (is.null(DT.priors)) {
        prior.assay <- lapply(levels(DT$Assay), function(a) list(V = 1, nu = 1, alpha.mu = 0, alpha.V = 25^2))
      } else {
        DT.priors[, Assay := factor(Assay, levels = levels(DT$Assay))]
        prior.assay <- lapply(levels(DT$Assay), function(a) list(V = (log(2) * DT.priors[Assay == a, s])^2, nu = DT.priors[Assay == a, df]))
      }
      names(prior.assay) <- paste0("Assay", levels(DT$Assay))
    } else {
      for (l in levels(DT$Assay)) DT[, paste0("Assay", l) := ifelse(Assay == l, 1, 0)]
      random.assay <- paste(paste0("idh(Assay", levels(DT$Assay), "):Component"), collapse = "+")
      if (is.null(DT.priors)) {
        prior.assay <- lapply(levels(DT$Assay), function(a) list(V = 1, nu = 1, alpha.mu = 0, alpha.V = 25^2))
      } else {
        DT.priors[, Assay := factor(Assay, levels = levels(DT$Assay))]
        prior.assay <- lapply(levels(DT$Assay), function(a) list(V = (log(2) * DT.priors[Assay == a, s])^2, nu = DT.priors[Assay == a, df]))
      }
      names(prior.assay) <- paste0("Assay", levels(DT$Assay))
    }

    # merge prior
    prior <- list(R = prior.rcov)
    if (!is.null(random.component)) {
      if (!is.null(random.assay)) {
        random <- as.formula(paste("~", random.component, "+", random.assay))
        prior$G <- c(prior.component, prior.assay)
      } else {
        random <- as.formula(paste("~", random.component))
        prior$G <- prior.component
      }
    } else {
      random <- as.formula(paste("~", random.assay))
      prior$G <- prior.assay
    }

    # family
    if (ctrl@error.model == "l" || ctrl@error.model == "lognormal") {
      DT$Count0 <- log(DT$Count0)
      if("Count1" %in% colnames(DT)) {
        DT[, Count1 := log(Count1)]
        family <- "cengaussian"
       } else {
         family <- "gaussian"
       }
    } else {
      DT[, Count0 := ceiling(Count0)]
      DT[, Count1 := ceiling(Count1)]
      if("Count1" %in% colnames(DT)) {
        family <- "cenpoisson"
      } else {
        family <- "poisson"
       }
    }

    ### RUN MODEL
    output <- list()

    # MCMCglmm fails if only one row... add a second
    if (nrow(DT) == 1) {
      DT <- rbind(DT, DT)
      DT[2, Count0 := NA]
      if ("Count1" %in% colnames(DT)) DT[2, Count1 := NA]
    }

    # MCMCglmm can very rarely fail on a dataset, try again with slightly perturbed values
    fit.model <- NULL
    attempt <- 1
    if ("summaries" %in% ctrl@keep) output$DT.summaries <- as.character(Sys.time())
    while (is.null(fit.model) && attempt <= 10) {
      if (attempt != 1) {
        if (family == "poisson" || family == "cenpoisson") {
          DT$Count0 <- DT$Count0 + 1
          if ("Count1" %in% names(DT)) DT$Count1 <- DT$Count1 + 1
        } else {
          rn <- rnorm(length(DT$Count0), sd = 1e-5)
          DT$Count0 <- DT$Count0 + rn
          if ("Count1" %in% names(DT)) DT$Count1 <- DT$Count1 + rn
        }
      }

      set.seed(ctrl@random.seed + (group - 1) * ctrl@nchain + (chain - 1))
      try(output$DT.timings <- system.time(fit.model <- MCMCglmm::MCMCglmm(
        fixed, random, rcov, family, data = DT, prior = prior,
        nitt = nitt, burnin = ctrl@nwarmup, thin = ctrl@thin, pr = T, verbose = F, singular.ok = T
      )))

      attempt <- attempt + 1
    }
    if (is.null(fit.model)) stop(paste0("[", Sys.time(), "] ERROR: MCMCglmm failed more than 10 times"))
    output$DT.timings <- data.table(Group = group, chain = chain, as.data.table(t(as.matrix(output$DT.timings))))

    if ("summaries" %in% ctrl@keep) {
      options(max.print = 99999)
      output$DT.summaries <- data.table(Group = group, chain = chain, Summary = paste(c(output$DT.summaries, capture.output(print(summary(fit.model))), as.character(Sys.time())), collapse = "\n"))
    }

    # create emmeans ref grid
    class(fit.model) <- "MCMCglmm_seaMass"
    frg <- emmeans::ref_grid(fit.model, data = DT, rg.limit = 10000000)

    ### MODEL0: EXTRACT ASSAY DEVIATIONS TO CALCULATE PRIORS
    if (input == "model0") {
      if (ctrl@assay.model == "measurement") {
        output$DT.assay.deviations <- as.data.table(fit.model$Sol[, grep("^Assay.+\\.Component:Measurement\\..+\\..+$", colnames(fit.model$Sol)), drop = F])
        output$DT.assay.deviations[, sample := 1:nrow(output$DT.assay.deviations)]
        output$DT.assay.deviations <- melt(output$DT.assay.deviations, variable.name = "Measurement", id.vars = "sample")
        output$DT.assay.deviations[, Assay := as.integer(sub("^Assay(.+)\\.Component:Measurement\\..+\\..+$", "\\1", Measurement))]
        output$DT.assay.deviations[, Component := as.integer(sub("^Assay.+\\.Component:Measurement\\.(.+)\\..+$", "\\1", Measurement))]
        output$DT.assay.deviations[, Measurement := as.integer(sub("^Assay.+\\.Component:Measurement\\..+\\.(.+)$", "\\1", Measurement))]
        output$DT.assay.deviations[, Group := group]
        output$DT.assay.deviations[, chain := chain]
        output$DT.assay.deviations[, value := value / log(2)]
        setcolorder(output$DT.assay.deviations, c("Assay", "Group", "Component", "Measurement", "chain", "sample"))
      } else if (ctrl@assay.model == "component") {
        output$DT.assay.deviations <- as.data.table(fit.model$Sol[, grep("^Assay.+\\.Component\\..+$", colnames(fit.model$Sol)), drop = F])
        output$DT.assay.deviations[, sample := 1:nrow(output$DT.assay.deviations)]
        output$DT.assay.deviations <- melt(output$DT.assay.deviations, variable.name = "Component", id.vars = "sample")
        output$DT.assay.deviations[, Assay := as.integer(sub("^Assay(.+)\\.Component\\..+$", "\\1", Component))]
        output$DT.assay.deviations[, Component := as.integer(sub("^Assay.+\\.Component\\.(.+)$", "\\1", Component))]
        output$DT.assay.deviations[, Group := group]
        output$DT.assay.deviations[, chain := chain]
        output$DT.assay.deviations[, value := value / log(2)]
        setcolorder(output$DT.assay.deviations, c("Assay", "Group", "Component", "chain", "sample"))
      }
    }

    ### EXTRACT FIXED EFFECT OUTPUT
    if (input != "model0" || "model0" %in% ctrl@keep) {
      # measurement means
      if  ("measurements" %in% ctrl@summarise || "measurement.means" %in% ctrl@keep) {
        if (nM == 1) {
          if (nA == 1) {
            output$DT.measurement.means <- as.data.table(fit.model$Sol[, 1, drop = F])
          } else {
            output$DT.measurement.means <- as.data.table(coda::as.mcmc(emmeans::emmeans(frg, "1")))
          }
          colnames(output$DT.measurement.means) <- "value"
          output$DT.measurement.means[, Component := as.integer(as.character(DT[1, Component]))]
          output$DT.measurement.means[, Measurement := as.integer(sub("^.+\\.(.+)$", "\\1", as.character(DT[1, Measurement])))]
          output$DT.measurement.means[, sample := 1:nrow(output$DT.measurement.means)]
        } else {
          output$DT.measurement.means <- as.data.table(coda::as.mcmc(emmeans::emmeans(frg, "Measurement")))
          output$DT.measurement.means[, sample := 1:nrow(output$DT.measurement.means)]
          output$DT.measurement.means <- melt(output$DT.measurement.means, variable.name = "Measurement", id.vars = "sample")
          output$DT.measurement.means[, Component := as.integer(sub("^Measurement (.+)\\..+$", "\\1", Measurement))]
          output$DT.measurement.means[, Measurement := as.integer(sub("^Measurement .+\\.(.+)$", "\\1", Measurement))]
        }
        output$DT.measurement.means[, chain := chain]
        output$DT.measurement.means[, value := value / log(2)]
        output$DT.measurement.means[, Group := group]
        setcolorder(output$DT.measurement.means, c("Group", "Component", "Measurement", "chain", "sample"))
      }

      # component means
      if ("components" %in% ctrl@summarise || "component.means" %in% ctrl@keep) {
        if (nM == 1) {
          if (nA == 1) {
            output$DT.component.means <- as.data.table(fit.model$Sol[, 1, drop = F])
          } else {
            output$DT.component.means <- as.data.table(coda::as.mcmc(emmeans::emmeans(frg, "1")))
          }
          colnames(output$DT.component.means) <- "value"
          output$DT.component.means[, Component := as.integer(as.character(DT[1, Component]))]
          output$DT.component.means[, sample := 1:nrow(output$DT.component.means)]
        } else {
          # do each level separately otherwise huge memory problem
          newlevs <- sub("\\..+$", "", levels(DT$Measurement))
          output$DT.component.means <- as.data.table(sapply(levels(DT$Component), function (l) {
            frg2 <- emmeans::add_grouping(frg, "Component", "Measurement", ifelse(newlevs == l, newlevs, NA))
            return(coda::as.mcmc(emmeans::emmeans(frg2, "Component")))
          }))
          output$DT.component.means[, sample := 1:nrow(output$DT.component.means)]
          output$DT.component.means <- melt(output$DT.component.means, variable.name = "Component", id.vars = "sample")
          output$DT.component.means[, Component := as.integer(as.character(Component))]
        }
        output$DT.component.means[, chain := chain]
        output$DT.component.means[, value := value / log(2)]
        output$DT.component.means[, Group := group]
        setcolorder(output$DT.component.means, c("Group", "Component", "chain", "sample"))
      }

      # component deviations
      if (ctrl@component.model == "independent" && ("components" %in% ctrl@summarise || "component.deviations" %in% ctrl@keep)) {
        if (nC == 1) {
          output$DT.component.deviations <- as.data.table(fit.model$Sol[, grep("^Component:Assay\\..+\\..+$", colnames(fit.model$Sol)), drop = F])
          output$DT.component.deviations[, sample := 1:nrow(output$DT.component.deviations)]
          output$DT.component.deviations <- melt(output$DT.component.deviations, variable.name = "Component", id.vars = "sample")
          output$DT.component.deviations[, Assay := as.integer(sub("^Component:Assay\\..+\\.(.+)$", "\\1", Component))]
          output$DT.component.deviations[, Component := as.integer(sub("^Component:Assay\\.(.+)\\..+$", "\\1", Component))]
        } else {
          output$DT.component.deviations <- as.data.table(fit.model$Sol[, grep("^Component.+\\.Assay\\..+$", colnames(fit.model$Sol)), drop = F])
          output$DT.component.deviations[, sample := 1:nrow(output$DT.component.deviations)]
          output$DT.component.deviations <- melt(output$DT.component.deviations, variable.name = "Component", id.vars = "sample")
          output$DT.component.deviations[, Assay := as.integer(sub("^Component.+\\.Assay\\.(.+)$", "\\1", Component))]
          output$DT.component.deviations[, Component := as.integer(sub("^Component(.+)\\.Assay\\..+$", "\\1", Component))]
        }
        output$DT.component.deviations[, Group := group]
        output$DT.component.deviations[, chain := chain]
        output$DT.component.deviations[, value := value / log(2)]
        setcolorder(output$DT.component.deviations, c("Group", "Component", "Assay", "chain", "sample"))
      }

      # group means
      if ("groups" %in% ctrl@summarise || "group.means" %in% ctrl@keep) {
        output$DT.group.means <- as.data.table(coda::as.mcmc(emmeans::emmeans(frg, "1")))
        colnames(output$DT.group.means) <- "value"
        output$DT.group.means[, sample := 1:nrow(output$DT.group.means)]
        output$DT.group.means[, chain := chain]
        output$DT.group.means[, value := value / log(2)]
        output$DT.group.means[, Group := group]
        setcolorder(output$DT.group.means, c("Group", "chain", "sample"))
      }

      # group quants
      if ("groups" %in% ctrl@summarise || "group.quants" %in% ctrl@keep) {
        if (nA == 1) {
          if (nM == 1) {
            output$DT.group.quants <- as.data.table(fit.model$Sol[, 1, drop = F])
          } else {
            output$DT.group.quants <- as.data.table(coda::as.mcmc(emmeans::emmeans(frg, "1")))
          }
          colnames(output$DT.group.quants) <- "value"
          output$DT.group.quants[, Assay := as.integer(DT[1, Assay])]
          output$DT.group.quants[, sample := 1:nrow(output$DT.group.quants)]
        } else {
          output$DT.group.quants <- as.data.table(coda::as.mcmc(emmeans::emmeans(frg, "Assay")))
          output$DT.group.quants[, sample := 1:nrow(output$DT.group.quants)]
          output$DT.group.quants <- melt(output$DT.group.quants, variable.name = "Assay", id.vars = "sample")
          output$DT.group.quants[, Assay := as.integer(sub("^Assay (.+)$", "\\1", Assay))]
        }
        output$DT.group.quants[, Group := group]
        output$DT.group.quants[, chain := chain]
        output$DT.group.quants[, value := value / log(2)]
        setcolorder(output$DT.group.quants, c("Group", "Assay", "chain", "sample"))
      }
    }

    # DELETE FIXED EFFECTS TO CONSERVE MEMORY
    fit.model$Sol <- NULL

    ### EXTRACT RANDOM EFFECT OUTPUT
    # measurement stdevs
    if (input == "model0" || "measurements" %in% ctrl@summarise || "measurement.stdevs" %in% ctrl@keep) {
      if (ctrl@measurement.model == "single" || nM == 1) {
        output$DT.measurement.stdevs <- as.data.table(fit.model$VCV[, "units", drop = F])
      } else {
        output$DT.measurement.stdevs <- as.data.table(fit.model$VCV[, grep("^Measurement.+\\..+\\.units$", colnames(fit.model$VCV)), drop = F])
      }
      output$DT.measurement.stdevs[, sample := 1:nrow(output$DT.measurement.stdevs)]
      output$DT.measurement.stdevs <- melt(output$DT.measurement.stdevs, id.vars = "sample")

      # component and measurement
      if (ctrl@measurement.model == "single") {
        output$DT.measurement.stdevs[, Component := group]
        output$DT.measurement.stdevs[, Measurement := group]
      } else if (nM == 1) {
        output$DT.measurement.stdevs[, Component := as.integer(as.character(DT[1, Component]))]
        output$DT.measurement.stdevs[, Measurement := as.integer(sub("^.+\\.(.+)$", "\\1", as.character(DT[1, Measurement])))]
      } else {
        output$DT.measurement.stdevs[, Component := as.integer(sub("^Measurement(.+)\\..+\\.units", "\\1", variable))]
        output$DT.measurement.stdevs[, Measurement := as.integer(sub("^Measurement.+\\.(.+)\\.units", "\\1", variable))]
      }

      # rest
      output$DT.measurement.stdevs[, Group := group]
      output$DT.measurement.stdevs[, chain := chain]
      output$DT.measurement.stdevs[, value := sqrt(value) / log(2)]
      output$DT.measurement.stdevs[, variable := NULL]
      setcolorder(output$DT.measurement.stdevs, c("Group", "Component", "Measurement", "chain", "sample"))
    }

    # component stdevs
    if (ctrl@component.model != "" && (input == "model0" || "components" %in% ctrl@summarise || "component.stdevs" %in% ctrl@keep)) {
      if (ctrl@component.model == "single" || nC == 1) {
        output$DT.component.stdevs <- as.data.table(fit.model$VCV[, "Component:Assay", drop = F])
      } else {
        output$DT.component.stdevs <- as.data.table(fit.model$VCV[, grep("^Component.+\\.Assay$", colnames(fit.model$VCV)), drop = F])
      }
      output$DT.component.stdevs[, sample := 1:nrow(output$DT.component.stdevs)]
      output$DT.component.stdevs <- melt(output$DT.component.stdevs, id.vars = "sample")

      # component
      if (ctrl@component.model == "single") {
        output$DT.component.stdevs[, Component := group]
      } else if (nC == 1) {
        output$DT.component.stdevs[, Component := as.integer(as.character(DT[1, Component]))]
      } else {
        output$DT.component.stdevs[, Component := as.integer(sub("^Component(.+)\\.Assay$", "\\1", variable))]
      }

      # rest
      output$DT.component.stdevs[, Group := group]
      output$DT.component.stdevs[, chain := chain]
      output$DT.component.stdevs[, value := sqrt(value) / log(2)]
      output$DT.component.stdevs[, variable := NULL]
      setcolorder(output$DT.component.stdevs, c("Group", "Component", "chain", "sample"))
    }

    ### WRITE OUT
    # if large enough write out group quants now to conserve memory, otherwise don't to conserve disk space
    if (object.size(output$DT.group.quants) > 2^18) {
      filename <- file.path("group.quants", paste0(chain, ".", group, ".fst"))
      fst::write.fst(output$DT.group.quants, file.path(filepath(object), input, filename))

      if (chain == 1) {
        # construct index
        output$DT.index.group.quants <- output$DT.group.quants[, .(
          from = .I[!duplicated(output$DT.group.quants, by = c("Group", "Assay"))],
          to = .I[!duplicated(output$DT.group.quants, fromLast = T, by = c("Group", "Assay"))]
        )]
        output$DT.index.group.quants <- cbind(
          output$DT.group.quants[output$DT.index.group.quants$from, .(Group, Assay)],
          data.table(file = factor(filename)),
          output$DT.index.group.quants
        )
      }

      output$DT.group.quants <- data.table()
    } else {
      if (chain == 1) output$DT.index.group.quants <- data.table()
    }

    # if large enough write out measurement means now to conserve memory, otherwise don't to conserve disk space
    if (object.size(output$DT.measurement.means) > 2^18) {
      filename <- file.path("measurement.means", paste0(chain, ".", group, ".fst"))
      fst::write.fst(output$DT.measurement.means, file.path(filepath(object), input, filename))

      if (chain == 1) {
        # construct index
        output$DT.index.measurement.means <- output$DT.measurement.means[, .(
          from = .I[!duplicated(output$DT.measurement.means, by = c("Group", "Component", "Measurement"))],
          to = .I[!duplicated(output$DT.measurement.means, fromLast = T, by = c("Group", "Component", "Measurement"))]
        )]
        output$DT.index.measurement.means <- cbind(
          output$DT.measurement.means[output$DT.index.measurement.means$from, .(Group, Component, Measurement)],
          data.table(file = factor(filename)),
          output$DT.index.measurement.means
        )
      }

      output$DT.measurement.means <- data.table()
    } else {
      if (chain == 1) output$DT.index.measurement.means <- data.table()
    }

    # if large enough write out component means now to conserve memory, otherwise don't to conserve disk space
    if (object.size(output$DT.component.means) > 2^18) {
      filename <- file.path("component.means", paste0(chain, ".", group, ".fst"))
      fst::write.fst(output$DT.component.means, file.path(filepath(object), input, filename))

      if (chain == 1) {
        # construct index
        output$DT.index.component.means <- output$DT.component.means[, .(
          from = .I[!duplicated(output$DT.component.means, by = c("Group", "Component"))],
          to = .I[!duplicated(output$DT.component.means, fromLast = T, by = c("Group", "Component"))]
        )]
        output$DT.index.component.means <- cbind(
          output$DT.component.means[output$DT.index.component.means$from, .(Group, Component)],
          data.table(file = factor(filename)),
          output$DT.index.component.means
        )
      }

      output$DT.component.means <- data.table()
    } else {
      if (chain == 1) output$DT.index.component.means <- data.table()
    }

    # if large enough write out group means now to conserve memory, otherwise don't to conserve disk space
    if (object.size(output$DT.group.means) > 2^18) {
      filename <- file.path("group.means", paste0(chain, ".", group, ".fst"))
      fst::write.fst(output$DT.group.means, file.path(filepath(object), input, filename))

      if (chain == 1) {
        # construct index
        output$DT.index.group.means <- output$DT.group.means[, .(
          from = .I[!duplicated(output$DT.group.means, by = "Group")],
          to = .I[!duplicated(output$DT.group.means, fromLast = T, by = "Group")]
        )]
        output$DT.index.group.means <- cbind(
          output$DT.group.means[output$DT.index.group.means$from, .(Group)],
          data.table(file = factor(filename)),
          output$DT.index.group.means
        )
      }

      output$DT.group.means <- data.table()
    } else {
      if (chain == 1) output$DT.index.group.means <- data.table()
    }

    # if large enough write out component deviations now to conserve memory, otherwise don't to conserve disk space
    if ("DT.component.deviations" %in% names(output)) {
      if (object.size(output$DT.component.deviations) > 2^18) {
        filename <- file.path("component.deviations", paste0(chain, ".", group, ".fst"))
        fst::write.fst(output$DT.component.deviations, file.path(filepath(object), input, filename))

        if (chain == 1) {
          # construct index
          output$DT.index.component.deviations <- output$DT.component.deviations[, .(
            from = .I[!duplicated(output$DT.component.deviations, by = c("Group", "Component", "Assay"))],
            to = .I[!duplicated(output$DT.component.deviations, fromLast = T, by = c("Group", "Component", "Assay"))]
          )]
          output$DT.index.component.deviations <- cbind(
            output$DT.component.deviations[output$DT.index.component.deviations$from, .(Group, Component, Assay)],
            data.table(file = factor(filename)),
            output$DT.index.component.deviations
          )
        }

        output$DT.component.deviations <- data.table()
      } else {
        if (chain == 1) output$DT.index.component.deviations <- data.table()
      }
    }

    # if large enough write out assay deviations now to conserve memory, otherwise don't to conserve disk space
    if ("DT.assay.deviations" %in% names(output)) {
      if (object.size(output$DT.assay.deviations) > 2^18) {
        filename <- file.path("assay.deviations", paste0(chain, ".", group, ".fst"))
        fst::write.fst(output$DT.assay.deviations, file.path(filepath(object), input, filename))

        if (chain == 1) {
          # construct index
          if (ctrl@assay.model == "measurement") {
            output$DT.index.assay.deviations <- output$DT.assay.deviations[, .(
              from = .I[!duplicated(output$DT.assay.deviations, by = c("Assay", "Group", "Component", "Measurement"))],
              to = .I[!duplicated(output$DT.assay.deviations, fromLast = T, by = c("Assay", "Group", "Component", "Measurement"))]
            )]
            output$DT.index.assay.deviations <- cbind(
              output$DT.assay.deviations[output$DT.index.assay.deviations$from, .(Assay, Group, Component, Measurement)],
              data.table(file = factor(filename)),
              output$DT.index.assay.deviations
            )
          } else {
            output$DT.index.assay.deviations <- output$DT.assay.deviations[, .(
              from = .I[!duplicated(output$DT.assay.deviations, by = c("Assay", "Group", "Component"))],
              to = .I[!duplicated(output$DT.assay.deviations, fromLast = T, by = c("Assay", "Group", "Component"))]
            )]
            output$DT.index.assay.deviations <- cbind(
              output$DT.assay.deviations[output$DT.index.assay.deviations$from, .(Assay, Group, Component)],
              data.table(file = factor(filename)),
              output$DT.index.assay.deviations
            )
          }
        }

        output$DT.assay.deviations <- data.table()
      } else {
        if (chain == 1) output$DT.index.assay.deviations <- data.table()
      }
    }

    # if large enough write out measurement stdevs now to conserve memory, otherwise don't to conserve disk space
    if (object.size(output$DT.measurement.stdevs) > 2^18) {
      filename <- file.path("measurement.stdevs", paste0(chain, ".", group, ".fst"))
      fst::write.fst(output$DT.measurement.stdevs, file.path(filepath(object), input, filename))

      if (chain == 1) {
        # construct index
        output$DT.index.measurement.stdevs <- output$DT.measurement.stdevs[, .(
          from = .I[!duplicated(output$DT.measurement.stdevs, by = c("Group", "Component", "Measurement"))],
          to = .I[!duplicated(output$DT.measurement.stdevs, fromLast = T, by = c("Group", "Component", "Measurement"))]
        )]
        output$DT.index.measurement.stdevs <- cbind(
          output$DT.measurement.stdevs[output$DT.index.measurement.stdevs$from, .(Group, Component, Measurement)],
          data.table(file = factor(filename)),
          output$DT.index.measurement.stdevs
        )
      }

      output$DT.measurement.stdevs <- data.table()
    } else {
      if (chain == 1) output$DT.index.measurement.stdevs <- data.table()
    }

    # if large enough write out component stdevs now to conserve memory, otherwise don't to conserve disk space
    if ("DT.component.stdevs" %in% names(output)) {
      if (object.size(output$DT.component.stdevs) > 2^18) {
        filename <- file.path("component.stdevs", paste0(chain, ".", group, ".fst"))
        fst::write.fst(output$DT.component.stdevs, file.path(filepath(object), input, filename))

        if (chain == 1) {
          # construct index
          output$DT.index.component.stdevs <- output$DT.component.stdevs[, .(
            from = .I[!duplicated(output$DT.component.stdevs, by = c("Group", "Component"))],
            to = .I[!duplicated(output$DT.component.stdevs, fromLast = T, by = c("Group", "Component"))]
          )]
          output$DT.index.component.stdevs <- cbind(
            output$DT.component.stdevs[output$DT.index.component.stdevs$from, .(Group, Component)],
            data.table(file = factor(filename)),
            output$DT.index.component.stdevs
          )
        }

        output$DT.component.stdevs <- data.table()
      } else {
        if (chain == 1) output$DT.index.component.stdevs <- data.table()
      }
    }

    return(output)
  }, nthread = ctrl@nthread, pred = DT.index[, pred.time]))

  # write out
  DT.groups <- groups(object, as.data.table = T)
  DT.components <- components(object, as.data.table = T)
  DT.measurements <- measurements(object, as.data.table = T)
  DT.design <- assay_design(object, as.data.table = T)

  # write out summaries
  if ("summaries" %in% ctrl@keep) {
    outputs$DT.summaries[, Group := factor(Group, levels = 1:nlevels(DT.groups[, Group]), labels = levels(DT.groups[, Group]))]
    setkey(outputs$DT.summaries, Group, chain)
    fst::write.fst(outputs$DT.summaries, file.path(filepath(object), input, file.path("summaries", paste0(chain, ".fst"))))
    outputs$DT.summaries <- NULL
  }

  # write out timings
  outputs$DT.timing[, Group := factor(Group, levels = 1:nlevels(DT.groups[, Group]), labels = levels(DT.groups[, Group]))]
  setkey(outputs$DT.timings, Group, chain)
  fst::write.fst(outputs$DT.timings, file.path(filepath(object), input, file.path("timings", paste0(chain, ".fst"))))
  outputs$DT.timings <- NULL

  # write out group quants
  if ("DT.group.quants" %in% names(outputs)) {
    if (nrow(outputs$DT.group.quants) > 0) {
      setorder(outputs$DT.group.quants, Group, Assay, chain, sample)
      filename <- file.path("group.quants", paste0(chain, ".fst"))
      fst::write.fst(outputs$DT.group.quants, file.path(filepath(object), input, filename))
      # finish index construction
      if (chain == 1) {
        DT.index.group.quants <- outputs$DT.group.quants[, .(
          from = .I[!duplicated(outputs$DT.group.quants, by = c("Group", "Assay"))],
          to = .I[!duplicated(outputs$DT.group.quants, fromLast = T, by = c("Group", "Assay"))]
        )]
        outputs$DT.index.group.quants <- rbind(outputs$DT.index.group.quants, cbind(
          outputs$DT.group.quants[DT.index.group.quants$from, .(Group, Assay)],
          data.table(file = factor(filename)),
          DT.index.group.quants
        ))
      }
      outputs$DT.group.quants <- NULL
    }
    # write index
    if (chain == 1) {
      outputs$DT.index.group.quants[, Group := factor(Group, levels = 1:nlevels(DT.groups[, Group]), labels = levels(DT.groups[, Group]))]
      outputs$DT.index.group.quants[, Assay := factor(Assay, levels = 1:nlevels(DT.design[, Assay]), labels = levels(DT.design[, Assay]))]
      setkey(outputs$DT.index.group.quants, Group, file, from)
      fst::write.fst(outputs$DT.index.group.quants, file.path(filepath(object), input, paste0("group.quants.index.fst")))
    }
  }

  # write out measurement means
  if ("DT.measurement.means" %in% names(outputs)) {
    if (nrow(outputs$DT.measurement.means) > 0) {
      setorder(outputs$DT.measurement.means, Group, Component, Measurement, chain, sample)
      filename <- file.path("measurement.means", paste0(chain, ".fst"))
      fst::write.fst(outputs$DT.measurement.means, file.path(filepath(object), input, filename))
      # finish index construction
      if (chain == 1) {
        DT.index.measurement.means <- outputs$DT.measurement.means[, .(
          from = .I[!duplicated(outputs$DT.measurement.means, by = c("Group", "Component", "Measurement"))],
          to = .I[!duplicated(outputs$DT.measurement.means, fromLast = T, by = c("Group", "Component", "Measurement"))]
        )]
        outputs$DT.index.measurement.means <- rbind(outputs$DT.index.measurement.means, cbind(
          outputs$DT.measurement.means[DT.index.measurement.means$from, .(Group, Component, Measurement)],
          data.table(file = factor(filename)),
          DT.index.measurement.means
        ))
      }
      outputs$DT.measurement.means <- NULL
    }
    # write index
    if (chain == 1) {
      outputs$DT.index.measurement.means[, Group := factor(Group, levels = 1:nlevels(DT.groups[, Group]), labels = levels(DT.groups[, Group]))]
      outputs$DT.index.measurement.means[, Component := factor(Component, levels = 1:nlevels(DT.components[, Component]), labels = levels(DT.components[, Component]))]
      outputs$DT.index.measurement.means[, Measurement := factor(Measurement, levels = 1:nlevels(DT.measurements[, Measurement]), labels = levels(DT.measurements[, Measurement]))]
      setkey(outputs$DT.index.measurement.means, Group, file, from)
      fst::write.fst(outputs$DT.index.measurement.means, file.path(filepath(object), input, paste0("measurement.means.index.fst")))
    }
  }

  # write out components means
  if ("DT.component.means" %in% names(outputs)) {
    if (nrow(outputs$DT.component.means) > 0) {
      setorder(outputs$DT.component.means, Group, Component, chain, sample)
      filename <- file.path("component.means", paste0(chain, ".fst"))
      fst::write.fst(outputs$DT.component.means, file.path(filepath(object), input, filename))
      # finish index construction
      if (chain == 1) {
        DT.index.component.means <- outputs$DT.component.means[, .(
          from = .I[!duplicated(outputs$DT.component.means, by = c("Group", "Component"))],
          to = .I[!duplicated(outputs$DT.component.means, fromLast = T, by = c("Group", "Component"))]
        )]
        outputs$DT.index.component.means <- rbind(outputs$DT.index.component.means, cbind(
          outputs$DT.component.means[DT.index.component.means$from, .(Group, Component)],
          data.table(file = factor(filename)),
          DT.index.component.means
        ))
      }
      outputs$DT.component.means <- NULL
    }
    # write index
    if (chain == 1) {
      outputs$DT.index.component.means[, Group := factor(Group, levels = 1:nlevels(DT.groups[, Group]), labels = levels(DT.groups[, Group]))]
      outputs$DT.index.component.means[, Component := factor(Component, levels = 1:nlevels(DT.components[, Component]), labels = levels(DT.components[, Component]))]
       setkey(outputs$DT.index.component.means, Group, file, from)
      fst::write.fst(outputs$DT.index.component.means, file.path(filepath(object), input, paste0("component.means.index.fst")))
    }
  }

  # write out group means
  if ("DT.group.means" %in% names(outputs)) {
    if (nrow(outputs$DT.group.means) > 0) {
      setorder(outputs$DT.group.means, Group, chain, sample)
      filename <- file.path("group.means", paste0(chain, ".fst"))
      fst::write.fst(outputs$DT.group.means, file.path(filepath(object), input, filename))
      # finish index construction
      if (chain == 1) {
        DT.index.group.means <- outputs$DT.group.means[, .(
          from = .I[!duplicated(outputs$DT.group.means, by = "Group")],
          to = .I[!duplicated(outputs$DT.group.means, fromLast = T, by = "Group")]
        )]
        outputs$DT.index.group.means <- rbind(outputs$DT.index.group.means, cbind(
          outputs$DT.group.means[DT.index.group.means$from, .(Group)],
          data.table(file = factor(filename)),
          DT.index.group.means
        ))
      }
      outputs$DT.group.means <- NULL
    }
    # write index
    if (chain == 1) {
      outputs$DT.index.group.means[, Group := factor(Group, levels = 1:nlevels(DT.groups[, Group]), labels = levels(DT.groups[, Group]))]
      setkey(outputs$DT.index.group.means, Group, file, from)
      fst::write.fst(outputs$DT.index.group.means, file.path(filepath(object), input, paste0("group.means.index.fst")))
    }
  }

  # write out component deviations
  if ("DT.component.deviations" %in% names(outputs)) {
    if (nrow(outputs$DT.component.deviations) > 0) {
      setorder(outputs$DT.component.deviations, Group, Component, Assay, chain, sample)
      filename <- file.path("component.deviations", paste0(chain, ".fst"))
      fst::write.fst(outputs$DT.component.deviations, file.path(filepath(object), input, filename))
      # finish index construction
      if (chain == 1) {
        DT.index.component.deviations <- outputs$DT.component.deviations[, .(
          from = .I[!duplicated(outputs$DT.component.deviations, by = c("Group", "Component", "Assay"))],
          to = .I[!duplicated(outputs$DT.component.deviations, fromLast = T, by = c("Group", "Component", "Assay"))]
        )]
        outputs$DT.index.component.deviations <- rbind(outputs$DT.index.component.deviations, cbind(
          outputs$DT.component.deviations[DT.index.component.deviations$from, .(Group, Component, Assay)],
          data.table(file = factor(filename)),
          DT.index.component.deviations
        ))
      }
      outputs$DT.component.deviations <- NULL
    }
    # write index
    if (chain == 1) {
      outputs$DT.index.component.deviations[, Group := factor(Group, levels = 1:nlevels(DT.groups[, Group]), labels = levels(DT.groups[, Group]))]
      outputs$DT.index.component.deviations[, Component := factor(Component, levels = 1:nlevels(DT.components[, Component]), labels = levels(DT.components[, Component]))]
      outputs$DT.index.component.deviations[, Assay := factor(Assay, levels = 1:nlevels(DT.design[, Assay]), labels = levels(DT.design[, Assay]))]
      setkey(outputs$DT.index.component.deviations, Group, file, from)
      fst::write.fst(outputs$DT.index.component.deviations, file.path(filepath(object), input, paste0("component.deviations.index.fst")))
    }
  }

  # write out assay deviations
  if ("DT.assay.deviations" %in% names(outputs)) {
    if (nrow(outputs$DT.assay.deviations) > 0) {
      if (ctrl@assay.model == "measurement") {
        setorder(outputs$DT.assay.deviations, Assay, Group, Component, Measurement, chain, sample)
        filename <- file.path("assay.deviations", paste0(chain, ".fst"))
        fst::write.fst(outputs$DT.assay.deviations, file.path(filepath(object), input, filename))
        # finish index construction
        if (chain == 1) {
          DT.index.assay.deviations <- outputs$DT.assay.deviations[, .(
            from = .I[!duplicated(outputs$DT.assay.deviations, by = c("Assay", "Group", "Component", "Measurement"))],
            to = .I[!duplicated(outputs$DT.assay.deviations, fromLast = T, by = c("Assay", "Group", "Component", "Measurement"))]
          )]
          outputs$DT.index.assay.deviations <- rbind(outputs$DT.index.assay.deviations, cbind(
            outputs$DT.assay.deviations[DT.index.assay.deviations$from, .(Assay, Group, Component, Measurement)],
            data.table(file = factor(filename)),
            DT.index.assay.deviations
          ))
        }
       } else {
        if (nrow(outputs$DT.assay.deviations) > 0) {
          setorder(outputs$DT.assay.deviations, Assay, Group, Component, chain, sample)
          filename <- file.path("assay.deviations", paste0(chain, ".fst"))
          fst::write.fst(outputs$DT.assay.deviations, file.path(filepath(object), input, filename))
          # finish index construction
          if (chain == 1) {
            DT.index.assay.deviations <- outputs$DT.assay.deviations[, .(
              from = .I[!duplicated(outputs$DT.assay.deviations, by = c("Assay", "Group", "Component"))],
              to = .I[!duplicated(outputs$DT.assay.deviations, fromLast = T, by = c("Assay", "Group", "Component"))]
            )]
            outputs$DT.index.assay.deviations <- rbind(outputs$DT.index.assay.deviations, cbind(
              outputs$DT.assay.deviations[DT.index.assay.deviations$from, .(Assay, Group, Component)],
              data.table(file = factor(filename)),
              DT.index.assay.deviations
            ))
          }
        }
      }
      outputs$DT.assay.deviations <- NULL
    }
    # write index
    if (chain == 1) {
      #print("writing index")
      if (ctrl@assay.model == "measurement") {
        outputs$DT.index.assay.deviations[, Assay := factor(Assay, levels = 1:nlevels(DT.design[, Assay]), labels = levels(DT.design[, Assay]))]
        outputs$DT.index.assay.deviations[, Group := factor(Group, levels = 1:nlevels(DT.groups[, Group]), labels = levels(DT.groups[, Group]))]
        outputs$DT.index.assay.deviations[, Component := factor(Component, levels = 1:nlevels(DT.components[, Component]), labels = levels(DT.components[, Component]))]
        outputs$DT.index.assay.deviations[, Measurement := factor(Measurement, levels = 1:nlevels(DT.measurements[, Measurement]), labels = levels(DT.measurements[, Measurement]))]
        setkey(outputs$DT.index.assay.deviations, Assay, file, from)
        fst::write.fst(outputs$DT.index.assay.deviations, file.path(filepath(object), input, paste0("assay.deviations.index.fst")))
      } else {
        outputs$DT.index.assay.deviations[, Assay := factor(Assay, levels = 1:nlevels(DT.design[, Assay]), labels = levels(DT.design[, Assay]))]
        outputs$DT.index.assay.deviations[, Group := factor(Group, levels = 1:nlevels(DT.groups[, Group]), labels = levels(DT.groups[, Group]))]
        outputs$DT.index.assay.deviations[, Component := factor(Component, levels = 1:nlevels(DT.components[, Component]), labels = levels(DT.components[, Component]))]
        setkey(outputs$DT.index.assay.deviations, Assay, file, from)
        fst::write.fst(outputs$DT.index.assay.deviations, file.path(filepath(object), input, paste0("assay.deviations.index.fst")))
      }
    }
  }

  # write out measurement stdevs
  if ("DT.measurement.stdevs" %in% names(outputs)) {
    if (nrow(outputs$DT.measurement.stdevs) > 0) {
      # write out remaining measurement stdevs
      if(ctrl@measurement.model == "independent") {
        setorder(outputs$DT.measurement.stdevs, Group, Component, Measurement, chain, sample)
        filename <- file.path("measurement.stdevs", paste0(chain, ".fst"))
        fst::write.fst(outputs$DT.measurement.stdevs, file.path(filepath(object), input, filename))
        # finish index construction
        if (chain == 1) {
          DT.index.measurement.stdevs <- outputs$DT.measurement.stdevs[, .(
            from = .I[!duplicated(outputs$DT.measurement.stdevs, by = c("Group", "Component", "Measurement"))],
            to = .I[!duplicated(outputs$DT.measurement.stdevs, fromLast = T, by = c("Group", "Component", "Measurement"))]
          )]
          outputs$DT.index.measurement.stdevs <- rbind(outputs$DT.index.measurement.stdevs, cbind(
            outputs$DT.measurement.stdevs[DT.index.measurement.stdevs$from, .(Group, Component, Measurement)],
            data.table(file = factor(filename)),
            DT.index.measurement.stdevs
          ))
        }
      } else {
        setorder(outputs$DT.measurement.stdevs, Group, chain, sample)
        filename <- file.path("measurement.stdevs", paste0(chain, ".fst"))
        fst::write.fst(outputs$DT.measurement.stdevs, file.path(filepath(object), input, filename))
        # finish index construction
        if (chain == 1) {
          DT.index.measurement.stdevs <- outputs$DT.measurement.stdevs[, .(
            from = .I[!duplicated(outputs$DT.measurement.stdevs, by = c("Group"))],
            to = .I[!duplicated(outputs$DT.measurement.stdevs, fromLast = T, by = c("Group"))]
          )]
          outputs$DT.index.measurement.stdevs <- rbind(outputs$DT.index.measurement.stdevs, cbind(
            outputs$DT.measurement.stdevs[DT.index.measurement.stdevs$from, .(Group)],
            data.table(file = factor(filename)),
            DT.index.measurement.stdevs
          ))
        }
      }
      outputs$DT.measurement.stdevs <- NULL
    }
    # write index
    if (chain == 1) {
      if(ctrl@measurement.model == "independent") {
        outputs$DT.index.measurement.stdevs[, Group := factor(Group, levels = 1:nlevels(DT.groups[, Group]), labels = levels(DT.groups[, Group]))]
        outputs$DT.index.measurement.stdevs[, Component := factor(Component, levels = 1:nlevels(DT.components[, Component]), labels = levels(DT.components[, Component]))]
        outputs$DT.index.measurement.stdevs[, Measurement := factor(Measurement, levels = 1:nlevels(DT.measurements[, Measurement]), labels = levels(DT.measurements[, Measurement]))]
        setkey(outputs$DT.index.measurement.stdevs, Group, file, from)
        fst::write.fst(outputs$DT.index.measurement.stdevs, file.path(filepath(object), input, paste0("measurement.stdevs.index.fst")))
      } else {
        outputs$DT.index.measurement.stdevs[, Group := factor(Group, levels = 1:nlevels(DT.groups[, Group]), labels = levels(DT.groups[, Group]))]
        setkey(outputs$DT.index.measurement.stdevs, Group, file, from)
        fst::write.fst(outputs$DT.index.measurement.stdevs, file.path(filepath(object), input, paste0("measurement.stdevs.index.fst")))
      }
    }
  }

  # write out component stdevs
  if ("DT.component.stdevs" %in% names(outputs)) {
    if (nrow(outputs$DT.component.stdevs) > 0) {
      setorder(outputs$DT.component.stdevs, Group, Component, chain, sample)
      filename <- file.path("component.stdevs", paste0(chain, ".fst"))
      fst::write.fst(outputs$DT.component.stdevs, file.path(filepath(object), input, filename))
      # finish index construction
      if (chain == 1) {
        DT.index.component.stdevs <- outputs$DT.component.stdevs[, .(
          from = .I[!duplicated(outputs$DT.component.stdevs, by = c("Group", "Component"))],
          to = .I[!duplicated(outputs$DT.component.stdevs, fromLast = T, by = c("Group", "Component"))]
        )]
        outputs$DT.index.component.stdevs <- rbind(outputs$DT.index.component.stdevs, cbind(
          outputs$DT.component.stdevs[DT.index.component.stdevs$from, .(Group, Component)],
          data.table(file = factor(filename)),
          DT.index.component.stdevs
        ))
      }
      outputs$DT.component.stdevs <- NULL
    }
    # write index
    if (chain == 1) {
      outputs$DT.index.component.stdevs[, Group := factor(Group, levels = 1:nlevels(DT.groups[, Group]), labels = levels(DT.groups[, Group]))]
      outputs$DT.index.component.stdevs[, Component := factor(Component, levels = 1:nlevels(DT.components[, Component]), labels = levels(DT.components[, Component]))]
      setkey(outputs$DT.index.component.stdevs, Group, file, from)
      fst::write.fst(outputs$DT.index.component.stdevs, file.path(filepath(object), input, paste0("component.stdevs.index.fst")))
    }
  }

  return(invisible(NULL))
})