R/effectsizes.R

In bramburger/colvano: Provide Intervals for Volcano Plots

#' Get effectsizes from limma lmFit object
#'
#' This function calculates effect sizes for lmFit object, assuming a single independent variable.
#'
#' @param fit the lmFit object
#' @param data data.frame, matrix or similar containing only the data to be used to calculate effect sizes. (no protein names, etc)
#' @param protein.names a vector containing the names of the proteins. Should be the same length as the number of rows in \code{data}
get.effects.limma <- function (fit, data, protein.names) {
    rbindlist(lapply(seq_len(nrow(data)), function (i) {
        eff <- effsize::cohen.d(unlist(data[i, ]),
                                factor(-fit$design, ordered=FALSE),
                                na.rm=TRUE)
        data.table(protein = protein.names[i],
                   effsize = eff$estimate,
                   eff.lo = eff$conf.int[1],
                   eff.hi = eff$conf.int[2])
    }))
}


#' Get effectsizes from two groups
#'
#' This function calculates effect sizes for two NOT-PAIRED groups.
#'
#' @param groupA,groupB two vectors or matrices containing the values for the two groups.
#' @param protein.names a vector containing the names of the proteins. Should be the same length as the number of rows in \code{groupA,groupB}.
get.effects.groups <- function (groupA, groupB, protein.names = NULL) {
    ## needs sanity checks for inputs
    retval <- NULL
    if (is.null(dim(groupA)) && is.null(dim(groupB))) {
        eff <- effsize::cohen.d(c(groupA, groupB) ~
                                    factor(c(rep("A", length(groupA)),
                                             rep("B", length(groupB))),
                                           ordered=FALSE),
                         na.rm=TRUE)
        retval <- data.table(protein = ifelse(is.null(protein.names), 1,
                                              protein.names),
                             effsize = eff$estimate,
                             eff.lo = eff$conf.int[1],
                             eff.hi = eff$conf.int[2])
    } else if (nrow(groupA) == nrow(groupB)) {
        groupA <- as.matrix(groupA)
        groupB <- as.matrix(groupB)
        retval <-
            rbindlist(lapply(seq_along(protein.names), function (i) {
                eff <- effsize::cohen.d(unlist(c(groupA[i,], groupB[i, ])) ~
                                            factor(c(rep("A", length(groupA[i,])),
                                                     rep("B", length(groupB[i,]))),
                                                   ordered=FALSE),
                                        na.rm=TRUE)
                data.table(protein = protein.names[i],
                           effsize = eff$estimate,
                           eff.lo = eff$conf.int[1],
                           eff.hi = eff$conf.int[2])
            }))
    }
    retval
}