AMR/R/pca.R

183 lines
7.8 KiB
R
Raw Normal View History

2020-03-07 21:48:21 +01:00
# ==================================================================== #
# TITLE #
2022-10-05 09:12:22 +02:00
# AMR: An R Package for Working with Antimicrobial Resistance Data #
2020-03-07 21:48:21 +01:00
# #
# SOURCE #
2020-07-08 14:48:06 +02:00
# https://github.com/msberends/AMR #
2020-03-07 21:48:21 +01:00
# #
2022-10-05 09:12:22 +02:00
# CITE AS #
# Berends MS, Luz CF, Friedrich AW, Sinha BNM, Albers CJ, Glasner C #
# (2022). AMR: An R Package for Working with Antimicrobial Resistance #
# Data. Journal of Statistical Software, 104(3), 1-31. #
# doi:10.18637/jss.v104.i03 #
# #
2022-12-27 15:16:15 +01:00
# Developed at the University of Groningen and the University Medical #
# Center Groningen in The Netherlands, in collaboration with many #
# colleagues from around the world, see our website. #
2020-03-07 21:48:21 +01:00
# #
# This R package is free software; you can freely use and distribute #
# it for both personal and commercial purposes under the terms of the #
# GNU General Public License version 2.0 (GNU GPL-2), as published by #
# the Free Software Foundation. #
# We created this package for both routine data analysis and academic #
# research and it was publicly released in the hope that it will be #
# useful, but it comes WITHOUT ANY WARRANTY OR LIABILITY. #
2020-10-08 11:16:03 +02:00
# #
# Visit our website for the full manual and a complete tutorial about #
# how to conduct AMR data analysis: https://msberends.github.io/AMR/ #
2020-03-07 21:48:21 +01:00
# ==================================================================== #
#' Principal Component Analysis (for AMR)
2022-08-28 10:31:50 +02:00
#'
2020-03-08 11:18:59 +01:00
#' Performs a principal component analysis (PCA) based on a data set with automatic determination for afterwards plotting the groups and labels, and automatic filtering on only suitable (i.e. non-empty and numeric) variables.
2021-05-12 18:15:03 +02:00
#' @param x a [data.frame] containing [numeric] columns
2020-05-16 13:05:47 +02:00
#' @param ... columns of `x` to be selected for PCA, can be unquoted since it supports quasiquotation.
2020-03-07 21:48:21 +01:00
#' @inheritParams stats::prcomp
2020-03-08 11:18:59 +01:00
#' @details The [pca()] function takes a [data.frame] as input and performs the actual PCA with the \R function [prcomp()].
2022-08-28 10:31:50 +02:00
#'
2021-05-12 18:15:03 +02:00
#' The result of the [pca()] function is a [prcomp] object, with an additional attribute `non_numeric_cols` which is a vector with the column names of all columns that do not contain [numeric] values. These are probably the groups and labels, and will be used by [ggplot_pca()].
2020-03-08 11:18:59 +01:00
#' @return An object of classes [pca] and [prcomp]
2020-03-08 09:12:11 +01:00
#' @importFrom stats prcomp
2020-03-07 21:48:21 +01:00
#' @export
2022-08-28 10:31:50 +02:00
#' @examples
#' # `example_isolates` is a data set available in the AMR package.
2020-03-07 21:48:21 +01:00
#' # See ?example_isolates.
#'
#' \donttest{
#' if (require("dplyr")) {
2022-08-28 10:31:50 +02:00
#' # calculate the resistance per group first
#' resistance_data <- example_isolates %>%
#' group_by(
#' order = mo_order(mo), # group on anything, like order
#' genus = mo_genus(mo)
#' ) %>% # and genus as we do here;
#' filter(n() >= 30) %>% # filter on only 30 results per group
2023-01-21 23:47:20 +01:00
#' summarise_if(is.sir, resistance) # then get resistance of all drugs
2022-08-28 10:31:50 +02:00
#'
2022-11-13 13:44:25 +01:00
#' # now conduct PCA for certain antimicrobial drugs
2022-08-28 10:31:50 +02:00
#' pca_result <- resistance_data %>%
#' pca(AMC, CXM, CTX, CAZ, GEN, TOB, TMP, SXT)
#'
#' pca_result
#' summary(pca_result)
2022-08-28 10:31:50 +02:00
#'
2022-08-21 16:37:20 +02:00
#' # old base R plotting method:
#' biplot(pca_result)
2022-08-21 16:37:20 +02:00
#' # new ggplot2 plotting method using this package:
#' if (require("ggplot2")) {
2022-11-05 12:06:40 +01:00
#' ggplot_pca(pca_result)
#'
2022-08-21 16:37:20 +02:00
#' ggplot_pca(pca_result) +
#' scale_colour_viridis_d() +
#' labs(title = "Title here")
#' }
#' }
2020-05-16 13:05:47 +02:00
#' }
2020-03-08 11:18:59 +01:00
pca <- function(x,
...,
retx = TRUE,
2022-08-28 10:31:50 +02:00
center = TRUE,
2020-03-08 11:18:59 +01:00
scale. = TRUE,
tol = NULL,
rank. = NULL) {
meet_criteria(x, allow_class = "data.frame")
meet_criteria(retx, allow_class = "logical", has_length = 1)
meet_criteria(center, allow_class = "logical", has_length = 1)
meet_criteria(scale., allow_class = "logical", has_length = 1)
meet_criteria(tol, allow_class = "numeric", has_length = 1, allow_NULL = TRUE)
meet_criteria(rank., allow_class = "numeric", has_length = 1, allow_NULL = TRUE)
2022-08-28 10:31:50 +02:00
2020-03-08 11:18:59 +01:00
# unset data.table, tibble, etc.
2020-03-07 21:48:21 +01:00
# also removes groups made by dplyr::group_by
x <- as.data.frame(x, stringsAsFactors = FALSE)
x.bak <- x
2022-08-28 10:31:50 +02:00
2020-05-16 13:05:47 +02:00
# defuse R expressions, this replaces rlang::enquos()
dots <- substitute(list(...))
if (length(dots) > 1) {
2020-03-07 21:48:21 +01:00
new_list <- list(0)
2020-05-16 13:05:47 +02:00
for (i in seq_len(length(dots) - 1)) {
new_list[[i]] <- tryCatch(eval(dots[[i + 1]], envir = x),
2022-08-28 10:31:50 +02:00
error = function(e) stop(e$message, call. = FALSE)
)
2020-03-07 21:48:21 +01:00
if (length(new_list[[i]]) == 1) {
2022-10-05 09:12:22 +02:00
if (is.character(new_list[[i]]) && new_list[[i]] %in% colnames(x)) {
# this is to support quoted variables: df %pm>% pca("mycol1", "mycol2")
2020-05-16 13:05:47 +02:00
new_list[[i]] <- x[, new_list[[i]]]
2020-03-07 21:48:21 +01:00
} else {
# remove item - it's an argument like `center`
2020-03-07 21:48:21 +01:00
new_list[[i]] <- NULL
}
}
}
2022-08-28 10:31:50 +02:00
2020-03-07 21:48:21 +01:00
x <- as.data.frame(new_list, stringsAsFactors = FALSE)
if (any(vapply(FUN.VALUE = logical(1), x, function(y) !is.numeric(y)))) {
warning_("in `pca()`: be sure to first calculate the resistance (or susceptibility) of variables with antimicrobial test results, since PCA works with numeric variables only. See Examples in ?pca.", call = FALSE)
2020-03-07 21:48:21 +01:00
}
2022-08-28 10:31:50 +02:00
2020-03-07 21:48:21 +01:00
# set column names
2020-05-16 13:05:47 +02:00
tryCatch(colnames(x) <- as.character(dots)[2:length(dots)],
2022-08-28 10:31:50 +02:00
error = function(e) warning("column names could not be set")
)
2021-05-18 11:29:31 +02:00
# keep only numeric columns
2022-08-27 20:49:37 +02:00
x <- x[, vapply(FUN.VALUE = logical(1), x, function(y) is.numeric(y)), drop = FALSE]
2020-03-07 21:48:21 +01:00
# bind the data set with the non-numeric columns
x <- cbind(x.bak[, vapply(FUN.VALUE = logical(1), x.bak, function(y) !is.numeric(y) & !all(is.na(y))), drop = FALSE], x)
2020-03-07 21:48:21 +01:00
}
2022-08-28 10:31:50 +02:00
x <- ungroup(x) # would otherwise select the grouping vars
2020-05-16 13:05:47 +02:00
x <- x[rowSums(is.na(x)) == 0, ] # remove columns containing NAs
2022-08-28 10:31:50 +02:00
2022-08-27 20:49:37 +02:00
pca_data <- x[, which(vapply(FUN.VALUE = logical(1), x, function(x) is.numeric(x))), drop = FALSE]
2022-08-28 10:31:50 +02:00
message_(
"Columns selected for PCA: ", vector_and(font_bold(colnames(pca_data), collapse = NULL), quotes = TRUE),
". Total observations available: ", nrow(pca_data), "."
)
2021-05-24 00:06:28 +02:00
if (getRversion() < "3.4.0") {
# stats::prcomp prior to 3.4.0 does not have the 'rank.' argument
pca_model <- prcomp(pca_data, retx = retx, center = center, scale. = scale., tol = tol)
} else {
pca_model <- prcomp(pca_data, retx = retx, center = center, scale. = scale., tol = tol, rank. = rank.)
}
groups <- x[, vapply(FUN.VALUE = logical(1), x, function(y) !is.numeric(y) & !all(is.na(y))), drop = FALSE]
rownames(groups) <- NULL
attr(pca_model, "non_numeric_cols") <- groups
2020-03-08 11:18:59 +01:00
class(pca_model) <- c("pca", class(pca_model))
pca_model
2020-03-07 21:48:21 +01:00
}
#' @method print pca
#' @export
#' @noRd
print.pca <- function(x, ...) {
a <- attributes(x)$non_numeric_cols
if (!is.null(a)) {
print_pca_group(a)
class(x) <- class(x)[class(x) != "pca"]
}
print(x, ...)
}
#' @method summary pca
#' @export
#' @noRd
summary.pca <- function(object, ...) {
a <- attributes(object)$non_numeric_cols
if (!is.null(a)) {
print_pca_group(a)
class(object) <- class(object)[class(object) != "pca"]
}
summary(object, ...)
}
print_pca_group <- function(a) {
grps <- sort(unique(a[, 1, drop = TRUE]))
cat("Groups (n=", length(grps), ", named as '", colnames(a)[1], "'):\n", sep = "")
print(grps)
cat("\n")
}