R/freq_pair_wtd_t_test.R

In philstraforelli/ggsigmark: An Implementation of Grammar of Graphics to Denote Statistically Significant Differences

#' Frequentist Test Comparisons Between Weighted Means
#'
#' @description Calculate weighted pairwise comparisons between group levels with corrections for
#' multiple testing
#'
#' @param x A numerical vector that is being compared across subgroups
#' @param subgroup A grouping vector or factor of the subgroups. It must be the same length as `x`.
#' @param weight A post-stratification weight vector. It must be the same length as `x`.
#' @param p.adjust.method The method for adjusting p-values.
#' @param ... Additional arguments passed on to \link[stats]{t.test}
#'
#' @details This is a wrapper for \link[stats]{pairwise.t.test} but allowing for means with post-stratification weights applied.
#'
#' @return Object of class "pairwise.htest"
#'
#' @seealso \link[stats]{pairwise.t.test}, \link[Hmisc]{wtd.mean}
#'
#' @importFrom stats p.adjust.methods
#'
#' @export
#'
#' @examples
#' #Without weights
#' Ozone <- airquality$Ozone
#' Month <- factor(airquality$Month, labels = month.abb[5:9])
#' freq_pair_wtd_t_test(Ozone, Month)
#' freq_pair_wtd_t_test(Ozone, Month, p.adj = "bonf")
#'
#' #With weights
#' Weights <- rnorm(153, mean = 1, sd = 0.3)
#' freq_pair_wtd_t_test(Ozone, Month, Weights)
#' freq_pair_wtd_t_test(Ozone, Month, Weights, p.adj = "bonf")

freq_pair_wtd_t_test <- function(x, subgroup, weight = NULL, p.adjust.method = p.adjust.methods, ...)
{
	if (is.null(weight)) {
		weight <- base::rep(1, length(x))
	}
	DNAME <- base::paste(base::deparse(base::substitute(x)), "and", base::deparse(base::substitute(subgroup)))
	x <- x[which(!is.na(subgroup))]
	subgroup <- subgroup[which(!is.na(subgroup))]
	weight <- weight[which(!is.na(subgroup))]
	subgroup <- base::factor(subgroup)
	p.adjust.method <- base::match.arg(p.adjust.method)

	xbar_df <- Hmisc::summarize(cbind(x, weight), Hmisc::llist(subgroup), function(y) Hmisc::wtd.mean(y[, 1], y[, 2]), stat.name = "wtd.mean") #Taken from the help page on Hmisc::wtd.stats
	xbar <- xbar_df$wtd.mean
	names(xbar) <- base::levels(subgroup)
	sd_df <- Hmisc::summarize(cbind(x, weight), Hmisc::llist(subgroup), function(y) sqrt(Hmisc::wtd.var(y[, 1], y[, 2])), stat.name = "wtd.sd") #Taken from the help page on Hmisc::wtd.stats
	sd <- sd_df$wtd.sd
	names(sd) <- base::levels(subgroup)
	n <- base::tapply(!is.na(x), subgroup, sum)
	degf <- n - 1
	total.degf <- base::sum(degf)
	pooled.sd <- base::sqrt(sum(sd^2 * degf) / total.degf)

	compare.levels <- function(i, j) {
		dif <- xbar[i] - xbar[j]
		se.dif <- pooled.sd * base::sqrt(1 / n[i] + 1 / n[j])
		t.val <- dif / se.dif
		2 * stats::pt(-base::abs(t.val), total.degf)
	}
	PVAL <- stats::pairwise.table(compare.levels, base::levels(subgroup), p.adjust.method)
	ans <- base::list(method = "Unpaired T test", data.name = DNAME, p.value = PVAL, p.adjust.method = p.adjust.method)
	class(ans) <- "pairwise.htest"
	return(ans)
}