R/test_counts.R

In michbur/dpcR: Digital PCR Analysis

#' Test counts
#'
#' The test for comparing counts from two or more digital PCR experiments.
#'
#' @aliases test_counts
#' @param input object of class \code{\linkS4class{adpcr}} or \code{\linkS4class{dpcr}}
#' with "nm" type.
#' @param model may have one of following values: \code{binomial}, \code{poisson},
#' \code{prop}, \code{ratio}. See Details.
#' @param conf.level confidence level of the intervals and groups.
#' @details
#' \code{test_counts} incorporates two different approaches to models: GLM (General Linear
#' Model) and multiple pair-wise tests. The GLM fits counts data from different
#' digital PCR experiments using quasibinomial or quasipoisson \code{\link[stats]{family}}.
#' Comparisons between single experiments utilize Tukey's contrast and multiple t-tests
#' (as provided by function \code{\link{glht}}).
#'
#' In case of pair-wise tests, (\code{\link[rateratio.test]{rateratio.test}} or
#' \code{\link[stats]{prop.test}}) are used to compare all pairs of experiments. The
#' p-values are adjusted using the Benjamini & Hochberg method (\code{\link[stats]{p.adjust}}).
#' Furthermore, confidence intervals are simultaneous.
#'
#' @note Mean number of template molecules per partition and its confidence intervals will
#' vary depending on input.
#' @export
#' @seealso
#' Functions used by \code{test_counts}:
#' \itemize{
#' \item \code{\link[stats]{glm}},
#' \item \code{\link[multcomp]{glht}},
#' \item \code{\link[multcomp]{cld}},
#' \item \code{\link[stats]{prop.test}},
#' \item \code{\link[rateratio.test]{rateratio.test}}
#' }
#'
#' GUI presenting capabilities of the test: \code{\link{test_counts_gui}}.
#'
#' @return an object of class \code{\linkS4class{count_test}}.
#' @author Michal Burdukiewicz, Stefan Roediger, Piotr Sobczyk.
#' @references Bretz F, Hothorn T, Westfall P, \emph{Multiple comparisons using R}.
#' Boca Raton, Florida, USA: Chapman & Hall/CRC Press (2010).
#' @examples
#' # be warned, the examples of test_counts are time-consuming
#' \dontrun{
#' adpcr1 <- sim_adpcr(m = 10, n = 765, times = 1000, pos_sums = FALSE, n_panels = 3)
#' adpcr2 <- sim_adpcr(m = 60, n = 550, times = 1000, pos_sums = FALSE, n_panels = 3)
#' adpcr2 <- rename_dpcr(adpcr2, exper = "Experiment2")
#' adpcr3 <- sim_adpcr(m = 10, n = 600, times = 1000, pos_sums = FALSE, n_panels = 3)
#' adpcr3 <- rename_dpcr(adpcr3, exper = "Experiment3")
#'
#' # compare experiments using binomial regression
#' two_groups_bin <- test_counts(bind_dpcr(adpcr1, adpcr2), model = "binomial")
#' summary(two_groups_bin)
#' plot(two_groups_bin)
#' # plot aggregated results
#' plot(two_groups_bin, aggregate = TRUE)
#' # get coefficients
#' coef(two_groups_bin)
#'
#' # this time use Poisson regression
#' two_groups_pois <- test_counts(bind_dpcr(adpcr1, adpcr2), model = "poisson")
#' summary(two_groups_pois)
#' plot(two_groups_pois)
#'
#' # see how test behaves when results aren't significantly different
#' one_group <- test_counts(bind_dpcr(adpcr1, adpcr3))
#' summary(one_group)
#' plot(one_group)
#' }
#'
test_counts <- function(input, model = "ratio", conf.level = 0.95) {
  if (!(model %in% c("binomial", "poisson", "prop", "ratio"))) {
    stop("Must must have one of following values: 'binomial', 'poisson', 'ratio' or 'prop'")
  }

  if (model %in% c("prop", "ratio")) {
    test_function <- if (model == "prop") prop.test else rateratio.test

    # extract number of positive partitions
    positives <- if (slot(input, "type") == "tnp") {
      slot(input, ".Data")[1, ]
    } else {
      colSums(input > 0, na.rm = TRUE)
    }
    total <- slot(input, "n")
    # change sequence of rows to create output similar to glm
    test_ids <- combn(1L:length(total), 2)[c(2, 1), ]

    # for only two experiments
    if (!is.matrix(test_ids)) {
      test_ids <- as.matrix(test_ids)
    }
    # perform one-against-one multiple proportion tests
    all_combns <- apply(test_ids, 2, function(i) {
      test_function(positives[i], total[i])
    })

    # adjust p-values
    p_vals <- p.adjust(sapply(all_combns, function(i) {
      i[["p.value"]]
    }), method = "BH")

    #     #for only two experiments
    #     if(!is.matrix(p_vals))
    #       p_vals <- as.matrix(p_vals)

    # values of chi square statistic
    statistics <- sapply(all_combns, function(i) {
      i[["statistic"]]
    })

    # what should be in res, depends on the model
    test_res <- if (model == "prop") {
      matrix(c(statistics, p_vals),
        ncol = 2,
        dimnames = list(
          apply(matrix(names(positives)[test_ids],
            ncol = 2,
            byrow = TRUE
          ), 1, function(i) {
            paste(i, collapse = " - ")
          }),
          c("X_squared", "p_value")
        )
      )
    } else {
      matrix(p_vals,
        ncol = 1,
        dimnames = list(
          apply(matrix(names(positives)[test_ids],
            ncol = 2,
            byrow = TRUE
          ), 1, function(i) {
            paste(i, collapse = " - ")
          }),
          c("p_value")
        )
      )
    }


    # split data in groups
    # calculate confidence intervals using Sidak's unequality
    group_vals <- fl(binom.confint(positives,
      total,
      conf.level = (conf.level)^(1 / ncol(input)),
      "wilson"
    )[, 4L:6])

    # only unsignif in reality
    only_signif <- test_ids[, p_vals > 1 - conf.level]
    if (!is.matrix(only_signif)) {
      only_signif <- as.matrix(only_signif)
    }
    # every experiment belongs to different group when dim(only_signif)[2] == 0
    groups <- if (dim(only_signif)[2] == 0) {
      as.list(1L:length(total))
    } else {
      unique(lapply(1L:length(total), function(i) {
        sort(unique(as.vector(only_signif[, as.logical(colSums(only_signif == i))])))
      }))
    }

    # detect single experiment groups - they are empty before this
    groups_length <- sapply(groups, length)
    if (0 %in% groups_length) {
      groups <- groups[groups_length != 0]
      groups <- c(groups, as.list((1L:length(total))[!(1L:length(total) %in%
        unlist(groups))]))
    }

    group_redundancy <- lapply(1L:length(groups), function(i) {
      which(sapply(groups, function(j) {
        all(groups[[i]] %in% j)
      }))
    })

    # group_redundancy <- do.call(rbind, lapply(1L:length(groups), function(single_gr) {
    #   groups_in <- sapply(groups, function(other_gr)
    #     all(groups[[single_gr]] %in% other_gr))
    #   data.frame(gr = single_gr, other_gr = 1L:length(groups), is_in = groups_in, n_gr = sum(groups_in))
    # }))
    #
    # largest_membership <- lapply(1L:length(groups), function(single_gr) {
    #   gr_subdf <- group_redundancy[group_redundancy[["other_gr"]] == single_gr, ]
    #   only_in <- gr_subdf[gr_subdf[["is_in"]], ]
    #   only_in[only_in[["n_gr"]] == max(only_in[["n_gr"]]), "gr"]
    # })
    #
    # lapply(unique(unlist(largest_membership[lengths(largest_membership) > 1])),
    #        function(single_doubtful)
    #          sum(sapply(largest_membership, function(other_gr)
    #            single_doubtful %in% other_gr)))

    group_matrix <- sapply(1L:length(total), function(experiment) {
      sapply(groups[lengths(group_redundancy) == 1], function(single_group) experiment %in% single_group)
    })

    gr_order <- order(sapply(groups[lengths(group_redundancy) == 1], function(single_group) {
      mean(group_vals[single_group, 1])
    }))

    # all experiments in one group
    if (is.null(dim(group_matrix))) {
      group_matrix <- matrix(group_matrix, nrow = 1)
    }
    # name groups using the abc convention and at the same time reorder them along to value
    dimnames(group_matrix) <- list(
      (1L:sum(lengths(group_redundancy) == 1))[gr_order],
      names(positives)
    )

    group_coef <- data.frame(
      apply(group_matrix, 2, function(i) {
        paste(sort(as.numeric(names(i[which(i)]))), collapse = ".")
      }),
      group_vals
    )
    colnames(group_coef) <- c("group", "lambda", "lambda.low", "lambda.up")
  } else {
    if (slot(input, "type") == "tnp") {
      stop("GLM does not work with 'tnp' type")
    }

    # choose proper family
    if (model == "binomial") {
      if (!(slot(input, "type") %in% c("tp"))) {
        # binarize input which isn't already binary
        input <- binarize(input)
      }
      # family for model
      fam <- quasibinomial(link = "log")
      # function transforming coefficients of model to lambdas
      trans_fun <- function(x) fl(exp(x))
    }

    if (model == "poisson") {
      # family for model
      fam <- quasipoisson(link = "log")
      # function transforming coefficients of model to lambdas
      trans_fun <- function(x) exp(x)
    }

    # dpcr version of melt
    n_vector <- slot(input, "n")
    m_dpcr <- do.call(rbind, lapply(1L:length(n_vector), function(i) {
      vals <- input[1L:n_vector[i], i]
      data.frame(experiment = rep(colnames(input)[i], length(vals)), values = vals)
    }))


    # fit model
    fit <- glm(values ~ experiment + 0, data = m_dpcr, family = fam)

    # do multiple comparision
    multi_comp <- glht(fit, linfct = mcp(experiment = "Tukey"))

    coefs <- summary(fit)[["coefficients"]][, 1:2]
    lambdas <- trans_fun(matrix(c(
      coefs[, 1],
      coefs[, 1] - coefs[, 2],
      coefs[, 1] + coefs[, 2]
    ), ncol = 3))

    summ_mc <- suppressWarnings(summary(multi_comp))
    groups <- suppressWarnings(cld(multi_comp)[["mcletters"]][["LetterMatrix"]])
    if (is.matrix(groups)) {
      # more than one group
      groups_vector <- apply(groups, 1, function(i) {
        paste0(colnames(groups)[i], collapse = "")
      })
    } else {
      # only one group
      groups_vector <- rep("a", ncol(input))
      names(groups_vector) <- colnames(input)
    }

    group_coef <- data.frame(groups_vector, lambdas)
    colnames(group_coef) <- c("group", "lambda", "lambda.low", "lambda.up")
    rownames(group_coef) <- colnames(input)
    test_res <- cbind(
      t = summ_mc[["test"]][["tstat"]],
      p_value = summ_mc[["test"]][["pvalues"]]
    )
  }

  summ <- summary(input, print = FALSE)[["summary"]]

  group_coef <- cbind(group_coef, summ[summ[["method"]] == "bhat", c("experiment", "replicate", "k", "n")])

  new("count_test",
    group_coef = group_coef, test_res = test_res,
    model = model
  )
}