diff --git a/R/ancova.R b/R/ancova.R index 2c6dcf57..9806f41a 100644 --- a/R/ancova.R +++ b/R/ancova.R @@ -298,7 +298,7 @@ AncovaInternal <- function(jaspResults, dataset = NULL, options) { model <- aov(model.formula, dataset, weights=WLS) - modelError <- try(silent = TRUE, lm(model.formula, dataset, weights=WLS, singular.ok = FALSE)) + modelError <- try(silent = TRUE, lm(model.formula, dataset, weights=WLS, singular.ok = TRUE)) # Make afex model for later use in contrasts if (!isTryError(modelError)) { @@ -335,7 +335,10 @@ AncovaInternal <- function(jaspResults, dataset = NULL, options) { } else if (.extractErrorMessage(model$modelError) == "residual sum of squares is 0 (within rounding error)") { anovaContainer$setError(gettext("Residual sum of squares is 0; this might be due to extremely low variance of your dependent variable")) return() - } else { + } else if (grepl(x = .extractErrorMessage(model$modelError), "Rank deficient model matrix")) { + anovaContainer[["model"]] <- createJaspState(object = model$model) + return() + }else { anovaContainer$setError(gettextf("An error occurred while computing the ANOVA: %s", .extractErrorMessage(model$modelError))) return() } @@ -667,9 +670,12 @@ AncovaInternal <- function(jaspResults, dataset = NULL, options) { anovaContainer[["contrastContainer"]] <- contrastContainer - if (!ready || anovaContainer$getError()) + if (!ready || anovaContainer$getError()) { return() - + } else if(is.null(anovaContainer[["afexModel"]]$object)) { + contrastContainer$setError(gettext("Insufficient data to estimate full model/contrasts. Likely empty cells in between-subjects design (i.e., bad data structure)")) + return() + } afexModel <- anovaContainer[["afexModel"]]$object ## Computation @@ -879,21 +885,18 @@ AncovaInternal <- function(jaspResults, dataset = NULL, options) { } if (options$postHocTypeStandard) - .anovaPostHocTable(postHocContainer, dataset, options, anovaContainer[["model"]]$object) + .anovaPostHocTable(postHocContainer, dataset, options, anovaContainer[["model"]]$object, anovaContainer[["afexModel"]]$object) if (options$postHocTypeGames) - gamesPostHoc <- .anovaGamesTable(postHocContainer, dataset, options, anovaContainer[["model"]]$object) + .anovaGamesTable(postHocContainer, dataset, options, anovaContainer[["model"]]$object) if (options$postHocTypeDunnet) - dunnettPostHoc <- .anovaDunnettTable(postHocContainer, dataset, options, anovaContainer[["model"]]$object) - - if (options$postHocTypeDunn) - dunnPostHoc <- .anovaDunnTable(postHocContainer, dataset, options, anovaContainer[["model"]]$object) + .anovaDunnettTable(postHocContainer, dataset, options, anovaContainer[["model"]]$object) return() } -.anovaPostHocTable <- function(postHocContainer, dataset, options, model) { +.anovaPostHocTable <- function(postHocContainer, dataset, options, model, afexModel = NULL) { if (!is.null(postHocContainer[["postHocStandardContainer"]])) return() @@ -902,7 +905,7 @@ AncovaInternal <- function(jaspResults, dataset = NULL, options) { "postHocCorrectionBonferroni", "postHocCorrectionHolm", "postHocCorrectionScheffe", "postHocCorrectionTukey", "postHocCorrectionSidak", "postHocSignificanceFlag", "postHocTypeStandardBootstrap", "postHocTypeStandardBootstrapSamples", - "postHocCi", "postHocCiLevel")) + "postHocCi", "postHocCiLevel", "postHocLetterTable", "postHocLetterAlpha")) postHocContainer[["postHocStandardContainer"]] <- postHocStandardContainer @@ -916,7 +919,20 @@ AncovaInternal <- function(jaspResults, dataset = NULL, options) { postHocStandardContainer[[thisVarName]] <- .createPostHocStandardTable(thisVarName, interactionTerm, options, options$postHocTypeStandardBootstrap) - } + if (options$postHocLetterTable) { + letterTable <- createJaspTable(title = paste0("Letter-Based Grouping - ", thisVarName)) + for (letterVar in postHocVariables[[postHocVarIndex]]) + letterTable$addColumnInfo(name=letterVar, type="string", combine = TRUE) + + letterTable$addColumnInfo(name="Letter", type="string") + letterTable$addFootnote("If two or more means share the same grouping symbol, + then we cannot show them to be different, but we also did not show them to be the same. ") + letterTable$showSpecifiedColumnsOnly <- TRUE + + postHocStandardContainer[[paste0(thisVarName, "LetterTable")]] <- letterTable + } + + postHocStandardContainer[[paste0(thisVarName, "LetterTable")]]} for (postHocVarIndex in 1:length(postHocVariables)) { @@ -937,10 +953,16 @@ AncovaInternal <- function(jaspResults, dataset = NULL, options) { bonfAdjustCIlevel <- .computeBonferroniConfidence(options$postHocCiLevel, numberOfLevels = numberOfLevels) - effectSizeResult <- as.data.frame(emmeans::eff_size(postHocRef[[postHocVarIndex]], - sigma = sqrt(mean(sigma(model)^2)), - edf = df.residual(model), - level = bonfAdjustCIlevel)) + if (isTRUE(options[["postHocLetterTable"]])) { + letterResult <- multcomp::cld(postHocRef[[postHocVarIndex]], + method = "pairwise", + Letters = letters, + alpha = options[["postHocLetterAlpha"]]) + letterResult <- letterResult[c(postHocVariables[[postHocVarIndex]], ".group")] + colnames(letterResult)[ncol(letterResult)] <- "Letter" + letterResult <- letterResult[order(as.numeric(rownames(letterResult))), ] + postHocStandardContainer[[paste0(thisVarName, "LetterTable")]]$setData(letterResult) + } allContrasts <- strsplit(as.character(resultPostHoc[[1]]$contrast), split = " - ") @@ -962,9 +984,17 @@ AncovaInternal <- function(jaspResults, dataset = NULL, options) { resultPostHoc[["contrast_A"]] <- lapply(allContrasts, `[[`, 1L) resultPostHoc[["contrast_B"]] <- lapply(allContrasts, `[[`, 2L) - resultPostHoc[["cohenD"]] <- effectSizeResult[["effect.size"]] - resultPostHoc[["cohenD_LowerCI"]] <- effectSizeResult[["lower.CL"]] - resultPostHoc[["cohenD_UpperCI"]] <- effectSizeResult[["upper.CL"]] + if (isFALSE(is.null(afexModel))) { + effectSizeResult <- as.data.frame(emmeans::eff_size(postHocRef[[postHocVarIndex]], + sigma = sqrt(mean(sigma(model)^2)), + edf = df.residual(model), + level = bonfAdjustCIlevel)) + resultPostHoc[["cohenD"]] <- effectSizeResult[["effect.size"]] + resultPostHoc[["cohenD_LowerCI"]] <- effectSizeResult[["lower.CL"]] + resultPostHoc[["cohenD_UpperCI"]] <- effectSizeResult[["upper.CL"]] + } else { + postHocStandardContainer[[thisVarName]]$addFootnote(message = gettext("Some parameters were not estimable due to missingness.")) + } if (options$postHocTypeStandardBootstrap) { @@ -1059,105 +1089,6 @@ AncovaInternal <- function(jaspResults, dataset = NULL, options) { } } -.anovaDunnTable <- function(postHocContainer, dataset, options, model) { - if (!is.null(postHocContainer[["postHocDunnContainer"]])) - return() - - postHocDunnContainer <- createJaspContainer(title = gettext("Dunn")) - postHocDunnContainer$dependOn(c("postHocTypeDunn", "postHocSignificanceFlag")) - postHocContainer[["postHocDunnContainer"]] <- postHocDunnContainer - - postHocVariables <- unlist(options$postHocTerms, recursive = FALSE) - postHocVariablesListV <- unname(lapply(postHocVariables, .v)) - - dunnVariables <- unique(unlist(options$postHocTerms)) - dependentVar <- options$dependent - - .createPostHocDunnTable <- function(myTitle) { - - postHocTable <- createJaspTable(title = gettextf("Dunn's Post Hoc Comparisons - %s", myTitle)) - - postHocTable$addColumnInfo(name="contrast", title=gettext("Comparison"), type="string") - postHocTable$addColumnInfo(name="z", type="number") - postHocTable$addColumnInfo(name="wA", title=gettext("Wi"), type="number") - postHocTable$addColumnInfo(name="wB", title=gettext("Wj"), type="number") - postHocTable$addColumnInfo(name="rbs", title=gettext("Rank-Biserial Correlation"), type="number") - postHocTable$addColumnInfo(name="pval", title=gettext("p"), type="pvalue") - postHocTable$addColumnInfo(name="bonferroni", title=gettext("pbonf"), type="pvalue") - postHocTable$addColumnInfo(name="holm", title=gettext("pholm"), type="pvalue") - - return(postHocTable) - } - - for (dunnVar in dunnVariables) { - - postHocDunnContainer[[dunnVar]] <- .createPostHocDunnTable(dunnVar) - dunnResult <- data.frame(contrast = character(), - z = numeric(), - wA = numeric(), - wB = numeric(), - rbs = numeric(), - pval = numeric(), - bonferroni = numeric(), - holm = numeric()) - - variableLevels <- levels(droplevels(dataset[[ .v(dunnVar) ]])) - nLevels <- length(variableLevels) - nPerGroup <- unname(unlist(table(dataset[[ .v(dunnVar) ]]))) - bigN <- sum(nPerGroup) - - fullRanks <- rank(dataset[[ dependentVar ]]) - ranksPerGroup <- by(fullRanks, dataset[[ dunnVar ]], list) - sumPerGroup <- unlist(lapply(ranksPerGroup, FUN = sum)) - meanPerGroup <- unname(sumPerGroup/nPerGroup) - - tab <- table(unlist(ranksPerGroup)) - nTies <- tab[tab > 1] - nTies <- sum(nTies^3 - nTies) - - for (i in 1:(nLevels - 1)) { - - for (j in (i + 1):nLevels) { - - contrast <- paste0(variableLevels[[i]], " - ", variableLevels[[j]]) - - sigmaAB <- sqrt( ( (bigN * (bigN + 1))/12 - nTies/(12 * (bigN - 1)) ) * (1/nPerGroup[i] + 1/nPerGroup[j] ) ) - zAB <- (meanPerGroup[i] - meanPerGroup[j]) / sigmaAB - pValAB <- 2 * pnorm(abs(zAB), lower.tail = FALSE) # make two-sided p-value - - a <- dataset[[ dependentVar ]][dataset[[ dunnVar ]] == variableLevels[[i]]] - b <- dataset[[ dependentVar ]][dataset[[ dunnVar ]] == variableLevels[[j]]] - u <- wilcox.test(a, b)$statistic - rbs <- abs(as.numeric(1-(2*u)/(nPerGroup[i]*nPerGroup[j]))) - - dunnResult <- rbind(dunnResult, data.frame(contrast = contrast, - z = zAB, - wA = meanPerGroup[i], - wB = meanPerGroup[j], - rbs = rbs, - pval = pValAB, - bonferroni = pValAB, - holm = pValAB)) - - } - - } - - allP <- dunnResult[["pval"]] - dunnResult[["bonferroni"]] <- p.adjust(allP, method = "bonferroni") - dunnResult[["holm"]] <- p.adjust(allP, method = "holm") - - postHocDunnContainer[[dunnVar]]$setData(dunnResult) - postHocDunnContainer[[dunnVar]]$addFootnote(message = gettext("Rank-biserial correlation based on individual Mann-Whitney tests.")) - - if (options$postHocSignificanceFlag) - .anovaAddSignificanceSigns(someTable = postHocDunnContainer[[dunnVar]], - allPvalues = cbind(dunnResult[, c("pval", "bonferroni", "holm")]), - resultRowNames = rownames(dunnResult)) - } - - return() -} .anovaGamesTable <- function(postHocContainer, dataset, options, model) { if (!is.null(postHocContainer[["postHocGamesContainer"]])) @@ -1661,7 +1592,7 @@ AncovaInternal <- function(jaspResults, dataset = NULL, options) { anovaContainer[["kruskalContainer"]] <- createJaspContainer(title = gettext("Kruskal-Wallis Test"), dependencies = c("kruskalWallisFactors", "kruskalCiLevel", - "kruskalEpsilon", "kruskalEta", + "kruskalEpsilon", "kruskalEta", "postHocTypeDunn", "kruskalEffectSizeEstimates")) kruskalTable <- createJaspTable(title = gettext("Kruskal-Wallis Test")) @@ -1687,7 +1618,8 @@ AncovaInternal <- function(jaspResults, dataset = NULL, options) { } kruskalTable$showSpecifiedColumnsOnly <- TRUE - + if (length(options[["kruskalWallisFactors"]]) > 1) + kruskalTable$addFootnote(gettext("Tests are conducted separately for each variable, without accounting for multivariate effects.")) if (!ready || anovaContainer$getError()) return() @@ -1714,6 +1646,107 @@ AncovaInternal <- function(jaspResults, dataset = NULL, options) { lowerEtaCI = sapply(kruskalEtaList, function(x) x$CI_low), upperEtaCI = sapply(kruskalEtaList, function(x) x$CI_high))) + if (options$postHocTypeDunn) + .anovaDunnTable(anovaContainer[["kruskalContainer"]], dataset, options, anovaContainer[["model"]]$object) + + + return() } +.anovaDunnTable <- function(kruskalContainer, dataset, options, model) { + if (!is.null(kruskalContainer[["dunnContainer"]])) + return() + + kruskalContainer[["dunnContainer"]] <- createJaspContainer(title = gettext("Dunn")) + + dunnVariables <- options[["kruskalWallisFactors"]] + dependentVar <- options[["dependent"]] + + .createPostHocDunnTable <- function(myTitle) { + + postHocTable <- createJaspTable(title = gettextf("Dunn's Post Hoc Comparisons - %s", myTitle)) + + postHocTable$addColumnInfo(name="contrast", title=gettext("Comparison"), type="string") + postHocTable$addColumnInfo(name="z", type="number") + postHocTable$addColumnInfo(name="wA", title=gettext("Wi"), type="number") + postHocTable$addColumnInfo(name="wB", title=gettext("Wj"), type="number") + postHocTable$addColumnInfo(name="rbs", title=gettext("Rank-Biserial Correlation"), type="number") + postHocTable$addColumnInfo(name="pval", title=gettext("p"), type="pvalue") + postHocTable$addColumnInfo(name="bonferroni", title=gettext("pbonf"), type="pvalue") + postHocTable$addColumnInfo(name="holm", title=gettext("pholm"), type="pvalue") + + return(postHocTable) + } + + for (dunnVar in dunnVariables) { + + kruskalContainer[["dunnContainer"]][[dunnVar]] <- .createPostHocDunnTable(dunnVar) + dunnResult <- data.frame(contrast = character(), + z = numeric(), + wA = numeric(), + wB = numeric(), + rbs = numeric(), + pval = numeric(), + bonferroni = numeric(), + holm = numeric()) + + variableLevels <- levels(droplevels(dataset[[ .v(dunnVar) ]])) + nLevels <- length(variableLevels) + nPerGroup <- unname(unlist(table(dataset[[ .v(dunnVar) ]]))) + bigN <- sum(nPerGroup) + + fullRanks <- rank(dataset[[ dependentVar ]]) + ranksPerGroup <- by(fullRanks, dataset[[ dunnVar ]], list) + sumPerGroup <- unlist(lapply(ranksPerGroup, FUN = sum)) + meanPerGroup <- unname(sumPerGroup/nPerGroup) + + tab <- table(unlist(ranksPerGroup)) + nTies <- tab[tab > 1] + nTies <- sum(nTies^3 - nTies) + + for (i in 1:(nLevels - 1)) { + + for (j in (i + 1):nLevels) { + + contrast <- paste0(variableLevels[[i]], " - ", variableLevels[[j]]) + + sigmaAB <- sqrt( ( (bigN * (bigN + 1))/12 - nTies/(12 * (bigN - 1)) ) * (1/nPerGroup[i] + 1/nPerGroup[j] ) ) + zAB <- (meanPerGroup[i] - meanPerGroup[j]) / sigmaAB + pValAB <- 2 * pnorm(abs(zAB), lower.tail = FALSE) # make two-sided p-value + + a <- dataset[[ dependentVar ]][dataset[[ dunnVar ]] == variableLevels[[i]]] + b <- dataset[[ dependentVar ]][dataset[[ dunnVar ]] == variableLevels[[j]]] + u <- wilcox.test(a, b)$statistic + rbs <- abs(as.numeric(1-(2*u)/(nPerGroup[i]*nPerGroup[j]))) + + dunnResult <- rbind(dunnResult, data.frame(contrast = contrast, + z = zAB, + wA = meanPerGroup[i], + wB = meanPerGroup[j], + rbs = rbs, + pval = pValAB, + bonferroni = pValAB, + holm = pValAB)) + + } + + } + + allP <- dunnResult[["pval"]] + dunnResult[["bonferroni"]] <- p.adjust(allP, method = "bonferroni") + dunnResult[["holm"]] <- p.adjust(allP, method = "holm") + + kruskalContainer[["dunnContainer"]][[dunnVar]]$setData(dunnResult) + kruskalContainer[["dunnContainer"]][[dunnVar]]$addFootnote(message = gettext("Rank-biserial correlation based on individual Mann-Whitney tests.")) + + if (options$postHocSignificanceFlag) + .anovaAddSignificanceSigns(someTable = kruskalContainer[["dunnContainer"]][[dunnVar]], + allPvalues = cbind(dunnResult[, c("pval", "bonferroni", "holm")]), + resultRowNames = rownames(dunnResult)) + } + + return() +} + + diff --git a/inst/help/Ancova.md b/inst/help/Ancova.md index fc06c512..7af75133 100755 --- a/inst/help/Ancova.md +++ b/inst/help/Ancova.md @@ -72,27 +72,29 @@ ANCOVA allows the user to analyze the difference between multiple group means, w ### Post Hoc Tests - To perform a post hoc test, drag one or more factor names to the right column. Several options are available: -- Type: Different types of post hoc tests can be selected. -- Standard (LSD): Pairwise t-tests are performed. All the corrections can be applied to this method. This option is selected by default. - - Confidence intervals: When this option is selected, the confidence interval for the mean difference is calculated. This is done for every post hoc method except for Dunn. By default this is set to 95% but this can be adjusted into the desired percentage. - - From `...` bootstraps: By selecting this option, the bootstrapped post hoc test is applied. By default, the number of replications is set to 1000. This can be changed into the desired number. - - Effect size: By selecting this option, the effect size (i.e., the magnitude of the observed effect) will be displayed. The used measure for the effect size is Cohen's d. The effect size will only be displayed for the post hoc type `Standard`. - - - Games-Howell: This method can be used when equal group/level variances are not assumed. The p-values are corrected with the Tukey method. - - Dunnett: When selecting this method, all the levels are compared to one specific level, for example to the control group. At the moment, it is not possible to manually specify to which level the others levels are compared, but it is based on the order of the levels. To change the order of the levels, the level labels can be adjusted. -
- GIF demonstration: Adjust level labels - -
- - - Dunn: This is a non-parametric test that can be used for testing small subsets of pairs. This post hoc test is a follow up for the Kruskal-Wallis test. The p-values are corrected with the Bonferroni and Holm methods. - -- Correction: To correct for multiple comparison testing and avoid Type I errors, different methods for correcting the p-value are available: - - Tukey: Compare all possible pairs of group means. This correction can be used when the groups of the independent variable have an equal sample size and variance. This method is commonly used and is selected by default. - - Scheffe: Adjusting significance levels in a linear regression, to account for multiple comparisons. This method is considered to be quite conservative. - - Bonferroni: This correction is considered conservative. The risk of Type I error is reduced, however the statistical power decreases as well. - - Holm: This method is also called sequential Bonferroni, and considered less conservative than the Bonferroni method. -- Flag Significant Comparisons: Add asterisks to the table to indicate 3 levels of significance. + - Type: Different types of post hoc tests can be selected. + - Standard (LSD): Pairwise t-tests are performed. All the corrections can be applied to this method. This option is selected by default. + - Confidence intervals: When this option is selected, the confidence interval for the mean difference is calculated. This is done for every post hoc method except for Dunn. By default this is set to 95% but this can be adjusted into the desired percentage. + - From `...` bootstraps: By selecting this option, the bootstrapped post hoc test is applied. By default, the number of replications is set to 1000. This can be changed into the desired number. + - Effect size: By selecting this option, the effect size (i.e., the magnitude of the observed effect) will be displayed. The used measure for the effect size is Cohen's d. The effect size will only be displayed for the post hoc type `Standard`. + + - Games-Howell: This method can be used when equal group/level variances are not assumed. The p-values are corrected with the Tukey method. + - Dunnett: When selecting this method, all the levels are compared to one specific level, for example to the control group. At the moment, it is not possible to manually specify to which level the others levels are compared, but it is based on the order of the levels. To change the order of the levels, the level labels can be adjusted. +
+ GIF demonstration: Adjust level labels + +
+ + - Dunn: This is a non-parametric test that can be used for testing small subsets of pairs. This post hoc test is a follow up for the Kruskal-Wallis test. The p-values are corrected with the Bonferroni and Holm methods. + + - Correction: To correct for multiple comparison testing and avoid Type I errors, different methods for correcting the p-value are available: + - Tukey: Compare all possible pairs of group means. This correction can be used when the groups of the independent variable have an equal sample size and variance. This method is commonly used and is selected by default. + - Scheffe: Adjusting significance levels in a linear regression, to account for multiple comparisons. This method is considered to be quite conservative. + - Bonferroni: This correction is considered conservative. The risk of Type I error is reduced, however the statistical power decreases as well. + - Holm: This method is also called sequential Bonferroni, and considered less conservative than the Bonferroni method. + - Flag significant comparisons: Add asterisks to the table to indicate 3 levels of significance. + - Letter-based grouping table: Set up a compact letter display of all pair-wise comparisons, based on 'multcomp::cld' and emmeans. + ### Descriptive Plots - To create a descriptive plot, select the independent variable to be placed on the horizontal axis. If there are more than one independent variable, the variables can be displayed in one plot by putting the other variable in the box Separate lines, or the variables can be displayed in separate plots by selecting the other variable in the box Separate plots. diff --git a/inst/help/Anova.md b/inst/help/Anova.md index 9367163d..9f223fde 100755 --- a/inst/help/Anova.md +++ b/inst/help/Anova.md @@ -89,6 +89,8 @@ ANOVA allows the user to analyze the difference between multiple group means. - Dunn: This is a non-parametric test that can be used for testing small subsets of pairs. This post hoc test is a follow up for the Kruskal-Wallis test. The p-values are corrected with the Bonferroni and Holm methods. + - Flag significant comparisons: Add asterisks to the table to indicate 3 levels of significance. + - Letter-based grouping table: Set up a compact letter display of all pair-wise comparisons, based on 'multcomp::cld' and emmeans. ### Descriptive Plots - To create a descriptive plot, select the independent variable to be placed on the horizontal axis. If there are more than one independent variable, the variables can be displayed in one plot by putting the other variable in the box Separate lines, or the variables can be displayed in separate plots by selecting the other variable in the box Separate plots. diff --git a/inst/qml/common/classical/PostHocDisplay.qml b/inst/qml/common/classical/PostHocDisplay.qml index 7c3efeb5..fb25cbbf 100644 --- a/inst/qml/common/classical/PostHocDisplay.qml +++ b/inst/qml/common/classical/PostHocDisplay.qml @@ -30,5 +30,19 @@ Group childrenOnSameRow: true CIField { name: "postHocCiLevel" } } - CheckBox { name: "postHocSignificanceFlag"; label: qsTr("Flag Significant Comparisons") } + CheckBox { name: "postHocSignificanceFlag"; label: qsTr("Flag significant comparisons") } + CheckBox + { + name: "postHocLetterTable" + label: qsTr("Letter-based grouping table") + childrenOnSameRow: true + DoubleField + { + name: "postHocLetterAlpha" + label: qsTr("α-level:") + min: 0 + max: 1 + defaultValue: 0.05} + } + }