Nothing
R/cf_correction.R
In DACF: Data Analysis with Ceiling and/or Floor Data
Defines functions
cf_correction
Documented in
cf_correction
#' cf_correction
#' @importFrom utils combn
#' @description Extends Tukey's HSD and Bonferroni procedure to account for ceiling and floor effects
#' @param x a dataframe of data with ceiling/floor effects and corresponding group variables or aov results
#' @param tests a character string specifying the desired multiple-comparison procedure: "all" (default, Tukey's HSD and Bonferroni), "tukey" (Tukey's HSD), or "bonf" (Bonferroni)
#' @param df.adjustment a character string specifying the desired method for adjusting the degree of freedom: "trunc" (default, Liu and Wang's truncated-normal corrections) or "unadj" (unadjusted)
#' @param gh.correction a character string specifying if the Welch & Games-Howell correction for heteroscedasticity should be included: "no_gh"(default) or "yes_gh"
#' @param alpha a (non-empty) numeric value of desired of alpha level
#' @param flr a (non-empty) numeric value of the floor threshold (e.g., the minimum score of the measurement scale)
#' @param ceil a (non-empty) numeric value of the ceiling threshold (e.g., the maximum score of the measurement scale)
#' @return a matrix containing pairwise comparison results. Columns depend on the
#' \code{tests} and \code{gh.correction} arguments:
#' \itemize{
#' \item \code{Comparison_i}: First group in the comparison
#' \item \code{Comparison_j}: Second group in the comparison
#' \item \code{mean_i}: The adjusted mean of the first group
#' \item \code{mean_j}: The adjusted mean of the second group
#' \item \code{diff_in_means}: The difference in adjusted means
#' \item \code{tukey.CI_lwr}: Lower bound of Tukey's HSD confidence interval (tests = "tukey" or "all")
#' \item \code{tukey.CI_upr}: Upper bound of Tukey's HSD confidence interval (tests = "tukey" or "all")
#' \item \code{hedges_g}: Hedges' g effect size (tests = "tukey" or "all")
#' \item \code{Q}: Adjusted Tukey's Q statistic (tests = "tukey" or "all")
#' \item \code{p}: p-value for Tukey's Q statistic (tests = "tukey" or "all")
#' \item \code{t}: t-test statistic (tests = "bonf" or "all")
#' \item \code{t_p}: p-value for the t-test statistic (tests = "bonf" or "all")
#' \item \code{p.bonferroni}: Bonferroni-adjusted p-value (tests = "bonf" or "all")
#' \item \code{bonf.CI_lwr}: Lower bound of Bonferroni confidence interval (tests = "bonf" or "all")
#' \item \code{bonf.CI_upr}: Upper bound of Bonferroni confidence interval (tests = "bonf" or "all")
#' \item \code{gh.CI_lwr}: Lower bound of Games-Howell confidence interval (gh.correction = "yes_gh", tests = "tukey" or "all")
#' \item \code{gh.CI_upr}: Upper bound of Games-Howell confidence interval (gh.correction = "yes_gh", tests = "tukey" or "all")
#' \item \code{games.howell.t}: Games-Howell t statistic (gh.correction = "yes_gh", tests = "tukey" or "all")
#' \item \code{games.howell.p}: Games-Howell p-value (gh.correction = "yes_gh", tests = "tukey" or "all")
#' \item \code{games.howell.g}: Games-Howell Hedges' g (gh.correction = "yes_gh", tests = "tukey" or "all")
#' \item \code{welch.satterthwaite.df}: Welch-Satterthwaite degrees of freedom (gh.correction = "yes_gh")
#' \item \code{welch.bonf.CI_lwr}: Lower bound of Welch-Bonferroni confidence interval (gh.correction = "yes_gh", tests = "bonf" or "all")
#' \item \code{welch.bonf.CI_upr}: Upper bound of Welch-Bonferroni confidence interval (gh.correction = "yes_gh", tests = "bonf" or "all")
#' \item \code{t_Welch}: Welch t statistic (gh.correction = "yes_gh", tests = "bonf" or "all")
#' \item \code{p_Welch}: Welch p-value (gh.correction = "yes_gh", tests = "bonf" or "all")
#' \item \code{p_Welch.bonferroni}: Welch Bonferroni-adjusted p-value (gh.correction = "yes_gh", tests = "bonf" or "all")
#' }
#'
#' @references
#'
#' Aitkin MA. Correlation in a Singly Truncated Bivariate Normal Distribution.
#' \emph{Psychometrika}. 1964;29(3):263-270.
#' \doi{10.1007/BF02289723}
#'
#' Cohen, A. C. (1959). Simplified Estimators for the Normal Distribution When Samples Are Singly Censored or Truncated.
#' \emph{AnnalsTechnometrics}, 1(3), 217–237.
#' \doi{10.1080/00401706.1959.10489859}
#'
#' Greene, W. H. (2002). Econometric Analysis.
#' \emph{In Econometric Analysis}.
#'
#' Liu, Q., Wang, L. (2020) t-Test and ANOVA for data with ceiling and/or floor effects.
#' \emph{Behav Res} 53, 264–277 .
#' \doi{10.3758/s13428-020-01407-2}
#'
#' Qi, H., Dong, Z., Wei, Q., Chen, X., & Luo, Y. (2025). Are Gamers Happier? Multidimensional Well-Being Differences in Risk Groups for Problematic Internet Use and Internet Gaming Disorder.
#' \emph{Personality and Individual Differences} 246, 113334.
#' \doi{10.2139/ssrn.5148047}
#'
#' @examples
#' data("Qi_data_males") # (Qi et al., 2025)
#' cf_correction(Qi_data_males, tests = "all", df.adjustment = "trunc",
#' gh.correction = "yes_gh", alpha = .05, flr = 0, ceil = 9)
#' cf_correction(Qi_data_males, tests = "tukey", df.adjustment = "unadj",
#' gh.correction = "no_gh", alpha = .05, flr = 0, ceil = 9)
#' cf_correction(Qi_data_males, tests = "all", df.adjustment = "unadj",
#' gh.correction = "no_gh", alpha = .05, flr = 0, ceil = 9)
#' @export cf_correction
cf_correction = function(x, tests="all", df.adjustment="trunc", gh.correction="no_gh", alpha=.05, flr, ceil) {
if (inherits(x, "aov")) {
input <- "aov"
} else if (is.data.frame(x) || is.matrix(x)) {
input <- "data"
} else {
stop("Input must be a data frame or an aov object.")
}
if (input == "data") {
data <- as.data.frame(x)
colnames(data) <- c("group_id", "value")
} else if (input == "aov") {
aov_obj <- x
data <- model.frame(aov_obj)
colnames(data) <- c("value", "group_id")
data <- data[, c("group_id", "value")]
}
if (flr >= ceil) stop("flr must be less than ceil.")
data <- data[order(data$group_id), ]
ns <- table(data$group_id)
kgroups <- length(ns)
unique_ids <- unique(data$group_id)
comparisons <- combn(unique_ids, 2)
rownames(comparisons) <- c("i", "j")
num_comparisons <- dim(comparisons)[2]
floor_counts <- tapply(data$value, data$group_id, function(x) sum(x == flr))
ceil_counts <- tapply(data$value, data$group_id, function(x) sum(x == ceil))
group_means_adj <- c()
group_var_adj <- c()
for(g in 1:kgroups){
y <- data[data$group_id == unique_ids[g], 2]
adj_vals <- rec.mean.var(y, flr=flr, ceil=ceil)
group_means_adj <- c(group_means_adj, adj_vals$est.mean)
group_var_adj <- c(group_var_adj, adj_vals$est.var)
}
names(group_means_adj) <- unique_ids
names(group_var_adj) <- unique_ids
group_sd_adj <- sqrt(group_var_adj)
SS_within_adj <- sum((ns - 1) * group_var_adj)
df_within_adj <- dim(data)[1] - kgroups
if (df.adjustment == "trunc") {
df <- sum(tapply(data$value, data$group_id,
function(x) length(x) - (sum(x == flr) - 1) - (sum(x == ceil) - 1)))
} else if (df.adjustment == "unadj") {
df <- sum(tapply(data$value, data$group_id, function(x) length(x) - 1))
}
MSE_adj <- SS_within_adj / df
results <- c()
for(c in 1:num_comparisons){
grp_i <- unname(comparisons[1, c])
grp_j <- unname(comparisons[2, c])
mean_i <- unname(group_means_adj[grp_i])
mean_j <- unname(group_means_adj[grp_j])
diff_in_means <- unname(mean_i - mean_j)
abs_diff <- unname(abs(mean_i - mean_j))
n_i <- unname(ns[grp_i])
n_j <- unname(ns[grp_j])
var_i <- unname(group_var_adj[grp_i])
var_j <- unname(group_var_adj[grp_j])
sqrt_term <- sqrt((MSE_adj/2) * (1/n_i + 1/n_j))
rmse <- sqrt(MSE_adj)
# Tukey
ci_term <- (qtukey(.95, kgroups, df=df) / sqrt(2)) * sqrt(MSE_adj * (1/n_i + 1/n_j))
lwr <- unname(diff_in_means - ci_term)
upr <- unname(diff_in_means + ci_term)
hedges_g <- unname(diff_in_means / rmse)
Q <- unname(abs_diff / sqrt_term)
p <- unname(ptukey(q=Q, nmeans=kgroups, df=df, lower.tail=FALSE))
# Bonferroni
t <- unname(diff_in_means / sqrt(MSE_adj * (1/n_i + 1/n_j)))
t_p <- unname(2 * pt(-abs(t), df, lower.tail=TRUE))
t_p_bonferroni <- unname(pmin(1, num_comparisons * t_p))
adjusted_alpha <- alpha / num_comparisons
ci_term_bonf <- qt(1 - adjusted_alpha/2, df) * sqrt(MSE_adj * (1/n_i + 1/n_j))
lwr_bonf <- unname(diff_in_means - ci_term_bonf)
upr_bonf <- unname(diff_in_means + ci_term_bonf)
# Games-Howell / Welch
vi <- unname(var_i / n_i)
vj <- unname(var_j / n_j)
df_i <- unname(n_i - floor_counts[grp_i] - ceil_counts[grp_i] - 1)
df_j <- unname(n_j - floor_counts[grp_j] - ceil_counts[grp_j] - 1)
se_ij <- unname(sqrt(vi + vj))
t_ij <- unname(diff_in_means / se_ij)
df_ij <- unname((vi + vj)^2 / (vi^2 / df_i + vj^2 / df_j))
df_ij <- max(df_ij, 1)
p_gh <- unname(ptukey(abs(t_ij) * sqrt(2), kgroups, df_ij, lower.tail=FALSE))
se_gh <- unname(sqrt(0.5 * (var_i / n_i + var_j / n_j)))
lwr_gh <- unname(diff_in_means - qtukey(p=0.95, nmeans=kgroups, df=df_ij) * se_gh)
upr_gh <- unname(diff_in_means + qtukey(p=0.95, nmeans=kgroups, df=df_ij) * se_gh)
hedges_g_gh <- unname(diff_in_means / sqrt(mean(group_var_adj)))
t_p_ij <- unname(2 * pt(-abs(t_ij), df_ij, lower.tail=TRUE))
t_p_bonf_ij <- unname(pmin(1, num_comparisons * t_p_ij))
ci_term_gh <- unname(qt(1 - (alpha / num_comparisons) / 2, df_ij) * se_ij)
lwr_gh_bonf <- unname(diff_in_means - ci_term_gh)
upr_gh_bonf <- unname(diff_in_means + ci_term_gh)
if (gh.correction == "no_gh") {
if (tests == "all") {
results <- rbind(results, c(
Comparison_i = grp_i,
Comparison_j = grp_j,
mean_i = mean_i,
mean_j = mean_j,
diff_in_means = diff_in_means,
tukey.CI_lwr = lwr,
tukey.CI_upr = upr,
hedges_g = hedges_g,
Q = Q,
p = p,
t = t,
t_p = t_p,
p.bonferroni = t_p_bonferroni,
bonf.CI_lwr = lwr_bonf,
bonf.CI_upr = upr_bonf
))
} else if (tests == "tukey") {
results <- rbind(results, c(
Comparison_i = grp_i,
Comparison_j = grp_j,
mean_i = mean_i,
mean_j = mean_j,
diff_in_means = diff_in_means,
tukey.CI_lwr = lwr,
tukey.CI_upr = upr,
hedges_g = hedges_g,
Q = Q,
p = p
))
} else if (tests == "bonf") {
results <- rbind(results, c(
Comparison_i = grp_i,
Comparison_j = grp_j,
mean_i = mean_i,
mean_j = mean_j,
diff_in_means = diff_in_means,
t = t,
t_p = t_p,
p.bonferroni = t_p_bonferroni,
bonf.CI_lwr = lwr_bonf,
bonf.CI_upr = upr_bonf
))
}
} else {
if (tests == "all") {
results <- rbind(results, c(
Comparison_i = grp_i,
Comparison_j = grp_j,
mean_i = mean_i,
mean_j = mean_j,
diff_in_means = diff_in_means,
tukey.CI_lwr = lwr,
tukey.CI_upr = upr,
hedges_g = hedges_g,
Q = Q,
p = p,
t = t,
t_p = t_p,
p.bonferroni = t_p_bonferroni,
bonf.CI_lwr = lwr_bonf,
bonf.CI_upr = upr_bonf,
gh.CI_lwr = lwr_gh,
gh.CI_upr = upr_gh,
games.howell.t = t_ij,
games.howell.p = p_gh,
games.howell.g = hedges_g_gh,
welch.satterthwaite.df = df_ij,
welch.bonf.CI_lwr = lwr_gh_bonf,
welch.bonf.CI_upr = upr_gh_bonf,
t_Welch = t_ij,
p_Welch = t_p_ij,
p_Welch.bonferroni = t_p_bonf_ij
))
} else if (tests == "tukey") {
results <- rbind(results, c(
Comparison_i = grp_i,
Comparison_j = grp_j,
mean_i = mean_i,
mean_j = mean_j,
diff_in_means = diff_in_means,
tukey.CI_lwr = lwr,
tukey.CI_upr = upr,
hedges_g = hedges_g,
Q = Q,
p = p,
gh.CI_lwr = lwr_gh,
gh.CI_upr = upr_gh,
games.howell.t = t_ij,
games.howell.p = p_gh,
games.howell.g = hedges_g_gh,
welch.satterthwaite.df = df_ij
))
} else if (tests == "bonf") {
results <- rbind(results, c(
Comparison_i = grp_i,
Comparison_j = grp_j,
mean_i = mean_i,
mean_j = mean_j,
diff_in_means = diff_in_means,
t = t,
t_p = t_p,
p.bonferroni = t_p_bonferroni,
bonf.CI_lwr = lwr_bonf,
bonf.CI_upr = upr_bonf,
welch.bonf.CI_lwr = lwr_gh_bonf,
welch.bonf.CI_upr = upr_gh_bonf,
t_Welch = t_ij,
p_Welch = t_p_ij,
p_Welch.bonferroni = t_p_bonf_ij
))
}
}
} # end comparisons loop
return(results)
}
Try the DACF package in your browser
Any scripts or data that you put into this service are public.
DACF documentation built on July 17, 2026, 1:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions cf_correction
Documented in cf_correction
#' cf_correction
#' @importFrom utils combn
#' @description Extends Tukey's HSD and Bonferroni procedure to account for ceiling and floor effects
#' @param x a dataframe of data with ceiling/floor effects and corresponding group variables or aov results
#' @param tests a character string specifying the desired multiple-comparison procedure: "all" (default, Tukey's HSD and Bonferroni), "tukey" (Tukey's HSD), or "bonf" (Bonferroni)
#' @param df.adjustment a character string specifying the desired method for adjusting the degree of freedom: "trunc" (default, Liu and Wang's truncated-normal corrections) or "unadj" (unadjusted)
#' @param gh.correction a character string specifying if the Welch & Games-Howell correction for heteroscedasticity should be included: "no_gh"(default) or "yes_gh"
#' @param alpha a (non-empty) numeric value of desired of alpha level
#' @param flr a (non-empty) numeric value of the floor threshold (e.g., the minimum score of the measurement scale)
#' @param ceil a (non-empty) numeric value of the ceiling threshold (e.g., the maximum score of the measurement scale)
#' @return a matrix containing pairwise comparison results. Columns depend on the
#' \code{tests} and \code{gh.correction} arguments:
#' \itemize{
#' \item \code{Comparison_i}: First group in the comparison
#' \item \code{Comparison_j}: Second group in the comparison
#' \item \code{mean_i}: The adjusted mean of the first group
#' \item \code{mean_j}: The adjusted mean of the second group
#' \item \code{diff_in_means}: The difference in adjusted means
#' \item \code{tukey.CI_lwr}: Lower bound of Tukey's HSD confidence interval (tests = "tukey" or "all")
#' \item \code{tukey.CI_upr}: Upper bound of Tukey's HSD confidence interval (tests = "tukey" or "all")
#' \item \code{hedges_g}: Hedges' g effect size (tests = "tukey" or "all")
#' \item \code{Q}: Adjusted Tukey's Q statistic (tests = "tukey" or "all")
#' \item \code{p}: p-value for Tukey's Q statistic (tests = "tukey" or "all")
#' \item \code{t}: t-test statistic (tests = "bonf" or "all")
#' \item \code{t_p}: p-value for the t-test statistic (tests = "bonf" or "all")
#' \item \code{p.bonferroni}: Bonferroni-adjusted p-value (tests = "bonf" or "all")
#' \item \code{bonf.CI_lwr}: Lower bound of Bonferroni confidence interval (tests = "bonf" or "all")
#' \item \code{bonf.CI_upr}: Upper bound of Bonferroni confidence interval (tests = "bonf" or "all")
#' \item \code{gh.CI_lwr}: Lower bound of Games-Howell confidence interval (gh.correction = "yes_gh", tests = "tukey" or "all")
#' \item \code{gh.CI_upr}: Upper bound of Games-Howell confidence interval (gh.correction = "yes_gh", tests = "tukey" or "all")
#' \item \code{games.howell.t}: Games-Howell t statistic (gh.correction = "yes_gh", tests = "tukey" or "all")
#' \item \code{games.howell.p}: Games-Howell p-value (gh.correction = "yes_gh", tests = "tukey" or "all")
#' \item \code{games.howell.g}: Games-Howell Hedges' g (gh.correction = "yes_gh", tests = "tukey" or "all")
#' \item \code{welch.satterthwaite.df}: Welch-Satterthwaite degrees of freedom (gh.correction = "yes_gh")
#' \item \code{welch.bonf.CI_lwr}: Lower bound of Welch-Bonferroni confidence interval (gh.correction = "yes_gh", tests = "bonf" or "all")
#' \item \code{welch.bonf.CI_upr}: Upper bound of Welch-Bonferroni confidence interval (gh.correction = "yes_gh", tests = "bonf" or "all")
#' \item \code{t_Welch}: Welch t statistic (gh.correction = "yes_gh", tests = "bonf" or "all")
#' \item \code{p_Welch}: Welch p-value (gh.correction = "yes_gh", tests = "bonf" or "all")
#' \item \code{p_Welch.bonferroni}: Welch Bonferroni-adjusted p-value (gh.correction = "yes_gh", tests = "bonf" or "all")
#' }
#'
#' @references
#'
#' Aitkin MA. Correlation in a Singly Truncated Bivariate Normal Distribution.
#' \emph{Psychometrika}. 1964;29(3):263-270.
#' \doi{10.1007/BF02289723}
#'
#' Cohen, A. C. (1959). Simplified Estimators for the Normal Distribution When Samples Are Singly Censored or Truncated.
#' \emph{AnnalsTechnometrics}, 1(3), 217–237.
#' \doi{10.1080/00401706.1959.10489859}
#'
#' Greene, W. H. (2002). Econometric Analysis.
#' \emph{In Econometric Analysis}.
#'
#' Liu, Q., Wang, L. (2020) t-Test and ANOVA for data with ceiling and/or floor effects.
#' \emph{Behav Res} 53, 264–277 .
#' \doi{10.3758/s13428-020-01407-2}
#'
#' Qi, H., Dong, Z., Wei, Q., Chen, X., & Luo, Y. (2025). Are Gamers Happier? Multidimensional Well-Being Differences in Risk Groups for Problematic Internet Use and Internet Gaming Disorder.
#' \emph{Personality and Individual Differences} 246, 113334.
#' \doi{10.2139/ssrn.5148047}
#'
#' @examples
#' data("Qi_data_males") # (Qi et al., 2025)
#' cf_correction(Qi_data_males, tests = "all", df.adjustment = "trunc",
#' gh.correction = "yes_gh", alpha = .05, flr = 0, ceil = 9)
#' cf_correction(Qi_data_males, tests = "tukey", df.adjustment = "unadj",
#' gh.correction = "no_gh", alpha = .05, flr = 0, ceil = 9)
#' cf_correction(Qi_data_males, tests = "all", df.adjustment = "unadj",
#' gh.correction = "no_gh", alpha = .05, flr = 0, ceil = 9)
#' @export cf_correction
cf_correction = function(x, tests="all", df.adjustment="trunc", gh.correction="no_gh", alpha=.05, flr, ceil) {
if (inherits(x, "aov")) {
input <- "aov"
} else if (is.data.frame(x) || is.matrix(x)) {
input <- "data"
} else {
stop("Input must be a data frame or an aov object.")
}
if (input == "data") {
data <- as.data.frame(x)
colnames(data) <- c("group_id", "value")
} else if (input == "aov") {
aov_obj <- x
data <- model.frame(aov_obj)
colnames(data) <- c("value", "group_id")
data <- data[, c("group_id", "value")]
}
if (flr >= ceil) stop("flr must be less than ceil.")
data <- data[order(data$group_id), ]
ns <- table(data$group_id)
kgroups <- length(ns)
unique_ids <- unique(data$group_id)
comparisons <- combn(unique_ids, 2)
rownames(comparisons) <- c("i", "j")
num_comparisons <- dim(comparisons)[2]
floor_counts <- tapply(data$value, data$group_id, function(x) sum(x == flr))
ceil_counts <- tapply(data$value, data$group_id, function(x) sum(x == ceil))
group_means_adj <- c()
group_var_adj <- c()
for(g in 1:kgroups){
y <- data[data$group_id == unique_ids[g], 2]
adj_vals <- rec.mean.var(y, flr=flr, ceil=ceil)
group_means_adj <- c(group_means_adj, adj_vals$est.mean)
group_var_adj <- c(group_var_adj, adj_vals$est.var)
}
names(group_means_adj) <- unique_ids
names(group_var_adj) <- unique_ids
group_sd_adj <- sqrt(group_var_adj)
SS_within_adj <- sum((ns - 1) * group_var_adj)
df_within_adj <- dim(data)[1] - kgroups
if (df.adjustment == "trunc") {
df <- sum(tapply(data$value, data$group_id,
function(x) length(x) - (sum(x == flr) - 1) - (sum(x == ceil) - 1)))
} else if (df.adjustment == "unadj") {
df <- sum(tapply(data$value, data$group_id, function(x) length(x) - 1))
}
MSE_adj <- SS_within_adj / df
results <- c()
for(c in 1:num_comparisons){
grp_i <- unname(comparisons[1, c])
grp_j <- unname(comparisons[2, c])
mean_i <- unname(group_means_adj[grp_i])
mean_j <- unname(group_means_adj[grp_j])
diff_in_means <- unname(mean_i - mean_j)
abs_diff <- unname(abs(mean_i - mean_j))
n_i <- unname(ns[grp_i])
n_j <- unname(ns[grp_j])
var_i <- unname(group_var_adj[grp_i])
var_j <- unname(group_var_adj[grp_j])
sqrt_term <- sqrt((MSE_adj/2) * (1/n_i + 1/n_j))
rmse <- sqrt(MSE_adj)
# Tukey
ci_term <- (qtukey(.95, kgroups, df=df) / sqrt(2)) * sqrt(MSE_adj * (1/n_i + 1/n_j))
lwr <- unname(diff_in_means - ci_term)
upr <- unname(diff_in_means + ci_term)
hedges_g <- unname(diff_in_means / rmse)
Q <- unname(abs_diff / sqrt_term)
p <- unname(ptukey(q=Q, nmeans=kgroups, df=df, lower.tail=FALSE))
# Bonferroni
t <- unname(diff_in_means / sqrt(MSE_adj * (1/n_i + 1/n_j)))
t_p <- unname(2 * pt(-abs(t), df, lower.tail=TRUE))
t_p_bonferroni <- unname(pmin(1, num_comparisons * t_p))
adjusted_alpha <- alpha / num_comparisons
ci_term_bonf <- qt(1 - adjusted_alpha/2, df) * sqrt(MSE_adj * (1/n_i + 1/n_j))
lwr_bonf <- unname(diff_in_means - ci_term_bonf)
upr_bonf <- unname(diff_in_means + ci_term_bonf)
# Games-Howell / Welch
vi <- unname(var_i / n_i)
vj <- unname(var_j / n_j)
df_i <- unname(n_i - floor_counts[grp_i] - ceil_counts[grp_i] - 1)
df_j <- unname(n_j - floor_counts[grp_j] - ceil_counts[grp_j] - 1)
se_ij <- unname(sqrt(vi + vj))
t_ij <- unname(diff_in_means / se_ij)
df_ij <- unname((vi + vj)^2 / (vi^2 / df_i + vj^2 / df_j))
df_ij <- max(df_ij, 1)
p_gh <- unname(ptukey(abs(t_ij) * sqrt(2), kgroups, df_ij, lower.tail=FALSE))
se_gh <- unname(sqrt(0.5 * (var_i / n_i + var_j / n_j)))
lwr_gh <- unname(diff_in_means - qtukey(p=0.95, nmeans=kgroups, df=df_ij) * se_gh)
upr_gh <- unname(diff_in_means + qtukey(p=0.95, nmeans=kgroups, df=df_ij) * se_gh)
hedges_g_gh <- unname(diff_in_means / sqrt(mean(group_var_adj)))
t_p_ij <- unname(2 * pt(-abs(t_ij), df_ij, lower.tail=TRUE))
t_p_bonf_ij <- unname(pmin(1, num_comparisons * t_p_ij))
ci_term_gh <- unname(qt(1 - (alpha / num_comparisons) / 2, df_ij) * se_ij)
lwr_gh_bonf <- unname(diff_in_means - ci_term_gh)
upr_gh_bonf <- unname(diff_in_means + ci_term_gh)
if (gh.correction == "no_gh") {
if (tests == "all") {
results <- rbind(results, c(
Comparison_i = grp_i,
Comparison_j = grp_j,
mean_i = mean_i,
mean_j = mean_j,
diff_in_means = diff_in_means,
tukey.CI_lwr = lwr,
tukey.CI_upr = upr,
hedges_g = hedges_g,
Q = Q,
p = p,
t = t,
t_p = t_p,
p.bonferroni = t_p_bonferroni,
bonf.CI_lwr = lwr_bonf,
bonf.CI_upr = upr_bonf
))
} else if (tests == "tukey") {
results <- rbind(results, c(
Comparison_i = grp_i,
Comparison_j = grp_j,
mean_i = mean_i,
mean_j = mean_j,
diff_in_means = diff_in_means,
tukey.CI_lwr = lwr,
tukey.CI_upr = upr,
hedges_g = hedges_g,
Q = Q,
p = p
))
} else if (tests == "bonf") {
results <- rbind(results, c(
Comparison_i = grp_i,
Comparison_j = grp_j,
mean_i = mean_i,
mean_j = mean_j,
diff_in_means = diff_in_means,
t = t,
t_p = t_p,
p.bonferroni = t_p_bonferroni,
bonf.CI_lwr = lwr_bonf,
bonf.CI_upr = upr_bonf
))
}
} else {
if (tests == "all") {
results <- rbind(results, c(
Comparison_i = grp_i,
Comparison_j = grp_j,
mean_i = mean_i,
mean_j = mean_j,
diff_in_means = diff_in_means,
tukey.CI_lwr = lwr,
tukey.CI_upr = upr,
hedges_g = hedges_g,
Q = Q,
p = p,
t = t,
t_p = t_p,
p.bonferroni = t_p_bonferroni,
bonf.CI_lwr = lwr_bonf,
bonf.CI_upr = upr_bonf,
gh.CI_lwr = lwr_gh,
gh.CI_upr = upr_gh,
games.howell.t = t_ij,
games.howell.p = p_gh,
games.howell.g = hedges_g_gh,
welch.satterthwaite.df = df_ij,
welch.bonf.CI_lwr = lwr_gh_bonf,
welch.bonf.CI_upr = upr_gh_bonf,
t_Welch = t_ij,
p_Welch = t_p_ij,
p_Welch.bonferroni = t_p_bonf_ij
))
} else if (tests == "tukey") {
results <- rbind(results, c(
Comparison_i = grp_i,
Comparison_j = grp_j,
mean_i = mean_i,
mean_j = mean_j,
diff_in_means = diff_in_means,
tukey.CI_lwr = lwr,
tukey.CI_upr = upr,
hedges_g = hedges_g,
Q = Q,
p = p,
gh.CI_lwr = lwr_gh,
gh.CI_upr = upr_gh,
games.howell.t = t_ij,
games.howell.p = p_gh,
games.howell.g = hedges_g_gh,
welch.satterthwaite.df = df_ij
))
} else if (tests == "bonf") {
results <- rbind(results, c(
Comparison_i = grp_i,
Comparison_j = grp_j,
mean_i = mean_i,
mean_j = mean_j,
diff_in_means = diff_in_means,
t = t,
t_p = t_p,
p.bonferroni = t_p_bonferroni,
bonf.CI_lwr = lwr_bonf,
bonf.CI_upr = upr_bonf,
welch.bonf.CI_lwr = lwr_gh_bonf,
welch.bonf.CI_upr = upr_gh_bonf,
t_Welch = t_ij,
p_Welch = t_p_ij,
p_Welch.bonferroni = t_p_bonf_ij
))
}
}
} # end comparisons loop
return(results)
}
Try the DACF package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.