R/commonAnovaBayesian.R
In jasp-stats/jaspAnova: ANOVA Module for JASP
Defines functions
.BANOVAestimateModels
.BANOVAmissingInteractionCells
.BANOVAerrorhandling
.BANOVArunAnalysis
#
# Copyright (C) 2013-2022 University of Amsterdam
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Abbreviations:
# HF = helper function
.BANOVArunAnalysis <- function(jaspResults, dataset, options, analysisType = c("ANOVA", "ANCOVA", "RM-ANOVA")) {
# the main workhorse of the Bayesian ANOVA, ANCOVA, and Repeated Measures.
analysisType <- match.arg(analysisType)
dataset <- .BANOVAreadData(dataset, options, analysisType)
if (isTryError(dataset))
.quitAnalysis(gettext("Error while loading data. Please verify your repeated measures observations."))
dataInfo <- .BANOVAerrorhandling(dataset, options, analysisType)
model <- .BANOVAestimateModels(jaspResults, dataset, options, dataInfo, analysisType)
model[["posteriors"]] <- .BANOVAestimatePosteriors(jaspResults, dataset, options, model)
.BANOVAeffectsTable (jaspResults, options, model)
.BANOVAestimatesTable(jaspResults, options, model)
.BANOVArsquaredTable (jaspResults, options, model)
# model averaged plots
.BANOVAposteriorPlot(jaspResults, dataset, options, model)
.BANOVAqqplot (jaspResults, options, model)
.BANOVArsqplot (jaspResults, options, model)
.BANOVAnullControlPostHocTable(jaspResults, dataset, options, model)
# single model plots and tables
.BANOVAsmi(jaspResults, dataset, options, model)
# descriptives
.BANOVAdescriptives(jaspResults, dataset, options, dataInfo, analysisType)
return()
}
.BANOVAerrorhandling <- function(dataset, options, analysisType) {
customChecks <- list()
if (analysisType != "RM-ANOVA") {
hasDV <- options$dependent != ""
hasIV <- any(lengths(options[c("fixedFactors", "covariates")]) != 0)
fixed <- options$fixedFactors
noVariables <- !(hasDV && hasIV)
target <- c(options$covariates, options$dependent)
# check if the last model has an interaction effect
mostComplexModel <- NULL
if (length(options[["modelTerms"]]) > 0L) {
mostComplexModel <- options[["modelTerms"]][[length(options[["modelTerms"]])]][["components"]]
mostComplexModel <- intersect(mostComplexModel, fixed)
}
if (length(mostComplexModel) > 1L)
customChecks[["missingInteractionCells"]] <- .BANOVAmissingInteractionCells
if (!noVariables) {
.hasErrors(
dataset = dataset,
type = c("infinity", "observations", "variance", "factorLevels", "duplicateColumns"),
infinity.target = target,
variance.target = target,
observations.target = target,
observations.amount = paste("<", length(options$modelTerms) + 1),
factorLevels.target = fixed,
factorLevels.amount = " < 2",
# custom checks
custom = customChecks,
missingInteractionCells.target = mostComplexModel,
exitAnalysisIfErrors = TRUE
)
}
} else {
hasDV <- !any(options$repeatedMeasuresCells == "")
hasIV <- any(lengths(options[c("betweenSubjectFactors", "covariates")]) != 0)
fixed <- unlist(options$betweenSubjectFactors)
covariates <- unlist(options$covariates)
noVariables <- !hasDV
target <- c(covariates, fixed)
# check if the last model has an interaction effect
mostComplexModel <- NULL
if (length(options[["modelTerms"]]) > 0L) {
mostComplexModel <- options[["modelTerms"]][[length(options[["modelTerms"]])]][["components"]]
mostComplexModel <- intersect(mostComplexModel, fixed)
}
if (length(mostComplexModel) > 1L)
customChecks[["missingInteractionCells"]] <- .BANOVAmissingInteractionCells
if (length(target) > 0L) {
.hasErrors(
dataset = dataset,
type = c("infinity", "observations", "variance", "factorLevels", "duplicateColumns"),
infinity.target = target,
variance.target = covariates,
observations.target = target,
observations.amount = paste("<", length(options$modelTerms) + 1),
factorLevels.target = fixed,
factorLevels.amount = " < 2",
duplicateColumns.target = target,
# custom checks
custom = customChecks,
missingInteractionCells.target = mostComplexModel,
exitAnalysisIfErrors = TRUE
)
}
if (!noVariables) {
# messages from commonmessages.R
customChecks <- list(
infinity = function() {
x <- dataset[, .BANOVAdependentName]
if (any(is.infinite(x)))
return(gettext("Infinity found in repeated measures cells."))
return(NULL)
},
observations = function() {
x <- dataset[, .BANOVAdependentName]
nObs <- length(options$modelTerms) + 1
if (length(na.omit(x)) <= nObs)
return(gettextf("Number of observations is < %s in repeated measures cells", nObs))
return(NULL)
},
variance = function() {
x <- dataset[, .BANOVAdependentName]
validValues <- x[is.finite(x)]
variance <- 0
if (length(validValues) > 1)
variance <- stats::var(validValues)
if (variance == 0)
return(gettext("The variance in the repeated measures cells is 0."))
return(NULL)
},
duplicateColumns = function() {
cols <- c(.BANOVAsubjectName, target)
if (length(cols) == 1L)
return()
datasetList <- as.list(dataset[, cols])
duplicatedCols <- duplicated(datasetList) | duplicated(datasetList, fromLast = TRUE)
if (any(duplicatedCols)) {
if (duplicatedCols[1L]) {
msg <- gettextf("Duplicate variables encountered in repeated measures cells, %s",
paste(target[duplicatedCols[-1L]], collapse = ", "))
} else {
msg <- gettextf("Duplicate variables encountered in %s",
paste(target[duplicatedCols[-1L]], collapse = ", "))
}
return(msg)
}
return(NULL)
},
redundantParticipantColumn = function() {
if (nlevels(dataset[[.BANOVAsubjectName]]) == nlevels(dataset[[target]])) {
return(gettextf("Duplicate participant column added to model terms, %s", target))
} else {
return()
}
}
)
.hasErrors(
dataset = dataset,
target = target,
custom = customChecks,
exitAnalysisIfErrors = TRUE
)
}
}
return(list(noVariables = noVariables, hasIV = hasIV, hasDV = hasDV))
}
.BANOVAmissingInteractionCells <- function(dataset, target) {
# sorts a data frame by all its columns
sortDataFrame <- function(df) {
df[do.call(order, df), ]
}
uniqueObservedCombinations <- sortDataFrame(unique(dataset[, target, drop = FALSE]))
uniqueExpectedCombinations <- sortDataFrame(do.call(expand.grid, c(lapply(uniqueObservedCombinations, unique), stringsAsFactors = FALSE)))
# if the rows match then it must be
if (nrow(uniqueObservedCombinations) == nrow(uniqueExpectedCombinations))
return(NULL)
# find missing combinations
observedStrings <- unname(apply(uniqueObservedCombinations, 1L, paste, collapse = jaspBase::interactionSymbol))
expectedStrings <- unname(apply(uniqueExpectedCombinations, 1L, paste, collapse = jaspBase::interactionSymbol))
missingStrings <- setdiff(expectedStrings, observedStrings)
missingStringsCollapsed <- paste(missingStrings[1:min(length(missingStrings), 10)], collapse = ", ")
interactionVariable <- paste(colnames(uniqueObservedCombinations), collapse = jaspBase::interactionSymbol)
if (length(missingStrings) <= 10) {
return(gettextf(
"The interaction effect of %1$s has no observations for the combination(s): %2$s. Please adjust the model and remove any interactions with missing cells",
interactionVariable, missingStringsCollapsed
))
} else {
return(gettextf(
"The interaction effect of %1$s has no observations for the combination(s): %2$s. Only the first 10 missing combinations are shown. Please adjust the model and remove any interactions with missing cells",
interactionVariable, missingStringsCollapsed
))
}
}
# model comparison ----
.BANOVAestimateModels <- function(jaspResults, dataset, options, errors,
analysisType = c("ANOVA", "ANCOVA", "RM-ANOVA")) {
# also makes the model comparison table
stateObj <- jaspResults[["tableModelComparisonState"]]$object
if (!is.null(jaspResults[["tableModelComparison"]]) && !is.null(stateObj)) {
stateObj$completelyReused <- TRUE # means that posteriors won't need to be resampled
return(stateObj)
} else if (errors$noVariables) {
modelTable <- .BANOVAinitModelComparisonTable(options)
jaspResults[["tableModelComparison"]] <- modelTable
return(list(analysisType = analysisType))
} else if (!is.null(stateObj) &&
.BANOVAmodelTermsUnchanged(stateObj, options) &&
.BANOVAmarginalityOptionsUnchanged(stateObj, options) &&
(.BANOVAmodelBFTypeOrOrderChanged(stateObj, options) || identical(stateObj[["shownTableSize"]], options[c("modelsShown", "numModelsShown")]))) {
# if the statement above is TRUE then no new variables were added (or seed changed)
# and the only change is in the Bayes factor type, the ordering, or the number of models shown
modelTable <- .BANOVAinitModelComparisonTable(options)
modelTable[["Models"]] <- .BANOVAgetModelTitlesWithAllTerms(stateObj[["models"]], stateObj[["model.list"]], analysisType, options[["hideNuisanceParameters"]])
internalTable <- stateObj$internalTableObj$internalTable
nmodels <- nrow(internalTable)
.BANOVAestimateModels
shownTableSize <- .BANOVAgetShownTableSizeAndAddFootnote(modelTable, nmodels, options[["modelsShown"]], options[["numModelsShown"]])
if (.BANOVAmodelPriorOptionsChanged(stateObj, options)) {
priorProbs <- .BANOVAcomputePriorModelProbs(stateObj$model.list, stateObj$nuisance, options)
internalTable[["P(M)"]] <- priorProbs
stateObj$completelyReused <- FALSE # model-averaged posteriors need to be resampled
} else {
stateObj$completelyReused <- TRUE # model-averaged posteriors do not need to be resampled
}
internalTableObj <- .BANOVAfinalizeInternalTable(options, internalTable)
stateObj$postProbs <- internalTableObj$internalTable[["P(M|data)"]]
stateObj$priorProbs <- internalTableObj$internalTable[["P(M)"]]
stateObj$internalTableObj <- internalTableObj
modelTable$setData(internalTableObj$table[seq_len(shownTableSize), ])
jaspResults[["tableModelComparison"]] <- modelTable
return(stateObj)
}
modelTerms <- options$modelTerms
dependent <- options$dependent
if (analysisType == "RM-ANOVA") {
modelTerms[[length(modelTerms) + 1L]] <- list(components = .BANOVAsubjectName, isNuisance = TRUE)
dependent <- .BANOVAdependentName
}
fixedFactors <- options$fixedFactors
randomFactors <- .BANOVAgetRandomFactors(options, analysisType)
tempRScale <- .BANOVAgetRScale(options, analysisType)
rscaleFixed <- tempRScale[["rscaleFixed"]]
rscaleRandom <- tempRScale[["rscaleRandom"]]
rscaleCont <- tempRScale[["rscaleCont"]]
rscaleEffects <- tempRScale[["rscaleEffects"]]
tmp <- .BANOVAcreateModelFormula(dependent, modelTerms)
model.formula <- tmp$model.formula
nuisance <- tmp$nuisance
effects <- tmp$effects
if (analysisType == "RM-ANOVA") {
legacy <- options[["legacyResults"]]
rmFactors <- vapply(options[["repeatedMeasuresFactors"]], `[[`, character(1L), "name")
} else {
legacy <- FALSE
rmFactors <- NULL
}
#Make a list of models to compare
temp <- .BANOVAgenerateAllModelFormulas(
formula = model.formula,
nuisance = nuisance,
analysisType = analysisType,
enforcePrincipleOfMarginalityFixedEffects = options[["enforcePrincipleOfMarginalityFixedEffects"]],
enforcePrincipleOfMarginalityRandomSlopes = options[["enforcePrincipleOfMarginalityRandomSlopes"]],
rmFactors = rmFactors,
legacy = legacy
)
model.list <- temp$modelList
nuisance <- temp$nuisance
nuisanceRandomSlopes <- temp$nuisanceRandomSlopes
if (length(model.list) == 1L) {
modelTable <- .BANOVAinitModelComparisonTable(options)
modelTable$setError(gettext("Bayes factor is undefined -- all effects are specified as nuisance."))
jaspResults[["tableModelComparison"]] <- modelTable
return(list(analysisType = analysisType))
}
# TODO: discuss if we want a hard limit if length(model.list) > ...
# that would avoid peoples computers from slowing down tremendously.
# JCG: A hard limit would be a bit cruel right? How about a warning instead? (in QML already?)
#Run all other models and store the results in the model.object
neffects <- length(effects)
nmodels <- length(model.list)
modelObject <- vector("list", nmodels)
reorderedEffects <- .BANOVAreorderTerms(effects)
reorderedNuisance <- if (is.null(nuisance)) NULL else .BANOVAreorderTerms(nuisance)
if (nmodels > 0L && neffects > 0L) {
# effects.matrix contains for each row (model), which predictors (columns) are in the model
# interactions.matrix contains for each predictor (column) what predictors (rows) is it composed of
interactions.matrix <- .BANOVAcomputeInteractionsMatrix(effects)
effects.matrix <- matrix(data = FALSE, nrow = nmodels, ncol = neffects)
colnames(effects.matrix) <- effects
for (m in seq_len(nmodels)) {
modelObject[[m]] <- list()
if (m == 1L) {
if (is.null(nuisance)) { # intercept only
modelObject[[m]]$title <- gettext("Null model")
next # all effects are FALSE anyway
} else {
if (analysisType == "RM-ANOVA" && !legacy) {
tempNuisance <- setdiff(nuisance, nuisanceRandomSlopes)
modelObject[[m]]$title <- sprintf(ngettext(
length(tempNuisance),
"Null model (incl. %s and random slopes)",
"Null model (incl. %s, and random slopes)"
), paste(.BANOVAdecodeNuisance(tempNuisance), collapse = ", "))
} else {
modelObject[[m]]$title <- gettextf("Null model (incl. %s)", paste(.BANOVAdecodeNuisance(nuisance), collapse = ", "))
}
}
}
model.effects <- .BANOVAgetFormulaComponents(model.list[[m]])
reorderedModelEffects <- .BANOVAreorderTerms(model.effects)
idx <- match(reorderedModelEffects, reorderedEffects, nomatch = 0L)
idx <- idx[!is.na(idx)]
effects.matrix[m, idx] <- TRUE
if (m > 1L) {
model.title <- model.effects[match(reorderedModelEffects, reorderedNuisance, 0L) == 0L]
modelObject[[m]]$title <- jaspBase::gsubInteractionSymbol(paste(model.title, collapse = " + "))
}
}
rownames(effects.matrix) <- sapply(modelObject, `[[`, "title")
}
modelTable <- .BANOVAinitModelComparisonTable(options)
modelNames <- .BANOVAgetModelTitlesWithAllTerms(modelObject, model.list, analysisType, options[["hideNuisanceParameters"]])
shownTableSize <- .BANOVAgetShownTableSizeAndAddFootnote(modelTable, nmodels, options[["modelsShown"]], options[["numModelsShown"]])
modelTable[["Models"]] <- modelNames[seq_len(shownTableSize)]
jaspResults[["tableModelComparison"]] <- modelTable
# internalTable is an internal representation of the model comparison table
# internalTable <- matrix(NA, nmodels, 6L, dimnames = list(modelNames, c("Models", "P(M)", "P(M|data)", "BFM", "BF10", "error %")))
internalTable <- data.frame(
Models = character(nmodels),
`P(M)` = numeric(nmodels),
`P(M|data)` = numeric(nmodels),
BFM = numeric(nmodels),
BF10 = numeric(nmodels),
`error %` = numeric(nmodels),
check.names = FALSE
)
# set model names, BF null model and p(M)
internalTable[["Models"]] <- modelNames
internalTable[1L, "BF10"] <- 0
internalTable[["P(M)"]] <- .BANOVAcomputePriorModelProbs(model.list, nuisance, options)
#Now compute Bayes Factors for each model in the list, and populate the tables accordingly
startProgressbar(nmodels, gettext("Computing Bayes factors"))
# check if any models can be resued from the state
if (!is.null(stateObj[["modelTerms"]])) {
oldFormulas <- sapply(stateObj$model.list, .BANOVAreorderFormulas)
newFormulas <- sapply(model.list, .BANOVAreorderFormulas)
reuseable <- match(newFormulas, oldFormulas)
} else {
reuseable <- rep(NA, nmodels)
}
if (analysisType == "RM-ANOVA") {
# the default null-model contains subject
anyNuisance <- TRUE
} else {
# the default null-model is intercept only
anyNuisance <- length(nuisance) > 0L
}
# without these there is no error
bfIterations <- if (options[["samplingMethodNumericAccuracy"]] == "auto") 1e4L else options[["samplesNumericAccuracy"]]
if (options[["integrationMethod"]] == "automatic") {
bfIntegrationMethod <- "auto"
} else {
bfIntegrationMethod <- "laplace"
modelTable$addFootnote(gettext("The Laplace approximation is faster but also less accurate than the automatic method. This is fine when exploring the data and obtaining a rough estimate. However, we recommend using the automatic integration method to obtain the final results."))
modelTable$addFootnote(gettext("The Laplace approximation does not always provide an error estimate."), colNames = "error %")
}
for (m in seq_len(nmodels)) {
# loop over all models, where the first is the null-model and the last the most complex model
if (is.na(reuseable[m])) {
if (!is.null(model.list[[m]])) {
.setSeedJASP(options)
bf <- try(BayesFactor::lmBF(
formula = model.list[[m]],
data = dataset,
whichRandom = randomFactors,
progress = FALSE,
posterior = FALSE,
rscaleFixed = rscaleFixed,
rscaleRandom = rscaleRandom,
rscaleCont = rscaleCont,
rscaleEffects = rscaleEffects,
iterations = bfIterations,
method = bfIntegrationMethod
))
if (isTryError(bf)) {
modelName <- strsplit(.BANOVAas.character.formula(model.list[[m]]), "~ ")[[1L]][2L]
.quitAnalysis(gettextf("Bayes factor could not be computed for model: %1$s.\nThe error message was: %2$s.",
.BANOVAdecodeSubject(modelName), .extractErrorMessage(bf)))
} else {
# delete the data object -- otherwise it gets saved in the state
bf@data <- data.frame()
}
} else {
bf <- NULL
}
} else {
bf <- stateObj$models[[reuseable[m]]]$bf
}
modelObject[[m]]$bf <- bf
if (!is.null(bf)) {
# if (isTryError(bf)) {
# message <- .extractErrorMessage(bf)
# modelObject[[m]]$error.message <- paste("Bayes factor computation failed with error:", message)
# } else {
# bfObj can have a modified denominator, but the saved objects are always compared against intercept only
bfObj <- modelObject[[m]]$bf
if (anyNuisance)
bfObj <- bfObj / modelObject[[1L]]$bf
internalTable[m, "BF10"] <- bfObj@bayesFactor[, "bf"] # always LogBF10!
internalTable[m, "error %"] <- bfObj@bayesFactor[, "error"] # always NA if bfIntegrationMethod == "laplace"
# }
# disable filling of Bayes factor column (see https://github.com/jasp-stats/jasp-test-release/issues/1018 for discussion)
# modelTable[["BF10"]] <- .recodeBFtype(internalTable[["BF10"]],
# newBFtype = options$bayesFactorType,
# oldBFtype = "LogBF10")
modelTable[["error %"]] <- internalTable[seq_len(shownTableSize), "error %"]
}
progressbarTick()
}
internalTableObj <- .BANOVAfinalizeInternalTable(options, internalTable)
modelTable$setData(internalTableObj$table[seq_len(shownTableSize), ])
if (length(internalTableObj[["footnotes"]]) > 0L) {
idxRow <- internalTableObj[["footnotes"]][["rows"]]
idxCol <- internalTableObj[["footnotes"]][["cols"]]
tb <- internalTableObj$table
modelTable$addFootnote(message = internalTableObj[["footnotes"]][["message"]],
symbol = " <sup>✝</sup>", # latin cross
rowNames = rownames(tb)[idxRow],
colNames = colnames(tb)[idxCol])
}
if (anyNuisance) {
if (analysisType == "RM-ANOVA" && !options[["legacyResults"]]) {
message <- gettextf("All models include %s, and random slopes for all repeated measures factors.",
paste0(.BANOVAdecodeNuisance(setdiff(nuisance, nuisanceRandomSlopes)), collapse = ", "))
} else {
if (analysisType == "RM-ANOVA")
modelTable$addFootnote(message = gettext("Random slopes for repeated measures factors are excluded."), symbol = .BANOVAGetWarningSymbol())
message <- gettextf("All models include %s.", paste0(.BANOVAdecodeNuisance(nuisance), collapse = ", "))
}
modelTable$addFootnote(message = message)
}
model <- list(
models = modelObject,
postProbs = internalTableObj$internalTable[["P(M|data)"]],
priorProbs = internalTableObj$internalTable[["P(M)"]],
internalTableObj = internalTableObj,
effects = effects.matrix,
interactions.matrix = interactions.matrix,
nuisance = nuisance,
nuisanceRandomSlopes = nuisanceRandomSlopes,
analysisType = analysisType,
# these are necessary for partial reusage of the state (e.g., when a fixedFactor is added/ removed)
model.list = model.list,
fixedFactors = fixedFactors,
randomFactors = if (analysisType == "RM-ANOVA") randomFactors else options[["randomFactors"]], # stored because they are modified in RM-ANOVA
modelTerms = modelTerms,
covariates = if (analysisType == "ANOVA") NULL else options[["covariates"]],
seed = options[["seed"]],
setSeed = options[["setSeed"]],
reuseable = reuseable,
RMFactors = options[["repeatedMeasuresFactors"]],
modelPriorOptions = options[.BANOVAmodelSpaceDependencies(options[["modelPrior"]])],
hideNuisanceEffects = options[["hideNuisanceEffects"]],
shownTableSize = options[c("modelsShown", "numModelsShown")]
)
# save state
stateObj <- createJaspState(object = model, dependencies = c(
# does NOT depend on any factors or covariates, to facilitate reusing previous models
"dependent", "repeatedMeasuresCells", "samplingMethodNumericAccuracy", "samplesNumericAccuracy", "seed", "setSeed",
.BANOVArscaleDependencies(options[["priorSpecificationMode"]])
))
jaspResults[["tableModelComparisonState"]] <- stateObj
return(model)
}
.BANOVAeffectsTable <- function(jaspResults, options, model) {
if (!is.null(jaspResults[["tableEffects"]]) || !options[["effects"]])
return()
# isTRUE should handle a state issue, see https://github.com/jasp-stats/jasp-test-release/issues/839
if (isTRUE(model[["analysisType"]] != "RM-ANOVA" && options[["dependent"]] != "")) {
title <- gettextf("Analysis of Effects - %s", options[["dependent"]])
} else {
title <- gettext("Analysis of Effects")
}
effectsTable <- createJaspTable(title = title)
effectsTable$position <- 2
effectsTable$dependOn(c(
.BANOVAdataDependencies(),
"effects", "effectsType",
"samplingMethodNumericAccuracy", "samplesNumericAccuracy", "bayesFactorType",
"integrationMethod",
.BANOVAmodelSpaceDependencies(options[["modelPrior"]]),
.BANOVArscaleDependencies(options[["priorSpecificationMode"]])
))
effectsTable$addCitation(.BANOVAcitations[1:2])
inclusion.title <- switch(
options$bayesFactorType,
"LogBF10" = gettext("Log(BF<sub>incl</sub>)"),
"BF01" = gettext("BF<sub>excl</sub>"),
"BF10" = gettext("BF<sub>incl</sub>")
)
effectsTable$addColumnInfo(name = "Effects", title = gettext("Effects"), type = "string")
effectsTable$addColumnInfo(name = "P(incl)", title = gettext("P(incl)"), type = "number")
effectsTable$addColumnInfo(name = "P(excl)", title = gettext("P(excl)"), type = "number")
effectsTable$addColumnInfo(name = "P(incl|data)", title = gettext("P(incl|data)"), type = "number")
effectsTable$addColumnInfo(name = "P(excl|data)", title = gettext("P(excl|data)"), type = "number")
effectsTable$addColumnInfo(name = "BFInclusion", title = inclusion.title, type = "number")
if (options$effectsType == "matchedModels") {
effectsTable$addFootnote(gettextf("Compares models that contain the effect to equivalent models stripped of the effect. Higher-order interactions are excluded. Analysis suggested by Sebastiaan Math%st.", "\u00F4"))
}
if (is.null(model$models)) {
jaspResults[["tableEffects"]] <- effectsTable
return()
}
effects.matrix <- model$effects
no.effects <- ncol(effects.matrix)
effectNames <- colnames(model$effects)
idxNotNan <- c(TRUE, !is.nan(model$postProbs[-1L]))
if (options$effectsType == "allModels") {
if (any(!idxNotNan)) {
# renormalize the prior and posterior probabilities
priorInclProb <- crossprod(effects.matrix[idxNotNan, , drop = FALSE], model[["priorProbs"]][idxNotNan] / sum(model[["priorProbs"]][idxNotNan]))
postInclProb <- crossprod(effects.matrix[idxNotNan, , drop = FALSE], model[["postProbs"]][idxNotNan] / sum(model[["postProbs"]][idxNotNan]))
} else {
# do not renormalize
priorInclProb <- crossprod(effects.matrix, model[["priorProbs"]])
postInclProb <- crossprod(effects.matrix, model[["postProbs"]])
}
# deal with numerical error
postInclProb[postInclProb > 1] <- 1
postInclProb[postInclProb < 0] <- 0
bfIncl <- (postInclProb / (1 - postInclProb)) / (priorInclProb / (1 - priorInclProb))
priorExclProb <- 1 - priorInclProb
postExclProb <- 1 - postInclProb
} else {
tmp <- BANOVAcomputMatchedInclusion(
effectNames, effects.matrix, model$interactions.matrix,
model$priorProbs, model$postProbs
)
priorInclProb <- tmp[["priorInclProb"]]
postInclProb <- tmp[["postInclProb"]]
bfIncl <- tmp[["bfIncl"]]
priorExclProb <- tmp[["priorExclProb"]]
postExclProb <- tmp[["postExclProb"]]
# show BFinclusion for nuisance predictors as 1, rather than NaN
priorInclIs1 <- is.nan(bfIncl) & abs(1 - priorInclProb) <= sqrt(.Machine$double.eps)
bfIncl[priorInclIs1] <- 1
}
if (sum(!idxNotNan) > 1L) { # null model is always omitted, so 2 or more omitted indicates some models failed
effectsTable$addFootnote(message = gettext("Some Bayes factors could not be calculated. Inclusion probabilities and Bayes factors are calculated while excluding these models. The results may be uninterpretable!"),
symbol = .BANOVAGetWarningSymbol()
)
}
effectsTable[["Effects"]] <- jaspBase::gsubInteractionSymbol(effectNames)
effectsTable[["P(incl)"]] <- priorInclProb
effectsTable[["P(incl|data)"]] <- postInclProb
effectsTable[["BFInclusion"]] <- switch(
options$bayesFactorType,
"LogBF10" = log(bfIncl),
"BF01" = 1 / bfIncl,
"BF10" = bfIncl
)
effectsTable[["P(excl)"]] <- priorExclProb
effectsTable[["P(excl|data)"]] <- postExclProb
jaspResults[["tableEffects"]] <- effectsTable
return()
}
BANOVAcomputMatchedInclusion <- function(effectNames, effects.matrix, interactions.matrix,
priorProbs, postProbs) {
# this method is inspired by this post: https://www.cogsci.nl/blog/interpreting-bayesian-repeated-measures-in-jasp
priorInclProb <- postInclProb <- bfIncl <- priorExclProb <- postExclProb <-
numeric(length(effectNames))
for (i in seq_along(effectNames)) {
effect <- effectNames[i]
# get all higher order interactions of which this effect is a component
# e.g., V1 is a component of V1:V2
idx1 <- interactions.matrix[effect, ]
# get all models that exclude the predictor, but that always include the lower order main effects
# e.g., V1:V2 is compared against models that always include V1 and V2
idx2 <- interactions.matrix[, effect] # if effect is V1:V2, idx2 contains c(V1, V2)
# idx3 is FALSE if a model contains higher order interactions of effect, TRUE otherwise
idx3 <- !matrixStats::rowAnys(effects.matrix[, idx1, drop = FALSE])
# all models that include the effect, without higher order interactions
idx4 <- idx3 & effects.matrix[, i]
priorInclProb[i] <- sum(idx4 * priorProbs)
postInclProb[i] <- sum(idx4 * postProbs)
# the models to consider for the prior/ posterior exclusion probability.
# idx5 includes models that have: all subcomponents & no higher order interaction & not the effect
idx5 <- matrixStats::rowAlls(effects.matrix[, idx2, drop = FALSE]) & idx3 & !effects.matrix[, i]
priorExclProb[i] <- sum(idx5 * priorProbs)
postExclProb[i] <- sum(idx5 * postProbs)
# compute inclusion BF
bfIncl[i] <- (postInclProb[i] / postExclProb[i]) / (priorInclProb[i] / priorExclProb[i])
}
return(list(
priorInclProb = priorInclProb,
postInclProb = postInclProb,
priorExclProb = priorExclProb,
postExclProb = postExclProb,
bfIncl = bfIncl
))
}
.BANOVAinitModelComparisonTable <- function(options) {
# function that creates an empty JASP table to be filled later
modelTable <- createJaspTable(title = gettext("Model Comparison"))
modelTable$position <- 1L
modelTable$addCitation(.BANOVAcitations[1:2])
modelTable$dependOn(c(
.BANOVAdataDependencies(),
"samplingMethodNumericAccuracy", "samplesNumericAccuracy", "integrationMethod",
"bayesFactorType", "bayesFactorOrder",
"hideNuisanceParameters", "legacyResults",
"modelsShown",
.BANOVAmodelSpaceDependencies(options[["modelPrior"]]),
.BANOVArscaleDependencies(options[["priorSpecificationMode"]])
))
if (options[["modelsShown"]] == "limited")
modelTable$dependOn("numModelsShown")
switch(options$bayesFactorType,
BF10 = {
bfm.title <- gettext("BF<sub>M</sub>")
bf.title <- gettext("BF<sub>10</sub>")
},
BF01 = {
bfm.title <- gettext("BF<sub>M</sub>")
bf.title <- gettext("BF<sub>01</sub>")
},
LogBF10 = {
bfm.title <- gettext("Log(BF<sub>M</sub>)")
bf.title <- gettext("Log(BF<sub>10</sub>)")
}
)
modelTable$addColumnInfo(name = "Models", title = gettext("Models"), type = "string")
modelTable$addColumnInfo(name = "P(M)", title = gettext("P(M)"), type = "number")
modelTable$addColumnInfo(name = "P(M|data)", title = gettext("P(M|data)"), type = "number")
modelTable$addColumnInfo(name = "BFM", title = bfm.title, type = "number")
modelTable$addColumnInfo(name = "BF10", title = bf.title, type = "number")
modelTable$addColumnInfo(name = "error %", title = gettextf("error %%"), type = "number")
return(modelTable)
}
.BANOVAfinalizeInternalTable <- function(options, internalTable) {
# function that actually fills in the table created by .BANOVAinitModelComparisonTable
footnotes <- NULL
logSumExp <- matrixStats::logSumExp
logbfs <- internalTable[["BF10"]]
logprior <- log(internalTable[["P(M)"]])
if (!anyNA(internalTable[["BF10"]])) {
# no errors, proceed normally and complete the table
logsumbfs <- logSumExp(logbfs + logprior)
internalTable[["P(M|data)"]] <- exp(logbfs + logprior - logsumbfs)
nmodels <- nrow(internalTable)
for (i in seq_len(nmodels)) {
internalTable[i, "BFM"] <- logbfs[i] - logSumExp(logbfs[-i]) + log(nmodels - 1L)
}
} else {
# create table excluding failed models
idxGood <- !is.na(logbfs)
widxGood <- which(idxGood)
logsumbfs <- logSumExp(logbfs[idxGood])
internalTable[["P(M|data)"]] <- exp(logbfs - logsumbfs)
nmodels <- sum(idxGood)
widxBad <- which(!idxGood)
for (i in widxGood) {
internalTable[["BFM"]][i] <- logbfs[i] - logSumExp(logbfs[-c(i, widxBad)]) + log(nmodels - 1L)
}
internalTable[widxBad, "P(M|data)"] <- NaN
internalTable[widxBad, "BFM"] <- NaN
internalTable[widxBad, "BF10"] <- NaN
footnotes <- list(
message = gettext("<b><em>Warning.</em></b> Some Bayes factors could not be calculated. Multi model inference is carried out while excluding these models. The results may be uninterpretable!")
)
}
# create the output table
table <- as.data.frame(internalTable)
if (options[["bayesFactorType"]] == "LogBF10") {
table[["BFM"]] <- internalTable[["BFM"]]
} else {
table[["BFM"]] <- exp(internalTable[["BFM"]])
}
o <- order(table[["BF10"]], decreasing = TRUE)
table <- table[o, ]
idxNull <- which(o == 1L)
if (options[["bayesFactorOrder"]] == "nullModelTop") {
table[idxNull, "error %"] <- NA
table <- table[c(idxNull, seq_len(nrow(table))[-idxNull]), ]
} else {
table[["BF10"]] <- table[["BF10"]] - table[1L, "BF10"]
# recompute error (see BayesFactor:::`.__T__/:base`$`BFBayesFactor#BFBayesFactor`)
table[idxNull, "error %"] <- 0
table[["error %"]] <- sqrt(table[["error %"]]^2 + table[["error %"]][1L]^2)
table[1L, "error %"] <- NA
}
table[["BF10"]] <- .recodeBFtype(table[["BF10"]], newBFtype = options[["bayesFactorType"]], oldBFtype = "LogBF10")
table[["error %"]] <- 100 * table[["error %"]]
if (!is.null(footnotes)) {
idxNan <- which(do.call(cbind, lapply(table[-ncol(table)], is.nan)), arr.ind = TRUE)
footnotes[["rows"]] <- idxNan[, "row"]
footnotes[["cols"]] <- idxNan[, "col"]
}
return(list(table = table, internalTable = internalTable, footnotes = footnotes))
}
# posterior inference ----
.BANOVAestimatePosteriors <- function(jaspResults, dataset, options, model) {
userNeedsPosteriorSamples <- options$posteriorEstimates || options$modelAveragedPosteriorPlot || options$qqPlot ||
options$rsqPlot || options$criTable
if (is.null(model$models) || !userNeedsPosteriorSamples)
return()
# model[["completelyReused"]] is needed because it can happen that some posterior samples can be reused (e.g.,
# when the modelTerms change)
posteriors <- jaspResults[["statePosteriors"]]$object
if (!is.null(posteriors) && isTRUE(model[["completelyReused"]])) { # can all posterior be reused?
posteriorsCRI <- jaspResults[["statePosteriorsCRI"]]$object
if (is.null(posteriorsCRI)) {# do we need to recompute credible intervals?
posteriorsCRI <- .BANOVAcomputePosteriorCRI(dataset, options, model, posteriors)
statePosteriorsCRI <- createJaspState(object = posteriorsCRI)
statePosteriorsCRI$dependOn(options = "credibleInterval",
optionsFromObject = jaspResults[["tableModelComparisonState"]])
jaspResults[["statePosteriorCRI"]] <- statePosteriorsCRI
}
} else { # compute posteriors and credible intervals
posteriors <- .BANOVAsamplePosteriors(jaspResults, dataset, options, model, posteriors)
posteriorsCRI <- .BANOVAcomputePosteriorCRI(dataset, options, model, posteriors)
statePosteriors <- createJaspState(object = posteriors)
statePosteriors$dependOn(
c(.BANOVAmodelSpaceDependencies(options[["modelPrior"]]), "samplingMethodMCMC", "samplesMCMC"),
optionsFromObject = jaspResults[["tableModelComparisonState"]]
)
jaspResults[["statePosteriors"]] <- statePosteriors
statePosteriorsCRI <- createJaspState(object = posteriorsCRI)
statePosteriorsCRI$dependOn(options = "credibleInterval",
optionsFromObject = jaspResults[["statePosteriors"]])
jaspResults[["statePosteriorCRI"]] <- statePosteriorsCRI
}
.BANOVASamplingIssuesTable(jaspResults, posteriors[["samplingIssues"]])
return(c(posteriors, posteriorsCRI))
}
.BANOVAestimatesTable <- function(jaspResults, options, model) {
if (!is.null(jaspResults[["tablePosteriorEstimates"]]) || !options[["posteriorEstimates"]])
return()
estsTable <- createJaspTable(title = gettext("Model Averaged Posterior Summary"))
estsTable$position <- 3
estsTable$dependOn(c(
.BANOVAdataDependencies(),
"posteriorEstimates",
"samplingMethodMCMC", "samplesMCMC", "credibleInterval",
.BANOVAmodelSpaceDependencies(options[["modelPrior"]]),
.BANOVArscaleDependencies(options[["priorSpecificationMode"]])
))
overTitle <- gettextf("%s%% Credible Interval", format(100 * options[["credibleInterval"]], digits = 3))
estsTable$addColumnInfo(name = "Variable", type = "string")
estsTable$addColumnInfo(name = "Level", type = "string")
estsTable$addColumnInfo(name = "Mean", type = "number")
estsTable$addColumnInfo(name = "SD", type = "number")
estsTable$addColumnInfo(name = "Lower", type = "number", overtitle = overTitle)
estsTable$addColumnInfo(name = "Upper", type = "number", overtitle = overTitle)
if (!is.null(model[["posteriors"]])) {
.BANOVAfillEstimatesTable(
jaspTable = estsTable,
mus = model$posteriors$weightedMeans,
sds = model$posteriors$weightedSds,
cri = model$posteriors$weightedCRIs,
hasNoLevels = model$posteriors$allContinuous,
isRandom = model$posteriors$isRandom
)
}
jaspResults[["tablePosteriorEstimates"]] <- estsTable
return()
}
.BANOVAfillEstimatesTable <- function(jaspTable, mus, sds, cri, hasNoLevels, isRandom = NULL) {
if (!is.null(isRandom) && any(isRandom)) {
# remove random effects
mus <- mus[!isRandom]
sds <- sds[!isRandom]
cri <- cri[!isRandom, ]
}
table <- .BANOVAgetLevelsFromParamNames(names(mus))
colnames(table) <- c("Variable", "Level")
# set repeated parameter names to ""
idxDup <- duplicated(table[, "Variable"])
table[idxDup, "Variable"] <- ""
# decode base64 variables
table[!idxDup, "Variable"] <- jaspBase::gsubInteractionSymbol(table[!idxDup, "Variable"])
# rename mu to Intercept
table[1L, "Variable"] <- "Intercept"
# attach posterior means and sds
table <- cbind(as.data.frame(table),
"Mean" = mus,
"SD" = sds,
"Lower" = cri[, 1L],
"Upper" = cri[, 2L])
if (hasNoLevels)
table <- table[, -2L]
jaspTable$setData(table)
return()
}
.BANOVArsquaredTable <- function(jaspResults, options, model) {
if (!is.null(jaspResults[["tableBMACRI"]]) || !options[["criTable"]])
return()
criTable <- createJaspTable(title = gettextf("Model Averaged R%s", "\u00B2"))
criTable$position <- 3.5
criTable$dependOn(c(
.BANOVAdataDependencies(),
"criTable",
"samplingMethodMCMC", "samplesMCMC", "bayesFactorType", "credibleInterval",
.BANOVAmodelSpaceDependencies(options[["modelPrior"]]),
.BANOVArscaleDependencies(options[["priorSpecificationMode"]])
))
overTitle <- gettextf("%s%% Credible Interval", format(100 * options[["credibleInterval"]], digits = 3))
criTable$addColumnInfo(name = "rsq", type = "string", title = "")
criTable$addColumnInfo(name = "Mean", type = "number")
criTable$addColumnInfo(name = "Lower", type = "number", overtitle = overTitle)
criTable$addColumnInfo(name = "Upper", type = "number", overtitle = overTitle)
if (!is.null(model[["posteriors"]])) {
cri <- model[["posteriors"]][["weightedRsqCri"]]
df <- data.frame(
rsq = "R\u00B2",
Mean = model[["posteriors"]][["weightedRsqMean"]],
Lower = cri[1L],
Upper = cri[2L]
)
criTable$setData(df)
} else {
criTable[["rsq"]] <- "R\u00B2"
}
jaspResults[["tableBMACRI"]] <- criTable
return()
}
.BANOVASamplingIssuesTable <- function(jaspResults, samplingIssues) {
if (is.null(samplingIssues) || length(samplingIssues) == 0L || !is.null(jaspResults[["tableSamplingIssues"]]))
return()
issuesTable <- createJaspTable(title = gettext("<b><em>Warning: sampling issues encountered!</em></b>"),
# below model and effects table (which may be unaffected) but above estimates table
position = 2.5)
issuesTable$addColumnInfo(name = "Models", title = gettext("Affected model"), type = "string")
issuesTable$addColumnInfo(name = "percentageSucces", title = gettextf("%% useful samples"), type = "number")
issuesTable$addColumnInfo(name = "noSuccess", title = gettext("Useful samples"), type = "integer")
issuesTable$addColumnInfo(name = "total", title = gettext("Total samples drawn"), type = "integer")
df <- data.frame(
Models = vapply(samplingIssues, `[[`, FUN.VALUE = character(1L), "model"),
noSuccess = vapply(samplingIssues, `[[`, FUN.VALUE = integer(1L), "remainingRows"),
total = vapply(samplingIssues, `[[`, FUN.VALUE = integer(1L), "originalRows")
)
df[["percentageSucces"]] <- df[["noSuccess"]] / df[["total"]]
if (any(df[["percentageSucces"]] < 0.25) || any(df[["remainingRows"]] < 1000L))
issuesTable$addFootnote(
gettextf("For some affected models, more than 75%% of the posterior samples failed, or fewer than 1000 samples remained for subsequent results. All model-averaged output may be biased and uninterpretable. Check the model specification and data for any odd patterns."),
symbol = .BANOVAGetWarningSymbol()
)
issuesTable$setData(df)
issuesTable$dependOn(
# these options correspond to userNeedsPosteriorSamples inside .BANOVAestimatePosteriors
options = c("posteriorEstimates", "modelAveragedPosteriorPlot", "qqPlot", "rsqPlot", "criTable", "modelTerms"),
optionsFromObject = jaspResults[["statePosteriors"]]
)
jaspResults[["tableSamplingIssues"]] <- issuesTable
}
.BANOVAGetWarningSymbol <- function() {
gettext("<b><em>Warning.</em></b>")
}
# Plots wrappers ----
.BANOVAposteriorPlot <- function(jaspResults, dataset, options, model) {
# meta wrapper for model averaged posterior plots
if (!is.null(jaspResults[["posteriorPlot"]]) || !options$modelAveragedPosteriorPlot)
return()
posteriorPlotContainer <- createJaspContainer(title = gettext("Model Averaged Posterior Distributions"))
jaspResults[["posteriorPlot"]] <- posteriorPlotContainer
posteriorPlotContainer$dependOn(c("modelAveragedPosteriorPlot", "modelTerms", "credibleInterval", "groupPosterior",
"repeatedMeasuresCells", "dependent"))
posteriorPlotContainer$position <- 4
if (is.null(model$models)) {
posteriorPlotContainer[["dummyplot"]] <- createJaspPlot(title = gettext("Posterior distribution"), width = 400, height = 400,
plot = NULL)
} else {
posteriorPlotContainer$dependOn(optionsFromObject = jaspResults[["tableModelComparisonState"]])
.BANOVAfillPosteriorPlotContainer(
container = posteriorPlotContainer,
densities = model$posterior$weightedDensities[, -1L, ], # omit intercept
cris = model$posterior$weightedCRIs[-1L, ], # omit intercept
isRandom = model$posteriors$isRandom[-1L], # omit intercept
groupParameters = identical(options[["groupPosterior"]], "grouped")
)
}
return()
}
.BANOVAfillPosteriorPlotContainer <- function(container, densities, cris, isRandom = NULL, groupParameters = FALSE) {
allParamNames <- colnames(densities)
tmp <- .BANOVAgetLevelsFromParamNames(allParamNames)
plotTitles <- jaspBase::gsubInteractionSymbol(tmp[, "parameter"])
xNames <- tmp[, "level"]
if (is.null(isRandom)) {
indices <- seq_along(allParamNames)
} else {
indices <- which(!isRandom)
}
if (groupParameters) {
indices <- split(indices, plotTitles[indices])[order(unique(plotTitles[indices]))]
nms <- names(indices)
for (i in seq_along(indices)) {
ind <- indices[[i]]
if (all(c(densities[, ind, "y"]) == 0)) { # only TRUE is never sampled and set to 0
plot <- createJaspPlot(title = nms[i], width = 400, height = 400)
plot$setError(gettext("Could not plot posterior distribution. Check if an error occurred for models including this parameter."))
} else {
# make prior posterior plot
dfLines <- data.frame(x = c(densities[, ind, "x"]),
y = c(densities[, ind, "y"]),
g = rep(xNames[ind], each = nrow(densities)))
xBreaks <- jaspGraphs::getPrettyAxisBreaks(dfLines$x)
yBreaks <- jaspGraphs::getPrettyAxisBreaks(c(0, dfLines$y))
breaksYmax <- yBreaks[length(yBreaks)]
obsYmax <- max(dfLines$y)
newymax <- max(1.25 * obsYmax, breaksYmax)
lInd <- length(ind)
if (lInd == 1L) {
plusmin <- 0
showLegend <- FALSE
} else {
# some jitter in CRI height
showLegend <- TRUE
plusmin <- rep(-1, lInd)
plusmin[seq(2, lInd, 2)] <- 0
if (lInd > 2)
plusmin[seq(3, lInd, 3)] <- 1
}
dfCri <- data.frame(
xmin = cris[ind, 1L],
xmax = cris[ind, 2L],
y = (newymax - obsYmax)/2 + obsYmax + plusmin * (newymax - obsYmax) / 10,
g = xNames[ind]
)
if (showLegend) {
# if two distributions are remarkably similar, i.e., have nearly identical credible intervals,
# we add a linetype aestethic
tol <- 1e-4
distMin <- which(stats::dist(dfCri$xmin) < tol)
distMax <- which(stats::dist(dfCri$xmin) < tol)
distBoth <- intersect(distMin, distMax)
if (length(distBoth) > 0) {
aesLine <- ggplot2::aes(x = x, y = y, g = g, color = g, linetype = g)
aesErrorbar <- ggplot2::aes(xmin = xmin, xmax = xmax, y = y, group = g, color = g, linetype = g)
} else {
aesLine <- ggplot2::aes(x = x, y = y, g = g, color = g)
aesErrorbar <- ggplot2::aes(xmin = xmin, xmax = xmax, y = y, group = g, color = g)
}
} else {
# don't use color for single density plots (covariates)
aesLine <- ggplot2::aes(x = x, y = y, g = g)
aesErrorbar <- ggplot2::aes(xmin = xmin, xmax = xmax, y = y, group = g)
}
maxheight <- min(newymax - dfCri$y[1:min(lInd, 3)])
xlab <- nms[i]
# ggplot doesn't like our fancy unicode * so we need to escape it
xlab <- stringi::stri_escape_unicode(xlab)
# change escaped version into x
xlab <- gsub(pattern = "\\u2009\\u273b\\u2009", replacement = " x ", x = xlab, fixed = TRUE)
ncolLegend <- ceiling(lInd / 14)
guideLegend <- ggplot2::guide_legend(title = gettext("Level"), keywidth = 0.25, keyheight = 0.1, default.unit = "inch",
ncol = ncolLegend)
p <- ggplot2::ggplot(data = dfLines, mapping = aesLine) +
ggplot2::geom_line(size = 1.1) +
ggplot2::geom_errorbarh(data = dfCri, mapping = aesErrorbar, height = maxheight, size = 1.1,
inherit.aes = FALSE) +
ggplot2::scale_x_continuous(name = xlab, breaks = xBreaks, limits = range(xBreaks)) +
ggplot2::scale_y_continuous(name = gettext("Density"), breaks = yBreaks, limits = c(0, newymax)) +
colorspace::scale_color_discrete_qualitative() +
ggplot2::scale_linetype() +
ggplot2::guides(color = guideLegend, linetype = guideLegend)
p <- jaspGraphs::themeJasp(p, legend.position = if (showLegend) "right" else "none")
plot <- createJaspPlot(title = nms[i], width = 400, height = 400, plot = p)
}
container[[allParamNames[i]]] <- plot
}
} else {
for (i in indices) {
# make prior posterior plot
df <- data.frame(x = densities[, i, "x"],
y = densities[, i, "y"])
p <- jaspGraphs::PlotPriorAndPosterior(
dfLines = df,
xName = xNames[i],
CRI = cris[i, ],
drawCRItxt = FALSE
)
plot <- createJaspPlot(title = plotTitles[i], width = 400, height = 400, plot = p)
container[[allParamNames[i]]] <- plot
}
}
return()
}
.BANOVAqqplot <- function(jaspResults, options, model) {
if (!is.null(jaspResults[["QQplot"]]) || !options[["qqPlot"]])
return()
plot <- createJaspPlot(
title = gettext("Model Averaged Q-Q Plot"),
width = 400,
height = 400,
aspectRatio = 1
)
if (!is.null(model[["models"]])) {
plot$plotObject <- jaspGraphs::plotQQnorm(
residuals = model[["posteriors"]][["weightedResidSumStats"]][,"mean"],
lower = model[["posteriors"]][["weightedResidSumStats"]][,"cri.2.5%"],
upper = model[["posteriors"]][["weightedResidSumStats"]][,"cri.97.5%"],
ablineColor = "darkred"
)
plot$dependOn(optionsFromObject = jaspResults[["tableModelComparisonState"]])
}
plot$dependOn(c("qqPlot", "modelTerms", "credibleInterval", "repeatedMeasuresCells", "dependent"))
plot$position <- 5
jaspResults[["QQplot"]] <- plot
return()
}
.BANOVArsqplot <- function(jaspResults, options, model) {
if (!is.null(jaspResults[["rsqplot"]]) || !options[["rsqPlot"]])
return()
plot <- createJaspPlot(
title = gettextf("Model Averaged Posterior R%s","\u00B2"),
width = 400,
height = 400,
aspectRatio = 1
)
if (!is.null(model[["models"]])) {
dd <- model[["posteriors"]][["weightedRsqDens"]]
rsqCri <- model[["posteriors"]][["weightedRsqCri"]]
df <- data.frame(x = dd$x, y = dd$y)
xName <- expression(R^2)
plot$plotObject <- jaspGraphs::PlotPriorAndPosterior(dfLines = df, xName = xName, CRI = rsqCri, drawCRItxt = FALSE)
plot$dependOn(optionsFromObject = jaspResults[["tableModelComparisonState"]])
}
plot$dependOn(c("rsqPlot", "modelTerms", "credibleInterval", "repeatedMeasuresCells", "dependent"))
plot$position <- 6
jaspResults[["rsqplot"]] <- plot
return()
}
# Post hoc comparison ----
.BANOVAnullControlPostHocTable <- function(jaspResults, dataset, options, model) {
if (length(options$postHocTerms) == 0L)
return()
postHocCollection <- jaspResults[["collectionPosthoc"]]
if (is.null(postHocCollection)) {
postHocCollection <- createJaspContainer(title = gettext("Post Hoc Tests"))
postHocCollection$position <- 8
postHocCollection$addCitation(.BANOVAcitations[3:4])
postHocCollection$dependOn(c("dependent", "repeatedMeasuresCells", "postHocNullControl", "bayesFactorType"))
jaspResults[["collectionPosthoc"]] <- postHocCollection
}
# the same footnote for all the tables
footnote <- gsub("[\r\n\t]", "",
gettext("The posterior odds have been corrected for multiple testing by fixing to 0.5 the prior probability that the null hypothesis holds across all comparisons (Westfall, Johnson, & Utts, 1997). Individual comparisons are based on the default t-test with a Cauchy (0, r = 1/sqrt(2)) prior. The \"U\" in the Bayes factor denotes that it is uncorrected."))
bfTxt <- if (options[["postHocNullControl"]]) ", U" else ""
switch(options[["bayesFactorType"]],
BF10 = { bf.title <- paste0("BF<sub>10", bfTxt, "</sub>") },
BF01 = { bf.title <- paste0("BF<sub>01", bfTxt, "</sub>") },
LogBF10 = { bf.title <- paste0("Log(BF<sub>10", bfTxt, "</sub>)") }
)
priorWidth <- 1 / sqrt(2)
posthoc.variables <- unlist(options[["postHocTerms"]])
if (model[["analysisType"]] == "RM-ANOVA") {
dependent <- .BANOVAdependentName
} else {
dependent <- options[["dependent"]]
}
for (posthoc.var in posthoc.variables) {
# does the table already exist?
if (!is.null(postHocCollection[[paste0("postHoc_", posthoc.var)]]))
next
postHocTable <- createJaspTable(title = paste0("Post Hoc Comparisons - ", posthoc.var))
postHocTable$addColumnInfo(name = "(I)", type = "string", title = "", combine=TRUE)
postHocTable$addColumnInfo(name = "(J)", type = "string", title = "")
postHocTable$addColumnInfo(name = "Prior Odds", type = "number", title = gettext("Prior Odds"))
postHocTable$addColumnInfo(name = "Posterior Odds", type = "number", title = gettext("Posterior Odds"))
postHocTable$addColumnInfo(name = "BF", type = "number", title = bf.title)
postHocTable$addColumnInfo(name = "error %", type = "number", title = gettextf("error %%"))
postHocTable$dependOn(optionContainsValue = list("postHocTerms" = posthoc.var))
if (options[["postHocNullControl"]])
postHocTable$addFootnote(footnote)
if (is.null(model$models)) { # only show empty table
postHocCollection[[paste0("postHoc_", posthoc.var)]] <- postHocTable
next
}
fixed <- unlist(c(options$fixedFactors, sapply(options$repeatedMeasuresFactors, `[[`, "name")))
if (model$analysisType == "RM-ANOVA" && posthoc.var %in% fixed && !posthoc.var %in% options$betweenSubjectFactors) {
variable.levels <- options$repeatedMeasuresFactors[[which(lapply(options$repeatedMeasuresFactors, function(x) x$name) == posthoc.var)]]$levels
paired <- TRUE
} else if (posthoc.var %in% c(options$fixedFactors, options$betweenSubjectFactors, options$randomFactors)) {
variable.levels <- levels(dataset[[posthoc.var]])
paired <- FALSE
} else {
next
}
if (length(variable.levels) < 2L)
next
pairs <- utils::combn(variable.levels, 2)
splitname <- if (model[["analysisType"]] == "RM-ANOVA") .BANOVAdependentName else dependent
allSplits <- split(dataset[[splitname]], dataset[[posthoc.var]])
errMessages <- NULL
for (i in 1:ncol(pairs)) {
if (!is.null(model$models)) {
x <- na.omit(allSplits[[pairs[1L, i]]])
y <- na.omit(allSplits[[pairs[2L, i]]])
ttest <- try(BayesFactor::ttestBF(x = x, y = y, rscale = priorWidth, paired = paired), silent=TRUE)
if (isTryError(ttest)) {
priorOdds <- postOdds <- logBF <- error <- "NaN"
errorMsg <- .extractErrorMessage(ttest)
if (errorMsg %in% names(errMessages)) {
errMessages[[errorMsg]][["row_names"]] <- c(errMessages[[errorMsg]][["row_names"]], i)
} else {
errMessages[[errorMsg]] <- list(
message = if (endsWith(errorMsg, ".")) errorMsg else paste0(errorMsg, "."),
row_names = i
)
}
} else {
if (options[["postHocNullControl"]]) {
pH0 <- 0.5^(2 / length(variable.levels))
} else {
pH0 <- 0.5
}
logBF <- ttest@bayesFactor$bf # <- log(BF10)
if (options[["bayesFactorType"]] == "BF01") {
priorOdds <- pH0 / (1 - pH0)
logBF <- -logBF
} else {
priorOdds <- (1 - pH0) / pH0
}
postOdds <- log(priorOdds) + logBF
postOdds <- exp(postOdds)
if (options[["bayesFactorType"]] != "LogBF10")
logBF <- exp(logBF)
error <- ttest@bayesFactor$error * 100
}
}
row <- list(
"(I)" = pairs[1L, i],
"(J)" = pairs[2L, i],
"Prior Odds" = priorOdds,
"Posterior Odds" = postOdds,
"BF" = logBF,
"error %" = error
)
postHocTable$addRows(row, rowNames = paste0("row", i))
}
if (!is.null(errMessages)) {
for (i in seq_along(errMessages))
postHocTable$addFootnote(message = errMessages[[i]][["message"]],
rowNames = paste0("row", errMessages[[i]][["row_names"]]))
}
postHocCollection[[paste0("postHoc_", posthoc.var)]] <- postHocTable
}
return()
}
# Data reading ----
.BANOVAreadData <- function(dataset, options, analysisType) {
if (is.null(dataset)) {
if (analysisType == "RM-ANOVA")
return(.BANOVAreadRManovaData(dataset, options))
numeric.vars <- c(unlist(options$covariates), unlist(options$dependent))
numeric.vars <- numeric.vars[numeric.vars != ""]
factor.vars <- c(unlist(options$fixedFactors), unlist(options$randomFactors))
factor.vars <- factor.vars[factor.vars != ""]
dataset <- .readDataSetToEnd(
columns.as.numeric = numeric.vars,
columns.as.factor = factor.vars,
exclude.na.listwise = c(numeric.vars, factor.vars)
)
}
return(dataset)
}
.BANOVAdependentName <- "JaspColumn_.dependent._Encoded"
.BANOVAsubjectName <- "JaspColumn_.subject._Encoded"
.BANOVAdecodeSubject <- function(string) {
return(gsub(.BANOVAsubjectName, "subject", string))
}
.BANOVAdecodeNuisance <- function(nuisance) {
# .BANOVAsubjectName needs to be handled separately
idx <- nuisance == .BANOVAsubjectName
nuisance[idx] <- "subject"
nuisance[!idx] <- jaspBase::gsubInteractionSymbol(nuisance[!idx])
return(nuisance)
}
.BANOVAreadRManovaData <- function(dataset, options) {
if (!("" %in% options$repeatedMeasuresCells)) {
rm.vars <- options$repeatedMeasuresCells
bs.factors <- options$betweenSubjectFactors
bs.covariates <- options$covariates
rm.factors <- options$repeatedMeasuresFactors
all.variables <- c (bs.factors, bs.covariates, rm.vars)
dataset <- .readDataSetToEnd(
columns.as.numeric = c(rm.vars, bs.covariates),
columns.as.factor = bs.factors,
exclude.na.listwise = all.variables
)
dataset <- try(
.shortToLong(dataset, rm.factors, rm.vars, c(bs.factors, bs.covariates), dependentName = .BANOVAdependentName, subjectName = .BANOVAsubjectName),
silent = TRUE
)
}
return(dataset)
}
# Descriptives ----
.BANOVAdescriptives <- function(jaspResults, dataset, options, errors, analysisType, ready = TRUE, position = 9001) {
if (!ready)
return()
# the main use of this function is that descriptives can now be reused for the frequentist ANOVAs
# without the container, the position could mess things up
descriptivesContainer <- jaspResults[["descriptivesContainer"]]
if (is.null(descriptivesContainer)) {
descriptivesContainer <- createJaspContainer(title = gettext("Descriptives"))
descriptivesContainer$dependOn(c("dependent", "repeatedMeasuresCells"))
descriptivesContainer$position <- position
jaspResults[["descriptivesContainer"]] <- descriptivesContainer
}
.BANOVAdescriptivesTable (descriptivesContainer, dataset, options, errors, analysisType)
.BANOVAdescriptivesPlots (descriptivesContainer, dataset, options, errors, analysisType)
.BANOVAbarPlots (descriptivesContainer, dataset, options, errors, analysisType)
.BANOVArainCloudPlots (descriptivesContainer, dataset, options, errors, analysisType)
return()
}
.BANOVAdescriptivesTable <- function(jaspContainer, dataset, options, errors, analysisType) {
if (!options[["descriptives"]] || !is.null(jaspContainer[["tableDescriptives"]]))
return()
if (analysisType == "RM-ANOVA") {
dependent <- .BANOVAdependentName
fixed <- unlist(c(lapply(options[["repeatedMeasuresFactors"]], `[[`, "name"), options[["betweenSubjectFactors"]]))
title <- gettext("Descriptives")
} else {
dependent <- options[["dependent"]]
fixed <- options[["fixedFactors"]]
if (!is.null(dependent) && dependent != "") {
title <- gettextf("Descriptives - %s", options$dependent)
} else {
title <- gettext("Descriptives")
}
}
descriptivesTable <- createJaspTable(title = title)
descriptivesTable$position <- 1
descriptivesTable$dependOn(c("dependent", "fixedFactors", "betweenSubjectFactors", "descriptives",
"credibleInterval"))
# internal names use base64 so they don't have " " which R changes into "." because R does that to dataframe names.
# Also adds a . in case the base64 is magically "Mean", "SD" or "N"
fixedDot <- paste0(fixed, ".")
for (i in seq_along(fixed))
descriptivesTable$addColumnInfo(name = fixedDot[i], type = "string", title = fixed[i], combine = TRUE)
descriptivesTable$addColumnInfo(name = "N", title = gettext("N"), type = "integer")
descriptivesTable$addColumnInfo(name = "Mean", title = gettext("Mean"), type = "number")
descriptivesTable$addColumnInfo(name = "SD", title = gettext("SD"), type = "number")
descriptivesTable$addColumnInfo(name = "SE", title = gettext("SE"), type = "number")
descriptivesTable$addColumnInfo(name = "coefOfVariation", title = gettext("Coefficient of variation"), type = "number")
if (!is.null(options[["credibleInterval"]])) {
overTitle <- gettextf("%s%% Credible Interval", format(100 * options[["credibleInterval"]], digits = 3))
descriptivesTable$addColumnInfo(name = gettext("Lower"), type = "number", overtitle = overTitle)
descriptivesTable$addColumnInfo(name = gettext("Upper"), type = "number", overtitle = overTitle)
}
descriptivesTable$showSpecifiedColumnsOnly <- TRUE
jaspContainer[["tableDescriptives"]] <- descriptivesTable
if (errors$noVariables)
return()
# order the data to show
dataset2 <- dataset[do.call(order, dataset[, fixed, drop = FALSE]), c(dependent, fixed)]
# by pasting the fixedFactors together we obtain the unique indices to group on. This excludes
# non-existent combinations. A "." is added to deal with the level "".
ind <- apply(dataset2[, fixed, drop = FALSE], 1L, paste0, ".", collapse = "")
# temporary function to calculate all descriptives
tmpFun <- function(data, fixed, dependent, cri) {
row <- list()
for (j in fixed)
row[[paste0(j, ".")]] <- as.character(data[1L, j])
N <- nrow(data)
row[["N"]] <- N
if (N == 0L) {
row[["Mean"]] <- row[["SD"]] <- row[["SE"]] <- row[["coefOfVariation"]] <- NA
} else if (N == 1L) {
row[["Mean"]] <- data[[dependent]]
row[["SD"]] <- row[["SE"]] <- row[["coefOfVariation"]] <- row[["Lower"]] <- row[["Upper"]] <- NA
} else {
row[["Mean"]] <- mean(data[[dependent]])
row[["SD"]] <- stats::sd(data[[dependent]])
row[["SE"]] <- stats::sd(data[[dependent]]) / sqrt(N)
row[["coefOfVariation"]] <- stats::sd(data[[dependent]]) / mean(data[[dependent]])
tmp <- jaspTTests::.posteriorSummaryGroupMean(data[[dependent]], cri)
row[["Lower"]] <- tmp[["ciLower"]]
row[["Upper"]] <- tmp[["ciUpper"]]
}
return(row)
}
# apply tempFun on each subset defined by ind
rows <- by(dataset2, ind, tmpFun, fixed = fixed,
dependent = dependent, cri = options[["credibleInterval"]])
# do.call(rbind, rows) turns rows into a data.frame (from a list) for jaspResults
data <- do.call(rbind.data.frame, rows)
# the results of `by` are in alphabetical order, so we reorder the data.frame to match the factor level order. Fixes https://github.com/jasp-stats/jasp-issues/issues/1741
newOrder <- match(unique(ind), rownames(data))
data <- data[newOrder, ]
# add footnote if there are unobserved combinations
nObserved <- nrow(data)
nPossible <- prod(sapply(dataset2[, fixed, drop = FALSE], nlevels))
if (nObserved != nPossible) {
descriptivesTable$addFootnote(
message = gettextf(
"Some combinations of factors are not observed and hence omitted (%1$g out of %2$g combinations are unobserved).",
nPossible - nObserved, nPossible
)
)
}
descriptivesTable$setData(data)
return()
}
.BANOVAdescriptivesPlots <- function(jaspContainer, dataset, options, errors, analysisType) {
if (length(options[["descriptivePlotHorizontalAxis"]]) == 0L
|| options[["descriptivePlotHorizontalAxis"]] == ""
|| !is.null(jaspContainer[["containerDescriptivesPlots"]]))
return()
descriptivesPlotContainer <- createJaspContainer(title = gettext("Descriptives plots"))
descriptivesPlotContainer$position <- 2
jaspContainer[["containerDescriptivesPlots"]] <- descriptivesPlotContainer
# either Bayesian or Frequentist anova
if (is.null(options$descriptivePlotErrorBarType)) { # TRUE implies Bayesian
plotErrorBars <- options$descriptivePlotCi
errorBarType <- "ci"
conf.interval <- options$descriptivePlotCiLevel
normalizeErrorBars <- TRUE
descriptivesPlotContainer$dependOn(c("dependent", "descriptivePlotCi", "descriptivePlotCiLevel"))
} else {
plotErrorBars <- options$descriptivePlotErrorBar
errorBarType <- options$descriptivePlotErrorBarType
conf.interval <- options$descriptivePlotCiLevel
normalizeErrorBars <- if (is.null(options[["normalizeErrorBarsDescriptives"]])) FALSE else options[["normalizeErrorBarsDescriptives"]]
descriptivesPlotContainer$dependOn(c("dependent", "descriptivePlotErrorBar", "descriptivePlotErrorBarType", "descriptivePlotCiLevel",
"descriptivePlotErrorBarPooled", "normalizeErrorBarsDescriptives"))
}
usePooledSE <- if (is.null(options[["descriptivePlotErrorBarPooled"]])) FALSE else options[["descriptivePlotErrorBarPooled"]]
descriptivesPlotContainer$dependOn(c("descriptivePlotHorizontalAxis", "descriptivePlotSeparateLines", "descriptivePlotSeparatePlot", "descriptivePlotYAxisLabel"))
if (errors$noVariables) {
descriptivesPlotContainer[["dummyplot"]] <- createJaspPlot(title = gettext("Descriptives Plot"))
return()
}
groupVars <- c(options$descriptivePlotHorizontalAxis, options$descriptivePlotSeparateLines, options$descriptivePlotSeparatePlot)
groupVars <- groupVars[groupVars != ""]
if (analysisType == "RM-ANOVA") {
dependent <- .BANOVAdependentName
yLabel <- options[["descriptivePlotYAxisLabel"]]
} else {
dependent<- options$dependent
yLabel <- options[["dependent"]]
}
betweenSubjectFactors <- groupVars[groupVars %in% options$betweenSubjectFactors]
repeatedMeasuresFactors <- groupVars[groupVars %in% sapply(options$repeatedMeasuresFactors, `[[`, "name")]
if (length(repeatedMeasuresFactors) == 0 || isFALSE(normalizeErrorBars)) {
summaryStat <- jaspTTests::.summarySE(as.data.frame(dataset), measurevar = dependent, groupvars = groupVars,
conf.interval = conf.interval, na.rm = TRUE, .drop = FALSE,
errorBarType = errorBarType, dependentName = .BANOVAdependentName,
subjectName = .BANOVAsubjectName)
} else {
summaryStat <- jaspTTests::.summarySEwithin(as.data.frame(dataset), measurevar= .BANOVAdependentName,
betweenvars = betweenSubjectFactors,
withinvars = repeatedMeasuresFactors,
idvar = .BANOVAsubjectName,
conf.interval = conf.interval,
na.rm=TRUE, .drop = FALSE, errorBarType = errorBarType,
usePooledSE = usePooledSE,
dependentName = .BANOVAdependentName,
subjectName = .BANOVAsubjectName)
}
if (options[["descriptivePlotHorizontalAxis"]] %in% options[["covariates"]]) {
splitScatterOptions <- options
splitScatterOptions[["colorPalette"]] <- "colorblind3"
splitScatterOptions[["scatterPlotLegend"]] <- TRUE
splitScatterOptions[["scatterPlotRegressionLine"]] <- TRUE
splitScatterOptions[["scatterPlotRegressionLineCi"]] <- plotErrorBars
splitScatterOptions[["scatterPlotRegressionLineType"]] <- "linear"
splitScatterOptions[["scatterPlotGraphTypeAbove"]] <- "none"
splitScatterOptions[["scatterPlotGraphTypeRight"]] <- "none"
splitScatterOptions[["scatterPlotRegressionLineCiLevel"]] <- options[["descriptivePlotCiLevel"]]
if (options$descriptivePlotSeparatePlot != "") {
for (thisLevel in levels(dataset[[options[["descriptivePlotSeparatePlot"]]]])) {
subData <- dataset[dataset[[options[["descriptivePlotSeparatePlot"]]]] == thisLevel, ]
thisPlotName <- paste0(options[["descriptivePlotHorizontalAxis"]], " - ", options[["dependent"]], ": ",
options[["descriptivePlotSeparatePlot"]], " = ", thisLevel)
jaspDescriptives::.descriptivesScatterPlots(descriptivesPlotContainer, subData, c(options[["descriptivePlotHorizontalAxis"]], options[["dependent"]]),
split = options[["descriptivePlotSeparateLines"]], options = splitScatterOptions, name = thisPlotName,
dependOnVariables = FALSE)
}
} else {
jaspDescriptives::.descriptivesScatterPlots(descriptivesPlotContainer, dataset, c(options[["descriptivePlotHorizontalAxis"]], options[["dependent"]]),
split = options[["descriptivePlotSeparateLines"]], options = splitScatterOptions, dependOnVariables = FALSE)
}
return()
}
colnames(summaryStat)[colnames(summaryStat) == dependent] <- "dependent"
if (options$descriptivePlotHorizontalAxis != "") {
colnames(summaryStat)[colnames(summaryStat) == options$descriptivePlotHorizontalAxis] <- "descriptivePlotHorizontalAxis"
}
if (options$descriptivePlotSeparateLines != "") {
colnames(summaryStat)[colnames(summaryStat) == options$descriptivePlotSeparateLines] <- "descriptivePlotSeparateLines"
}
if (options$descriptivePlotSeparatePlot != "") {
colnames(summaryStat)[colnames(summaryStat) == options$descriptivePlotSeparatePlot] <- "descriptivePlotSeparatePlot"
}
base_breaks_x <- function(x){
b <- unique(as.numeric(x))
d <- data.frame(y=-Inf, yend=-Inf, x=min(b), xend=max(b))
list(ggplot2::geom_segment(data=d, ggplot2::aes(x=x, y=y, xend=xend, yend=yend), inherit.aes=FALSE, size = 1))
}
base_breaks_y <- function(x, plotErrorBars){
if (plotErrorBars) {
ci.pos <- c(x[,"dependent"], x[,"dependent"]-x[,"ci"],x[,"dependent"]+x[,"ci"])
b <- pretty(ci.pos)
d <- data.frame(x=-Inf, xend=-Inf, y=min(b), yend=max(b))
list(ggplot2::geom_segment(data=d, ggplot2::aes(x=x, y=y, xend=xend, yend=yend), inherit.aes=FALSE, size = 1),
ggplot2::scale_y_continuous(breaks=c(min(b),max(b))))
} else {
b <- pretty(x[,"dependent"])
d <- data.frame(x=-Inf, xend=-Inf, y=min(b), yend=max(b))
list(ggplot2::geom_segment(data=d, ggplot2::aes(x=x, y=y, xend=xend, yend=yend), inherit.aes=FALSE, size = 1),
ggplot2::scale_y_continuous(breaks=c(min(b),max(b))))
}
}
if (options$descriptivePlotSeparatePlot != "") {
subsetPlots <- levels(summaryStat[,"descriptivePlotSeparatePlot"])
nPlots <- length(subsetPlots)
} else {
nPlots <- 1
}
for (i in seq_len(nPlots)) {
if (nPlots > 1L) {
title <- paste(options$descriptivePlotSeparatePlot,": ",subsetPlots[i], sep = "")
} else {
title <- ""
}
descriptivesPlot <- createJaspPlot(title = title)
descriptivesPlotContainer[[title]] <- descriptivesPlot
descriptivesPlot$height <- 300
if (options$descriptivePlotSeparateLines != "") {
descriptivesPlot$width <- 430
} else {
descriptivesPlot$width <- 300
}
if (options$descriptivePlotSeparatePlot != "") {
summaryStatSubset <- subset(summaryStat,summaryStat[,"descriptivePlotSeparatePlot"] == subsetPlots[i])
} else {
summaryStatSubset <- summaryStat
}
if (options$descriptivePlotSeparateLines == "") {
p <- ggplot2::ggplot(summaryStatSubset, ggplot2::aes(x=descriptivePlotHorizontalAxis,
y=dependent,
group=1))
} else {
p <- ggplot2::ggplot(summaryStatSubset, ggplot2::aes(x=descriptivePlotHorizontalAxis,
y=dependent,
group=descriptivePlotSeparateLines,
shape=descriptivePlotSeparateLines,
fill=descriptivePlotSeparateLines))
}
if (plotErrorBars && !(options[["descriptivePlotHorizontalAxis"]] %in% options[["covariates"]])) {
pd <- ggplot2::position_dodge(.2)
p = p + ggplot2::geom_errorbar(ggplot2::aes(ymin=ciLower,
ymax=ciUpper),
colour="black", width=.2, position=pd)
} else {
pd <- ggplot2::position_dodge(0)
}
guideLegend <- ggplot2::guide_legend(nrow = min(10, nlevels(summaryStatSubset$descriptivePlotSeparateLines)),
title = options$descriptivePlotSeparateLines, keywidth = 0.1, keyheight = 0.3,
default.unit = "inch")
if (options[["descriptivePlotHorizontalAxis"]] %in% options[["covariates"]]) {
line <- ggplot2::geom_smooth(method = "lm", size = .7, color = "black", se = FALSE)
addHorizontalVar <- summaryStatSubset[,"descriptivePlotHorizontalAxis"]
} else {
line <- ggplot2::geom_line(position=pd, size = .7)
}
if (plotErrorBars) {
ci.pos <- c(summaryStatSubset[,"dependent"],
summaryStatSubset[,"dependent"]-summaryStatSubset[,"ci"],
summaryStatSubset[,"dependent"]+summaryStatSubset[,"ci"],
min(summaryStatSubset[,"dependent"])*1.1,
max(summaryStatSubset[,"dependent"])*1.1)
yBreaks <- jaspGraphs::getPrettyAxisBreaks(ci.pos)
} else {
yBreaks <- jaspGraphs::getPrettyAxisBreaks(c(summaryStatSubset[,"dependent"],
min(summaryStatSubset[,"dependent"])*1.1,
max(summaryStatSubset[,"dependent"])*1.1))
}
if (options[["descriptivePlotHorizontalAxis"]] %in% options[["covariates"]]) {
ggXaxis <- ggplot2::scale_x_continuous(breaks = jaspGraphs::getPrettyAxisBreaks(summaryStatSubset[,"descriptivePlotHorizontalAxis"]))
} else {
ggXaxis <- ggplot2::scale_x_discrete(breaks = jaspGraphs::getPrettyAxisBreaks(summaryStatSubset[,"descriptivePlotHorizontalAxis"]))
}
p <- p + line +
ggplot2::geom_point(position=pd, size=4) +
ggplot2::scale_fill_manual(values = c(rep(c("white","black"),5),rep("grey",100)), guide=guideLegend) +
ggplot2::scale_shape_manual(values = c(rep(c(21:25),each=2),21:25,7:14,33:112), guide=guideLegend) +
ggplot2::scale_color_manual(values = rep("black",200),guide=guideLegend) +
ggplot2::labs(y = yLabel, x = options[["descriptivePlotHorizontalAxis"]]) +
ggplot2::scale_y_continuous(breaks = yBreaks, limits = range(yBreaks)) +
ggXaxis +
jaspGraphs::geom_rangeframe() +
jaspGraphs::themeJaspRaw(legend.position = "right")
descriptivesPlot$plotObject <- p
}
return()
}
.BANOVAbarPlots <- function(jaspContainer, dataset, options, errors, analysisType) {
if (length(options[["barPlotHorizontalAxis"]]) == 0L
|| options[["barPlotHorizontalAxis"]] == ""
|| !is.null(jaspContainer[["containerBarPlots"]]))
return()
barPlotContainer <- createJaspContainer(title = gettext("Bar plots"))
barPlotContainer$position <- 3
jaspContainer[["containerBarPlots"]] <- barPlotContainer
plotErrorBars <- options[["barPlotErrorBars"]]
errorBarType <- options[["barPlotErrorBarType"]]
confInterval <- options[["barPlotCiInterval"]]
barPlotContainer$dependOn(c("dependent", "barPlotErrorBars", "barPlotErrorBarType", "applyMoreyCorrectionErrorBarsBarplot",
"barPlotHorizontalZeroFix", "barPlotCiInterval", "usePooledStandErrorCITwo"))
usePooledSE <- if (is.null(options[["usePooledStandErrorCITwo"]])) FALSE else options[["usePooledStandErrorCITwo"]]
normalizeErrorBars <- if (is.null(options[["normalizeErrorBarsBarplot"]])) TRUE else options[["normalizeErrorBarsBarplot"]]
barPlotHorizontalZeroFix <- options[["barPlotHorizontalZeroFix"]]
barPlotContainer$dependOn(c("barPlotHorizontalAxis", "barPlotSeparatePlots", "labelYAxisTwo"))
if (errors$noVariables) {
barPlotContainer[["dummyplot"]] <- createJaspPlot(title = gettext("Bar Plot"))
return()
}
groupVars <- c(options[["barPlotHorizontalAxis"]], options[["barPlotSeparatePlots"]])
groupVars <- groupVars[groupVars != ""]
if (analysisType == "RM-ANOVA") {
dependent <- .BANOVAdependentName
yLabel <- options[["labelYAxisTwo"]]
} else {
dependent<- options[["dependent"]]
yLabel <- options[["dependent"]]
}
betweenSubjectFactors <- groupVars[groupVars %in% options[["betweenSubjectFactors"]]]
repeatedMeasuresFactors <- groupVars[groupVars %in% sapply(options[["repeatedMeasuresFactors"]], `[[`, "name")]
if (length(repeatedMeasuresFactors) == 0 || isFALSE(normalizeErrorBars)) {
summaryStat <- jaspTTests::.summarySE(as.data.frame(dataset), measurevar = dependent, groupvars = groupVars,
conf.interval = confInterval, na.rm = TRUE, .drop = FALSE,
errorBarType = errorBarType, dependentName = .BANOVAdependentName,
subjectName = .BANOVAsubjectName)
} else {
summaryStat <- jaspTTests::.summarySEwithin(as.data.frame(dataset), measurevar = .BANOVAdependentName,
betweenvars = betweenSubjectFactors,
withinvars = repeatedMeasuresFactors,
idvar = .BANOVAsubjectName,
conf.interval = confInterval,
na.rm=TRUE, .drop = FALSE, errorBarType = errorBarType,
usePooledSE = usePooledSE,
dependentName = .BANOVAdependentName,
subjectName = .BANOVAsubjectName)
}
colnames(summaryStat)[colnames(summaryStat) == dependent] <- "dependent"
if (options[["barPlotHorizontalAxis"]] != "") {
colnames(summaryStat)[colnames(summaryStat) == options[["barPlotHorizontalAxis"]]] <- "barPlotHorizontalAxis"
}
if (options[["barPlotSeparatePlots"]] != "") {
colnames(summaryStat)[colnames(summaryStat) == options[["barPlotSeparatePlots"]]] <- "barPlotSeparatePlots"
subsetPlots <- levels(summaryStat[, "barPlotSeparatePlots"])
nPlots <- length(subsetPlots)
} else {
nPlots <- 1L
}
for (i in seq_len(nPlots)) {
title <- if (nPlots > 1L) {
paste0(options[["barPlotSeparatePlots"]], ": ", subsetPlots[i])
} else {
""
}
barPlot <- createJaspPlot(title = title)
barPlotContainer[[title]] <- barPlot
barPlot$height <- 500
barPlot$width <- 500
if (options[["barPlotSeparatePlots"]] != "") {
summaryStatSubset <- subset(summaryStat, summaryStat[, "barPlotSeparatePlots"] == subsetPlots[i])
} else {
summaryStatSubset <- summaryStat
}
error <- NULL
if (plotErrorBars) {
pd <- ggplot2::position_dodge(.2)
error <- ggplot2::geom_errorbar(ggplot2::aes(ymin = ciLower, ymax = ciUpper),
colour = "black", width = .2, position = pd)
}
values <- 1.1 * range(summaryStat[["dependent"]])
if (barPlotHorizontalZeroFix)
values <- c(0, values)
if (plotErrorBars) {
ci.pos <- c(summaryStat[,"dependent"],
summaryStat[,"dependent"]-summaryStat[,"ci"],
summaryStat[,"dependent"]+summaryStat[,"ci"],
min(summaryStat[,"dependent"])*1.1,
max(summaryStat[,"dependent"])*1.1)
yBreaks <- jaspGraphs::getPrettyAxisBreaks(c(values, ci.pos))
} else {
yBreaks <- jaspGraphs::getPrettyAxisBreaks(values)
}
pd2 <- ggplot2::position_dodge2(preserve = "single")
p <- ggplot2::ggplot(summaryStatSubset, ggplot2::aes(x = barPlotHorizontalAxis, y = dependent, group = 1)) +
ggplot2::geom_hline(yintercept = 0, color = "#858585", size = 0.3) +
ggplot2::geom_bar(stat = "identity", fill = "grey", col = "black", width = .6, position = pd2) +
error +
ggplot2::labs(y = yLabel, x = options[["barPlotHorizontalAxis"]]) +
ggplot2::scale_y_continuous(breaks = yBreaks, limits = range(yBreaks), oob = scales::rescale_none) +
ggplot2::scale_x_discrete(breaks = jaspGraphs::getPrettyAxisBreaks(summaryStatSubset[, "barPlotHorizontalAxis"])) +
jaspGraphs::geom_rangeframe(sides = "l") +
jaspGraphs::themeJaspRaw()
barPlot$plotObject <- p
}
return()
}
.BANOVArainCloudPlots <- function(jaspContainer, dataset, options, errors, analysisType) {
if (length(options[["rainCloudHorizontalAxis"]]) == 0L
|| options[["rainCloudHorizontalAxis"]] == ""
|| !is.null(jaspContainer[["containerRainCloudPlots"]]))
return()
rainCloudPlotsContainer <- createJaspContainer(title = gettext("Raincloud plots"))
rainCloudPlotsContainer$position <- 4
jaspContainer[["containerRainCloudPlots"]] <- rainCloudPlotsContainer
rainCloudPlotsContainer$dependOn(c("dependent", "rainCloudHorizontalAxis", "rainCloudSeparatePlots",
"rainCloudYAxisLabel", "rainCloudHorizontalDisplay"))
if (errors$noVariables) {
rainCloudPlotsContainer[["dummyplot"]] <- createJaspPlot(title = "")
return()
}
groupVar <- options[["rainCloudHorizontalAxis"]]
if (analysisType == "RM-ANOVA") {
addLines <- !(groupVar %in% unlist(options[["betweenSubjectFactors"]]))
dependentV <- .BANOVAdependentName
yLabel <- options[["rainCloudYAxisLabel"]]
if (trimws(yLabel) == "") {
title <- gettext("Dependent")
yLabel <- NULL
} else {
title <- yLabel
}
} else {
addLines <- FALSE
dependentV <- options[["dependent"]]
yLabel <- options[["dependent"]]
title <- options[["dependent"]]
}
if (!is.null(options$rainCloudHorizontalDisplay) && options$rainCloudHorizontalDisplay)
horiz <- TRUE
else
horiz <- FALSE
if (options$rainCloudSeparatePlots != "") {
for (thisLevel in levels(dataset[[options[["rainCloudSeparatePlots"]]]])) {
subData <- dataset[dataset[[options[["rainCloudSeparatePlots"]]]] == thisLevel, ]
thisPlotName <- paste0(title, ": ", options[["rainCloudSeparatePlots"]], ": ", thisLevel)
subPlot <- createJaspPlot(title = thisPlotName, width = 480, height = 320)
rainCloudPlotsContainer[[thisLevel]] <- subPlot
p <- try(jaspTTests::.descriptivesPlotsRainCloudFill(subData, dependentV, groupVar, yLabel, groupVar, addLines, horiz, NULL))
if(isTryError(p))
subPlot$setError(.extractErrorMessage(p))
else
subPlot$plotObject <- p
}
} else {
singlePlot <- createJaspPlot(title = title, width = 480, height = 320)
rainCloudPlotsContainer[["rainCloudPlotSingle"]] <- singlePlot
p <- try(jaspTTests::.descriptivesPlotsRainCloudFill(dataset, dependentV, groupVar, yLabel, groupVar, addLines, horiz, NULL))
if(isTryError(p))
singlePlot$setError(.extractErrorMessage(p))
else
singlePlot$plotObject <- p
}
return()
}
# Sample posteriors ----
.BANOVAsamplePosteriors <- function(jaspResults, dataset, options, model, state) {
# TODO: the density approximation can become more efficient with a fast parametric density approximation
# TODO: Bayes factor samples unobserved interaction levels, what to do?
# if the most complex model is retrieved from the state?
nIter <- if (options[["samplingMethodMCMC"]] == "auto") 1e4L else options[["samplesMCMC"]]
nmodels <- length(model[["models"]])
postProbs <- model[["postProbs"]]
statistics <- vector("list", nmodels)
randomFactors <- .BANOVAgetRandomFactors(options, model[["analysisType"]])
dataTypes <- .BANOVAgetDataTypes(dataset, model[["model.list"]][[nmodels]], randomFactors)
levelInfo <- .BANOVAgetLevelInfo(dataset, model[["model.list"]][[nmodels]], dataTypes)
renameFrom <- renameTo <- NULL
if (!is.null(state)) { # can we reuse some posteriors?
reuseable <- model[["reuseable"]]
# it's possible that a user just renames a level of a repeated measures factor
# in that case everything can be reused, but we have to rename some parameters.
# the same may happen when a variable becomes a nuisance parameter, because then
# BayesFactor changes its order in the names.
oldLevelInfo <- state[["levelInfo"]]$levelNames
newLevelInfo <- levelInfo$levelNames
sortTerms <- function(x) {
# split b:a to c("b", "a"), sort it, and then paste it back
# otherwise posteriors samples between different runs are not correctly retrieved from the state.
sapply(strsplit(x, ":", fixed = TRUE), function(x) paste(sort(x), collapse = ":"))
}
tmp_old <- sortTerms(names(oldLevelInfo))
tmp_new <- sortTerms(names(newLevelInfo))
matches <- match(tmp_new, tmp_old)
for (i in seq_along(matches)) {
if (is.na(matches[i]))
next
j <- matches[i]
if (!identical(oldLevelInfo[[j]], newLevelInfo[[i]])) {
if (length(oldLevelInfo[[j]]) == length(newLevelInfo[[i]])) {
renameFrom <- c(renameFrom, oldLevelInfo[[j]])
renameTo <- c(renameTo, newLevelInfo[[i]])
}
}
}
} else {
reuseable <- rep(NA, nmodels)
}
# NOTE: some code checks require saving all samples. To do so, change all samples to samples[[i]] and uncomment:
# samples <- vector("list", nmodels)
allParamNames <- c("mu", unlist(levelInfo$levelNames))
nparam <- length(allParamNames)
weightedMeans <- weights <- numeric(nparam)
names(weights) <- names(weightedMeans) <- allParamNames
allContinuous <- TRUE
h <- (1 - options[["credibleInterval"]]) / 2
probs <- c(h, 1-h)
originalFun <- BayesFactor:::makeChainNeater
on.exit(jaspBase:::assignFunctionInPackage(originalFun, "makeChainNeater", "BayesFactor"))
jaspBase:::assignFunctionInPackage(.BANOVAmakeChainNeater, "makeChainNeater", "BayesFactor")
samplingIssues <- list()
startProgressbar(nmodels, gettext("Sampling posteriors"))
for (i in seq_len(nmodels)) {
# check if the state is reuseable and if it's not NULL, which means it didn't crash
if (is.na(reuseable[i]) || is.null(state$statistics[[reuseable[i]]])) {
if (i == 1L && is.null(model[["models"]][[i]][["bf"]])) {
# NULL model only contains an intercept, use custom sampler
# NOTE: RM-ANOVA never enters here (and would crash if it did)
.setSeedJASP(options)
samples <- .BANOVAsampleNullModel(dataset[[options[["dependent"]]]], nsamples = nIter)
types <- NULL
} else {
# NOTE: we have to sample the random effects, otherwise we cant make predictions (needed for residuals and R^2)
# put the dataset back in
.setSeedJASP(options)
bfObj <- model[["models"]][[i]][["bf"]]
bfObj@data <- dataset
samples <- try(BayesFactor::posterior(bfObj, iterations = nIter))
if (isTryError(samples)) {
# this may error whenever the Bayes factor couldn't be calculated
next
}
types <- samples@model@dataTypes
# BayesFactor stores an internal copy of the dataset, so it can still have old names
if (length(renameFrom) > 0) {
cnms <- colnames(samples)
cnms <- plyr::mapvalues(cnms, renameFrom, renameTo, warn_missing = FALSE) # NOTE: dependency could be removed
colnames(samples) <- cnms
}
# keep only relevant columns, drop sig2, g_xxx, ...
idx <- match(allParamNames, colnames(samples), nomatch = 0L)
samples <- samples[, idx, drop = FALSE]
# for some odd reason, Bayesfactor uses as column name contcor1-contcor1
# if there is a covariate AND fixed factors, but only contcor1 if all variables are continuous...
if (all(types == "continuous")) {
cnms <- colnames(samples)[-1L] # omit the intercept (mu) which is not changed by Bayesfactor
colnames(samples)[-1L] <- paste0(cnms, "-", cnms)
} else {
allContinuous <- FALSE
}
# so some models may yield a bunch of NAs, for example,
# BayesFactor::lmBF(contNormal ~ facGender + contGamma + facGender * contGamma, data = "debug.csv", whichRandom = "facGender")
# we add a footnote and try to not to crash.
if (anyNA(samples)) {
originalRows <- nrow(samples)
samples <- samples[complete.cases(samples), , drop = FALSE]
remainingRows <- nrow(samples)
samplingIssues[[length(samplingIssues) + 1L]] <- list(
model = model[["models"]][[i]][["title"]],
originalRows = originalRows,
remainingRows = remainingRows
)
next
}
}
# although matrixStats::colMeans2 is faster than .colMeans the cost of matrixStats:: is not worth it.
nms <- colnames(samples)
statistics[[i]]$names <- nms # <- these are the names for all objects within one sublist
statistics[[i]]$mean <- .colMeans (samples, m = nIter, n = NCOL(samples))
statistics[[i]]$var <- matrixStats::colVars (samples)
statistics[[i]]$cri <- matrixStats::colQuantiles(samples, probs = probs)
statistics[[i]]$approx <- .BANOVAfitDensity(samples = samples)
statistics[[i]]$types <- types
} else { # reuse state
statistics[[i]] <- state$statistics[[reuseable[i]]]
nms <- statistics[[i]]$names
types <- statistics[[i]]$types
if (length(renameFrom) > 0) {
nms <- plyr::mapvalues(nms, renameFrom, renameTo, warn_missing = FALSE) # NOTE: dependency could be removed
statistics[[i]]$names <- nms
}
if (allContinuous)
allContinuous <- all(types == "continuous")
}
# compute model averaged posterior means
weightedMeans[nms] <- weightedMeans[nms] + postProbs[i] * statistics[[i]]$mean
weights [nms] <- weights [nms] + postProbs[i]
progressbarTick()
}
# find out which parameters are random -- this uses types from the last iteration above
isRandom <- logical(nparam)
idx <- which(types == "random")
for (i in idx)
isRandom <- isRandom | startsWith(names(weights), names(types)[i])
# the weights used above don't sum to 1 because we consider a subset of the models.
# Now we renormalize to ensure the weights used in each weighted mean do sum to 1.
weightedMeans <- weightedMeans / weights
# given the model averaged posterior means calculate the weighted posterior standard deviations
weightedSds <- 0 * weightedMeans # keeps the names
r <- nIter / (nIter - 1)
# statistics has length 0 iff the posterior sampling errored
widxGoodStatistics <- which(lengths(statistics) > 0L)
for (i in widxGoodStatistics) {
nms <- statistics[[i]]$names
var <- statistics[[i]]$var # ~ sum((x - mean(x))^2
mu <- statistics[[i]]$mean
# cc = sum((x - y)^2) - sum((x - mean(x))^2), where y is the weighted mean
# hence, we can get the weighted sum of squares from the individual posterior variances
cc <- weightedMeans[nms]^2 + mu^2 - 2 * mu * weightedMeans[nms]
weightedSds[nms] <- weightedSds[nms] + postProbs[i] * (var + cc * r)
}
weightedSds <- sqrt(weightedSds / weights)
# the loops above are optimized versions of the code below that also calculates the
# weighted mean and weighted sd, but needs to store all samples at once.
# xx <- lapply(samples, function(x, y) {
# if (y %in% colnames(x)) {dd <- model$posteriors$weightedRsqDens
# x[, y]
# } else {
# NULL
# }
# }, y = "XY29udEJpbm9t-1") # XY29udEJpbm9t = contBinom
# idx <- lengths(xx) > 0
# xx <- unlist(xx[idx])
# w <- postProbs[idx]
# sum(w) # equals weights[2]
# weighted.mean(xx, rep(w / sum(w), each = nIter)) # equals weightedMeans[2]
# sqrt(Hmisc::wtd.var(xx, rep(w / sum(w), each = nIter), method = "unbiased")) # equals weightedSds[2]
# compute model averaged densities
steps <- 2^9 # grid size for densities is 2^steps
weightedDensities <- array(0, dim = c(steps, nparam, 2), dimnames = list(NULL, names(weightedMeans), c("x", "y")))
ranges <- matrix(0, nparam, 2L, dimnames = list(allParamNames, NULL))
# get the outermost x-values for each densities this could be vectorized with matrixStats::rowRanges if the
# data are stored as a matrix but that is memory inefficient
for (i in widxGoodStatistics) {
nms <- statistics[[i]]$names
indices <- which(nms %in% allParamNames)
for (j in indices) {
ranges[j, ] <- range(ranges[j, ], statistics[[i]]$approx$xRanges[j, ])
}
}
# construct one common grid for each observed density
for (i in seq_len(nparam)) {
weightedDensities[, i, 1L] <- seq(ranges[i, 1], ranges[i, 2], length.out = steps)
}
for (i in widxGoodStatistics) {
nms <- statistics[[i]]$names
ind <- match(nms, allParamNames)
for (j in seq_along(ind)) {
# approximate all distributions on a common grid
ap <- approx(x = statistics[[i]]$approx$fit[, j, 1L],
y = statistics[[i]]$approx$fit[, j, 2L],
xout = weightedDensities [, ind[j], 1L],
# not observed is approximated to 0
yleft = 0, yright = 0)
weightedDensities[, ind[j], 2L] <- weightedDensities[, ind[j], 2L] + ap$y * postProbs[i]
}
}
# postProbs don't sum to one so we renormalize the densities
weightedDensities[, , 2L] <- sweep(weightedDensities[, , 2L], 2, weights, FUN = `/`)
# if all model with one parameter failed, the line above introduces NaNs since 0 / 0 is NaN
weightedDensities[is.nan(weightedDensities)] <- 0
# # compute weighted CRIs
# weightedCRIs <- matrix(NA, nparam, 2L, dimnames = list(names(weights), NULL))
# cri <- options[["credibleInterval"]]
# for (i in seq_len(nparam)) {
# weightedCRIs[i, ] <- .BANOVAapproxCRI(weightedDensities[, i, 1L], weightedDensities[, i, 2L], cri)
# }
#
# # compute residuals and r-squared
# # sample from the joint posterior over models and parameters
# tmp <- .BANOVAgetSMIResidRsq(weightedDensities, dataset, model$model.list[[nmodels]], nIter, weights)
# means <- rowMeans(tmp$resid)
# h <- (1 - cri) / 2
# quants <- matrixStats::rowQuantiles(tmp$resid, probs = c(h, 1 - h))
#
# # the code above is equivalent to the code below, but the code below needs to keep all posterior samples of
# # all models in memory.
# # weights <- rep(postProbs, each = nIter)
# # independentVariable <- all.vars(.BANOVAgetModelFormulaFromBFobj(model$models[[2L]]))[1L]
# # resids <- matrix(NA, nrow(dataset), 0)
# # rsq <- vector("list", nmodels)
# # # get residuals of all models individually
# # for (i in seq_len(nmodels)) {
# # if (is.null(model$models[[i]]$bf)) {
# # tmp2 <- .BANOVAresidualsNullModel(nIter, dataset[[independentVariable]])
# # } else {
# # tmp2 <- .BANOVAgetSMIResidRsq(
# # posterior = samples[[i]],
# # dataset = dataset,
# # formula = .BANOVAgetModelFormulaFromBFobj(model$models[[i]])
# # )
# # }
# # resids <- cbind(resids, tmp2$resid)
# # rsq[[i]] <- tmp2$rsq
# # }
# # # compute weighted mean for each row
# # means2 <- tcrossprod(weights / nIter, resids)
# # plot(means, means2); abline(0, 1)
# # quants2 <- apply(resids, 1L, Hmisc::wtd.quantile, weights = weights, probs = c(0.025, 0.975))
# # plot(quants[, 1], quants2[1, ]); abline(0, 1)
# # plot(quants[, 2], quants2[2, ]); abline(0, 1)
#
# # all information for q-q plot of residuals
# weightedResidSumStats <- matrix(c(means, quants), nrow = length(means), ncol = 3L,
# dimnames = list(NULL, c("mean", "cri.2.5%", "cri.97.5%")))
#
# # all information for r-squared density plot
# weightedRsqDens <- density(tmp$rsq, n = 2^11, from = 0, to = 1)
# weightedRsqCri <- quantile(tmp$rsq, probs = c(h, 1 - h))
return(list(
statistics = statistics, weights = weights, weightedMeans = weightedMeans, weightedSds = weightedSds,
weightedDensities = weightedDensities,
#weightedCRIs = weightedCRIs, weightedResidSumStats = weightedResidSumStats,
#weightedRsqDens = weightedRsqDens, weightedRsqCri = weightedRsqCri,
allContinuous = allContinuous, isRandom = isRandom, levelInfo = levelInfo,
samplingIssues = samplingIssues
))
}
.BANOVAcomputePosteriorCRI <- function(dataset, options, model, posterior) {
nIter <- if (options[["samplingMethodMCMC"]] == "auto") 1e4L else options[["samplesMCMC"]]
weightedDensities <- posterior[["weightedDensities"]]
weights <- posterior[["weights"]]
nmodels <- length(model[["models"]])
nparam <- length(weights)
# compute weighted CRIs
weightedCRIs <- matrix(0, nparam, 2L, dimnames = list(names(weights), NULL))
cri <- options[["credibleInterval"]]
for (i in seq_len(nparam)) {
if (!all(weightedDensities[, i, 2L] == 0)) { # only TRUE if we explicitly set it to 0
weightedCRIs[i, ] <- .BANOVAapproxCRI(weightedDensities[, i, 1L], weightedDensities[, i, 2L], cri)
}
}
# compute residuals and r-squared
# sample from the joint posterior over models and parameters
tmp <- .BANOVAgetSMIResidRsq(weightedDensities, dataset, model$model.list[[nmodels]], nIter, model, options)
means <- rowMeans(tmp$resid)
h <- (1 - cri) / 2
quants <- matrixStats::rowQuantiles(tmp$resid, probs = c(h, 1 - h))
# the code above is equivalent to the code below, but the code below needs to keep all posterior samples of
# all models in memory. plug this body of this function into the call inside .BANOVAsamplePosteriors() to check
# weights <- rep(postProbs, each = nIter)
# independentVariable <- all.vars(.BANOVAgetModelFormulaFromBFobj(model$models[[2L]]))[1L]
# resids <- matrix(NA, nrow(dataset), 0)
# rsq <- vector("list", nmodels)
# # get residuals of all models individually
# for (i in seq_len(nmodels)) {
# if (is.null(model$models[[i]]$bf)) {
# tmp2 <- .BANOVAresidualsNullModel(nIter, dataset[[independentVariable]])
# } else {
# tmp2 <- .BANOVAgetSMIResidRsq(
# posterior = samples[[i]],
# dataset = dataset,
# formula = .BANOVAgetModelFormulaFromBFobj(model$models[[i]])
# )
# }
# resids <- cbind(resids, tmp2$resid)
# rsq[[i]] <- tmp2$rsq
# }
# # compute weighted mean for each row
# means2 <- tcrossprod(weights / nIter, resids)
# plot(means, means2); abline(0, 1)
# quants2 <- apply(resids, 1L, Hmisc::wtd.quantile, weights = weights, probs = c(0.025, 0.975))
# plot(quants[, 1], quants2[1, ]); abline(0, 1)
# plot(quants[, 2], quants2[2, ]); abline(0, 1)
# all information for q-q plot of residuals
weightedResidSumStats <- matrix(c(means, quants), nrow = length(means), ncol = 3L,
dimnames = list(NULL, c("mean", "cri.2.5%", "cri.97.5%")))
# all information for r-squared density plot
weightedRsqDens <- density(tmp$rsq, n = 2^11, from = 0, to = 1)
weightedRsqCri <- quantile(tmp$rsq, probs = c(h, 1 - h))
weightedRsqMean <- mean(tmp$rsq)
return(list(
weightedCRIs = weightedCRIs, weightedResidSumStats = weightedResidSumStats,
weightedRsqDens = weightedRsqDens, weightedRsqCri = weightedRsqCri, weightedRsqMean = weightedRsqMean
))
}
.BANOVAsampleNullModel <- function(nsamples, dependent) {
# sample from posterior under NULL model, needed to compute model averaged residuals.
rt.scaled <- function (n, df, mean = 0, sd = 1, ncp) {
mean + sd * stats::rt(n, df, ncp = ncp)
}
# sample from the marginal posterior distribution of the mean t-distribution, based on Murphy (2007)
n <- length(dependent)
muObs <- mean(dependent)
varObs <- var(dependent)
# uninformative priors
k0 <- 1e-6
a0 <- 1e-5
b0 <- 1e-5
an <- a0 + n / 2
bn <- b0 + 0.5 * (n - 1L) * varObs + k0 * n * muObs / (2 * (k0 + n))
kn <- k0 + n
mun <- n * muObs / kn
samples <- matrix(rt.scaled(n = nsamples, df = 2 * an, mean = mun, sd = bn / (an * kn), ncp = 0), nsamples,
1L, dimnames = list(NULL, "mu"))
return(samples)
}
.BANOVAresidualsNullModel <- function(nsamples, dependent) {
# sample from posterior under NULL model, needed to compute model averaged residuals.
rt.scaled <- function (n, df, mean = 0, sd = 1, ncp) {
mean + sd * stats::rt(n, df, ncp = ncp)
}
# sample from the marginal posterior distribution of the mean t-distribution, based on Murphy (2007)
n <- length(dependent)
muObs <- mean(dependent)
varObs <- var(dependent)
# uninformative priors
k0 <- 1e-6
a0 <- 1e-5
b0 <- 1e-5
an <- a0 + n / 2
bn <- b0 + 0.5 * (n - 1L) * varObs + k0 * n * muObs / (2 * (k0 + n))
kn <- k0 + n
mun <- n * muObs / kn
samples <- rt.scaled(n = nsamples, df = 2 * an, mean = mun, sd = bn / (an * kn), ncp = 0)
# compute all pairwise differences (residuals)
preds <- tcrossprod(rep(1, n), samples)
resids <- dependent - preds
rsq <- .BANOVAcomputeRsq(dependent, preds)
return(list(resids = resids, rsq = rsq))
}
.BANOVAfitDensity <- function(samples, gridsize = 100L, getXRange = TRUE) {
# ideally we don't do kernel density estimation but instead use some parametric approximation that is suited
# for unimodal distributions. It should also be fast, e.g., fit using method of moments (or another analytic method)
nc <- ncol(samples)
fits <- array(
data = NA,
dim = c(gridsize, nc, 2),
dimnames = list(NULL, colnames(samples), c("x", "y"))
)
n <- nrow(samples)
del0 <- (1/(4 * pi))^(1/10) # normal kernel
for (i in seq_len(nc)) {
# bandwidth definition from KernSmooth::bkde, but adjusted when sqrt(var(samples[, i])) is infinite
vx <- sqrt(var(samples[, i]))
if (is.infinite(vx)) # contWide can cause the variance to be infinite
vx <- sqrt(.Machine$double.xmax) # some big number as a fallback
bandwidth <- del0 * (243/(35 * n))^(1/5) * vx
fits[, i, ] <- do.call(cbind, KernSmooth::bkde(samples[, i], gridsize = gridsize, bandwidth = bandwidth))
}
if (getXRange) {
if (nc == 1L) {
xRanges <- matrix(range(fits[, , "x"]), 1L, 2L)
} else {
xRanges <- matrixStats::colRanges(fits[, , "x"])
}
rownames(xRanges) <- colnames(samples)
} else {
xRanges <- NULL
}
return(list(fit = fits, xRanges = xRanges))
}
.BANOVAapproxCRI <- function(x, y, cri = 0.95) {
# approximate cdf
y2 <- cumsum(y)
y2 <- y2 / y2[length(y2)]
h <- (1 - cri) / 2
# find closest observed match
idx <- c(
which.min(abs(y2 - h)),
which.min(abs(y2 - 1 + h))
)
return(x[idx])
}
# plot posteriors ----
.BANOVAgetBMAdensity <- function(samples, weights, fromTo = NULL, n = 2^10) {
# @param samples, list of samples
# @param weigths, vector of numeric weights
# @param fromTo vector of length 2 specifying lower and upper bound (optional).
# @return a list with $x the x-coordinates and $y the y-coordinates.
if (length(samples) != length(weights))
stop(gettext("length of samples must be equal to length of weights!"))
# remove NULL indices
idxNonNull <- lengths(samples) > 0
weights <- weights[idxNonNull]
samples <- samples[idxNonNull]
# renormalize the weights so output is a proper density function (TODO: do we want this, or BAS style?)
weights <- weights / sum(weights)
# create x-grid for density
if (is.null(fromTo))
fromTo <- range(sapply(samples, range))
xs <- seq(fromTo[1L], fromTo[2L], length.out = 2^10)
# compute weighted density
ys <- numeric(n)
for (i in seq_along(samples)) {
ys <- ys + weights[i] * density(samples[[i]], from = fromTo[1L], to = fromTo[2L], n = n)$y
}
return(list(x = xs, y = ys))
}
.BANOVAreSample <- function(n, x, y, prop0 = NULL) {
cdf <- cumsum(y)
max <- cdf[length(cdf)]
if (max == 0) # can't salvage this
return(numeric(n))
cdf <- cdf / max
if (!is.null(prop0) && prop0 <= (1 - sqrt(.Machine$double.eps))) {
i1 <- runif(n) <= prop0
samples <- numeric(n)
samples[i1] <- approx(cdf, x, runif(sum(i1)), rule = 2)$y
} else {
samples <- approx(cdf, x, runif(n), rule = 2)$y
}
return(samples)
}
.BANOVAgetSMIResidRsq <- function(posterior, dataset, formula, nIter, model, options) {
# @param posterior object from Bayesfactor package, or SxPx2 array of weighted densities
# @param dataset dataset
# @param formula the formula for this specific model. Supply the formula of the most complex model for BMA inference.
# @param nIter number of posterior samples
# @param model model estimates
# @param levelInfo for each variable, how many levels there are.
#
# @return a list with residuals, predictions, and r-squared
# @details the matrix multiplication in this function allocates an array of nobs * nsamples, which can be enormous.
# if is more memory efficient to use for loops (but slower in R) to calculate predictions for each observation
# and posterior sample. This could be done in the future if performance is an issue. However, this likely cannot be
# done efficiently in R. For modelAveraged inference, model and levelInfo need to be given.
# here we need the data in a one-hot encoded way
# first one is dependent variable, rest are independent variables
dvs <- all.vars(formula)[1L]
ivs <- all.vars(formula)[-1L]
# idx contains variables that need to be one-hot encoded
idx <- sapply(dataset[ivs], is.factor)
# idx <- isFactor & names(isFactor) %in% dvs
# do one-hot encoding (thank base R for calling a function `contrast` and giving it the argument `contrast`...)
datOneHot <- model.matrix(formula, data = dataset[c(dvs, ivs)],
contrasts.arg = lapply(dataset[ivs][idx], contrasts, contrasts = FALSE))
nobs <- nrow(datOneHot)
# possible switch for memory efficient alternatives:
# size of 1 numeric element in bytes
# doubleSize <- as.numeric(utils::object.size(double(2)) - utils::object.size(double(1)))
# 1073741824 is 1 GB in bytes
# 536870912 is 512 MB in bytes
# the size of an empty matrix is ignored here
# useMatmult <- ((doubleSize * nobs * nIter) / 536870912.0) <= 1
# if (useMatmult) {
if (length(dim(posterior)) == 3L) {
# we're doing model averaged inference
# recompute the levelInfo, this is relatively cheap and we cannot be sure that the levelInfo in the
# posterior was not retrieved from state (and contains too many or too few variables)
nmodels <- length(model[["model.list"]])
randomFactors <- .BANOVAgetRandomFactors(options, model[["analysisType"]])
dataTypes <- .BANOVAgetDataTypes(dataset, model[["model.list"]][[nmodels]], randomFactors)
levelInfo <- .BANOVAgetLevelInfo(dataset, model[["model.list"]][[nmodels]], dataTypes)
# which parameters are nuisance?
parameterNames <- c("mu", names(levelInfo[["levelNames"]]))
# the variables that are nuisance, e.g., contBinom
isNuisanceVar <- c(TRUE, parameterNames[-1L] %in% model[["nuisance"]])
names(isNuisanceVar) <- parameterNames
# the parameters that are nuisance, e.g., contBinom-0 and contBinom-1
isNuisanceParam <- rep(isNuisanceVar, c(1, lengths(levelInfo[["levelNames"]])))
names(isNuisanceParam) <- c("mu", unlist(levelInfo[["levelNames"]]))
if (!all(names(isNuisanceParam) %in% colnames(posterior)))
stop("!all(names(isNuisance) %in% colnames(posterior)) was not true.")
# we cannot use model[["effects"]] because that exludes nuisance parameters, so we rebuild one that doesn't exclude these
effects <- matrix(FALSE, nmodels, length(isNuisanceParam), dimnames = list(NULL, names(isNuisanceParam)))
for (i in seq_along(model[["model.list"]])) {
if (is.null(model[["model.list"]][[i]])) { # implies an intercept-only null model
cnms <- "mu"
} else {
idx <- names(levelInfo[["levelNames"]]) %in% attr(terms(model[["model.list"]][[i]]), "term.labels")
cnms <- c("mu", unlist(levelInfo[["levelNames"]][idx]))
}
effects[i, cnms] <- TRUE#c(TRUE, names(levelInfo[["levelNames"]]) %in% attr(terms(model[["model.list"]][[i]]), "term.labels"))
}
# sample from the BMA posterior
postProbs <- model[["postProbs"]]
postProbs[is.na(postProbs)] <- 0
.setSeedJASP(options)
modelIdx <- sample(length(postProbs), nIter, TRUE, postProbs) # resample the models according to posterior prob
samples <- matrix(0, nrow = nIter, ncol = ncol(datOneHot), dimnames = list(NULL, names(isNuisanceParam)))
# resample nuisance parameters
for (i in which(isNuisanceParam)) {
.setSeedJASP(options)
samples[, i] <- .BANOVAreSample(n = nIter, x = posterior[, i, "x"], y = posterior[, i, "y"])
}
# resample non nuisance parameters
for (i in which(!isNuisanceParam)) {
idx <- names(isNuisanceParam)[i]
mult <- effects[modelIdx, idx] # how often models with this variable are sampled
nTemp <- sum(mult)
if (nTemp > 0L){
.setSeedJASP(options)
samples[mult, i] <- .BANOVAreSample(n = sum(mult), x = posterior[, i, "x"], y = posterior[, i, "y"])
}
}
preds <- tcrossprod(datOneHot, samples)
} else {
# TODO: ensure that column order is always correct for both!
# otherwise we're doing inference for a single model
preds <- tcrossprod(datOneHot, posterior)
}
# calculate residuals (correctly recycles dat[[dvs]])
resid <- dataset[[dvs]] - preds
rsq <- .BANOVAcomputeRsq(dataset[[dvs]], preds)
# eta <- .BANOVAcomputeEtasq(dat[[dvs]], preds)
# } else {
# slow but memory efficient stuff
# }
return(list(resid = resid, preds = preds, rsq = rsq))
}
# HF computation ----
.BANOVAinitBayesFactor <- function() {
defaults <- list(
BFMaxModels = 50000,
BFpretestIterations = 100,
BFapproxOptimizer = "optim",
BFapproxLimits = c(-15, 15),
BFprogress = interactive(),
BFfactorsMax = 5
)
idx <- setdiff(names(defaults), names(options()))
options(defaults[idx])
}
.BANOVAcomputeRsq <- function(obs, preds) {
# NOTE: R^2 != cor(obs, predict) because the predictions from the posterior samples are not OLS estimates.
if (is.null(dim(preds)))
preds <- matrix(preds, ncol = 1L)
# definition from http://www.stat.columbia.edu/~gelman/research/unpublished/bayes_R2.pdf
colVars2 <- function(x) {matrixStats::colVars(sweep(x, 2L, matrixStats::colMeans2(x)))}
ee <- colVars2(preds)
ff <- colVars2(obs - preds)
gg <- ee / (ee + ff)
if (anyNA(gg)) {
idx <- is.na(gg)
if (all(idx)) {
# this means that either the posterior sampling failed or the input data is bogus.
.quitAnalysis(gettext("All predictions had zero variance. Check the predictors/DV for anomalies or increase the number of posterior samples."))
} else {
# this shouldn't happen, but is necessary because otherwise stuff will crash downstream.
gg[is.na(gg)] <- mean(gg, na.rm = TRUE)
}
}
return(gg)
}
.BANOVAcomputeEtasq <- function(obs, preds) {
# partial eta^2
if (is.null(dim(preds)))
preds <- matrix(preds, ncol = 1L)
# return(1 - matrixStats::colVars(preds - obs) / var(obs))
# return(1 - matrixStats::rowVars(preds - obs) / var(obs))
return(matrixStats::colVars(preds) / var(obs))
}
.BANOVAmakeChainNeater <- function(chains, Xnames, formula, data, dataTypes, gMap, unreduce, continuous, columnFilter) {
# so BayesFactor:::makeChainNeater doesn't always make things neater
# see also https://github.com/jasp-stats/jaspAnova/pull/43#discussion_r656430090
# example where BayesFactor goes wrong:
# dd <- structure(list(
# rm_factor = structure(c(1L, 2L, 1L, 2L, 1L, 2L), .Label = c("Level 1", "Level 2"), class = "factor"),
# dependent = c(6, 6, 5, 10, 6, 9),
# age = c(59.2161793866222, 59.2161793866222, 53.4457078686802, 53.4457078686802, 50.3248450015473, 50.3248450015473),
# subject = structure(c(1L, 1L, 2L, 2L, 3L, 3L), .Label = c("1", "2", "3"), class = "factor")),
# row.names = c(NA, 6L), class = "data.frame"
# )
#
# bf_obj_1 <- BayesFactor::lmBF(dependent ~ age + subject, data = dd, whichRandom = "subject")
# bf_obj_2 <- BayesFactor::lmBF(dependent ~ age + subject + rm_factor, data = dd, whichRandom = "subject")
#
# print(BayesFactor::posterior(bf_obj_1, iterations = 3)) # <- ugly
# print(BayesFactor::posterior(bf_obj_2, iterations = 3)) # <- good
#
# originalFun <- BayesFactor:::makeChainNeater
# jaspBase:::assignFunctionInPackage(.BANOVAmakeChainNeater, "makeChainNeater", "BayesFactor")
# print(BayesFactor::posterior(bf_obj_1, iterations = 3)) # <- good
#
# jaspBase:::assignFunctionInPackage(originalFun, "makeChainNeater", "BayesFactor")
# print(BayesFactor::posterior(bf_obj_1, iterations = 3)) # <- ugly
# this part is identical to BayesFactor:::makeChainNeater
P = length(gMap)
nGs = max(gMap) + 1
factors = BayesFactor:::fmlaFactors(formula, data)[-1]
dataTypes = dataTypes[names(dataTypes) %in% factors]
types = BayesFactor:::termTypes(formula, data, dataTypes)
lastPars = ncol(chains) + (-nGs):0
if (any(continuous)) {
gNames = paste("g", c(names(types[types != "continuous"]), "continuous"), sep = "_")
} else {
gNames = paste("g", names(types), sep = "_")
}
if (is.null(columnFilter)) {
ignoreCols = ignoreCols = rep(0, P)
} else {
ignoreCols = BayesFactor:::filterVectorLogical(columnFilter, names(gMap))
}
# here we're also checking here if any dataTypes are "random" instead of only checking for fixed dataTypes
if (!unreduce | !any(dataTypes %in% c("fixed", "random"))) {
labels = c("mu", Xnames[!ignoreCols], "sig2", gNames)
colnames(chains) = labels
return(chains)
}
parLabels = BayesFactor:::makeLabelList(formula, data, dataTypes, unreduce,
columnFilter)
labels = c("mu", parLabels)
betaChains = chains[, 1:(ncol(chains) - 2 - nGs) + 1, drop = FALSE]
betaChains = BayesFactor:::ureduceChains(betaChains, formula, data, dataTypes, gMap, ignoreCols)
newChains = cbind(chains[, 1], betaChains, chains[, lastPars])
labels = c(labels, "sig2", gNames)
colnames(newChains) = labels
return(newChains)
}
.BANOVAcomputeInteractionsMatrix <- function(effects) {
neffects <- length(effects)
interactions.matrix <- matrix(FALSE, nrow = neffects, ncol = neffects)
rownames(interactions.matrix) <- colnames(interactions.matrix) <- effects
if (neffects > 1L) {
effect.components <- sapply(effects, strsplit, split = ":", fixed = TRUE)
for (e in seq_len(neffects)) {
interactions.matrix[e, ] <- sapply(1:neffects, function(ee) {
(sum(effect.components[[e]] %in% effect.components[[ee]]) == length(effect.components[[e]]))
})
}
diag(interactions.matrix) <- FALSE
}
return(interactions.matrix)
}
# Model prior ----
.BANOVAcomputePriorModelProbs <- function(models, nuisance, options) {
if (options[["modelPrior"]] == "uniform") {
modelprobs <- rep(1 / length(models), length(models))
} else if (options[["modelPrior"]] == "custom") {
inclusionProbabilities <- vapply(options[["customPriorSpecification"]], `[[`, FUN.VALUE = numeric(1L), "inclusionProbability")
modelprobs <- .BANOVAcustomInclusionProbabilitiesToModelProbabilities(models, nuisance, inclusionProbabilities, enforceMarginality = options[["enforcePrincipleOfMarginalityFixedEffects"]])
} else {
noPredictorsPerModel <- vapply(models, function(x) if (is.null(x)) 0L else length(attr(terms(x), "term.labels")), FUN.VALUE = integer(1L))
noNuisancePredictors <- noPredictorsPerModel[1L]
noPredictorsPerModel <- noPredictorsPerModel - noNuisancePredictors
totalNoPredictors <- noPredictorsPerModel[length(noPredictorsPerModel)]
if (options[["modelPrior"]] %in% c("betaBinomial", "Wilson", "Castillo")) {
switch (options[["modelPrior"]],
"betaBinomial" = {alpha = options[["betaBinomialParameterA"]]; beta = options[["betaBinomialParameterB"]] },
"Wilson" = {alpha = 1.0; beta = totalNoPredictors * options[["wilsonParameterLambda"]] },
"Castillo" = {alpha = 1.0; beta = totalNoPredictors ^ options[["castilloParameterU"]] }
)
modelprobs <- dBetaBinomialModelPrior(noPredictorsPerModel, totalNoPredictors, alpha, beta)
} else if (options[["modelPrior"]] == "Bernoulli") {
modelprobs <- dBernoulliModelPrior(noPredictorsPerModel, totalNoPredictors, options[["bernoulliParameter"]])
}
# both the betabinomial and bernoulli model priors do not sum to 1 when marginality is respected
modelprobs <- modelprobs / sum(modelprobs)
}
return(modelprobs)
}
dBetaBinomialModelPrior <- function(k, n, alpha = 1.0, beta = 1.0, log = FALSE) {
# NOTE: this is not the pdf of the betabinomial, but the pdf divided by choose(n, k), so that it's the probabiltiy of a particular model
logprobability <- lbeta(k + alpha, n - k + beta) - lbeta(alpha, beta)
if (log)
return(logprobability)
return(exp(logprobability))
}
dBernoulliModelPrior <- function(k, n, prob = 0.5, log = FALSE) {
logprobability <- k * log(prob) + (n - k) * log(1 - prob)
if (log)
return(logprobability)
return(exp(logprobability))
}
.BANOVAcustomInclusionProbabilitiesToModelProbabilities <- function(modelList, nuisance, inclusionProbabilities, enforceMarginality = TRUE) {
# TODO: discuss with the team:
# when marginality is enforced, how should the prior on interaction effects be interpreted when a user wants matched models?
formulaFullModel <- modelList[[length(modelList)]]
termsFullModel <- terms(formulaFullModel)
termLabels <- attr(termsFullModel, "term.labels")
fullNonNuisance <- !(termLabels %in% nuisance)
termLabels <- termLabels[fullNonNuisance]
# the principle of marginality is enforced (or not) through the function getPresentInteractions
if (enforceMarginality) {
# in this case, interaction effects do not appear without the main effects, so the user provided probability
# is interpreted as a conditional probability given that the main effects are included.
# thus, p(interaction) = 0 if the main effects are missing, and the user provided probability otherwise.
#
# for example, consider p(A) = P(B) = .5, P(A:B) = .2
#
# model P(model)
# A P(A) * (1 - P(B)) # <- here we do NOT do '* (1 - P(A:B))'
# A, B p(A) * P(B) * (1 - P(A:B))
# A, B, A:B p(A) * P(B) * P(A:B)
termFactors <- attr(termsFullModel, "factors")[, termLabels, drop = FALSE]
termFactors <- termFactors[rownames(termFactors) %in% termLabels, , drop = FALSE]
rnmsTermFactors <- rownames(termFactors)
termOrders <- attr(termsFullModel, "order")[fullNonNuisance]
# create a list where the names refer to the interactions and the elements refer to all children of that interaction
listOfDescendants <- setNames(lapply(termLabels, function(x) {
idx <- which(termFactors[, x] == 1L)
rnmsTermFactors[idx]
}), termLabels)
termIsNotInteraction <- termOrders == 1L
getPresentInteractions <- function(currentTermLabels) {
vapply(listOfDescendants, function(descendants, currentTermLabels) {
all(descendants %in% currentTermLabels)
}, logical(1L), currentTermLabels) | termIsNotInteraction
}
# NOTE: the loop below fails for some combinations of user priors, but I can see how this behavior is desirable
# the users provide the raw p(interaction) however, the code below assumes that this is
# p(interaction | main effects), so to compensate we divide each interaction by the product of the components.
# however, for some inclusion probabilities of the main effects you get impossible values, for example:
# p(A_) = 0.1, p(B) = 0.1, P(A*B) = 0.9 => 0.9 / 0.1^2 = 90, which is bigger than 1 and should be impossible.
# for (i in which(!termIsNotInteraction)) { # loop over all interaction terms
# idxSubterms <- which(termLabels %in% listOfDescendants[[termLabels[i]]])
# inclusionProbabilities[i] <- inclusionProbabilities[i] / prod(inclusionProbabilities[idxSubterms])
# }
} else {
# in this case, interaction effects do appear without the main effects, so the user provided probability
# is interpreted in the same way for interaction effects and fixed effects.
# it's interpreted as an unconditional inclusion probability
getPresentInteractions <- function(...) {
rep(TRUE, length(termLabels))
}
}
exclusionProbabilties <- 1 - inclusionProbabilities
# so terms(y ~ b + a:b) reorders a:b to b:a but terms(y ~ a + a:b) does not...
termLabelOrder <- strsplit(termLabels, ":", fixed = TRUE)
modelProbs <- numeric(length(modelList))
for (i in seq_along(modelList)) {
if (is.null(modelList[[i]])) {
excludedTerms <- rep(TRUE, length(termLabels))
excludedTerms <- excludedTerms & getPresentInteractions(NULL)
modelProbs[i] <- prod(exclusionProbabilties)
} else {
currentTerms <- terms(modelList[[i]])
currentTermLabels <- setdiff(labels(currentTerms), nuisance)
currentTermLabels <- .BANOVAorderTermsByKnownOrder(currentTerms, currentTermLabels, termLabelOrder)
presentInteractions <- getPresentInteractions(currentTermLabels)
includedTerms <- termLabels %in% currentTermLabels
excludedTerms <- !includedTerms
includedTerms <- includedTerms & presentInteractions
excludedTerms <- excludedTerms & presentInteractions
modelProbs[i] <- prod(inclusionProbabilities[includedTerms], exclusionProbabilties[excludedTerms])
}
}
return(modelProbs)
}
# Other ----
.BANOVAgetRScale <- function(options, analysisType) {
if (options[["priorSpecificationMode"]] == "acrossParameters") {
rscaleFixed <- options[["cauchyPriorScaleFixedEffects"]]
rscaleRandom <- options[["cauchyPriorScaleRandomEffects"]]
rscaleCont <- if (analysisType == "ANOVA") "medium" else options[["cauchyPriorScaleCovariates"]]
rscaleEffects <- NULL
} else {
# NOTE: this still need a value otherwise BayesFactor returns NA for the bf
# default values of lmBF
rscaleFixed <- "medium"
rscaleRandom <- "nuisance"
rscaleCont <- if (analysisType == "ANOVA") "medium" else options[["cauchyPriorScaleCovariates"]]
rscaleEffectsNames <- vapply(options[["customPriorSpecification"]], FUN.VALUE = character(1L), function(x) {
paste(x[["components"]], collapse = ":")
})
rscaleEffects <- vapply(options[["customPriorSpecification"]], FUN.VALUE = numeric(1L), `[[`, "scaleFixedEffects")
rscaleEffectsKeep <- if (analysisType == "ANOVA") {
rep(TRUE, length(rscaleEffects))
} else {
vapply(options[["customPriorSpecification"]], FUN.VALUE = logical(1L), function(x) {
is.null(options[["covariates"]]) || !any(x[["components"]] %in% options[["covariates"]])
})
}
rscaleEffects <- rscaleEffects[rscaleEffectsKeep]
names(rscaleEffects) <- rscaleEffectsNames[rscaleEffectsKeep]
}
return(list(rscaleFixed = rscaleFixed, rscaleRandom = rscaleRandom, rscaleCont = rscaleCont, rscaleEffects = rscaleEffects))
}
.BANOVAgetShownTableSizeAndAddFootnote <- function(modelTable, nmodels, modelsShown, numModelsShown) {
if (modelsShown == "unlimited") {
shownTableSize <- nmodels
} else {
shownTableSize <- min(nmodels, numModelsShown)
if (shownTableSize < nmodels)
modelTable$addFootnote(gettextf("Showing the best %1$d out of %2$d models.", numModelsShown, nmodels))
}
return(shownTableSize)
}
# # .BANOVAupdateRscales <- function() {
#
# # for some reason, BayesFactor supports custom rscales for but throws a error if you try to set them for continuous parameters
# # this temporarily overwrite the stop call in this BayesFactor:::createRscales
# # the function will throw an error if BayesFactor:::createRscales is updated
#
# # originalCreateRscales <- BayesFactor:::createRscales
# # originalBody <- body(originalCreateRscales)
# # if (!identical(originalBody[[17]][[3]][[2]][[1]], quote(invisible))) {
# # if (!identical(originalBody[[17]][[3]][[2]][[1]], quote(stop)))
# # stop("Bayesfactor got an update, and createRscales/ .BANOVAupdateRscales needs to be fixed!", domain = NA)
# # originalBody[[17L]][[3L]][[2L]][[1L]] <- substitute(invisible)
# # newCreateRscales <- originalCreateRscales
# # body(newCreateRscales) <- originalBody
# #
# # jaspBase::assignFunctionInPackage(newCreateRscales, "createRscales", "BayesFactor")
# # }
# #
# # tmp <- methods::getMethod(BayesFactor::compare, list(numerator="BFlinearModel", denominator="missing", data="data.frame"))
# # fun <- tmp@.Data
# # newfun <- fun
# # body <- body(fun)
# # if (!identical(body[[19L]][[2L]][[2L]][[4L]][[2L]], quote(FALSE))) {
# # if (!identical(body[[19L]][[2L]][[2L]][[4L]][[2L]], quote(all(relevantDataTypes == "continuous"))))
# # stop("Bayesfactor got an update, and BayesFactor::`.__T__compare:BayesFactor`$`BFlinearModel#missing#data.frame`/ .BANOVAupdateRscales needs to be fixed!", domain = NA)
# # body[[19L]][[2L]][[2L]][[4L]][[2L]] <- substitute(FALSE) # do not take the fast path
# # body(newfun) <- body
# #
# # # horrible but necessary
# # env <- environment()#getNamespace(jaspAnova)
# # assign("compare", BayesFactor::compare, envir = env)
# # setMethod("compare", signature(numerator = "BFlinearModel", denominator = "missing", data = "data.frame"), newfun, where = env)
#
# # env <- getNamespace("BayesFactor")
# # unlockBinding("compare", env)
#
# # debugonce(setMethod)
# # setMethod("compare", signature(numerator = "BFlinearModel", denominator = "missing", data = "data.frame"), newfun, where = env)
#
# # tmp@.Data <- newfun
# # s4env <- BayesFactor:::`.__T__compare:BayesFactor`
# # unlockBinding("BFlinearModel#missing#data.frame", s4env)
# # s4env$`BFlinearModel#missing#data.frame` <- tmp
# # lockBinding("BFlinearModel#missing#data.frame", s4env)
# # jaspBase::assignFunctionInPackage(s4env, ".__T__compare:BayesFactor", "BayesFactor")
# # }
#
# # return()
#
# # we could reset the changes with on.exit, but I'm not sure if it's necessary and doing that cleanly requires withr::defer and adds a dependency
# # do.call(on.exit({
# # jaspBase::assignFunctionInPackage(tempRscales[["originalCreateRscales"]], "createRscales", "BayesFactor")
# #
# # s4env <- BayesFactor:::`.__T__compare:BayesFactor`
# # unlockBinding("BFlinearModel#missing#data.frame", s4env)
# # s4env$`BFlinearModel#missing#data.frame` <- newfun
# # lockBinding("BFlinearModel#missing#data.frame", s4env)
# # jaspBase::assignFunctionInPackage(s4env, ".__T__compare:BayesFactor", "BayesFactor")
# # }, add = TRUE, after = FALSE)
# }
# Dependencies ----
.BANOVAdataDependencies <- function() {
# seed is in here because basically everything depends on the seed
c("dependent", "randomFactors", "covariates", "fixedFactors", "betweenSubjectFactors", "repeatedMeasuresFactors", "repeatedMeasuresCells", "modelTerms",
"seed", "setSeed")
}
.BANOVAmodelSpaceDependencies <- function(modelPrior) {
c("modelPrior", "betaBinomialParameterA", "betaBinomialParameterB", "wilsonParameterLambda", "castilloParamerU", "enforcePrincipleOfMarginalityFixedEffects", "enforcePrincipleOfMarginalityRandomSlopes",
if (modelPrior == "custom") "customPriorSpecification"
)
}
.BANOVArscaleDependencies <- function(priorSpecificationMode) {
c("cauchyPriorScaleFixedEffects", "cauchyPriorScaleRandomEffects", "cauchyPriorScaleCovariates", "priorSpecificationMode",
if (priorSpecificationMode == "perTerm") "customPriorSpecification"
)
}
.BANOVAmodelPriorOptionsChanged <- function(state, options) {
!identical(state[["modelPriorOptions"]][.BANOVAmodelSpaceDependencies(options[["modelPrior"]])], options[.BANOVAmodelSpaceDependencies(options[["modelPrior"]])])
}
.BANOVAmarginalityOptionsUnchanged <- function(state, options) {
stateOpts <- state[["modelPriorOptions"]]
stateOpts[["enforcePrincipleOfMarginalityFixedEffects"]] == options[["enforcePrincipleOfMarginalityFixedEffects"]] &&
stateOpts[["enforcePrincipleOfMarginalityRandomSlopes"]] == options[["enforcePrincipleOfMarginalityRandomSlopes"]]
}
.BANOVAmodelBFTypeOrOrderChanged <- function(state, options) {
nms <- c("fixedFactors", "modelTerms", "randomFactors", "covariates", "seed", "setSeed")
nms <- intersect(names(options), nms) # excludes covariates for ANOVA
identical(state[nms], options[nms])
}
.BANOVAmodelTermsUnchanged <- function(state, options) {
identical(state[["modelTerms"]], options[["modelTerms"]])
}
# HF formulas ----
.BANOVAgetFormulaComponents <- function(x, what = c("components", "variables")) {
what <- match.arg(what)
if (what == "components") {
return(colnames(attr(terms(x), "factors")))
} else {
return(all.vars(x)[-1L])
}
}
.BANOVAgetModelTitlesWithAllTerms <- function(modelObjects, modelList, analysisType, hideNuisance) {
if (hideNuisance)
return(vapply(modelObjects, FUN = `[[`, FUN.VALUE = character(1L), "title"))
nullModelName <- gettext("Null model")
res <- gsub("(.*)~\\s+", "", vapply(modelList, function(x) if (is.null(x)) nullModelName else .BANOVAas.character.formula(x), character(1L)))
if (analysisType == "RM-ANOVA")
res <- gsub(.BANOVAsubjectName, "subject", res)
res[res == ""] <- nullModelName
return(jaspBase::gsubInteractionSymbol(res))
}
.BANOVAgenerateAllModelFormulas <- function(formula, nuisance = NULL, analysisType = "RM-ANOVA",
enforcePrincipleOfMarginalityFixedEffects = TRUE,
enforcePrincipleOfMarginalityRandomSlopes = FALSE,
rmFactors = NULL, legacy = FALSE
) {
# TODO: it might be nicer to represent the NULL model as
# y ~ 1, rather than NULL. Right now the null model is a formula
# when there are nuisance variables and otherwise NULL, which leads to a bunch
# of annoying if (is.null(model[i])) exceptions.
neverExclude <- paste0("^", nuisance, "$")
if (!legacy && analysisType == "RM-ANOVA" && is.null(rmFactors))
stop(".BANOVAgenerateAllModelFormulas called with invalid arguments: analysisType = \"RM-ANOVA\", legacy = FALSE, rmFactors = NULL", domain = NA)
modelSpace <- try(BayesFactor::enumerateGeneralModels(
formula,
whichModels = if (enforcePrincipleOfMarginalityFixedEffects) "withmain" else "all",
neverExclude = neverExclude)
)
nuisanceRandomSlopes <- NULL
if (isTryError(modelSpace))
.quitAnalysis(gettextf("The following error occured in BayesFactor::enumerateGeneralModels: %s",
.extractErrorMessage(modelSpace)))
if (analysisType == "RM-ANOVA" && !legacy) {
# adjust the nuisance to include the random slopes, only done here because BayesFactor::enumerateGeneralModels
# does not handle this properly
allPossibleSlopes <- labels(stats::terms(stats::as.formula(paste("y~", paste(rmFactors, collapse = "*")))))
# drop the most complex interaction
allPossibleSlopes <- allPossibleSlopes[-length(allPossibleSlopes)]
# if at least one interaction among the repeated measures is considered
if (length(allPossibleSlopes) > 0L) {
# add interaction with subject
nuisanceRandomSlopes <- paste0(allPossibleSlopes, ":", .BANOVAsubjectName)
if (enforcePrincipleOfMarginalityRandomSlopes) {
for (i in seq_along(modelSpace)) {
presentLabels <- labels(stats::terms(modelSpace[[i]]))
termsToAddChar <- intersect(allPossibleSlopes, presentLabels)
if (length(termsToAddChar) > 0L) {
termsToAddChar <- paste0(termsToAddChar, ":", .BANOVAsubjectName)
termsToAdd <- as.formula(paste("~ . +", paste0(termsToAddChar, collapse = " + ")))
modelSpace[[i]] <- update.formula(modelSpace[[i]], new = termsToAdd)
}
}
} else {
termsToAdd <- as.formula(paste("~ . +", paste0(nuisanceRandomSlopes, collapse = " + ")))
modelSpace <- lapply(modelSpace, update.formula, new = termsToAdd)
}
# for the reordering done below. termsToAdd always contains the most complex random effects
formula <- update.formula(formula, new = termsToAdd)
# update nuisance
nuisance <- c(nuisance, nuisanceRandomSlopes)
}
}
if (is.null(nuisance)) {
return(list(modelList = c(list(NULL), modelSpace), nuisance = nuisance, nuisanceRandomSlopes = nuisanceRandomSlopes))
} else {
# put the null-model first
i <- length(modelSpace)
modelSpace <- c(modelSpace[[i]], modelSpace[-i])
# BayesFactor has the nasty habit of changing the order of interactions whenever one of the components is
# a nuisance variable. Here we ensure the order matches that of the input formula (which matches the order a user entered the factors).
# see also https://github.com/jasp-stats/jasp-test-release/issues/1181
originalTerms <- terms(formula)
originalOrd <- attr(originalTerms, "order")
originalIdxNonInteraction <- which(originalOrd == 1L)
originalLabels <- attr(originalTerms, "term.labels")
originalLabelsPieces <- strsplit(originalLabels[originalOrd > 1L], ":")
originalLabelsPiecesSorted <- lapply(originalLabelsPieces, sort)
dependent <- all.vars(modelSpace[[1]])[1L]
for (i in seq_along(modelSpace)) {
term <- terms(modelSpace[[i]])
ord <- attr(term, "order") # interaction order
idxNonInteraction <- which(ord == 1L)
termLabels <- attr(term, "term.labels")
termLabels[idxNonInteraction] <- originalLabels[originalIdxNonInteraction][match(termLabels[idxNonInteraction], originalLabels[originalIdxNonInteraction])]
for (j in which(ord > 1L)) {
labelPieces <- strsplit(termLabels[j], ":")[[1L]]
sortedLabelPieces <- sort(labelPieces)
for (k in seq_along(originalLabelsPieces))
if (identical(sortedLabelPieces, originalLabelsPiecesSorted[[k]]))
termLabels[j] <- paste(originalLabelsPieces[[k]], collapse = ":")
}
modelSpace[[i]] <- as.formula(paste(dependent, "~", paste(termLabels, collapse = " + ")))
}
return(list(modelList = modelSpace, nuisance = nuisance, nuisanceRandomSlopes = nuisanceRandomSlopes))
}
}
.BANOVAcreateModelFormula <- function(dependent, modelTerms) {
rhs <- "" # right hand side of the formula
nuisance <- NULL
effects <- NULL
for (term in modelTerms) {
comp <- term[["components"]]
newValue <- paste(comp, collapse = ":")
rhs <- if (rhs != "") paste0(rhs, " + ", newValue) else newValue
if (!is.null(term[["isNuisance"]]) && term[["isNuisance"]]) {
nuisance <- c(nuisance, newValue)
} else {
effects <- c(effects, newValue)
}
}
model.formula <- formula(paste(dependent, "~", rhs))
# this would be cleaner ideal if BayesFactor::enumerateGeneralModels would handle the nuisance properly.
# if (isRMANOVA && !legacy) {
# randomSlopes <- paste0(rmFactors, ":", .BANOVAsubjectName)
# termsToAdd <- as.formula(paste("~ . +", paste0(randomSlopes, collapse = " + ")))
# model.formula <- update.formula(model.formula, termsToAdd)
# nuisance <- c(nuisance, randomSlopes)
# }
return(list(model.formula = model.formula, nuisance = nuisance, effects = effects))
}
.BANOVAgetModelFormulaFromBFobj <- function(BayesFactorObj, asCharacter = FALSE) {
out <- BayesFactorObj$bf@numerator[[1L]]@identifier$formula
if (asCharacter) {
return(out)
} else {
return(as.formula(out))
}
}
.BANOVAgetLevelInfo <- function(dataset, formula, dataTypes) {
levelNames <- .BANOVAmakeLabelList(formula, dataset, dataTypes)
levelCounts <- lengths(levelNames) # counts of the factors including interaction terms
return(list(levelCounts = levelCounts, levelNames = levelNames))
}
.BANOVAmakeLabelList <- function(formula, data, dataTypes, unreduce = TRUE, columnFilter = NULL) {
# from BayesFactor 0.9.12.4.2
# identical to BayesFactor:::makeLabelList, except that we return a named list
# with
# names: the variables that make up a particular (interaction) effect
# values: the parameter names as returned by BayesFactor::posterior(...)
terms = attr(terms(formula, data = data), "term.labels")
if (!is.null(columnFilter))
terms = terms[!BayesFactor:::filterVectorLogical(columnFilter, terms)]
if (unreduce)
dataTypes[dataTypes == "fixed"] = "random"
labelList = lapply(terms, function(term, data, dataTypes) {
effects = strsplit(term, ":", fixed = TRUE)[[1]]
my.names = BayesFactor:::design.names.intList(effects, data, dataTypes)
return(paste(term, "-", my.names, sep = ""))
}, data = data, dataTypes = dataTypes)
# this part is different from BayesFactor
names(labelList) <- terms
return(labelList)
}
.BANOVAgetDataTypes <- function(dataset, formula, whichRandom = NULL) {
# from BayesFactor:::lmBF
BayesFactor:::createDataTypes(
formula,
whichRandom = whichRandom,
data = dataset,
analysis = "lm"
)
}
.BANOVAgetRandomFactors <- function(options, analysisType) {
if (analysisType == "RM-ANOVA")
return(.BANOVAsubjectName)
else
return(unlist(options[["randomFactors"]]))
}
.BANOVAgetLevelsFromParamNames <- function(names) {
# NOTE: this works because base64 does not contain "-"
# split on first "-"; 2 implies output of length 2, i.e., only split once
out <- do.call(rbind, stringr::str_split(names, "-", 2L))
# continous variables have as level name the variable name
idx <- out[, 1L] == out[, 2L]
out[idx, 2L] <- ""
# change dots into spaces for aesthetic purposes
out[, 2L] <- gsub(".", " ", out[, 2L], fixed = TRUE)
colnames(out) <- c("parameter", "level")
return(out)
}
.BANOVAas.character.formula <- function(x, ...) {
# we could also extend the S3 function as.character
Reduce(paste, trimws(deparse(x)))
}
.BANOVAreorderFormulas <- function(x) {
# This function reorders the terms of a formula such that they are alphabetical
# e.g., a ~ c + b + c:b becomes a ~ b + b:c
# This is necessary because BayesFactor::enumerateGeneralModels always appends the nuisance terms
# and a ~ b + c != a ~ c + b
# so without this function the state does not get reused whenever a user modifies the nuisance terms
if (is.null(x))
return("NULL") # as a string so it can become one big character vector for match()
s <- strsplit(attr(stats::terms.formula(x), "term.labels"), ":")
for (i in which(lengths(s) > 1L))
s[[i]] <- paste0(sort(s[[i]]), collapse = ":")
return(paste(all.vars(x)[1L], "~", paste(sort(unlist(s)), collapse = " + ")))
}
.BANOVAreorderTerms <- function(x) {
# This function reorders the terms of a formula such that they are alphabetical
# it essentially does the same as .BANOVAreorderFormulas but expects x to be character
# x should be a character string of the form c("a" "a:b"), so not a ~ b + c:d.
# For example, output from `all.vars(formula)` or `.BANOVAgetFormulaComponents(formula)`.
s <- strsplit(x, ":")
for (i in which(lengths(s) > 1L))
s[[i]] <- paste0(sort(s[[i]]), collapse = ":")
return(unlist(s))
}
.BANOVAorderTermsByKnownOrder <- function(currentTerms, currentTermLabels, termLabelOrder) {
# This function reorders the terms of a formula such that they follow the order of termLabelOrder
# currentTerms: output of stats::terms(formula)
# currentTermLabels: output of labels(stats::formula(formula)), possibly after filtering out nuisance
# termLabelOrder: character vector of term labels in the desired order
currentFactors <- attr(currentTerms, "factors")
if (any(currentFactors > 1L)) { # TRUE implies may need to reorder some terms according to termLabelOrder
for (j in grep(":", currentTermLabels, fixed = TRUE)) {
tmp <- strsplit(currentTermLabels[j], ":", fixed = TRUE)[[1L]]
for (order in termLabelOrder) {
if (all(tmp %in% order) && length(tmp) == length(order)) {
currentTermLabels[j] <- paste(order, collapse = ":")
break
}
}
}
}
return(currentTermLabels)
}
# Single Model Inference (SMI) ----
.BANOVAsmi <- function(jaspResults, dataset, options, model) {
userWantsSMI <- any(unlist(options[c(
"singleModelPosteriorPlot", "singleModelQqPlot", "singleModelRsqPlot", "singleModelEstimates", "singleModelCriTable"
)]))
if (!userWantsSMI)
return()
if (!is.null(jaspResults[["containerSingleModel"]])) {
singleModelContainer <- jaspResults[["containerSingleModel"]]
} else {
singleModelContainer <- createJaspContainer(title = gettext("Single Model Inference"))
singleModelContainer$dependOn(c(
"singleModelTerms", "dependent", "samplingMethodMCMC", "samplesMCMC",
"repeatedMeasuresCells", "seed", "setSeed",
.BANOVArscaleDependencies(options[["priorSpecificationMode"]])
))
jaspResults[["containerSingleModel"]] <- singleModelContainer
singleModelContainer$position <- 7
}
singleModel <- jaspResults[["singleModelState"]]$object
if (!is.null(model[["models"]]) && is.null(singleModel) && length(options$singleModelTerms) > 0L && userWantsSMI) {
singleModel <- try(.BANOVAsmiSamplePosterior(dataset, options, model[["analysisType"]]))
if (isTryError(singleModel)) {
singleModelContainer$setError(gettextf("Error in single model inference:\n%s", singleModel))
singleModel <- NULL
} else {
singleModelState <- createJaspState(object = singleModel)
singleModelState$dependOn(optionsFromObject = singleModelContainer)
jaspResults[["singleModelState"]] <- singleModelState
}
}
.BANOVAsmiEstimates(singleModelContainer, options, singleModel)
.BANOVAsmiRsquared (singleModelContainer, options, singleModel)
.BANOVAsmiQqPlot (singleModelContainer, options, singleModel)
.BANOVAsmiRsqPlot (singleModelContainer, options, singleModel)
.BANOVAsmiPosteriorPlot(singleModelContainer, dataset, options, singleModel)
return()
}
.BANOVAsmiSamplePosterior <- function(dataset, options, analysisType) {
nIter <- if (options[["samplingMethodMCMC"]] == "auto") 1e3L else options[["samplesMCMC"]]
modelTerms <- options$singleModelTerms
dependent <- options$dependent
if (analysisType == "RM-ANOVA") {
modelTerms[[length(modelTerms) + 1L]] <- list(components = .BANOVAsubjectName, isNuisance = TRUE)
dependent <- .BANOVAdependentName
}
tempRScale <- .BANOVAgetRScale(options, analysisType)
rscaleFixed <- tempRScale[["rscaleFixed"]]
rscaleRandom <- tempRScale[["rscaleRandom"]]
rscaleCont <- tempRScale[["rscaleCont"]]
rscaleEffects <- tempRScale[["rscaleEffects"]]
formula <- .BANOVAcreateModelFormula(dependent, modelTerms)$model.formula
randomFactors <- .BANOVAgetRandomFactors(options, analysisType)
dataTypes <- .BANOVAgetDataTypes(dataset, formula, randomFactors)
levelInfo <- .BANOVAgetLevelInfo(dataset, formula, dataTypes)
# if all variables are continuous, then continuous variables keep their name. E.g., "contGamma" stays "contGamma".
# if any variables are factors, then continuous variables names are doubled. E.g., "contGamma" becomes "contGamma-contGamma".
allParamNames <- c(
"mu",
if (all(dataTypes == "continuous")) names(levelInfo$levelNames) else unlist(levelInfo$levelNames)
)
.setSeedJASP(options)
samples <- BayesFactor::lmBF(
formula = formula,
data = dataset,
whichRandom = unlist(randomFactors),
progress = TRUE,
posterior = TRUE,
rscaleFixed = rscaleFixed,
rscaleRandom = rscaleRandom,
rscaleCont = rscaleCont,
rscaleEffects = rscaleEffects,
iterations = nIter
)
types <- samples@model@dataTypes
# keep only relevant columns, drop sig2, g_xxx, ...
idx <- match(allParamNames, colnames(samples), nomatch = 0L)
samples <- samples[, idx, drop = FALSE]
# for some odd reason, Bayesfactor uses as column name contcor1-contcor1
# if there is a covariate AND fixed factors, but only contcor1 if all variables are continuous...
allContinuous <- all(types == "continuous")
if (allContinuous) {
cnms <- colnames(samples)[-1L] # omit the intercept (mu) which is not changed by Bayesfactor
colnames(samples)[-1L] <- paste0(cnms, "-", cnms)
}
# find out which parameters are random
isRandom <- logical(ncol(samples))
idx <- which(types == "random")
for (i in idx)
isRandom <- isRandom | startsWith(colnames(samples), names(types)[i])
means <- colMeans(samples)
sds <- matrixStats::colSds(samples)
names(means) <- names(sds) <- colnames(samples)
h <- (1 - options[["credibleInterval"]]) / 2
probs <- c(h, 1-h)
cri <- matrixStats::colQuantiles(samples, probs = probs)
densities <- .BANOVAfitDensity(samples, 2^9, FALSE)
tmp <- .BANOVAgetSMIResidRsq(samples, dataset, formula, nIter, options = options)
residmeans <- rowMeans(tmp$resid)
residquants <- matrixStats::rowQuantiles(tmp$resid, probs = c(h, 1-h))
# all information for q-q plot of residuals
residSumStats <- matrix(c(residmeans, residquants), nrow = length(residmeans), ncol = 3L,
dimnames = list(NULL, c("mean", "cri.2.5%", "cri.97.5%")))
# all information for r-squared density plot
rsqDens <- density(tmp$rsq, n = 2^11, from = 0, to = 1)
rsqCri <- quantile(tmp$rsq, probs = probs)
rsqMean <- mean(tmp$rsq)
return(list(
means = means, sds = sds, CRIs = cri, densities = densities$fit,
residSumStats = residSumStats, rsqDens = rsqDens, rsqCri = rsqCri, rsqMean = rsqMean,
allContinuous = allContinuous, isRandom = isRandom
))
}
.BANOVAsmiEstimates <- function(jaspContainer, options, model) {
if (!is.null(jaspContainer[["SMItablePosteriorEstimates"]]) || !options[["singleModelEstimates"]])
return()
estsTable <- createJaspTable(title = gettext("Single Model Posterior Summary"))
estsTable$position <- 1
estsTable$dependOn(c("singleModelEstimates", "credibleInterval"))
overTitle <- gettextf("%s%% Credible Interval", format(100 * options[["credibleInterval"]], digits = 3))
estsTable$addColumnInfo(name = "Variable", title = gettext("Variable"), type = "string")
estsTable$addColumnInfo(name = "Level", title = gettext("Level"), type = "string")
estsTable$addColumnInfo(name = "Mean", title = gettext("Mean"), type = "number")
estsTable$addColumnInfo(name = "SD", title = gettext("SD"), type = "number")
estsTable$addColumnInfo(name = "Lower", title = gettext("Lower"), type = "number", overtitle = overTitle)
estsTable$addColumnInfo(name = "Upper", title = gettext("Upper"), type = "number", overtitle = overTitle)
if (!(is.null(model) || estsTable$getError())) {
.BANOVAfillEstimatesTable(
jaspTable = estsTable,
mus = model$means,
sds = model$sds,
cri = model$CRIs,
hasNoLevels = model$allContinuous,
isRandom = model$isRandom
)
}
jaspContainer[["SMItablePosteriorEstimates"]] <- estsTable
return()
}
.BANOVAsmiPosteriorPlot <- function(jaspContainer, dataset, options, model) {
# meta wrapper for model averaged posterior plots, single model posterior plots, and Q-Q plots
if (!is.null(jaspContainer[["SMIposteriorPlot"]]) || !options$singleModelPosteriorPlot)
return()
posteriorPlotContainer <- createJaspContainer(title = gettext("Single Model Posterior Distributions"))
jaspContainer[["SMIposteriorPlot"]] <- posteriorPlotContainer
posteriorPlotContainer$position <- 2
posteriorPlotContainer$dependOn(c("singleModelPosteriorPlot", "singleModelGroupPosterior"))
if (is.null(model) || posteriorPlotContainer$getError()) {
posteriorPlotContainer[["dummyplot"]] <- createJaspPlot(title = gettext("Posterior distribution"), width = 400, height = 400,
plot = NULL)
} else {
.BANOVAfillPosteriorPlotContainer(
container = posteriorPlotContainer,
densities = model$densities[, -1L, ], # omit intercept
cris = model$CRIs[-1L, ], # omit intercept
isRandom = model$isRandom[-1L], # omit intercept
groupParameters = identical(options[["groupPosterior"]], "grouped")
)
}
return()
}
.BANOVAsmiQqPlot <- function(jaspContainer, options, model) {
if (!is.null(jaspContainer[["QQplot"]]) || !options$singleModelQqPlot)
return()
if (is.null(model) || jaspContainer$getError()) {
p <- NULL
} else {
p <- jaspGraphs::plotQQnorm(
residuals = model$residSumStats[,"mean"],
lower = model$residSumStats[,"cri.2.5%"],
upper = model$residSumStats[,"cri.97.5%"],
ablineColor = "darkred"
)
}
plot <- createJaspPlot(
title = gettext("Q-Q Plot"),
width = 400,
height = 400,
plot = p,
aspectRatio = 1
)
plot$dependOn("singleModelQqPlot")
plot$position <- 3
jaspContainer[["QQplot"]] <- plot
return()
}
.BANOVAsmiRsqPlot <- function(jaspContainer, options, model) {
if (!is.null(jaspContainer[["smirsqplot"]]) || !options$singleModelRsqPlot)
return()
if (is.null(model) || jaspContainer$getError()) {
p <- NULL
} else {
dd <- model$rsqDens
rsqCri <- model$rsqCri
df <- data.frame(x = dd$x, y = dd$y)
xName <- expression(R^2)
p <- jaspGraphs::PlotPriorAndPosterior(dfLines = df, xName = xName, CRI = rsqCri, drawCRItxt = FALSE)
}
plot <- createJaspPlot(
title = gettextf("Posterior R%s", "\u00B2"),
width = 400,
height = 400,
plot = p,
aspectRatio = 1
)
plot$dependOn("singleModelRsqPlot")
plot$position <- 4
jaspContainer[["smirsqplot"]] <- plot
return()
}
.BANOVAsmiRsquared <- function(jaspResults, options, model) {
if (!is.null(jaspResults[["tableSMICRI"]]) || !options[["singleModelCriTable"]])
return()
criTable <- createJaspTable(title = gettextf("Single Model R%s", "\u00B2"))
criTable$position <- 3.5
criTable$dependOn(c("singleModelCriTable", "credibleInterval"))
overTitle <- gettextf("%s%% Credible Interval", format(100 * options[["credibleInterval"]], digits = 3))
criTable$addColumnInfo(name = "rsq", title = "", type = "string")
criTable$addColumnInfo(name = "Mean", title = gettext("Mean"), type = "number")
criTable$addColumnInfo(name = "Lower", title = gettext("Lower"), type = "number", overtitle = overTitle)
criTable$addColumnInfo(name = "Upper", title = gettext("Upper"), type = "number", overtitle = overTitle)
if (!is.null(model)) {
cri <- model[["rsqCri"]]
df <- data.frame(
rsq = "R\u00B2",
Mean = model[["rsqMean"]],
Lower = cri[1L],
Upper = cri[2L],
row.names = NULL
)
criTable$setData(df)
} else {
criTable[["rsq"]] <- "R\u00B2"
}
jaspResults[["tableSMICRI"]] <- criTable
return()
}
# Citations ----
.BANOVAcitations <- c(
"MoreyEtal2015" = "Morey, R. D. & Rouder, J. N. (2015). BayesFactor (Version 0.9.10-2)[Computer software].",
"RouderEtal2012" = "Rouder, J. N., Morey, R. D., Speckman, P. L., Province, J. M., (2012) Default Bayes Factors for ANOVA Designs. Journal of Mathematical Psychology. 56. p. 356-374.",
"Jeffreys1938" = "Jeffreys, H. (1938). Significance tests when several degrees of freedom arise simultaneously. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 165, 161-198.",
"WestfallEtal1997" = "Westfall, P. H., Johnson, W. O., & Utts, J. M. (1997). A Bayesian perspective on the Bonferroni adjustment. Biometrika, 84, 419-427."
)
jasp-stats/jaspAnova documentation built on April 5, 2025, 3:44 p.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 .BANOVAestimateModels .BANOVAmissingInteractionCells .BANOVAerrorhandling .BANOVArunAnalysis
#
# Copyright (C) 2013-2022 University of Amsterdam
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Abbreviations:
# HF = helper function
.BANOVArunAnalysis <- function(jaspResults, dataset, options, analysisType = c("ANOVA", "ANCOVA", "RM-ANOVA")) {
# the main workhorse of the Bayesian ANOVA, ANCOVA, and Repeated Measures.
analysisType <- match.arg(analysisType)
dataset <- .BANOVAreadData(dataset, options, analysisType)
if (isTryError(dataset))
.quitAnalysis(gettext("Error while loading data. Please verify your repeated measures observations."))
dataInfo <- .BANOVAerrorhandling(dataset, options, analysisType)
model <- .BANOVAestimateModels(jaspResults, dataset, options, dataInfo, analysisType)
model[["posteriors"]] <- .BANOVAestimatePosteriors(jaspResults, dataset, options, model)
.BANOVAeffectsTable (jaspResults, options, model)
.BANOVAestimatesTable(jaspResults, options, model)
.BANOVArsquaredTable (jaspResults, options, model)
# model averaged plots
.BANOVAposteriorPlot(jaspResults, dataset, options, model)
.BANOVAqqplot (jaspResults, options, model)
.BANOVArsqplot (jaspResults, options, model)
.BANOVAnullControlPostHocTable(jaspResults, dataset, options, model)
# single model plots and tables
.BANOVAsmi(jaspResults, dataset, options, model)
# descriptives
.BANOVAdescriptives(jaspResults, dataset, options, dataInfo, analysisType)
return()
}
.BANOVAerrorhandling <- function(dataset, options, analysisType) {
customChecks <- list()
if (analysisType != "RM-ANOVA") {
hasDV <- options$dependent != ""
hasIV <- any(lengths(options[c("fixedFactors", "covariates")]) != 0)
fixed <- options$fixedFactors
noVariables <- !(hasDV && hasIV)
target <- c(options$covariates, options$dependent)
# check if the last model has an interaction effect
mostComplexModel <- NULL
if (length(options[["modelTerms"]]) > 0L) {
mostComplexModel <- options[["modelTerms"]][[length(options[["modelTerms"]])]][["components"]]
mostComplexModel <- intersect(mostComplexModel, fixed)
}
if (length(mostComplexModel) > 1L)
customChecks[["missingInteractionCells"]] <- .BANOVAmissingInteractionCells
if (!noVariables) {
.hasErrors(
dataset = dataset,
type = c("infinity", "observations", "variance", "factorLevels", "duplicateColumns"),
infinity.target = target,
variance.target = target,
observations.target = target,
observations.amount = paste("<", length(options$modelTerms) + 1),
factorLevels.target = fixed,
factorLevels.amount = " < 2",
# custom checks
custom = customChecks,
missingInteractionCells.target = mostComplexModel,
exitAnalysisIfErrors = TRUE
)
}
} else {
hasDV <- !any(options$repeatedMeasuresCells == "")
hasIV <- any(lengths(options[c("betweenSubjectFactors", "covariates")]) != 0)
fixed <- unlist(options$betweenSubjectFactors)
covariates <- unlist(options$covariates)
noVariables <- !hasDV
target <- c(covariates, fixed)
# check if the last model has an interaction effect
mostComplexModel <- NULL
if (length(options[["modelTerms"]]) > 0L) {
mostComplexModel <- options[["modelTerms"]][[length(options[["modelTerms"]])]][["components"]]
mostComplexModel <- intersect(mostComplexModel, fixed)
}
if (length(mostComplexModel) > 1L)
customChecks[["missingInteractionCells"]] <- .BANOVAmissingInteractionCells
if (length(target) > 0L) {
.hasErrors(
dataset = dataset,
type = c("infinity", "observations", "variance", "factorLevels", "duplicateColumns"),
infinity.target = target,
variance.target = covariates,
observations.target = target,
observations.amount = paste("<", length(options$modelTerms) + 1),
factorLevels.target = fixed,
factorLevels.amount = " < 2",
duplicateColumns.target = target,
# custom checks
custom = customChecks,
missingInteractionCells.target = mostComplexModel,
exitAnalysisIfErrors = TRUE
)
}
if (!noVariables) {
# messages from commonmessages.R
customChecks <- list(
infinity = function() {
x <- dataset[, .BANOVAdependentName]
if (any(is.infinite(x)))
return(gettext("Infinity found in repeated measures cells."))
return(NULL)
},
observations = function() {
x <- dataset[, .BANOVAdependentName]
nObs <- length(options$modelTerms) + 1
if (length(na.omit(x)) <= nObs)
return(gettextf("Number of observations is < %s in repeated measures cells", nObs))
return(NULL)
},
variance = function() {
x <- dataset[, .BANOVAdependentName]
validValues <- x[is.finite(x)]
variance <- 0
if (length(validValues) > 1)
variance <- stats::var(validValues)
if (variance == 0)
return(gettext("The variance in the repeated measures cells is 0."))
return(NULL)
},
duplicateColumns = function() {
cols <- c(.BANOVAsubjectName, target)
if (length(cols) == 1L)
return()
datasetList <- as.list(dataset[, cols])
duplicatedCols <- duplicated(datasetList) | duplicated(datasetList, fromLast = TRUE)
if (any(duplicatedCols)) {
if (duplicatedCols[1L]) {
msg <- gettextf("Duplicate variables encountered in repeated measures cells, %s",
paste(target[duplicatedCols[-1L]], collapse = ", "))
} else {
msg <- gettextf("Duplicate variables encountered in %s",
paste(target[duplicatedCols[-1L]], collapse = ", "))
}
return(msg)
}
return(NULL)
},
redundantParticipantColumn = function() {
if (nlevels(dataset[[.BANOVAsubjectName]]) == nlevels(dataset[[target]])) {
return(gettextf("Duplicate participant column added to model terms, %s", target))
} else {
return()
}
}
)
.hasErrors(
dataset = dataset,
target = target,
custom = customChecks,
exitAnalysisIfErrors = TRUE
)
}
}
return(list(noVariables = noVariables, hasIV = hasIV, hasDV = hasDV))
}
.BANOVAmissingInteractionCells <- function(dataset, target) {
# sorts a data frame by all its columns
sortDataFrame <- function(df) {
df[do.call(order, df), ]
}
uniqueObservedCombinations <- sortDataFrame(unique(dataset[, target, drop = FALSE]))
uniqueExpectedCombinations <- sortDataFrame(do.call(expand.grid, c(lapply(uniqueObservedCombinations, unique), stringsAsFactors = FALSE)))
# if the rows match then it must be
if (nrow(uniqueObservedCombinations) == nrow(uniqueExpectedCombinations))
return(NULL)
# find missing combinations
observedStrings <- unname(apply(uniqueObservedCombinations, 1L, paste, collapse = jaspBase::interactionSymbol))
expectedStrings <- unname(apply(uniqueExpectedCombinations, 1L, paste, collapse = jaspBase::interactionSymbol))
missingStrings <- setdiff(expectedStrings, observedStrings)
missingStringsCollapsed <- paste(missingStrings[1:min(length(missingStrings), 10)], collapse = ", ")
interactionVariable <- paste(colnames(uniqueObservedCombinations), collapse = jaspBase::interactionSymbol)
if (length(missingStrings) <= 10) {
return(gettextf(
"The interaction effect of %1$s has no observations for the combination(s): %2$s. Please adjust the model and remove any interactions with missing cells",
interactionVariable, missingStringsCollapsed
))
} else {
return(gettextf(
"The interaction effect of %1$s has no observations for the combination(s): %2$s. Only the first 10 missing combinations are shown. Please adjust the model and remove any interactions with missing cells",
interactionVariable, missingStringsCollapsed
))
}
}
# model comparison ----
.BANOVAestimateModels <- function(jaspResults, dataset, options, errors,
analysisType = c("ANOVA", "ANCOVA", "RM-ANOVA")) {
# also makes the model comparison table
stateObj <- jaspResults[["tableModelComparisonState"]]$object
if (!is.null(jaspResults[["tableModelComparison"]]) && !is.null(stateObj)) {
stateObj$completelyReused <- TRUE # means that posteriors won't need to be resampled
return(stateObj)
} else if (errors$noVariables) {
modelTable <- .BANOVAinitModelComparisonTable(options)
jaspResults[["tableModelComparison"]] <- modelTable
return(list(analysisType = analysisType))
} else if (!is.null(stateObj) &&
.BANOVAmodelTermsUnchanged(stateObj, options) &&
.BANOVAmarginalityOptionsUnchanged(stateObj, options) &&
(.BANOVAmodelBFTypeOrOrderChanged(stateObj, options) || identical(stateObj[["shownTableSize"]], options[c("modelsShown", "numModelsShown")]))) {
# if the statement above is TRUE then no new variables were added (or seed changed)
# and the only change is in the Bayes factor type, the ordering, or the number of models shown
modelTable <- .BANOVAinitModelComparisonTable(options)
modelTable[["Models"]] <- .BANOVAgetModelTitlesWithAllTerms(stateObj[["models"]], stateObj[["model.list"]], analysisType, options[["hideNuisanceParameters"]])
internalTable <- stateObj$internalTableObj$internalTable
nmodels <- nrow(internalTable)
.BANOVAestimateModels
shownTableSize <- .BANOVAgetShownTableSizeAndAddFootnote(modelTable, nmodels, options[["modelsShown"]], options[["numModelsShown"]])
if (.BANOVAmodelPriorOptionsChanged(stateObj, options)) {
priorProbs <- .BANOVAcomputePriorModelProbs(stateObj$model.list, stateObj$nuisance, options)
internalTable[["P(M)"]] <- priorProbs
stateObj$completelyReused <- FALSE # model-averaged posteriors need to be resampled
} else {
stateObj$completelyReused <- TRUE # model-averaged posteriors do not need to be resampled
}
internalTableObj <- .BANOVAfinalizeInternalTable(options, internalTable)
stateObj$postProbs <- internalTableObj$internalTable[["P(M|data)"]]
stateObj$priorProbs <- internalTableObj$internalTable[["P(M)"]]
stateObj$internalTableObj <- internalTableObj
modelTable$setData(internalTableObj$table[seq_len(shownTableSize), ])
jaspResults[["tableModelComparison"]] <- modelTable
return(stateObj)
}
modelTerms <- options$modelTerms
dependent <- options$dependent
if (analysisType == "RM-ANOVA") {
modelTerms[[length(modelTerms) + 1L]] <- list(components = .BANOVAsubjectName, isNuisance = TRUE)
dependent <- .BANOVAdependentName
}
fixedFactors <- options$fixedFactors
randomFactors <- .BANOVAgetRandomFactors(options, analysisType)
tempRScale <- .BANOVAgetRScale(options, analysisType)
rscaleFixed <- tempRScale[["rscaleFixed"]]
rscaleRandom <- tempRScale[["rscaleRandom"]]
rscaleCont <- tempRScale[["rscaleCont"]]
rscaleEffects <- tempRScale[["rscaleEffects"]]
tmp <- .BANOVAcreateModelFormula(dependent, modelTerms)
model.formula <- tmp$model.formula
nuisance <- tmp$nuisance
effects <- tmp$effects
if (analysisType == "RM-ANOVA") {
legacy <- options[["legacyResults"]]
rmFactors <- vapply(options[["repeatedMeasuresFactors"]], `[[`, character(1L), "name")
} else {
legacy <- FALSE
rmFactors <- NULL
}
#Make a list of models to compare
temp <- .BANOVAgenerateAllModelFormulas(
formula = model.formula,
nuisance = nuisance,
analysisType = analysisType,
enforcePrincipleOfMarginalityFixedEffects = options[["enforcePrincipleOfMarginalityFixedEffects"]],
enforcePrincipleOfMarginalityRandomSlopes = options[["enforcePrincipleOfMarginalityRandomSlopes"]],
rmFactors = rmFactors,
legacy = legacy
)
model.list <- temp$modelList
nuisance <- temp$nuisance
nuisanceRandomSlopes <- temp$nuisanceRandomSlopes
if (length(model.list) == 1L) {
modelTable <- .BANOVAinitModelComparisonTable(options)
modelTable$setError(gettext("Bayes factor is undefined -- all effects are specified as nuisance."))
jaspResults[["tableModelComparison"]] <- modelTable
return(list(analysisType = analysisType))
}
# TODO: discuss if we want a hard limit if length(model.list) > ...
# that would avoid peoples computers from slowing down tremendously.
# JCG: A hard limit would be a bit cruel right? How about a warning instead? (in QML already?)
#Run all other models and store the results in the model.object
neffects <- length(effects)
nmodels <- length(model.list)
modelObject <- vector("list", nmodels)
reorderedEffects <- .BANOVAreorderTerms(effects)
reorderedNuisance <- if (is.null(nuisance)) NULL else .BANOVAreorderTerms(nuisance)
if (nmodels > 0L && neffects > 0L) {
# effects.matrix contains for each row (model), which predictors (columns) are in the model
# interactions.matrix contains for each predictor (column) what predictors (rows) is it composed of
interactions.matrix <- .BANOVAcomputeInteractionsMatrix(effects)
effects.matrix <- matrix(data = FALSE, nrow = nmodels, ncol = neffects)
colnames(effects.matrix) <- effects
for (m in seq_len(nmodels)) {
modelObject[[m]] <- list()
if (m == 1L) {
if (is.null(nuisance)) { # intercept only
modelObject[[m]]$title <- gettext("Null model")
next # all effects are FALSE anyway
} else {
if (analysisType == "RM-ANOVA" && !legacy) {
tempNuisance <- setdiff(nuisance, nuisanceRandomSlopes)
modelObject[[m]]$title <- sprintf(ngettext(
length(tempNuisance),
"Null model (incl. %s and random slopes)",
"Null model (incl. %s, and random slopes)"
), paste(.BANOVAdecodeNuisance(tempNuisance), collapse = ", "))
} else {
modelObject[[m]]$title <- gettextf("Null model (incl. %s)", paste(.BANOVAdecodeNuisance(nuisance), collapse = ", "))
}
}
}
model.effects <- .BANOVAgetFormulaComponents(model.list[[m]])
reorderedModelEffects <- .BANOVAreorderTerms(model.effects)
idx <- match(reorderedModelEffects, reorderedEffects, nomatch = 0L)
idx <- idx[!is.na(idx)]
effects.matrix[m, idx] <- TRUE
if (m > 1L) {
model.title <- model.effects[match(reorderedModelEffects, reorderedNuisance, 0L) == 0L]
modelObject[[m]]$title <- jaspBase::gsubInteractionSymbol(paste(model.title, collapse = " + "))
}
}
rownames(effects.matrix) <- sapply(modelObject, `[[`, "title")
}
modelTable <- .BANOVAinitModelComparisonTable(options)
modelNames <- .BANOVAgetModelTitlesWithAllTerms(modelObject, model.list, analysisType, options[["hideNuisanceParameters"]])
shownTableSize <- .BANOVAgetShownTableSizeAndAddFootnote(modelTable, nmodels, options[["modelsShown"]], options[["numModelsShown"]])
modelTable[["Models"]] <- modelNames[seq_len(shownTableSize)]
jaspResults[["tableModelComparison"]] <- modelTable
# internalTable is an internal representation of the model comparison table
# internalTable <- matrix(NA, nmodels, 6L, dimnames = list(modelNames, c("Models", "P(M)", "P(M|data)", "BFM", "BF10", "error %")))
internalTable <- data.frame(
Models = character(nmodels),
`P(M)` = numeric(nmodels),
`P(M|data)` = numeric(nmodels),
BFM = numeric(nmodels),
BF10 = numeric(nmodels),
`error %` = numeric(nmodels),
check.names = FALSE
)
# set model names, BF null model and p(M)
internalTable[["Models"]] <- modelNames
internalTable[1L, "BF10"] <- 0
internalTable[["P(M)"]] <- .BANOVAcomputePriorModelProbs(model.list, nuisance, options)
#Now compute Bayes Factors for each model in the list, and populate the tables accordingly
startProgressbar(nmodels, gettext("Computing Bayes factors"))
# check if any models can be resued from the state
if (!is.null(stateObj[["modelTerms"]])) {
oldFormulas <- sapply(stateObj$model.list, .BANOVAreorderFormulas)
newFormulas <- sapply(model.list, .BANOVAreorderFormulas)
reuseable <- match(newFormulas, oldFormulas)
} else {
reuseable <- rep(NA, nmodels)
}
if (analysisType == "RM-ANOVA") {
# the default null-model contains subject
anyNuisance <- TRUE
} else {
# the default null-model is intercept only
anyNuisance <- length(nuisance) > 0L
}
# without these there is no error
bfIterations <- if (options[["samplingMethodNumericAccuracy"]] == "auto") 1e4L else options[["samplesNumericAccuracy"]]
if (options[["integrationMethod"]] == "automatic") {
bfIntegrationMethod <- "auto"
} else {
bfIntegrationMethod <- "laplace"
modelTable$addFootnote(gettext("The Laplace approximation is faster but also less accurate than the automatic method. This is fine when exploring the data and obtaining a rough estimate. However, we recommend using the automatic integration method to obtain the final results."))
modelTable$addFootnote(gettext("The Laplace approximation does not always provide an error estimate."), colNames = "error %")
}
for (m in seq_len(nmodels)) {
# loop over all models, where the first is the null-model and the last the most complex model
if (is.na(reuseable[m])) {
if (!is.null(model.list[[m]])) {
.setSeedJASP(options)
bf <- try(BayesFactor::lmBF(
formula = model.list[[m]],
data = dataset,
whichRandom = randomFactors,
progress = FALSE,
posterior = FALSE,
rscaleFixed = rscaleFixed,
rscaleRandom = rscaleRandom,
rscaleCont = rscaleCont,
rscaleEffects = rscaleEffects,
iterations = bfIterations,
method = bfIntegrationMethod
))
if (isTryError(bf)) {
modelName <- strsplit(.BANOVAas.character.formula(model.list[[m]]), "~ ")[[1L]][2L]
.quitAnalysis(gettextf("Bayes factor could not be computed for model: %1$s.\nThe error message was: %2$s.",
.BANOVAdecodeSubject(modelName), .extractErrorMessage(bf)))
} else {
# delete the data object -- otherwise it gets saved in the state
bf@data <- data.frame()
}
} else {
bf <- NULL
}
} else {
bf <- stateObj$models[[reuseable[m]]]$bf
}
modelObject[[m]]$bf <- bf
if (!is.null(bf)) {
# if (isTryError(bf)) {
# message <- .extractErrorMessage(bf)
# modelObject[[m]]$error.message <- paste("Bayes factor computation failed with error:", message)
# } else {
# bfObj can have a modified denominator, but the saved objects are always compared against intercept only
bfObj <- modelObject[[m]]$bf
if (anyNuisance)
bfObj <- bfObj / modelObject[[1L]]$bf
internalTable[m, "BF10"] <- bfObj@bayesFactor[, "bf"] # always LogBF10!
internalTable[m, "error %"] <- bfObj@bayesFactor[, "error"] # always NA if bfIntegrationMethod == "laplace"
# }
# disable filling of Bayes factor column (see https://github.com/jasp-stats/jasp-test-release/issues/1018 for discussion)
# modelTable[["BF10"]] <- .recodeBFtype(internalTable[["BF10"]],
# newBFtype = options$bayesFactorType,
# oldBFtype = "LogBF10")
modelTable[["error %"]] <- internalTable[seq_len(shownTableSize), "error %"]
}
progressbarTick()
}
internalTableObj <- .BANOVAfinalizeInternalTable(options, internalTable)
modelTable$setData(internalTableObj$table[seq_len(shownTableSize), ])
if (length(internalTableObj[["footnotes"]]) > 0L) {
idxRow <- internalTableObj[["footnotes"]][["rows"]]
idxCol <- internalTableObj[["footnotes"]][["cols"]]
tb <- internalTableObj$table
modelTable$addFootnote(message = internalTableObj[["footnotes"]][["message"]],
symbol = " <sup>✝</sup>", # latin cross
rowNames = rownames(tb)[idxRow],
colNames = colnames(tb)[idxCol])
}
if (anyNuisance) {
if (analysisType == "RM-ANOVA" && !options[["legacyResults"]]) {
message <- gettextf("All models include %s, and random slopes for all repeated measures factors.",
paste0(.BANOVAdecodeNuisance(setdiff(nuisance, nuisanceRandomSlopes)), collapse = ", "))
} else {
if (analysisType == "RM-ANOVA")
modelTable$addFootnote(message = gettext("Random slopes for repeated measures factors are excluded."), symbol = .BANOVAGetWarningSymbol())
message <- gettextf("All models include %s.", paste0(.BANOVAdecodeNuisance(nuisance), collapse = ", "))
}
modelTable$addFootnote(message = message)
}
model <- list(
models = modelObject,
postProbs = internalTableObj$internalTable[["P(M|data)"]],
priorProbs = internalTableObj$internalTable[["P(M)"]],
internalTableObj = internalTableObj,
effects = effects.matrix,
interactions.matrix = interactions.matrix,
nuisance = nuisance,
nuisanceRandomSlopes = nuisanceRandomSlopes,
analysisType = analysisType,
# these are necessary for partial reusage of the state (e.g., when a fixedFactor is added/ removed)
model.list = model.list,
fixedFactors = fixedFactors,
randomFactors = if (analysisType == "RM-ANOVA") randomFactors else options[["randomFactors"]], # stored because they are modified in RM-ANOVA
modelTerms = modelTerms,
covariates = if (analysisType == "ANOVA") NULL else options[["covariates"]],
seed = options[["seed"]],
setSeed = options[["setSeed"]],
reuseable = reuseable,
RMFactors = options[["repeatedMeasuresFactors"]],
modelPriorOptions = options[.BANOVAmodelSpaceDependencies(options[["modelPrior"]])],
hideNuisanceEffects = options[["hideNuisanceEffects"]],
shownTableSize = options[c("modelsShown", "numModelsShown")]
)
# save state
stateObj <- createJaspState(object = model, dependencies = c(
# does NOT depend on any factors or covariates, to facilitate reusing previous models
"dependent", "repeatedMeasuresCells", "samplingMethodNumericAccuracy", "samplesNumericAccuracy", "seed", "setSeed",
.BANOVArscaleDependencies(options[["priorSpecificationMode"]])
))
jaspResults[["tableModelComparisonState"]] <- stateObj
return(model)
}
.BANOVAeffectsTable <- function(jaspResults, options, model) {
if (!is.null(jaspResults[["tableEffects"]]) || !options[["effects"]])
return()
# isTRUE should handle a state issue, see https://github.com/jasp-stats/jasp-test-release/issues/839
if (isTRUE(model[["analysisType"]] != "RM-ANOVA" && options[["dependent"]] != "")) {
title <- gettextf("Analysis of Effects - %s", options[["dependent"]])
} else {
title <- gettext("Analysis of Effects")
}
effectsTable <- createJaspTable(title = title)
effectsTable$position <- 2
effectsTable$dependOn(c(
.BANOVAdataDependencies(),
"effects", "effectsType",
"samplingMethodNumericAccuracy", "samplesNumericAccuracy", "bayesFactorType",
"integrationMethod",
.BANOVAmodelSpaceDependencies(options[["modelPrior"]]),
.BANOVArscaleDependencies(options[["priorSpecificationMode"]])
))
effectsTable$addCitation(.BANOVAcitations[1:2])
inclusion.title <- switch(
options$bayesFactorType,
"LogBF10" = gettext("Log(BF<sub>incl</sub>)"),
"BF01" = gettext("BF<sub>excl</sub>"),
"BF10" = gettext("BF<sub>incl</sub>")
)
effectsTable$addColumnInfo(name = "Effects", title = gettext("Effects"), type = "string")
effectsTable$addColumnInfo(name = "P(incl)", title = gettext("P(incl)"), type = "number")
effectsTable$addColumnInfo(name = "P(excl)", title = gettext("P(excl)"), type = "number")
effectsTable$addColumnInfo(name = "P(incl|data)", title = gettext("P(incl|data)"), type = "number")
effectsTable$addColumnInfo(name = "P(excl|data)", title = gettext("P(excl|data)"), type = "number")
effectsTable$addColumnInfo(name = "BFInclusion", title = inclusion.title, type = "number")
if (options$effectsType == "matchedModels") {
effectsTable$addFootnote(gettextf("Compares models that contain the effect to equivalent models stripped of the effect. Higher-order interactions are excluded. Analysis suggested by Sebastiaan Math%st.", "\u00F4"))
}
if (is.null(model$models)) {
jaspResults[["tableEffects"]] <- effectsTable
return()
}
effects.matrix <- model$effects
no.effects <- ncol(effects.matrix)
effectNames <- colnames(model$effects)
idxNotNan <- c(TRUE, !is.nan(model$postProbs[-1L]))
if (options$effectsType == "allModels") {
if (any(!idxNotNan)) {
# renormalize the prior and posterior probabilities
priorInclProb <- crossprod(effects.matrix[idxNotNan, , drop = FALSE], model[["priorProbs"]][idxNotNan] / sum(model[["priorProbs"]][idxNotNan]))
postInclProb <- crossprod(effects.matrix[idxNotNan, , drop = FALSE], model[["postProbs"]][idxNotNan] / sum(model[["postProbs"]][idxNotNan]))
} else {
# do not renormalize
priorInclProb <- crossprod(effects.matrix, model[["priorProbs"]])
postInclProb <- crossprod(effects.matrix, model[["postProbs"]])
}
# deal with numerical error
postInclProb[postInclProb > 1] <- 1
postInclProb[postInclProb < 0] <- 0
bfIncl <- (postInclProb / (1 - postInclProb)) / (priorInclProb / (1 - priorInclProb))
priorExclProb <- 1 - priorInclProb
postExclProb <- 1 - postInclProb
} else {
tmp <- BANOVAcomputMatchedInclusion(
effectNames, effects.matrix, model$interactions.matrix,
model$priorProbs, model$postProbs
)
priorInclProb <- tmp[["priorInclProb"]]
postInclProb <- tmp[["postInclProb"]]
bfIncl <- tmp[["bfIncl"]]
priorExclProb <- tmp[["priorExclProb"]]
postExclProb <- tmp[["postExclProb"]]
# show BFinclusion for nuisance predictors as 1, rather than NaN
priorInclIs1 <- is.nan(bfIncl) & abs(1 - priorInclProb) <= sqrt(.Machine$double.eps)
bfIncl[priorInclIs1] <- 1
}
if (sum(!idxNotNan) > 1L) { # null model is always omitted, so 2 or more omitted indicates some models failed
effectsTable$addFootnote(message = gettext("Some Bayes factors could not be calculated. Inclusion probabilities and Bayes factors are calculated while excluding these models. The results may be uninterpretable!"),
symbol = .BANOVAGetWarningSymbol()
)
}
effectsTable[["Effects"]] <- jaspBase::gsubInteractionSymbol(effectNames)
effectsTable[["P(incl)"]] <- priorInclProb
effectsTable[["P(incl|data)"]] <- postInclProb
effectsTable[["BFInclusion"]] <- switch(
options$bayesFactorType,
"LogBF10" = log(bfIncl),
"BF01" = 1 / bfIncl,
"BF10" = bfIncl
)
effectsTable[["P(excl)"]] <- priorExclProb
effectsTable[["P(excl|data)"]] <- postExclProb
jaspResults[["tableEffects"]] <- effectsTable
return()
}
BANOVAcomputMatchedInclusion <- function(effectNames, effects.matrix, interactions.matrix,
priorProbs, postProbs) {
# this method is inspired by this post: https://www.cogsci.nl/blog/interpreting-bayesian-repeated-measures-in-jasp
priorInclProb <- postInclProb <- bfIncl <- priorExclProb <- postExclProb <-
numeric(length(effectNames))
for (i in seq_along(effectNames)) {
effect <- effectNames[i]
# get all higher order interactions of which this effect is a component
# e.g., V1 is a component of V1:V2
idx1 <- interactions.matrix[effect, ]
# get all models that exclude the predictor, but that always include the lower order main effects
# e.g., V1:V2 is compared against models that always include V1 and V2
idx2 <- interactions.matrix[, effect] # if effect is V1:V2, idx2 contains c(V1, V2)
# idx3 is FALSE if a model contains higher order interactions of effect, TRUE otherwise
idx3 <- !matrixStats::rowAnys(effects.matrix[, idx1, drop = FALSE])
# all models that include the effect, without higher order interactions
idx4 <- idx3 & effects.matrix[, i]
priorInclProb[i] <- sum(idx4 * priorProbs)
postInclProb[i] <- sum(idx4 * postProbs)
# the models to consider for the prior/ posterior exclusion probability.
# idx5 includes models that have: all subcomponents & no higher order interaction & not the effect
idx5 <- matrixStats::rowAlls(effects.matrix[, idx2, drop = FALSE]) & idx3 & !effects.matrix[, i]
priorExclProb[i] <- sum(idx5 * priorProbs)
postExclProb[i] <- sum(idx5 * postProbs)
# compute inclusion BF
bfIncl[i] <- (postInclProb[i] / postExclProb[i]) / (priorInclProb[i] / priorExclProb[i])
}
return(list(
priorInclProb = priorInclProb,
postInclProb = postInclProb,
priorExclProb = priorExclProb,
postExclProb = postExclProb,
bfIncl = bfIncl
))
}
.BANOVAinitModelComparisonTable <- function(options) {
# function that creates an empty JASP table to be filled later
modelTable <- createJaspTable(title = gettext("Model Comparison"))
modelTable$position <- 1L
modelTable$addCitation(.BANOVAcitations[1:2])
modelTable$dependOn(c(
.BANOVAdataDependencies(),
"samplingMethodNumericAccuracy", "samplesNumericAccuracy", "integrationMethod",
"bayesFactorType", "bayesFactorOrder",
"hideNuisanceParameters", "legacyResults",
"modelsShown",
.BANOVAmodelSpaceDependencies(options[["modelPrior"]]),
.BANOVArscaleDependencies(options[["priorSpecificationMode"]])
))
if (options[["modelsShown"]] == "limited")
modelTable$dependOn("numModelsShown")
switch(options$bayesFactorType,
BF10 = {
bfm.title <- gettext("BF<sub>M</sub>")
bf.title <- gettext("BF<sub>10</sub>")
},
BF01 = {
bfm.title <- gettext("BF<sub>M</sub>")
bf.title <- gettext("BF<sub>01</sub>")
},
LogBF10 = {
bfm.title <- gettext("Log(BF<sub>M</sub>)")
bf.title <- gettext("Log(BF<sub>10</sub>)")
}
)
modelTable$addColumnInfo(name = "Models", title = gettext("Models"), type = "string")
modelTable$addColumnInfo(name = "P(M)", title = gettext("P(M)"), type = "number")
modelTable$addColumnInfo(name = "P(M|data)", title = gettext("P(M|data)"), type = "number")
modelTable$addColumnInfo(name = "BFM", title = bfm.title, type = "number")
modelTable$addColumnInfo(name = "BF10", title = bf.title, type = "number")
modelTable$addColumnInfo(name = "error %", title = gettextf("error %%"), type = "number")
return(modelTable)
}
.BANOVAfinalizeInternalTable <- function(options, internalTable) {
# function that actually fills in the table created by .BANOVAinitModelComparisonTable
footnotes <- NULL
logSumExp <- matrixStats::logSumExp
logbfs <- internalTable[["BF10"]]
logprior <- log(internalTable[["P(M)"]])
if (!anyNA(internalTable[["BF10"]])) {
# no errors, proceed normally and complete the table
logsumbfs <- logSumExp(logbfs + logprior)
internalTable[["P(M|data)"]] <- exp(logbfs + logprior - logsumbfs)
nmodels <- nrow(internalTable)
for (i in seq_len(nmodels)) {
internalTable[i, "BFM"] <- logbfs[i] - logSumExp(logbfs[-i]) + log(nmodels - 1L)
}
} else {
# create table excluding failed models
idxGood <- !is.na(logbfs)
widxGood <- which(idxGood)
logsumbfs <- logSumExp(logbfs[idxGood])
internalTable[["P(M|data)"]] <- exp(logbfs - logsumbfs)
nmodels <- sum(idxGood)
widxBad <- which(!idxGood)
for (i in widxGood) {
internalTable[["BFM"]][i] <- logbfs[i] - logSumExp(logbfs[-c(i, widxBad)]) + log(nmodels - 1L)
}
internalTable[widxBad, "P(M|data)"] <- NaN
internalTable[widxBad, "BFM"] <- NaN
internalTable[widxBad, "BF10"] <- NaN
footnotes <- list(
message = gettext("<b><em>Warning.</em></b> Some Bayes factors could not be calculated. Multi model inference is carried out while excluding these models. The results may be uninterpretable!")
)
}
# create the output table
table <- as.data.frame(internalTable)
if (options[["bayesFactorType"]] == "LogBF10") {
table[["BFM"]] <- internalTable[["BFM"]]
} else {
table[["BFM"]] <- exp(internalTable[["BFM"]])
}
o <- order(table[["BF10"]], decreasing = TRUE)
table <- table[o, ]
idxNull <- which(o == 1L)
if (options[["bayesFactorOrder"]] == "nullModelTop") {
table[idxNull, "error %"] <- NA
table <- table[c(idxNull, seq_len(nrow(table))[-idxNull]), ]
} else {
table[["BF10"]] <- table[["BF10"]] - table[1L, "BF10"]
# recompute error (see BayesFactor:::`.__T__/:base`$`BFBayesFactor#BFBayesFactor`)
table[idxNull, "error %"] <- 0
table[["error %"]] <- sqrt(table[["error %"]]^2 + table[["error %"]][1L]^2)
table[1L, "error %"] <- NA
}
table[["BF10"]] <- .recodeBFtype(table[["BF10"]], newBFtype = options[["bayesFactorType"]], oldBFtype = "LogBF10")
table[["error %"]] <- 100 * table[["error %"]]
if (!is.null(footnotes)) {
idxNan <- which(do.call(cbind, lapply(table[-ncol(table)], is.nan)), arr.ind = TRUE)
footnotes[["rows"]] <- idxNan[, "row"]
footnotes[["cols"]] <- idxNan[, "col"]
}
return(list(table = table, internalTable = internalTable, footnotes = footnotes))
}
# posterior inference ----
.BANOVAestimatePosteriors <- function(jaspResults, dataset, options, model) {
userNeedsPosteriorSamples <- options$posteriorEstimates || options$modelAveragedPosteriorPlot || options$qqPlot ||
options$rsqPlot || options$criTable
if (is.null(model$models) || !userNeedsPosteriorSamples)
return()
# model[["completelyReused"]] is needed because it can happen that some posterior samples can be reused (e.g.,
# when the modelTerms change)
posteriors <- jaspResults[["statePosteriors"]]$object
if (!is.null(posteriors) && isTRUE(model[["completelyReused"]])) { # can all posterior be reused?
posteriorsCRI <- jaspResults[["statePosteriorsCRI"]]$object
if (is.null(posteriorsCRI)) {# do we need to recompute credible intervals?
posteriorsCRI <- .BANOVAcomputePosteriorCRI(dataset, options, model, posteriors)
statePosteriorsCRI <- createJaspState(object = posteriorsCRI)
statePosteriorsCRI$dependOn(options = "credibleInterval",
optionsFromObject = jaspResults[["tableModelComparisonState"]])
jaspResults[["statePosteriorCRI"]] <- statePosteriorsCRI
}
} else { # compute posteriors and credible intervals
posteriors <- .BANOVAsamplePosteriors(jaspResults, dataset, options, model, posteriors)
posteriorsCRI <- .BANOVAcomputePosteriorCRI(dataset, options, model, posteriors)
statePosteriors <- createJaspState(object = posteriors)
statePosteriors$dependOn(
c(.BANOVAmodelSpaceDependencies(options[["modelPrior"]]), "samplingMethodMCMC", "samplesMCMC"),
optionsFromObject = jaspResults[["tableModelComparisonState"]]
)
jaspResults[["statePosteriors"]] <- statePosteriors
statePosteriorsCRI <- createJaspState(object = posteriorsCRI)
statePosteriorsCRI$dependOn(options = "credibleInterval",
optionsFromObject = jaspResults[["statePosteriors"]])
jaspResults[["statePosteriorCRI"]] <- statePosteriorsCRI
}
.BANOVASamplingIssuesTable(jaspResults, posteriors[["samplingIssues"]])
return(c(posteriors, posteriorsCRI))
}
.BANOVAestimatesTable <- function(jaspResults, options, model) {
if (!is.null(jaspResults[["tablePosteriorEstimates"]]) || !options[["posteriorEstimates"]])
return()
estsTable <- createJaspTable(title = gettext("Model Averaged Posterior Summary"))
estsTable$position <- 3
estsTable$dependOn(c(
.BANOVAdataDependencies(),
"posteriorEstimates",
"samplingMethodMCMC", "samplesMCMC", "credibleInterval",
.BANOVAmodelSpaceDependencies(options[["modelPrior"]]),
.BANOVArscaleDependencies(options[["priorSpecificationMode"]])
))
overTitle <- gettextf("%s%% Credible Interval", format(100 * options[["credibleInterval"]], digits = 3))
estsTable$addColumnInfo(name = "Variable", type = "string")
estsTable$addColumnInfo(name = "Level", type = "string")
estsTable$addColumnInfo(name = "Mean", type = "number")
estsTable$addColumnInfo(name = "SD", type = "number")
estsTable$addColumnInfo(name = "Lower", type = "number", overtitle = overTitle)
estsTable$addColumnInfo(name = "Upper", type = "number", overtitle = overTitle)
if (!is.null(model[["posteriors"]])) {
.BANOVAfillEstimatesTable(
jaspTable = estsTable,
mus = model$posteriors$weightedMeans,
sds = model$posteriors$weightedSds,
cri = model$posteriors$weightedCRIs,
hasNoLevels = model$posteriors$allContinuous,
isRandom = model$posteriors$isRandom
)
}
jaspResults[["tablePosteriorEstimates"]] <- estsTable
return()
}
.BANOVAfillEstimatesTable <- function(jaspTable, mus, sds, cri, hasNoLevels, isRandom = NULL) {
if (!is.null(isRandom) && any(isRandom)) {
# remove random effects
mus <- mus[!isRandom]
sds <- sds[!isRandom]
cri <- cri[!isRandom, ]
}
table <- .BANOVAgetLevelsFromParamNames(names(mus))
colnames(table) <- c("Variable", "Level")
# set repeated parameter names to ""
idxDup <- duplicated(table[, "Variable"])
table[idxDup, "Variable"] <- ""
# decode base64 variables
table[!idxDup, "Variable"] <- jaspBase::gsubInteractionSymbol(table[!idxDup, "Variable"])
# rename mu to Intercept
table[1L, "Variable"] <- "Intercept"
# attach posterior means and sds
table <- cbind(as.data.frame(table),
"Mean" = mus,
"SD" = sds,
"Lower" = cri[, 1L],
"Upper" = cri[, 2L])
if (hasNoLevels)
table <- table[, -2L]
jaspTable$setData(table)
return()
}
.BANOVArsquaredTable <- function(jaspResults, options, model) {
if (!is.null(jaspResults[["tableBMACRI"]]) || !options[["criTable"]])
return()
criTable <- createJaspTable(title = gettextf("Model Averaged R%s", "\u00B2"))
criTable$position <- 3.5
criTable$dependOn(c(
.BANOVAdataDependencies(),
"criTable",
"samplingMethodMCMC", "samplesMCMC", "bayesFactorType", "credibleInterval",
.BANOVAmodelSpaceDependencies(options[["modelPrior"]]),
.BANOVArscaleDependencies(options[["priorSpecificationMode"]])
))
overTitle <- gettextf("%s%% Credible Interval", format(100 * options[["credibleInterval"]], digits = 3))
criTable$addColumnInfo(name = "rsq", type = "string", title = "")
criTable$addColumnInfo(name = "Mean", type = "number")
criTable$addColumnInfo(name = "Lower", type = "number", overtitle = overTitle)
criTable$addColumnInfo(name = "Upper", type = "number", overtitle = overTitle)
if (!is.null(model[["posteriors"]])) {
cri <- model[["posteriors"]][["weightedRsqCri"]]
df <- data.frame(
rsq = "R\u00B2",
Mean = model[["posteriors"]][["weightedRsqMean"]],
Lower = cri[1L],
Upper = cri[2L]
)
criTable$setData(df)
} else {
criTable[["rsq"]] <- "R\u00B2"
}
jaspResults[["tableBMACRI"]] <- criTable
return()
}
.BANOVASamplingIssuesTable <- function(jaspResults, samplingIssues) {
if (is.null(samplingIssues) || length(samplingIssues) == 0L || !is.null(jaspResults[["tableSamplingIssues"]]))
return()
issuesTable <- createJaspTable(title = gettext("<b><em>Warning: sampling issues encountered!</em></b>"),
# below model and effects table (which may be unaffected) but above estimates table
position = 2.5)
issuesTable$addColumnInfo(name = "Models", title = gettext("Affected model"), type = "string")
issuesTable$addColumnInfo(name = "percentageSucces", title = gettextf("%% useful samples"), type = "number")
issuesTable$addColumnInfo(name = "noSuccess", title = gettext("Useful samples"), type = "integer")
issuesTable$addColumnInfo(name = "total", title = gettext("Total samples drawn"), type = "integer")
df <- data.frame(
Models = vapply(samplingIssues, `[[`, FUN.VALUE = character(1L), "model"),
noSuccess = vapply(samplingIssues, `[[`, FUN.VALUE = integer(1L), "remainingRows"),
total = vapply(samplingIssues, `[[`, FUN.VALUE = integer(1L), "originalRows")
)
df[["percentageSucces"]] <- df[["noSuccess"]] / df[["total"]]
if (any(df[["percentageSucces"]] < 0.25) || any(df[["remainingRows"]] < 1000L))
issuesTable$addFootnote(
gettextf("For some affected models, more than 75%% of the posterior samples failed, or fewer than 1000 samples remained for subsequent results. All model-averaged output may be biased and uninterpretable. Check the model specification and data for any odd patterns."),
symbol = .BANOVAGetWarningSymbol()
)
issuesTable$setData(df)
issuesTable$dependOn(
# these options correspond to userNeedsPosteriorSamples inside .BANOVAestimatePosteriors
options = c("posteriorEstimates", "modelAveragedPosteriorPlot", "qqPlot", "rsqPlot", "criTable", "modelTerms"),
optionsFromObject = jaspResults[["statePosteriors"]]
)
jaspResults[["tableSamplingIssues"]] <- issuesTable
}
.BANOVAGetWarningSymbol <- function() {
gettext("<b><em>Warning.</em></b>")
}
# Plots wrappers ----
.BANOVAposteriorPlot <- function(jaspResults, dataset, options, model) {
# meta wrapper for model averaged posterior plots
if (!is.null(jaspResults[["posteriorPlot"]]) || !options$modelAveragedPosteriorPlot)
return()
posteriorPlotContainer <- createJaspContainer(title = gettext("Model Averaged Posterior Distributions"))
jaspResults[["posteriorPlot"]] <- posteriorPlotContainer
posteriorPlotContainer$dependOn(c("modelAveragedPosteriorPlot", "modelTerms", "credibleInterval", "groupPosterior",
"repeatedMeasuresCells", "dependent"))
posteriorPlotContainer$position <- 4
if (is.null(model$models)) {
posteriorPlotContainer[["dummyplot"]] <- createJaspPlot(title = gettext("Posterior distribution"), width = 400, height = 400,
plot = NULL)
} else {
posteriorPlotContainer$dependOn(optionsFromObject = jaspResults[["tableModelComparisonState"]])
.BANOVAfillPosteriorPlotContainer(
container = posteriorPlotContainer,
densities = model$posterior$weightedDensities[, -1L, ], # omit intercept
cris = model$posterior$weightedCRIs[-1L, ], # omit intercept
isRandom = model$posteriors$isRandom[-1L], # omit intercept
groupParameters = identical(options[["groupPosterior"]], "grouped")
)
}
return()
}
.BANOVAfillPosteriorPlotContainer <- function(container, densities, cris, isRandom = NULL, groupParameters = FALSE) {
allParamNames <- colnames(densities)
tmp <- .BANOVAgetLevelsFromParamNames(allParamNames)
plotTitles <- jaspBase::gsubInteractionSymbol(tmp[, "parameter"])
xNames <- tmp[, "level"]
if (is.null(isRandom)) {
indices <- seq_along(allParamNames)
} else {
indices <- which(!isRandom)
}
if (groupParameters) {
indices <- split(indices, plotTitles[indices])[order(unique(plotTitles[indices]))]
nms <- names(indices)
for (i in seq_along(indices)) {
ind <- indices[[i]]
if (all(c(densities[, ind, "y"]) == 0)) { # only TRUE is never sampled and set to 0
plot <- createJaspPlot(title = nms[i], width = 400, height = 400)
plot$setError(gettext("Could not plot posterior distribution. Check if an error occurred for models including this parameter."))
} else {
# make prior posterior plot
dfLines <- data.frame(x = c(densities[, ind, "x"]),
y = c(densities[, ind, "y"]),
g = rep(xNames[ind], each = nrow(densities)))
xBreaks <- jaspGraphs::getPrettyAxisBreaks(dfLines$x)
yBreaks <- jaspGraphs::getPrettyAxisBreaks(c(0, dfLines$y))
breaksYmax <- yBreaks[length(yBreaks)]
obsYmax <- max(dfLines$y)
newymax <- max(1.25 * obsYmax, breaksYmax)
lInd <- length(ind)
if (lInd == 1L) {
plusmin <- 0
showLegend <- FALSE
} else {
# some jitter in CRI height
showLegend <- TRUE
plusmin <- rep(-1, lInd)
plusmin[seq(2, lInd, 2)] <- 0
if (lInd > 2)
plusmin[seq(3, lInd, 3)] <- 1
}
dfCri <- data.frame(
xmin = cris[ind, 1L],
xmax = cris[ind, 2L],
y = (newymax - obsYmax)/2 + obsYmax + plusmin * (newymax - obsYmax) / 10,
g = xNames[ind]
)
if (showLegend) {
# if two distributions are remarkably similar, i.e., have nearly identical credible intervals,
# we add a linetype aestethic
tol <- 1e-4
distMin <- which(stats::dist(dfCri$xmin) < tol)
distMax <- which(stats::dist(dfCri$xmin) < tol)
distBoth <- intersect(distMin, distMax)
if (length(distBoth) > 0) {
aesLine <- ggplot2::aes(x = x, y = y, g = g, color = g, linetype = g)
aesErrorbar <- ggplot2::aes(xmin = xmin, xmax = xmax, y = y, group = g, color = g, linetype = g)
} else {
aesLine <- ggplot2::aes(x = x, y = y, g = g, color = g)
aesErrorbar <- ggplot2::aes(xmin = xmin, xmax = xmax, y = y, group = g, color = g)
}
} else {
# don't use color for single density plots (covariates)
aesLine <- ggplot2::aes(x = x, y = y, g = g)
aesErrorbar <- ggplot2::aes(xmin = xmin, xmax = xmax, y = y, group = g)
}
maxheight <- min(newymax - dfCri$y[1:min(lInd, 3)])
xlab <- nms[i]
# ggplot doesn't like our fancy unicode * so we need to escape it
xlab <- stringi::stri_escape_unicode(xlab)
# change escaped version into x
xlab <- gsub(pattern = "\\u2009\\u273b\\u2009", replacement = " x ", x = xlab, fixed = TRUE)
ncolLegend <- ceiling(lInd / 14)
guideLegend <- ggplot2::guide_legend(title = gettext("Level"), keywidth = 0.25, keyheight = 0.1, default.unit = "inch",
ncol = ncolLegend)
p <- ggplot2::ggplot(data = dfLines, mapping = aesLine) +
ggplot2::geom_line(size = 1.1) +
ggplot2::geom_errorbarh(data = dfCri, mapping = aesErrorbar, height = maxheight, size = 1.1,
inherit.aes = FALSE) +
ggplot2::scale_x_continuous(name = xlab, breaks = xBreaks, limits = range(xBreaks)) +
ggplot2::scale_y_continuous(name = gettext("Density"), breaks = yBreaks, limits = c(0, newymax)) +
colorspace::scale_color_discrete_qualitative() +
ggplot2::scale_linetype() +
ggplot2::guides(color = guideLegend, linetype = guideLegend)
p <- jaspGraphs::themeJasp(p, legend.position = if (showLegend) "right" else "none")
plot <- createJaspPlot(title = nms[i], width = 400, height = 400, plot = p)
}
container[[allParamNames[i]]] <- plot
}
} else {
for (i in indices) {
# make prior posterior plot
df <- data.frame(x = densities[, i, "x"],
y = densities[, i, "y"])
p <- jaspGraphs::PlotPriorAndPosterior(
dfLines = df,
xName = xNames[i],
CRI = cris[i, ],
drawCRItxt = FALSE
)
plot <- createJaspPlot(title = plotTitles[i], width = 400, height = 400, plot = p)
container[[allParamNames[i]]] <- plot
}
}
return()
}
.BANOVAqqplot <- function(jaspResults, options, model) {
if (!is.null(jaspResults[["QQplot"]]) || !options[["qqPlot"]])
return()
plot <- createJaspPlot(
title = gettext("Model Averaged Q-Q Plot"),
width = 400,
height = 400,
aspectRatio = 1
)
if (!is.null(model[["models"]])) {
plot$plotObject <- jaspGraphs::plotQQnorm(
residuals = model[["posteriors"]][["weightedResidSumStats"]][,"mean"],
lower = model[["posteriors"]][["weightedResidSumStats"]][,"cri.2.5%"],
upper = model[["posteriors"]][["weightedResidSumStats"]][,"cri.97.5%"],
ablineColor = "darkred"
)
plot$dependOn(optionsFromObject = jaspResults[["tableModelComparisonState"]])
}
plot$dependOn(c("qqPlot", "modelTerms", "credibleInterval", "repeatedMeasuresCells", "dependent"))
plot$position <- 5
jaspResults[["QQplot"]] <- plot
return()
}
.BANOVArsqplot <- function(jaspResults, options, model) {
if (!is.null(jaspResults[["rsqplot"]]) || !options[["rsqPlot"]])
return()
plot <- createJaspPlot(
title = gettextf("Model Averaged Posterior R%s","\u00B2"),
width = 400,
height = 400,
aspectRatio = 1
)
if (!is.null(model[["models"]])) {
dd <- model[["posteriors"]][["weightedRsqDens"]]
rsqCri <- model[["posteriors"]][["weightedRsqCri"]]
df <- data.frame(x = dd$x, y = dd$y)
xName <- expression(R^2)
plot$plotObject <- jaspGraphs::PlotPriorAndPosterior(dfLines = df, xName = xName, CRI = rsqCri, drawCRItxt = FALSE)
plot$dependOn(optionsFromObject = jaspResults[["tableModelComparisonState"]])
}
plot$dependOn(c("rsqPlot", "modelTerms", "credibleInterval", "repeatedMeasuresCells", "dependent"))
plot$position <- 6
jaspResults[["rsqplot"]] <- plot
return()
}
# Post hoc comparison ----
.BANOVAnullControlPostHocTable <- function(jaspResults, dataset, options, model) {
if (length(options$postHocTerms) == 0L)
return()
postHocCollection <- jaspResults[["collectionPosthoc"]]
if (is.null(postHocCollection)) {
postHocCollection <- createJaspContainer(title = gettext("Post Hoc Tests"))
postHocCollection$position <- 8
postHocCollection$addCitation(.BANOVAcitations[3:4])
postHocCollection$dependOn(c("dependent", "repeatedMeasuresCells", "postHocNullControl", "bayesFactorType"))
jaspResults[["collectionPosthoc"]] <- postHocCollection
}
# the same footnote for all the tables
footnote <- gsub("[\r\n\t]", "",
gettext("The posterior odds have been corrected for multiple testing by fixing to 0.5 the prior probability that the null hypothesis holds across all comparisons (Westfall, Johnson, & Utts, 1997). Individual comparisons are based on the default t-test with a Cauchy (0, r = 1/sqrt(2)) prior. The \"U\" in the Bayes factor denotes that it is uncorrected."))
bfTxt <- if (options[["postHocNullControl"]]) ", U" else ""
switch(options[["bayesFactorType"]],
BF10 = { bf.title <- paste0("BF<sub>10", bfTxt, "</sub>") },
BF01 = { bf.title <- paste0("BF<sub>01", bfTxt, "</sub>") },
LogBF10 = { bf.title <- paste0("Log(BF<sub>10", bfTxt, "</sub>)") }
)
priorWidth <- 1 / sqrt(2)
posthoc.variables <- unlist(options[["postHocTerms"]])
if (model[["analysisType"]] == "RM-ANOVA") {
dependent <- .BANOVAdependentName
} else {
dependent <- options[["dependent"]]
}
for (posthoc.var in posthoc.variables) {
# does the table already exist?
if (!is.null(postHocCollection[[paste0("postHoc_", posthoc.var)]]))
next
postHocTable <- createJaspTable(title = paste0("Post Hoc Comparisons - ", posthoc.var))
postHocTable$addColumnInfo(name = "(I)", type = "string", title = "", combine=TRUE)
postHocTable$addColumnInfo(name = "(J)", type = "string", title = "")
postHocTable$addColumnInfo(name = "Prior Odds", type = "number", title = gettext("Prior Odds"))
postHocTable$addColumnInfo(name = "Posterior Odds", type = "number", title = gettext("Posterior Odds"))
postHocTable$addColumnInfo(name = "BF", type = "number", title = bf.title)
postHocTable$addColumnInfo(name = "error %", type = "number", title = gettextf("error %%"))
postHocTable$dependOn(optionContainsValue = list("postHocTerms" = posthoc.var))
if (options[["postHocNullControl"]])
postHocTable$addFootnote(footnote)
if (is.null(model$models)) { # only show empty table
postHocCollection[[paste0("postHoc_", posthoc.var)]] <- postHocTable
next
}
fixed <- unlist(c(options$fixedFactors, sapply(options$repeatedMeasuresFactors, `[[`, "name")))
if (model$analysisType == "RM-ANOVA" && posthoc.var %in% fixed && !posthoc.var %in% options$betweenSubjectFactors) {
variable.levels <- options$repeatedMeasuresFactors[[which(lapply(options$repeatedMeasuresFactors, function(x) x$name) == posthoc.var)]]$levels
paired <- TRUE
} else if (posthoc.var %in% c(options$fixedFactors, options$betweenSubjectFactors, options$randomFactors)) {
variable.levels <- levels(dataset[[posthoc.var]])
paired <- FALSE
} else {
next
}
if (length(variable.levels) < 2L)
next
pairs <- utils::combn(variable.levels, 2)
splitname <- if (model[["analysisType"]] == "RM-ANOVA") .BANOVAdependentName else dependent
allSplits <- split(dataset[[splitname]], dataset[[posthoc.var]])
errMessages <- NULL
for (i in 1:ncol(pairs)) {
if (!is.null(model$models)) {
x <- na.omit(allSplits[[pairs[1L, i]]])
y <- na.omit(allSplits[[pairs[2L, i]]])
ttest <- try(BayesFactor::ttestBF(x = x, y = y, rscale = priorWidth, paired = paired), silent=TRUE)
if (isTryError(ttest)) {
priorOdds <- postOdds <- logBF <- error <- "NaN"
errorMsg <- .extractErrorMessage(ttest)
if (errorMsg %in% names(errMessages)) {
errMessages[[errorMsg]][["row_names"]] <- c(errMessages[[errorMsg]][["row_names"]], i)
} else {
errMessages[[errorMsg]] <- list(
message = if (endsWith(errorMsg, ".")) errorMsg else paste0(errorMsg, "."),
row_names = i
)
}
} else {
if (options[["postHocNullControl"]]) {
pH0 <- 0.5^(2 / length(variable.levels))
} else {
pH0 <- 0.5
}
logBF <- ttest@bayesFactor$bf # <- log(BF10)
if (options[["bayesFactorType"]] == "BF01") {
priorOdds <- pH0 / (1 - pH0)
logBF <- -logBF
} else {
priorOdds <- (1 - pH0) / pH0
}
postOdds <- log(priorOdds) + logBF
postOdds <- exp(postOdds)
if (options[["bayesFactorType"]] != "LogBF10")
logBF <- exp(logBF)
error <- ttest@bayesFactor$error * 100
}
}
row <- list(
"(I)" = pairs[1L, i],
"(J)" = pairs[2L, i],
"Prior Odds" = priorOdds,
"Posterior Odds" = postOdds,
"BF" = logBF,
"error %" = error
)
postHocTable$addRows(row, rowNames = paste0("row", i))
}
if (!is.null(errMessages)) {
for (i in seq_along(errMessages))
postHocTable$addFootnote(message = errMessages[[i]][["message"]],
rowNames = paste0("row", errMessages[[i]][["row_names"]]))
}
postHocCollection[[paste0("postHoc_", posthoc.var)]] <- postHocTable
}
return()
}
# Data reading ----
.BANOVAreadData <- function(dataset, options, analysisType) {
if (is.null(dataset)) {
if (analysisType == "RM-ANOVA")
return(.BANOVAreadRManovaData(dataset, options))
numeric.vars <- c(unlist(options$covariates), unlist(options$dependent))
numeric.vars <- numeric.vars[numeric.vars != ""]
factor.vars <- c(unlist(options$fixedFactors), unlist(options$randomFactors))
factor.vars <- factor.vars[factor.vars != ""]
dataset <- .readDataSetToEnd(
columns.as.numeric = numeric.vars,
columns.as.factor = factor.vars,
exclude.na.listwise = c(numeric.vars, factor.vars)
)
}
return(dataset)
}
.BANOVAdependentName <- "JaspColumn_.dependent._Encoded"
.BANOVAsubjectName <- "JaspColumn_.subject._Encoded"
.BANOVAdecodeSubject <- function(string) {
return(gsub(.BANOVAsubjectName, "subject", string))
}
.BANOVAdecodeNuisance <- function(nuisance) {
# .BANOVAsubjectName needs to be handled separately
idx <- nuisance == .BANOVAsubjectName
nuisance[idx] <- "subject"
nuisance[!idx] <- jaspBase::gsubInteractionSymbol(nuisance[!idx])
return(nuisance)
}
.BANOVAreadRManovaData <- function(dataset, options) {
if (!("" %in% options$repeatedMeasuresCells)) {
rm.vars <- options$repeatedMeasuresCells
bs.factors <- options$betweenSubjectFactors
bs.covariates <- options$covariates
rm.factors <- options$repeatedMeasuresFactors
all.variables <- c (bs.factors, bs.covariates, rm.vars)
dataset <- .readDataSetToEnd(
columns.as.numeric = c(rm.vars, bs.covariates),
columns.as.factor = bs.factors,
exclude.na.listwise = all.variables
)
dataset <- try(
.shortToLong(dataset, rm.factors, rm.vars, c(bs.factors, bs.covariates), dependentName = .BANOVAdependentName, subjectName = .BANOVAsubjectName),
silent = TRUE
)
}
return(dataset)
}
# Descriptives ----
.BANOVAdescriptives <- function(jaspResults, dataset, options, errors, analysisType, ready = TRUE, position = 9001) {
if (!ready)
return()
# the main use of this function is that descriptives can now be reused for the frequentist ANOVAs
# without the container, the position could mess things up
descriptivesContainer <- jaspResults[["descriptivesContainer"]]
if (is.null(descriptivesContainer)) {
descriptivesContainer <- createJaspContainer(title = gettext("Descriptives"))
descriptivesContainer$dependOn(c("dependent", "repeatedMeasuresCells"))
descriptivesContainer$position <- position
jaspResults[["descriptivesContainer"]] <- descriptivesContainer
}
.BANOVAdescriptivesTable (descriptivesContainer, dataset, options, errors, analysisType)
.BANOVAdescriptivesPlots (descriptivesContainer, dataset, options, errors, analysisType)
.BANOVAbarPlots (descriptivesContainer, dataset, options, errors, analysisType)
.BANOVArainCloudPlots (descriptivesContainer, dataset, options, errors, analysisType)
return()
}
.BANOVAdescriptivesTable <- function(jaspContainer, dataset, options, errors, analysisType) {
if (!options[["descriptives"]] || !is.null(jaspContainer[["tableDescriptives"]]))
return()
if (analysisType == "RM-ANOVA") {
dependent <- .BANOVAdependentName
fixed <- unlist(c(lapply(options[["repeatedMeasuresFactors"]], `[[`, "name"), options[["betweenSubjectFactors"]]))
title <- gettext("Descriptives")
} else {
dependent <- options[["dependent"]]
fixed <- options[["fixedFactors"]]
if (!is.null(dependent) && dependent != "") {
title <- gettextf("Descriptives - %s", options$dependent)
} else {
title <- gettext("Descriptives")
}
}
descriptivesTable <- createJaspTable(title = title)
descriptivesTable$position <- 1
descriptivesTable$dependOn(c("dependent", "fixedFactors", "betweenSubjectFactors", "descriptives",
"credibleInterval"))
# internal names use base64 so they don't have " " which R changes into "." because R does that to dataframe names.
# Also adds a . in case the base64 is magically "Mean", "SD" or "N"
fixedDot <- paste0(fixed, ".")
for (i in seq_along(fixed))
descriptivesTable$addColumnInfo(name = fixedDot[i], type = "string", title = fixed[i], combine = TRUE)
descriptivesTable$addColumnInfo(name = "N", title = gettext("N"), type = "integer")
descriptivesTable$addColumnInfo(name = "Mean", title = gettext("Mean"), type = "number")
descriptivesTable$addColumnInfo(name = "SD", title = gettext("SD"), type = "number")
descriptivesTable$addColumnInfo(name = "SE", title = gettext("SE"), type = "number")
descriptivesTable$addColumnInfo(name = "coefOfVariation", title = gettext("Coefficient of variation"), type = "number")
if (!is.null(options[["credibleInterval"]])) {
overTitle <- gettextf("%s%% Credible Interval", format(100 * options[["credibleInterval"]], digits = 3))
descriptivesTable$addColumnInfo(name = gettext("Lower"), type = "number", overtitle = overTitle)
descriptivesTable$addColumnInfo(name = gettext("Upper"), type = "number", overtitle = overTitle)
}
descriptivesTable$showSpecifiedColumnsOnly <- TRUE
jaspContainer[["tableDescriptives"]] <- descriptivesTable
if (errors$noVariables)
return()
# order the data to show
dataset2 <- dataset[do.call(order, dataset[, fixed, drop = FALSE]), c(dependent, fixed)]
# by pasting the fixedFactors together we obtain the unique indices to group on. This excludes
# non-existent combinations. A "." is added to deal with the level "".
ind <- apply(dataset2[, fixed, drop = FALSE], 1L, paste0, ".", collapse = "")
# temporary function to calculate all descriptives
tmpFun <- function(data, fixed, dependent, cri) {
row <- list()
for (j in fixed)
row[[paste0(j, ".")]] <- as.character(data[1L, j])
N <- nrow(data)
row[["N"]] <- N
if (N == 0L) {
row[["Mean"]] <- row[["SD"]] <- row[["SE"]] <- row[["coefOfVariation"]] <- NA
} else if (N == 1L) {
row[["Mean"]] <- data[[dependent]]
row[["SD"]] <- row[["SE"]] <- row[["coefOfVariation"]] <- row[["Lower"]] <- row[["Upper"]] <- NA
} else {
row[["Mean"]] <- mean(data[[dependent]])
row[["SD"]] <- stats::sd(data[[dependent]])
row[["SE"]] <- stats::sd(data[[dependent]]) / sqrt(N)
row[["coefOfVariation"]] <- stats::sd(data[[dependent]]) / mean(data[[dependent]])
tmp <- jaspTTests::.posteriorSummaryGroupMean(data[[dependent]], cri)
row[["Lower"]] <- tmp[["ciLower"]]
row[["Upper"]] <- tmp[["ciUpper"]]
}
return(row)
}
# apply tempFun on each subset defined by ind
rows <- by(dataset2, ind, tmpFun, fixed = fixed,
dependent = dependent, cri = options[["credibleInterval"]])
# do.call(rbind, rows) turns rows into a data.frame (from a list) for jaspResults
data <- do.call(rbind.data.frame, rows)
# the results of `by` are in alphabetical order, so we reorder the data.frame to match the factor level order. Fixes https://github.com/jasp-stats/jasp-issues/issues/1741
newOrder <- match(unique(ind), rownames(data))
data <- data[newOrder, ]
# add footnote if there are unobserved combinations
nObserved <- nrow(data)
nPossible <- prod(sapply(dataset2[, fixed, drop = FALSE], nlevels))
if (nObserved != nPossible) {
descriptivesTable$addFootnote(
message = gettextf(
"Some combinations of factors are not observed and hence omitted (%1$g out of %2$g combinations are unobserved).",
nPossible - nObserved, nPossible
)
)
}
descriptivesTable$setData(data)
return()
}
.BANOVAdescriptivesPlots <- function(jaspContainer, dataset, options, errors, analysisType) {
if (length(options[["descriptivePlotHorizontalAxis"]]) == 0L
|| options[["descriptivePlotHorizontalAxis"]] == ""
|| !is.null(jaspContainer[["containerDescriptivesPlots"]]))
return()
descriptivesPlotContainer <- createJaspContainer(title = gettext("Descriptives plots"))
descriptivesPlotContainer$position <- 2
jaspContainer[["containerDescriptivesPlots"]] <- descriptivesPlotContainer
# either Bayesian or Frequentist anova
if (is.null(options$descriptivePlotErrorBarType)) { # TRUE implies Bayesian
plotErrorBars <- options$descriptivePlotCi
errorBarType <- "ci"
conf.interval <- options$descriptivePlotCiLevel
normalizeErrorBars <- TRUE
descriptivesPlotContainer$dependOn(c("dependent", "descriptivePlotCi", "descriptivePlotCiLevel"))
} else {
plotErrorBars <- options$descriptivePlotErrorBar
errorBarType <- options$descriptivePlotErrorBarType
conf.interval <- options$descriptivePlotCiLevel
normalizeErrorBars <- if (is.null(options[["normalizeErrorBarsDescriptives"]])) FALSE else options[["normalizeErrorBarsDescriptives"]]
descriptivesPlotContainer$dependOn(c("dependent", "descriptivePlotErrorBar", "descriptivePlotErrorBarType", "descriptivePlotCiLevel",
"descriptivePlotErrorBarPooled", "normalizeErrorBarsDescriptives"))
}
usePooledSE <- if (is.null(options[["descriptivePlotErrorBarPooled"]])) FALSE else options[["descriptivePlotErrorBarPooled"]]
descriptivesPlotContainer$dependOn(c("descriptivePlotHorizontalAxis", "descriptivePlotSeparateLines", "descriptivePlotSeparatePlot", "descriptivePlotYAxisLabel"))
if (errors$noVariables) {
descriptivesPlotContainer[["dummyplot"]] <- createJaspPlot(title = gettext("Descriptives Plot"))
return()
}
groupVars <- c(options$descriptivePlotHorizontalAxis, options$descriptivePlotSeparateLines, options$descriptivePlotSeparatePlot)
groupVars <- groupVars[groupVars != ""]
if (analysisType == "RM-ANOVA") {
dependent <- .BANOVAdependentName
yLabel <- options[["descriptivePlotYAxisLabel"]]
} else {
dependent<- options$dependent
yLabel <- options[["dependent"]]
}
betweenSubjectFactors <- groupVars[groupVars %in% options$betweenSubjectFactors]
repeatedMeasuresFactors <- groupVars[groupVars %in% sapply(options$repeatedMeasuresFactors, `[[`, "name")]
if (length(repeatedMeasuresFactors) == 0 || isFALSE(normalizeErrorBars)) {
summaryStat <- jaspTTests::.summarySE(as.data.frame(dataset), measurevar = dependent, groupvars = groupVars,
conf.interval = conf.interval, na.rm = TRUE, .drop = FALSE,
errorBarType = errorBarType, dependentName = .BANOVAdependentName,
subjectName = .BANOVAsubjectName)
} else {
summaryStat <- jaspTTests::.summarySEwithin(as.data.frame(dataset), measurevar= .BANOVAdependentName,
betweenvars = betweenSubjectFactors,
withinvars = repeatedMeasuresFactors,
idvar = .BANOVAsubjectName,
conf.interval = conf.interval,
na.rm=TRUE, .drop = FALSE, errorBarType = errorBarType,
usePooledSE = usePooledSE,
dependentName = .BANOVAdependentName,
subjectName = .BANOVAsubjectName)
}
if (options[["descriptivePlotHorizontalAxis"]] %in% options[["covariates"]]) {
splitScatterOptions <- options
splitScatterOptions[["colorPalette"]] <- "colorblind3"
splitScatterOptions[["scatterPlotLegend"]] <- TRUE
splitScatterOptions[["scatterPlotRegressionLine"]] <- TRUE
splitScatterOptions[["scatterPlotRegressionLineCi"]] <- plotErrorBars
splitScatterOptions[["scatterPlotRegressionLineType"]] <- "linear"
splitScatterOptions[["scatterPlotGraphTypeAbove"]] <- "none"
splitScatterOptions[["scatterPlotGraphTypeRight"]] <- "none"
splitScatterOptions[["scatterPlotRegressionLineCiLevel"]] <- options[["descriptivePlotCiLevel"]]
if (options$descriptivePlotSeparatePlot != "") {
for (thisLevel in levels(dataset[[options[["descriptivePlotSeparatePlot"]]]])) {
subData <- dataset[dataset[[options[["descriptivePlotSeparatePlot"]]]] == thisLevel, ]
thisPlotName <- paste0(options[["descriptivePlotHorizontalAxis"]], " - ", options[["dependent"]], ": ",
options[["descriptivePlotSeparatePlot"]], " = ", thisLevel)
jaspDescriptives::.descriptivesScatterPlots(descriptivesPlotContainer, subData, c(options[["descriptivePlotHorizontalAxis"]], options[["dependent"]]),
split = options[["descriptivePlotSeparateLines"]], options = splitScatterOptions, name = thisPlotName,
dependOnVariables = FALSE)
}
} else {
jaspDescriptives::.descriptivesScatterPlots(descriptivesPlotContainer, dataset, c(options[["descriptivePlotHorizontalAxis"]], options[["dependent"]]),
split = options[["descriptivePlotSeparateLines"]], options = splitScatterOptions, dependOnVariables = FALSE)
}
return()
}
colnames(summaryStat)[colnames(summaryStat) == dependent] <- "dependent"
if (options$descriptivePlotHorizontalAxis != "") {
colnames(summaryStat)[colnames(summaryStat) == options$descriptivePlotHorizontalAxis] <- "descriptivePlotHorizontalAxis"
}
if (options$descriptivePlotSeparateLines != "") {
colnames(summaryStat)[colnames(summaryStat) == options$descriptivePlotSeparateLines] <- "descriptivePlotSeparateLines"
}
if (options$descriptivePlotSeparatePlot != "") {
colnames(summaryStat)[colnames(summaryStat) == options$descriptivePlotSeparatePlot] <- "descriptivePlotSeparatePlot"
}
base_breaks_x <- function(x){
b <- unique(as.numeric(x))
d <- data.frame(y=-Inf, yend=-Inf, x=min(b), xend=max(b))
list(ggplot2::geom_segment(data=d, ggplot2::aes(x=x, y=y, xend=xend, yend=yend), inherit.aes=FALSE, size = 1))
}
base_breaks_y <- function(x, plotErrorBars){
if (plotErrorBars) {
ci.pos <- c(x[,"dependent"], x[,"dependent"]-x[,"ci"],x[,"dependent"]+x[,"ci"])
b <- pretty(ci.pos)
d <- data.frame(x=-Inf, xend=-Inf, y=min(b), yend=max(b))
list(ggplot2::geom_segment(data=d, ggplot2::aes(x=x, y=y, xend=xend, yend=yend), inherit.aes=FALSE, size = 1),
ggplot2::scale_y_continuous(breaks=c(min(b),max(b))))
} else {
b <- pretty(x[,"dependent"])
d <- data.frame(x=-Inf, xend=-Inf, y=min(b), yend=max(b))
list(ggplot2::geom_segment(data=d, ggplot2::aes(x=x, y=y, xend=xend, yend=yend), inherit.aes=FALSE, size = 1),
ggplot2::scale_y_continuous(breaks=c(min(b),max(b))))
}
}
if (options$descriptivePlotSeparatePlot != "") {
subsetPlots <- levels(summaryStat[,"descriptivePlotSeparatePlot"])
nPlots <- length(subsetPlots)
} else {
nPlots <- 1
}
for (i in seq_len(nPlots)) {
if (nPlots > 1L) {
title <- paste(options$descriptivePlotSeparatePlot,": ",subsetPlots[i], sep = "")
} else {
title <- ""
}
descriptivesPlot <- createJaspPlot(title = title)
descriptivesPlotContainer[[title]] <- descriptivesPlot
descriptivesPlot$height <- 300
if (options$descriptivePlotSeparateLines != "") {
descriptivesPlot$width <- 430
} else {
descriptivesPlot$width <- 300
}
if (options$descriptivePlotSeparatePlot != "") {
summaryStatSubset <- subset(summaryStat,summaryStat[,"descriptivePlotSeparatePlot"] == subsetPlots[i])
} else {
summaryStatSubset <- summaryStat
}
if (options$descriptivePlotSeparateLines == "") {
p <- ggplot2::ggplot(summaryStatSubset, ggplot2::aes(x=descriptivePlotHorizontalAxis,
y=dependent,
group=1))
} else {
p <- ggplot2::ggplot(summaryStatSubset, ggplot2::aes(x=descriptivePlotHorizontalAxis,
y=dependent,
group=descriptivePlotSeparateLines,
shape=descriptivePlotSeparateLines,
fill=descriptivePlotSeparateLines))
}
if (plotErrorBars && !(options[["descriptivePlotHorizontalAxis"]] %in% options[["covariates"]])) {
pd <- ggplot2::position_dodge(.2)
p = p + ggplot2::geom_errorbar(ggplot2::aes(ymin=ciLower,
ymax=ciUpper),
colour="black", width=.2, position=pd)
} else {
pd <- ggplot2::position_dodge(0)
}
guideLegend <- ggplot2::guide_legend(nrow = min(10, nlevels(summaryStatSubset$descriptivePlotSeparateLines)),
title = options$descriptivePlotSeparateLines, keywidth = 0.1, keyheight = 0.3,
default.unit = "inch")
if (options[["descriptivePlotHorizontalAxis"]] %in% options[["covariates"]]) {
line <- ggplot2::geom_smooth(method = "lm", size = .7, color = "black", se = FALSE)
addHorizontalVar <- summaryStatSubset[,"descriptivePlotHorizontalAxis"]
} else {
line <- ggplot2::geom_line(position=pd, size = .7)
}
if (plotErrorBars) {
ci.pos <- c(summaryStatSubset[,"dependent"],
summaryStatSubset[,"dependent"]-summaryStatSubset[,"ci"],
summaryStatSubset[,"dependent"]+summaryStatSubset[,"ci"],
min(summaryStatSubset[,"dependent"])*1.1,
max(summaryStatSubset[,"dependent"])*1.1)
yBreaks <- jaspGraphs::getPrettyAxisBreaks(ci.pos)
} else {
yBreaks <- jaspGraphs::getPrettyAxisBreaks(c(summaryStatSubset[,"dependent"],
min(summaryStatSubset[,"dependent"])*1.1,
max(summaryStatSubset[,"dependent"])*1.1))
}
if (options[["descriptivePlotHorizontalAxis"]] %in% options[["covariates"]]) {
ggXaxis <- ggplot2::scale_x_continuous(breaks = jaspGraphs::getPrettyAxisBreaks(summaryStatSubset[,"descriptivePlotHorizontalAxis"]))
} else {
ggXaxis <- ggplot2::scale_x_discrete(breaks = jaspGraphs::getPrettyAxisBreaks(summaryStatSubset[,"descriptivePlotHorizontalAxis"]))
}
p <- p + line +
ggplot2::geom_point(position=pd, size=4) +
ggplot2::scale_fill_manual(values = c(rep(c("white","black"),5),rep("grey",100)), guide=guideLegend) +
ggplot2::scale_shape_manual(values = c(rep(c(21:25),each=2),21:25,7:14,33:112), guide=guideLegend) +
ggplot2::scale_color_manual(values = rep("black",200),guide=guideLegend) +
ggplot2::labs(y = yLabel, x = options[["descriptivePlotHorizontalAxis"]]) +
ggplot2::scale_y_continuous(breaks = yBreaks, limits = range(yBreaks)) +
ggXaxis +
jaspGraphs::geom_rangeframe() +
jaspGraphs::themeJaspRaw(legend.position = "right")
descriptivesPlot$plotObject <- p
}
return()
}
.BANOVAbarPlots <- function(jaspContainer, dataset, options, errors, analysisType) {
if (length(options[["barPlotHorizontalAxis"]]) == 0L
|| options[["barPlotHorizontalAxis"]] == ""
|| !is.null(jaspContainer[["containerBarPlots"]]))
return()
barPlotContainer <- createJaspContainer(title = gettext("Bar plots"))
barPlotContainer$position <- 3
jaspContainer[["containerBarPlots"]] <- barPlotContainer
plotErrorBars <- options[["barPlotErrorBars"]]
errorBarType <- options[["barPlotErrorBarType"]]
confInterval <- options[["barPlotCiInterval"]]
barPlotContainer$dependOn(c("dependent", "barPlotErrorBars", "barPlotErrorBarType", "applyMoreyCorrectionErrorBarsBarplot",
"barPlotHorizontalZeroFix", "barPlotCiInterval", "usePooledStandErrorCITwo"))
usePooledSE <- if (is.null(options[["usePooledStandErrorCITwo"]])) FALSE else options[["usePooledStandErrorCITwo"]]
normalizeErrorBars <- if (is.null(options[["normalizeErrorBarsBarplot"]])) TRUE else options[["normalizeErrorBarsBarplot"]]
barPlotHorizontalZeroFix <- options[["barPlotHorizontalZeroFix"]]
barPlotContainer$dependOn(c("barPlotHorizontalAxis", "barPlotSeparatePlots", "labelYAxisTwo"))
if (errors$noVariables) {
barPlotContainer[["dummyplot"]] <- createJaspPlot(title = gettext("Bar Plot"))
return()
}
groupVars <- c(options[["barPlotHorizontalAxis"]], options[["barPlotSeparatePlots"]])
groupVars <- groupVars[groupVars != ""]
if (analysisType == "RM-ANOVA") {
dependent <- .BANOVAdependentName
yLabel <- options[["labelYAxisTwo"]]
} else {
dependent<- options[["dependent"]]
yLabel <- options[["dependent"]]
}
betweenSubjectFactors <- groupVars[groupVars %in% options[["betweenSubjectFactors"]]]
repeatedMeasuresFactors <- groupVars[groupVars %in% sapply(options[["repeatedMeasuresFactors"]], `[[`, "name")]
if (length(repeatedMeasuresFactors) == 0 || isFALSE(normalizeErrorBars)) {
summaryStat <- jaspTTests::.summarySE(as.data.frame(dataset), measurevar = dependent, groupvars = groupVars,
conf.interval = confInterval, na.rm = TRUE, .drop = FALSE,
errorBarType = errorBarType, dependentName = .BANOVAdependentName,
subjectName = .BANOVAsubjectName)
} else {
summaryStat <- jaspTTests::.summarySEwithin(as.data.frame(dataset), measurevar = .BANOVAdependentName,
betweenvars = betweenSubjectFactors,
withinvars = repeatedMeasuresFactors,
idvar = .BANOVAsubjectName,
conf.interval = confInterval,
na.rm=TRUE, .drop = FALSE, errorBarType = errorBarType,
usePooledSE = usePooledSE,
dependentName = .BANOVAdependentName,
subjectName = .BANOVAsubjectName)
}
colnames(summaryStat)[colnames(summaryStat) == dependent] <- "dependent"
if (options[["barPlotHorizontalAxis"]] != "") {
colnames(summaryStat)[colnames(summaryStat) == options[["barPlotHorizontalAxis"]]] <- "barPlotHorizontalAxis"
}
if (options[["barPlotSeparatePlots"]] != "") {
colnames(summaryStat)[colnames(summaryStat) == options[["barPlotSeparatePlots"]]] <- "barPlotSeparatePlots"
subsetPlots <- levels(summaryStat[, "barPlotSeparatePlots"])
nPlots <- length(subsetPlots)
} else {
nPlots <- 1L
}
for (i in seq_len(nPlots)) {
title <- if (nPlots > 1L) {
paste0(options[["barPlotSeparatePlots"]], ": ", subsetPlots[i])
} else {
""
}
barPlot <- createJaspPlot(title = title)
barPlotContainer[[title]] <- barPlot
barPlot$height <- 500
barPlot$width <- 500
if (options[["barPlotSeparatePlots"]] != "") {
summaryStatSubset <- subset(summaryStat, summaryStat[, "barPlotSeparatePlots"] == subsetPlots[i])
} else {
summaryStatSubset <- summaryStat
}
error <- NULL
if (plotErrorBars) {
pd <- ggplot2::position_dodge(.2)
error <- ggplot2::geom_errorbar(ggplot2::aes(ymin = ciLower, ymax = ciUpper),
colour = "black", width = .2, position = pd)
}
values <- 1.1 * range(summaryStat[["dependent"]])
if (barPlotHorizontalZeroFix)
values <- c(0, values)
if (plotErrorBars) {
ci.pos <- c(summaryStat[,"dependent"],
summaryStat[,"dependent"]-summaryStat[,"ci"],
summaryStat[,"dependent"]+summaryStat[,"ci"],
min(summaryStat[,"dependent"])*1.1,
max(summaryStat[,"dependent"])*1.1)
yBreaks <- jaspGraphs::getPrettyAxisBreaks(c(values, ci.pos))
} else {
yBreaks <- jaspGraphs::getPrettyAxisBreaks(values)
}
pd2 <- ggplot2::position_dodge2(preserve = "single")
p <- ggplot2::ggplot(summaryStatSubset, ggplot2::aes(x = barPlotHorizontalAxis, y = dependent, group = 1)) +
ggplot2::geom_hline(yintercept = 0, color = "#858585", size = 0.3) +
ggplot2::geom_bar(stat = "identity", fill = "grey", col = "black", width = .6, position = pd2) +
error +
ggplot2::labs(y = yLabel, x = options[["barPlotHorizontalAxis"]]) +
ggplot2::scale_y_continuous(breaks = yBreaks, limits = range(yBreaks), oob = scales::rescale_none) +
ggplot2::scale_x_discrete(breaks = jaspGraphs::getPrettyAxisBreaks(summaryStatSubset[, "barPlotHorizontalAxis"])) +
jaspGraphs::geom_rangeframe(sides = "l") +
jaspGraphs::themeJaspRaw()
barPlot$plotObject <- p
}
return()
}
.BANOVArainCloudPlots <- function(jaspContainer, dataset, options, errors, analysisType) {
if (length(options[["rainCloudHorizontalAxis"]]) == 0L
|| options[["rainCloudHorizontalAxis"]] == ""
|| !is.null(jaspContainer[["containerRainCloudPlots"]]))
return()
rainCloudPlotsContainer <- createJaspContainer(title = gettext("Raincloud plots"))
rainCloudPlotsContainer$position <- 4
jaspContainer[["containerRainCloudPlots"]] <- rainCloudPlotsContainer
rainCloudPlotsContainer$dependOn(c("dependent", "rainCloudHorizontalAxis", "rainCloudSeparatePlots",
"rainCloudYAxisLabel", "rainCloudHorizontalDisplay"))
if (errors$noVariables) {
rainCloudPlotsContainer[["dummyplot"]] <- createJaspPlot(title = "")
return()
}
groupVar <- options[["rainCloudHorizontalAxis"]]
if (analysisType == "RM-ANOVA") {
addLines <- !(groupVar %in% unlist(options[["betweenSubjectFactors"]]))
dependentV <- .BANOVAdependentName
yLabel <- options[["rainCloudYAxisLabel"]]
if (trimws(yLabel) == "") {
title <- gettext("Dependent")
yLabel <- NULL
} else {
title <- yLabel
}
} else {
addLines <- FALSE
dependentV <- options[["dependent"]]
yLabel <- options[["dependent"]]
title <- options[["dependent"]]
}
if (!is.null(options$rainCloudHorizontalDisplay) && options$rainCloudHorizontalDisplay)
horiz <- TRUE
else
horiz <- FALSE
if (options$rainCloudSeparatePlots != "") {
for (thisLevel in levels(dataset[[options[["rainCloudSeparatePlots"]]]])) {
subData <- dataset[dataset[[options[["rainCloudSeparatePlots"]]]] == thisLevel, ]
thisPlotName <- paste0(title, ": ", options[["rainCloudSeparatePlots"]], ": ", thisLevel)
subPlot <- createJaspPlot(title = thisPlotName, width = 480, height = 320)
rainCloudPlotsContainer[[thisLevel]] <- subPlot
p <- try(jaspTTests::.descriptivesPlotsRainCloudFill(subData, dependentV, groupVar, yLabel, groupVar, addLines, horiz, NULL))
if(isTryError(p))
subPlot$setError(.extractErrorMessage(p))
else
subPlot$plotObject <- p
}
} else {
singlePlot <- createJaspPlot(title = title, width = 480, height = 320)
rainCloudPlotsContainer[["rainCloudPlotSingle"]] <- singlePlot
p <- try(jaspTTests::.descriptivesPlotsRainCloudFill(dataset, dependentV, groupVar, yLabel, groupVar, addLines, horiz, NULL))
if(isTryError(p))
singlePlot$setError(.extractErrorMessage(p))
else
singlePlot$plotObject <- p
}
return()
}
# Sample posteriors ----
.BANOVAsamplePosteriors <- function(jaspResults, dataset, options, model, state) {
# TODO: the density approximation can become more efficient with a fast parametric density approximation
# TODO: Bayes factor samples unobserved interaction levels, what to do?
# if the most complex model is retrieved from the state?
nIter <- if (options[["samplingMethodMCMC"]] == "auto") 1e4L else options[["samplesMCMC"]]
nmodels <- length(model[["models"]])
postProbs <- model[["postProbs"]]
statistics <- vector("list", nmodels)
randomFactors <- .BANOVAgetRandomFactors(options, model[["analysisType"]])
dataTypes <- .BANOVAgetDataTypes(dataset, model[["model.list"]][[nmodels]], randomFactors)
levelInfo <- .BANOVAgetLevelInfo(dataset, model[["model.list"]][[nmodels]], dataTypes)
renameFrom <- renameTo <- NULL
if (!is.null(state)) { # can we reuse some posteriors?
reuseable <- model[["reuseable"]]
# it's possible that a user just renames a level of a repeated measures factor
# in that case everything can be reused, but we have to rename some parameters.
# the same may happen when a variable becomes a nuisance parameter, because then
# BayesFactor changes its order in the names.
oldLevelInfo <- state[["levelInfo"]]$levelNames
newLevelInfo <- levelInfo$levelNames
sortTerms <- function(x) {
# split b:a to c("b", "a"), sort it, and then paste it back
# otherwise posteriors samples between different runs are not correctly retrieved from the state.
sapply(strsplit(x, ":", fixed = TRUE), function(x) paste(sort(x), collapse = ":"))
}
tmp_old <- sortTerms(names(oldLevelInfo))
tmp_new <- sortTerms(names(newLevelInfo))
matches <- match(tmp_new, tmp_old)
for (i in seq_along(matches)) {
if (is.na(matches[i]))
next
j <- matches[i]
if (!identical(oldLevelInfo[[j]], newLevelInfo[[i]])) {
if (length(oldLevelInfo[[j]]) == length(newLevelInfo[[i]])) {
renameFrom <- c(renameFrom, oldLevelInfo[[j]])
renameTo <- c(renameTo, newLevelInfo[[i]])
}
}
}
} else {
reuseable <- rep(NA, nmodels)
}
# NOTE: some code checks require saving all samples. To do so, change all samples to samples[[i]] and uncomment:
# samples <- vector("list", nmodels)
allParamNames <- c("mu", unlist(levelInfo$levelNames))
nparam <- length(allParamNames)
weightedMeans <- weights <- numeric(nparam)
names(weights) <- names(weightedMeans) <- allParamNames
allContinuous <- TRUE
h <- (1 - options[["credibleInterval"]]) / 2
probs <- c(h, 1-h)
originalFun <- BayesFactor:::makeChainNeater
on.exit(jaspBase:::assignFunctionInPackage(originalFun, "makeChainNeater", "BayesFactor"))
jaspBase:::assignFunctionInPackage(.BANOVAmakeChainNeater, "makeChainNeater", "BayesFactor")
samplingIssues <- list()
startProgressbar(nmodels, gettext("Sampling posteriors"))
for (i in seq_len(nmodels)) {
# check if the state is reuseable and if it's not NULL, which means it didn't crash
if (is.na(reuseable[i]) || is.null(state$statistics[[reuseable[i]]])) {
if (i == 1L && is.null(model[["models"]][[i]][["bf"]])) {
# NULL model only contains an intercept, use custom sampler
# NOTE: RM-ANOVA never enters here (and would crash if it did)
.setSeedJASP(options)
samples <- .BANOVAsampleNullModel(dataset[[options[["dependent"]]]], nsamples = nIter)
types <- NULL
} else {
# NOTE: we have to sample the random effects, otherwise we cant make predictions (needed for residuals and R^2)
# put the dataset back in
.setSeedJASP(options)
bfObj <- model[["models"]][[i]][["bf"]]
bfObj@data <- dataset
samples <- try(BayesFactor::posterior(bfObj, iterations = nIter))
if (isTryError(samples)) {
# this may error whenever the Bayes factor couldn't be calculated
next
}
types <- samples@model@dataTypes
# BayesFactor stores an internal copy of the dataset, so it can still have old names
if (length(renameFrom) > 0) {
cnms <- colnames(samples)
cnms <- plyr::mapvalues(cnms, renameFrom, renameTo, warn_missing = FALSE) # NOTE: dependency could be removed
colnames(samples) <- cnms
}
# keep only relevant columns, drop sig2, g_xxx, ...
idx <- match(allParamNames, colnames(samples), nomatch = 0L)
samples <- samples[, idx, drop = FALSE]
# for some odd reason, Bayesfactor uses as column name contcor1-contcor1
# if there is a covariate AND fixed factors, but only contcor1 if all variables are continuous...
if (all(types == "continuous")) {
cnms <- colnames(samples)[-1L] # omit the intercept (mu) which is not changed by Bayesfactor
colnames(samples)[-1L] <- paste0(cnms, "-", cnms)
} else {
allContinuous <- FALSE
}
# so some models may yield a bunch of NAs, for example,
# BayesFactor::lmBF(contNormal ~ facGender + contGamma + facGender * contGamma, data = "debug.csv", whichRandom = "facGender")
# we add a footnote and try to not to crash.
if (anyNA(samples)) {
originalRows <- nrow(samples)
samples <- samples[complete.cases(samples), , drop = FALSE]
remainingRows <- nrow(samples)
samplingIssues[[length(samplingIssues) + 1L]] <- list(
model = model[["models"]][[i]][["title"]],
originalRows = originalRows,
remainingRows = remainingRows
)
next
}
}
# although matrixStats::colMeans2 is faster than .colMeans the cost of matrixStats:: is not worth it.
nms <- colnames(samples)
statistics[[i]]$names <- nms # <- these are the names for all objects within one sublist
statistics[[i]]$mean <- .colMeans (samples, m = nIter, n = NCOL(samples))
statistics[[i]]$var <- matrixStats::colVars (samples)
statistics[[i]]$cri <- matrixStats::colQuantiles(samples, probs = probs)
statistics[[i]]$approx <- .BANOVAfitDensity(samples = samples)
statistics[[i]]$types <- types
} else { # reuse state
statistics[[i]] <- state$statistics[[reuseable[i]]]
nms <- statistics[[i]]$names
types <- statistics[[i]]$types
if (length(renameFrom) > 0) {
nms <- plyr::mapvalues(nms, renameFrom, renameTo, warn_missing = FALSE) # NOTE: dependency could be removed
statistics[[i]]$names <- nms
}
if (allContinuous)
allContinuous <- all(types == "continuous")
}
# compute model averaged posterior means
weightedMeans[nms] <- weightedMeans[nms] + postProbs[i] * statistics[[i]]$mean
weights [nms] <- weights [nms] + postProbs[i]
progressbarTick()
}
# find out which parameters are random -- this uses types from the last iteration above
isRandom <- logical(nparam)
idx <- which(types == "random")
for (i in idx)
isRandom <- isRandom | startsWith(names(weights), names(types)[i])
# the weights used above don't sum to 1 because we consider a subset of the models.
# Now we renormalize to ensure the weights used in each weighted mean do sum to 1.
weightedMeans <- weightedMeans / weights
# given the model averaged posterior means calculate the weighted posterior standard deviations
weightedSds <- 0 * weightedMeans # keeps the names
r <- nIter / (nIter - 1)
# statistics has length 0 iff the posterior sampling errored
widxGoodStatistics <- which(lengths(statistics) > 0L)
for (i in widxGoodStatistics) {
nms <- statistics[[i]]$names
var <- statistics[[i]]$var # ~ sum((x - mean(x))^2
mu <- statistics[[i]]$mean
# cc = sum((x - y)^2) - sum((x - mean(x))^2), where y is the weighted mean
# hence, we can get the weighted sum of squares from the individual posterior variances
cc <- weightedMeans[nms]^2 + mu^2 - 2 * mu * weightedMeans[nms]
weightedSds[nms] <- weightedSds[nms] + postProbs[i] * (var + cc * r)
}
weightedSds <- sqrt(weightedSds / weights)
# the loops above are optimized versions of the code below that also calculates the
# weighted mean and weighted sd, but needs to store all samples at once.
# xx <- lapply(samples, function(x, y) {
# if (y %in% colnames(x)) {dd <- model$posteriors$weightedRsqDens
# x[, y]
# } else {
# NULL
# }
# }, y = "XY29udEJpbm9t-1") # XY29udEJpbm9t = contBinom
# idx <- lengths(xx) > 0
# xx <- unlist(xx[idx])
# w <- postProbs[idx]
# sum(w) # equals weights[2]
# weighted.mean(xx, rep(w / sum(w), each = nIter)) # equals weightedMeans[2]
# sqrt(Hmisc::wtd.var(xx, rep(w / sum(w), each = nIter), method = "unbiased")) # equals weightedSds[2]
# compute model averaged densities
steps <- 2^9 # grid size for densities is 2^steps
weightedDensities <- array(0, dim = c(steps, nparam, 2), dimnames = list(NULL, names(weightedMeans), c("x", "y")))
ranges <- matrix(0, nparam, 2L, dimnames = list(allParamNames, NULL))
# get the outermost x-values for each densities this could be vectorized with matrixStats::rowRanges if the
# data are stored as a matrix but that is memory inefficient
for (i in widxGoodStatistics) {
nms <- statistics[[i]]$names
indices <- which(nms %in% allParamNames)
for (j in indices) {
ranges[j, ] <- range(ranges[j, ], statistics[[i]]$approx$xRanges[j, ])
}
}
# construct one common grid for each observed density
for (i in seq_len(nparam)) {
weightedDensities[, i, 1L] <- seq(ranges[i, 1], ranges[i, 2], length.out = steps)
}
for (i in widxGoodStatistics) {
nms <- statistics[[i]]$names
ind <- match(nms, allParamNames)
for (j in seq_along(ind)) {
# approximate all distributions on a common grid
ap <- approx(x = statistics[[i]]$approx$fit[, j, 1L],
y = statistics[[i]]$approx$fit[, j, 2L],
xout = weightedDensities [, ind[j], 1L],
# not observed is approximated to 0
yleft = 0, yright = 0)
weightedDensities[, ind[j], 2L] <- weightedDensities[, ind[j], 2L] + ap$y * postProbs[i]
}
}
# postProbs don't sum to one so we renormalize the densities
weightedDensities[, , 2L] <- sweep(weightedDensities[, , 2L], 2, weights, FUN = `/`)
# if all model with one parameter failed, the line above introduces NaNs since 0 / 0 is NaN
weightedDensities[is.nan(weightedDensities)] <- 0
# # compute weighted CRIs
# weightedCRIs <- matrix(NA, nparam, 2L, dimnames = list(names(weights), NULL))
# cri <- options[["credibleInterval"]]
# for (i in seq_len(nparam)) {
# weightedCRIs[i, ] <- .BANOVAapproxCRI(weightedDensities[, i, 1L], weightedDensities[, i, 2L], cri)
# }
#
# # compute residuals and r-squared
# # sample from the joint posterior over models and parameters
# tmp <- .BANOVAgetSMIResidRsq(weightedDensities, dataset, model$model.list[[nmodels]], nIter, weights)
# means <- rowMeans(tmp$resid)
# h <- (1 - cri) / 2
# quants <- matrixStats::rowQuantiles(tmp$resid, probs = c(h, 1 - h))
#
# # the code above is equivalent to the code below, but the code below needs to keep all posterior samples of
# # all models in memory.
# # weights <- rep(postProbs, each = nIter)
# # independentVariable <- all.vars(.BANOVAgetModelFormulaFromBFobj(model$models[[2L]]))[1L]
# # resids <- matrix(NA, nrow(dataset), 0)
# # rsq <- vector("list", nmodels)
# # # get residuals of all models individually
# # for (i in seq_len(nmodels)) {
# # if (is.null(model$models[[i]]$bf)) {
# # tmp2 <- .BANOVAresidualsNullModel(nIter, dataset[[independentVariable]])
# # } else {
# # tmp2 <- .BANOVAgetSMIResidRsq(
# # posterior = samples[[i]],
# # dataset = dataset,
# # formula = .BANOVAgetModelFormulaFromBFobj(model$models[[i]])
# # )
# # }
# # resids <- cbind(resids, tmp2$resid)
# # rsq[[i]] <- tmp2$rsq
# # }
# # # compute weighted mean for each row
# # means2 <- tcrossprod(weights / nIter, resids)
# # plot(means, means2); abline(0, 1)
# # quants2 <- apply(resids, 1L, Hmisc::wtd.quantile, weights = weights, probs = c(0.025, 0.975))
# # plot(quants[, 1], quants2[1, ]); abline(0, 1)
# # plot(quants[, 2], quants2[2, ]); abline(0, 1)
#
# # all information for q-q plot of residuals
# weightedResidSumStats <- matrix(c(means, quants), nrow = length(means), ncol = 3L,
# dimnames = list(NULL, c("mean", "cri.2.5%", "cri.97.5%")))
#
# # all information for r-squared density plot
# weightedRsqDens <- density(tmp$rsq, n = 2^11, from = 0, to = 1)
# weightedRsqCri <- quantile(tmp$rsq, probs = c(h, 1 - h))
return(list(
statistics = statistics, weights = weights, weightedMeans = weightedMeans, weightedSds = weightedSds,
weightedDensities = weightedDensities,
#weightedCRIs = weightedCRIs, weightedResidSumStats = weightedResidSumStats,
#weightedRsqDens = weightedRsqDens, weightedRsqCri = weightedRsqCri,
allContinuous = allContinuous, isRandom = isRandom, levelInfo = levelInfo,
samplingIssues = samplingIssues
))
}
.BANOVAcomputePosteriorCRI <- function(dataset, options, model, posterior) {
nIter <- if (options[["samplingMethodMCMC"]] == "auto") 1e4L else options[["samplesMCMC"]]
weightedDensities <- posterior[["weightedDensities"]]
weights <- posterior[["weights"]]
nmodels <- length(model[["models"]])
nparam <- length(weights)
# compute weighted CRIs
weightedCRIs <- matrix(0, nparam, 2L, dimnames = list(names(weights), NULL))
cri <- options[["credibleInterval"]]
for (i in seq_len(nparam)) {
if (!all(weightedDensities[, i, 2L] == 0)) { # only TRUE if we explicitly set it to 0
weightedCRIs[i, ] <- .BANOVAapproxCRI(weightedDensities[, i, 1L], weightedDensities[, i, 2L], cri)
}
}
# compute residuals and r-squared
# sample from the joint posterior over models and parameters
tmp <- .BANOVAgetSMIResidRsq(weightedDensities, dataset, model$model.list[[nmodels]], nIter, model, options)
means <- rowMeans(tmp$resid)
h <- (1 - cri) / 2
quants <- matrixStats::rowQuantiles(tmp$resid, probs = c(h, 1 - h))
# the code above is equivalent to the code below, but the code below needs to keep all posterior samples of
# all models in memory. plug this body of this function into the call inside .BANOVAsamplePosteriors() to check
# weights <- rep(postProbs, each = nIter)
# independentVariable <- all.vars(.BANOVAgetModelFormulaFromBFobj(model$models[[2L]]))[1L]
# resids <- matrix(NA, nrow(dataset), 0)
# rsq <- vector("list", nmodels)
# # get residuals of all models individually
# for (i in seq_len(nmodels)) {
# if (is.null(model$models[[i]]$bf)) {
# tmp2 <- .BANOVAresidualsNullModel(nIter, dataset[[independentVariable]])
# } else {
# tmp2 <- .BANOVAgetSMIResidRsq(
# posterior = samples[[i]],
# dataset = dataset,
# formula = .BANOVAgetModelFormulaFromBFobj(model$models[[i]])
# )
# }
# resids <- cbind(resids, tmp2$resid)
# rsq[[i]] <- tmp2$rsq
# }
# # compute weighted mean for each row
# means2 <- tcrossprod(weights / nIter, resids)
# plot(means, means2); abline(0, 1)
# quants2 <- apply(resids, 1L, Hmisc::wtd.quantile, weights = weights, probs = c(0.025, 0.975))
# plot(quants[, 1], quants2[1, ]); abline(0, 1)
# plot(quants[, 2], quants2[2, ]); abline(0, 1)
# all information for q-q plot of residuals
weightedResidSumStats <- matrix(c(means, quants), nrow = length(means), ncol = 3L,
dimnames = list(NULL, c("mean", "cri.2.5%", "cri.97.5%")))
# all information for r-squared density plot
weightedRsqDens <- density(tmp$rsq, n = 2^11, from = 0, to = 1)
weightedRsqCri <- quantile(tmp$rsq, probs = c(h, 1 - h))
weightedRsqMean <- mean(tmp$rsq)
return(list(
weightedCRIs = weightedCRIs, weightedResidSumStats = weightedResidSumStats,
weightedRsqDens = weightedRsqDens, weightedRsqCri = weightedRsqCri, weightedRsqMean = weightedRsqMean
))
}
.BANOVAsampleNullModel <- function(nsamples, dependent) {
# sample from posterior under NULL model, needed to compute model averaged residuals.
rt.scaled <- function (n, df, mean = 0, sd = 1, ncp) {
mean + sd * stats::rt(n, df, ncp = ncp)
}
# sample from the marginal posterior distribution of the mean t-distribution, based on Murphy (2007)
n <- length(dependent)
muObs <- mean(dependent)
varObs <- var(dependent)
# uninformative priors
k0 <- 1e-6
a0 <- 1e-5
b0 <- 1e-5
an <- a0 + n / 2
bn <- b0 + 0.5 * (n - 1L) * varObs + k0 * n * muObs / (2 * (k0 + n))
kn <- k0 + n
mun <- n * muObs / kn
samples <- matrix(rt.scaled(n = nsamples, df = 2 * an, mean = mun, sd = bn / (an * kn), ncp = 0), nsamples,
1L, dimnames = list(NULL, "mu"))
return(samples)
}
.BANOVAresidualsNullModel <- function(nsamples, dependent) {
# sample from posterior under NULL model, needed to compute model averaged residuals.
rt.scaled <- function (n, df, mean = 0, sd = 1, ncp) {
mean + sd * stats::rt(n, df, ncp = ncp)
}
# sample from the marginal posterior distribution of the mean t-distribution, based on Murphy (2007)
n <- length(dependent)
muObs <- mean(dependent)
varObs <- var(dependent)
# uninformative priors
k0 <- 1e-6
a0 <- 1e-5
b0 <- 1e-5
an <- a0 + n / 2
bn <- b0 + 0.5 * (n - 1L) * varObs + k0 * n * muObs / (2 * (k0 + n))
kn <- k0 + n
mun <- n * muObs / kn
samples <- rt.scaled(n = nsamples, df = 2 * an, mean = mun, sd = bn / (an * kn), ncp = 0)
# compute all pairwise differences (residuals)
preds <- tcrossprod(rep(1, n), samples)
resids <- dependent - preds
rsq <- .BANOVAcomputeRsq(dependent, preds)
return(list(resids = resids, rsq = rsq))
}
.BANOVAfitDensity <- function(samples, gridsize = 100L, getXRange = TRUE) {
# ideally we don't do kernel density estimation but instead use some parametric approximation that is suited
# for unimodal distributions. It should also be fast, e.g., fit using method of moments (or another analytic method)
nc <- ncol(samples)
fits <- array(
data = NA,
dim = c(gridsize, nc, 2),
dimnames = list(NULL, colnames(samples), c("x", "y"))
)
n <- nrow(samples)
del0 <- (1/(4 * pi))^(1/10) # normal kernel
for (i in seq_len(nc)) {
# bandwidth definition from KernSmooth::bkde, but adjusted when sqrt(var(samples[, i])) is infinite
vx <- sqrt(var(samples[, i]))
if (is.infinite(vx)) # contWide can cause the variance to be infinite
vx <- sqrt(.Machine$double.xmax) # some big number as a fallback
bandwidth <- del0 * (243/(35 * n))^(1/5) * vx
fits[, i, ] <- do.call(cbind, KernSmooth::bkde(samples[, i], gridsize = gridsize, bandwidth = bandwidth))
}
if (getXRange) {
if (nc == 1L) {
xRanges <- matrix(range(fits[, , "x"]), 1L, 2L)
} else {
xRanges <- matrixStats::colRanges(fits[, , "x"])
}
rownames(xRanges) <- colnames(samples)
} else {
xRanges <- NULL
}
return(list(fit = fits, xRanges = xRanges))
}
.BANOVAapproxCRI <- function(x, y, cri = 0.95) {
# approximate cdf
y2 <- cumsum(y)
y2 <- y2 / y2[length(y2)]
h <- (1 - cri) / 2
# find closest observed match
idx <- c(
which.min(abs(y2 - h)),
which.min(abs(y2 - 1 + h))
)
return(x[idx])
}
# plot posteriors ----
.BANOVAgetBMAdensity <- function(samples, weights, fromTo = NULL, n = 2^10) {
# @param samples, list of samples
# @param weigths, vector of numeric weights
# @param fromTo vector of length 2 specifying lower and upper bound (optional).
# @return a list with $x the x-coordinates and $y the y-coordinates.
if (length(samples) != length(weights))
stop(gettext("length of samples must be equal to length of weights!"))
# remove NULL indices
idxNonNull <- lengths(samples) > 0
weights <- weights[idxNonNull]
samples <- samples[idxNonNull]
# renormalize the weights so output is a proper density function (TODO: do we want this, or BAS style?)
weights <- weights / sum(weights)
# create x-grid for density
if (is.null(fromTo))
fromTo <- range(sapply(samples, range))
xs <- seq(fromTo[1L], fromTo[2L], length.out = 2^10)
# compute weighted density
ys <- numeric(n)
for (i in seq_along(samples)) {
ys <- ys + weights[i] * density(samples[[i]], from = fromTo[1L], to = fromTo[2L], n = n)$y
}
return(list(x = xs, y = ys))
}
.BANOVAreSample <- function(n, x, y, prop0 = NULL) {
cdf <- cumsum(y)
max <- cdf[length(cdf)]
if (max == 0) # can't salvage this
return(numeric(n))
cdf <- cdf / max
if (!is.null(prop0) && prop0 <= (1 - sqrt(.Machine$double.eps))) {
i1 <- runif(n) <= prop0
samples <- numeric(n)
samples[i1] <- approx(cdf, x, runif(sum(i1)), rule = 2)$y
} else {
samples <- approx(cdf, x, runif(n), rule = 2)$y
}
return(samples)
}
.BANOVAgetSMIResidRsq <- function(posterior, dataset, formula, nIter, model, options) {
# @param posterior object from Bayesfactor package, or SxPx2 array of weighted densities
# @param dataset dataset
# @param formula the formula for this specific model. Supply the formula of the most complex model for BMA inference.
# @param nIter number of posterior samples
# @param model model estimates
# @param levelInfo for each variable, how many levels there are.
#
# @return a list with residuals, predictions, and r-squared
# @details the matrix multiplication in this function allocates an array of nobs * nsamples, which can be enormous.
# if is more memory efficient to use for loops (but slower in R) to calculate predictions for each observation
# and posterior sample. This could be done in the future if performance is an issue. However, this likely cannot be
# done efficiently in R. For modelAveraged inference, model and levelInfo need to be given.
# here we need the data in a one-hot encoded way
# first one is dependent variable, rest are independent variables
dvs <- all.vars(formula)[1L]
ivs <- all.vars(formula)[-1L]
# idx contains variables that need to be one-hot encoded
idx <- sapply(dataset[ivs], is.factor)
# idx <- isFactor & names(isFactor) %in% dvs
# do one-hot encoding (thank base R for calling a function `contrast` and giving it the argument `contrast`...)
datOneHot <- model.matrix(formula, data = dataset[c(dvs, ivs)],
contrasts.arg = lapply(dataset[ivs][idx], contrasts, contrasts = FALSE))
nobs <- nrow(datOneHot)
# possible switch for memory efficient alternatives:
# size of 1 numeric element in bytes
# doubleSize <- as.numeric(utils::object.size(double(2)) - utils::object.size(double(1)))
# 1073741824 is 1 GB in bytes
# 536870912 is 512 MB in bytes
# the size of an empty matrix is ignored here
# useMatmult <- ((doubleSize * nobs * nIter) / 536870912.0) <= 1
# if (useMatmult) {
if (length(dim(posterior)) == 3L) {
# we're doing model averaged inference
# recompute the levelInfo, this is relatively cheap and we cannot be sure that the levelInfo in the
# posterior was not retrieved from state (and contains too many or too few variables)
nmodels <- length(model[["model.list"]])
randomFactors <- .BANOVAgetRandomFactors(options, model[["analysisType"]])
dataTypes <- .BANOVAgetDataTypes(dataset, model[["model.list"]][[nmodels]], randomFactors)
levelInfo <- .BANOVAgetLevelInfo(dataset, model[["model.list"]][[nmodels]], dataTypes)
# which parameters are nuisance?
parameterNames <- c("mu", names(levelInfo[["levelNames"]]))
# the variables that are nuisance, e.g., contBinom
isNuisanceVar <- c(TRUE, parameterNames[-1L] %in% model[["nuisance"]])
names(isNuisanceVar) <- parameterNames
# the parameters that are nuisance, e.g., contBinom-0 and contBinom-1
isNuisanceParam <- rep(isNuisanceVar, c(1, lengths(levelInfo[["levelNames"]])))
names(isNuisanceParam) <- c("mu", unlist(levelInfo[["levelNames"]]))
if (!all(names(isNuisanceParam) %in% colnames(posterior)))
stop("!all(names(isNuisance) %in% colnames(posterior)) was not true.")
# we cannot use model[["effects"]] because that exludes nuisance parameters, so we rebuild one that doesn't exclude these
effects <- matrix(FALSE, nmodels, length(isNuisanceParam), dimnames = list(NULL, names(isNuisanceParam)))
for (i in seq_along(model[["model.list"]])) {
if (is.null(model[["model.list"]][[i]])) { # implies an intercept-only null model
cnms <- "mu"
} else {
idx <- names(levelInfo[["levelNames"]]) %in% attr(terms(model[["model.list"]][[i]]), "term.labels")
cnms <- c("mu", unlist(levelInfo[["levelNames"]][idx]))
}
effects[i, cnms] <- TRUE#c(TRUE, names(levelInfo[["levelNames"]]) %in% attr(terms(model[["model.list"]][[i]]), "term.labels"))
}
# sample from the BMA posterior
postProbs <- model[["postProbs"]]
postProbs[is.na(postProbs)] <- 0
.setSeedJASP(options)
modelIdx <- sample(length(postProbs), nIter, TRUE, postProbs) # resample the models according to posterior prob
samples <- matrix(0, nrow = nIter, ncol = ncol(datOneHot), dimnames = list(NULL, names(isNuisanceParam)))
# resample nuisance parameters
for (i in which(isNuisanceParam)) {
.setSeedJASP(options)
samples[, i] <- .BANOVAreSample(n = nIter, x = posterior[, i, "x"], y = posterior[, i, "y"])
}
# resample non nuisance parameters
for (i in which(!isNuisanceParam)) {
idx <- names(isNuisanceParam)[i]
mult <- effects[modelIdx, idx] # how often models with this variable are sampled
nTemp <- sum(mult)
if (nTemp > 0L){
.setSeedJASP(options)
samples[mult, i] <- .BANOVAreSample(n = sum(mult), x = posterior[, i, "x"], y = posterior[, i, "y"])
}
}
preds <- tcrossprod(datOneHot, samples)
} else {
# TODO: ensure that column order is always correct for both!
# otherwise we're doing inference for a single model
preds <- tcrossprod(datOneHot, posterior)
}
# calculate residuals (correctly recycles dat[[dvs]])
resid <- dataset[[dvs]] - preds
rsq <- .BANOVAcomputeRsq(dataset[[dvs]], preds)
# eta <- .BANOVAcomputeEtasq(dat[[dvs]], preds)
# } else {
# slow but memory efficient stuff
# }
return(list(resid = resid, preds = preds, rsq = rsq))
}
# HF computation ----
.BANOVAinitBayesFactor <- function() {
defaults <- list(
BFMaxModels = 50000,
BFpretestIterations = 100,
BFapproxOptimizer = "optim",
BFapproxLimits = c(-15, 15),
BFprogress = interactive(),
BFfactorsMax = 5
)
idx <- setdiff(names(defaults), names(options()))
options(defaults[idx])
}
.BANOVAcomputeRsq <- function(obs, preds) {
# NOTE: R^2 != cor(obs, predict) because the predictions from the posterior samples are not OLS estimates.
if (is.null(dim(preds)))
preds <- matrix(preds, ncol = 1L)
# definition from http://www.stat.columbia.edu/~gelman/research/unpublished/bayes_R2.pdf
colVars2 <- function(x) {matrixStats::colVars(sweep(x, 2L, matrixStats::colMeans2(x)))}
ee <- colVars2(preds)
ff <- colVars2(obs - preds)
gg <- ee / (ee + ff)
if (anyNA(gg)) {
idx <- is.na(gg)
if (all(idx)) {
# this means that either the posterior sampling failed or the input data is bogus.
.quitAnalysis(gettext("All predictions had zero variance. Check the predictors/DV for anomalies or increase the number of posterior samples."))
} else {
# this shouldn't happen, but is necessary because otherwise stuff will crash downstream.
gg[is.na(gg)] <- mean(gg, na.rm = TRUE)
}
}
return(gg)
}
.BANOVAcomputeEtasq <- function(obs, preds) {
# partial eta^2
if (is.null(dim(preds)))
preds <- matrix(preds, ncol = 1L)
# return(1 - matrixStats::colVars(preds - obs) / var(obs))
# return(1 - matrixStats::rowVars(preds - obs) / var(obs))
return(matrixStats::colVars(preds) / var(obs))
}
.BANOVAmakeChainNeater <- function(chains, Xnames, formula, data, dataTypes, gMap, unreduce, continuous, columnFilter) {
# so BayesFactor:::makeChainNeater doesn't always make things neater
# see also https://github.com/jasp-stats/jaspAnova/pull/43#discussion_r656430090
# example where BayesFactor goes wrong:
# dd <- structure(list(
# rm_factor = structure(c(1L, 2L, 1L, 2L, 1L, 2L), .Label = c("Level 1", "Level 2"), class = "factor"),
# dependent = c(6, 6, 5, 10, 6, 9),
# age = c(59.2161793866222, 59.2161793866222, 53.4457078686802, 53.4457078686802, 50.3248450015473, 50.3248450015473),
# subject = structure(c(1L, 1L, 2L, 2L, 3L, 3L), .Label = c("1", "2", "3"), class = "factor")),
# row.names = c(NA, 6L), class = "data.frame"
# )
#
# bf_obj_1 <- BayesFactor::lmBF(dependent ~ age + subject, data = dd, whichRandom = "subject")
# bf_obj_2 <- BayesFactor::lmBF(dependent ~ age + subject + rm_factor, data = dd, whichRandom = "subject")
#
# print(BayesFactor::posterior(bf_obj_1, iterations = 3)) # <- ugly
# print(BayesFactor::posterior(bf_obj_2, iterations = 3)) # <- good
#
# originalFun <- BayesFactor:::makeChainNeater
# jaspBase:::assignFunctionInPackage(.BANOVAmakeChainNeater, "makeChainNeater", "BayesFactor")
# print(BayesFactor::posterior(bf_obj_1, iterations = 3)) # <- good
#
# jaspBase:::assignFunctionInPackage(originalFun, "makeChainNeater", "BayesFactor")
# print(BayesFactor::posterior(bf_obj_1, iterations = 3)) # <- ugly
# this part is identical to BayesFactor:::makeChainNeater
P = length(gMap)
nGs = max(gMap) + 1
factors = BayesFactor:::fmlaFactors(formula, data)[-1]
dataTypes = dataTypes[names(dataTypes) %in% factors]
types = BayesFactor:::termTypes(formula, data, dataTypes)
lastPars = ncol(chains) + (-nGs):0
if (any(continuous)) {
gNames = paste("g", c(names(types[types != "continuous"]), "continuous"), sep = "_")
} else {
gNames = paste("g", names(types), sep = "_")
}
if (is.null(columnFilter)) {
ignoreCols = ignoreCols = rep(0, P)
} else {
ignoreCols = BayesFactor:::filterVectorLogical(columnFilter, names(gMap))
}
# here we're also checking here if any dataTypes are "random" instead of only checking for fixed dataTypes
if (!unreduce | !any(dataTypes %in% c("fixed", "random"))) {
labels = c("mu", Xnames[!ignoreCols], "sig2", gNames)
colnames(chains) = labels
return(chains)
}
parLabels = BayesFactor:::makeLabelList(formula, data, dataTypes, unreduce,
columnFilter)
labels = c("mu", parLabels)
betaChains = chains[, 1:(ncol(chains) - 2 - nGs) + 1, drop = FALSE]
betaChains = BayesFactor:::ureduceChains(betaChains, formula, data, dataTypes, gMap, ignoreCols)
newChains = cbind(chains[, 1], betaChains, chains[, lastPars])
labels = c(labels, "sig2", gNames)
colnames(newChains) = labels
return(newChains)
}
.BANOVAcomputeInteractionsMatrix <- function(effects) {
neffects <- length(effects)
interactions.matrix <- matrix(FALSE, nrow = neffects, ncol = neffects)
rownames(interactions.matrix) <- colnames(interactions.matrix) <- effects
if (neffects > 1L) {
effect.components <- sapply(effects, strsplit, split = ":", fixed = TRUE)
for (e in seq_len(neffects)) {
interactions.matrix[e, ] <- sapply(1:neffects, function(ee) {
(sum(effect.components[[e]] %in% effect.components[[ee]]) == length(effect.components[[e]]))
})
}
diag(interactions.matrix) <- FALSE
}
return(interactions.matrix)
}
# Model prior ----
.BANOVAcomputePriorModelProbs <- function(models, nuisance, options) {
if (options[["modelPrior"]] == "uniform") {
modelprobs <- rep(1 / length(models), length(models))
} else if (options[["modelPrior"]] == "custom") {
inclusionProbabilities <- vapply(options[["customPriorSpecification"]], `[[`, FUN.VALUE = numeric(1L), "inclusionProbability")
modelprobs <- .BANOVAcustomInclusionProbabilitiesToModelProbabilities(models, nuisance, inclusionProbabilities, enforceMarginality = options[["enforcePrincipleOfMarginalityFixedEffects"]])
} else {
noPredictorsPerModel <- vapply(models, function(x) if (is.null(x)) 0L else length(attr(terms(x), "term.labels")), FUN.VALUE = integer(1L))
noNuisancePredictors <- noPredictorsPerModel[1L]
noPredictorsPerModel <- noPredictorsPerModel - noNuisancePredictors
totalNoPredictors <- noPredictorsPerModel[length(noPredictorsPerModel)]
if (options[["modelPrior"]] %in% c("betaBinomial", "Wilson", "Castillo")) {
switch (options[["modelPrior"]],
"betaBinomial" = {alpha = options[["betaBinomialParameterA"]]; beta = options[["betaBinomialParameterB"]] },
"Wilson" = {alpha = 1.0; beta = totalNoPredictors * options[["wilsonParameterLambda"]] },
"Castillo" = {alpha = 1.0; beta = totalNoPredictors ^ options[["castilloParameterU"]] }
)
modelprobs <- dBetaBinomialModelPrior(noPredictorsPerModel, totalNoPredictors, alpha, beta)
} else if (options[["modelPrior"]] == "Bernoulli") {
modelprobs <- dBernoulliModelPrior(noPredictorsPerModel, totalNoPredictors, options[["bernoulliParameter"]])
}
# both the betabinomial and bernoulli model priors do not sum to 1 when marginality is respected
modelprobs <- modelprobs / sum(modelprobs)
}
return(modelprobs)
}
dBetaBinomialModelPrior <- function(k, n, alpha = 1.0, beta = 1.0, log = FALSE) {
# NOTE: this is not the pdf of the betabinomial, but the pdf divided by choose(n, k), so that it's the probabiltiy of a particular model
logprobability <- lbeta(k + alpha, n - k + beta) - lbeta(alpha, beta)
if (log)
return(logprobability)
return(exp(logprobability))
}
dBernoulliModelPrior <- function(k, n, prob = 0.5, log = FALSE) {
logprobability <- k * log(prob) + (n - k) * log(1 - prob)
if (log)
return(logprobability)
return(exp(logprobability))
}
.BANOVAcustomInclusionProbabilitiesToModelProbabilities <- function(modelList, nuisance, inclusionProbabilities, enforceMarginality = TRUE) {
# TODO: discuss with the team:
# when marginality is enforced, how should the prior on interaction effects be interpreted when a user wants matched models?
formulaFullModel <- modelList[[length(modelList)]]
termsFullModel <- terms(formulaFullModel)
termLabels <- attr(termsFullModel, "term.labels")
fullNonNuisance <- !(termLabels %in% nuisance)
termLabels <- termLabels[fullNonNuisance]
# the principle of marginality is enforced (or not) through the function getPresentInteractions
if (enforceMarginality) {
# in this case, interaction effects do not appear without the main effects, so the user provided probability
# is interpreted as a conditional probability given that the main effects are included.
# thus, p(interaction) = 0 if the main effects are missing, and the user provided probability otherwise.
#
# for example, consider p(A) = P(B) = .5, P(A:B) = .2
#
# model P(model)
# A P(A) * (1 - P(B)) # <- here we do NOT do '* (1 - P(A:B))'
# A, B p(A) * P(B) * (1 - P(A:B))
# A, B, A:B p(A) * P(B) * P(A:B)
termFactors <- attr(termsFullModel, "factors")[, termLabels, drop = FALSE]
termFactors <- termFactors[rownames(termFactors) %in% termLabels, , drop = FALSE]
rnmsTermFactors <- rownames(termFactors)
termOrders <- attr(termsFullModel, "order")[fullNonNuisance]
# create a list where the names refer to the interactions and the elements refer to all children of that interaction
listOfDescendants <- setNames(lapply(termLabels, function(x) {
idx <- which(termFactors[, x] == 1L)
rnmsTermFactors[idx]
}), termLabels)
termIsNotInteraction <- termOrders == 1L
getPresentInteractions <- function(currentTermLabels) {
vapply(listOfDescendants, function(descendants, currentTermLabels) {
all(descendants %in% currentTermLabels)
}, logical(1L), currentTermLabels) | termIsNotInteraction
}
# NOTE: the loop below fails for some combinations of user priors, but I can see how this behavior is desirable
# the users provide the raw p(interaction) however, the code below assumes that this is
# p(interaction | main effects), so to compensate we divide each interaction by the product of the components.
# however, for some inclusion probabilities of the main effects you get impossible values, for example:
# p(A_) = 0.1, p(B) = 0.1, P(A*B) = 0.9 => 0.9 / 0.1^2 = 90, which is bigger than 1 and should be impossible.
# for (i in which(!termIsNotInteraction)) { # loop over all interaction terms
# idxSubterms <- which(termLabels %in% listOfDescendants[[termLabels[i]]])
# inclusionProbabilities[i] <- inclusionProbabilities[i] / prod(inclusionProbabilities[idxSubterms])
# }
} else {
# in this case, interaction effects do appear without the main effects, so the user provided probability
# is interpreted in the same way for interaction effects and fixed effects.
# it's interpreted as an unconditional inclusion probability
getPresentInteractions <- function(...) {
rep(TRUE, length(termLabels))
}
}
exclusionProbabilties <- 1 - inclusionProbabilities
# so terms(y ~ b + a:b) reorders a:b to b:a but terms(y ~ a + a:b) does not...
termLabelOrder <- strsplit(termLabels, ":", fixed = TRUE)
modelProbs <- numeric(length(modelList))
for (i in seq_along(modelList)) {
if (is.null(modelList[[i]])) {
excludedTerms <- rep(TRUE, length(termLabels))
excludedTerms <- excludedTerms & getPresentInteractions(NULL)
modelProbs[i] <- prod(exclusionProbabilties)
} else {
currentTerms <- terms(modelList[[i]])
currentTermLabels <- setdiff(labels(currentTerms), nuisance)
currentTermLabels <- .BANOVAorderTermsByKnownOrder(currentTerms, currentTermLabels, termLabelOrder)
presentInteractions <- getPresentInteractions(currentTermLabels)
includedTerms <- termLabels %in% currentTermLabels
excludedTerms <- !includedTerms
includedTerms <- includedTerms & presentInteractions
excludedTerms <- excludedTerms & presentInteractions
modelProbs[i] <- prod(inclusionProbabilities[includedTerms], exclusionProbabilties[excludedTerms])
}
}
return(modelProbs)
}
# Other ----
.BANOVAgetRScale <- function(options, analysisType) {
if (options[["priorSpecificationMode"]] == "acrossParameters") {
rscaleFixed <- options[["cauchyPriorScaleFixedEffects"]]
rscaleRandom <- options[["cauchyPriorScaleRandomEffects"]]
rscaleCont <- if (analysisType == "ANOVA") "medium" else options[["cauchyPriorScaleCovariates"]]
rscaleEffects <- NULL
} else {
# NOTE: this still need a value otherwise BayesFactor returns NA for the bf
# default values of lmBF
rscaleFixed <- "medium"
rscaleRandom <- "nuisance"
rscaleCont <- if (analysisType == "ANOVA") "medium" else options[["cauchyPriorScaleCovariates"]]
rscaleEffectsNames <- vapply(options[["customPriorSpecification"]], FUN.VALUE = character(1L), function(x) {
paste(x[["components"]], collapse = ":")
})
rscaleEffects <- vapply(options[["customPriorSpecification"]], FUN.VALUE = numeric(1L), `[[`, "scaleFixedEffects")
rscaleEffectsKeep <- if (analysisType == "ANOVA") {
rep(TRUE, length(rscaleEffects))
} else {
vapply(options[["customPriorSpecification"]], FUN.VALUE = logical(1L), function(x) {
is.null(options[["covariates"]]) || !any(x[["components"]] %in% options[["covariates"]])
})
}
rscaleEffects <- rscaleEffects[rscaleEffectsKeep]
names(rscaleEffects) <- rscaleEffectsNames[rscaleEffectsKeep]
}
return(list(rscaleFixed = rscaleFixed, rscaleRandom = rscaleRandom, rscaleCont = rscaleCont, rscaleEffects = rscaleEffects))
}
.BANOVAgetShownTableSizeAndAddFootnote <- function(modelTable, nmodels, modelsShown, numModelsShown) {
if (modelsShown == "unlimited") {
shownTableSize <- nmodels
} else {
shownTableSize <- min(nmodels, numModelsShown)
if (shownTableSize < nmodels)
modelTable$addFootnote(gettextf("Showing the best %1$d out of %2$d models.", numModelsShown, nmodels))
}
return(shownTableSize)
}
# # .BANOVAupdateRscales <- function() {
#
# # for some reason, BayesFactor supports custom rscales for but throws a error if you try to set them for continuous parameters
# # this temporarily overwrite the stop call in this BayesFactor:::createRscales
# # the function will throw an error if BayesFactor:::createRscales is updated
#
# # originalCreateRscales <- BayesFactor:::createRscales
# # originalBody <- body(originalCreateRscales)
# # if (!identical(originalBody[[17]][[3]][[2]][[1]], quote(invisible))) {
# # if (!identical(originalBody[[17]][[3]][[2]][[1]], quote(stop)))
# # stop("Bayesfactor got an update, and createRscales/ .BANOVAupdateRscales needs to be fixed!", domain = NA)
# # originalBody[[17L]][[3L]][[2L]][[1L]] <- substitute(invisible)
# # newCreateRscales <- originalCreateRscales
# # body(newCreateRscales) <- originalBody
# #
# # jaspBase::assignFunctionInPackage(newCreateRscales, "createRscales", "BayesFactor")
# # }
# #
# # tmp <- methods::getMethod(BayesFactor::compare, list(numerator="BFlinearModel", denominator="missing", data="data.frame"))
# # fun <- tmp@.Data
# # newfun <- fun
# # body <- body(fun)
# # if (!identical(body[[19L]][[2L]][[2L]][[4L]][[2L]], quote(FALSE))) {
# # if (!identical(body[[19L]][[2L]][[2L]][[4L]][[2L]], quote(all(relevantDataTypes == "continuous"))))
# # stop("Bayesfactor got an update, and BayesFactor::`.__T__compare:BayesFactor`$`BFlinearModel#missing#data.frame`/ .BANOVAupdateRscales needs to be fixed!", domain = NA)
# # body[[19L]][[2L]][[2L]][[4L]][[2L]] <- substitute(FALSE) # do not take the fast path
# # body(newfun) <- body
# #
# # # horrible but necessary
# # env <- environment()#getNamespace(jaspAnova)
# # assign("compare", BayesFactor::compare, envir = env)
# # setMethod("compare", signature(numerator = "BFlinearModel", denominator = "missing", data = "data.frame"), newfun, where = env)
#
# # env <- getNamespace("BayesFactor")
# # unlockBinding("compare", env)
#
# # debugonce(setMethod)
# # setMethod("compare", signature(numerator = "BFlinearModel", denominator = "missing", data = "data.frame"), newfun, where = env)
#
# # tmp@.Data <- newfun
# # s4env <- BayesFactor:::`.__T__compare:BayesFactor`
# # unlockBinding("BFlinearModel#missing#data.frame", s4env)
# # s4env$`BFlinearModel#missing#data.frame` <- tmp
# # lockBinding("BFlinearModel#missing#data.frame", s4env)
# # jaspBase::assignFunctionInPackage(s4env, ".__T__compare:BayesFactor", "BayesFactor")
# # }
#
# # return()
#
# # we could reset the changes with on.exit, but I'm not sure if it's necessary and doing that cleanly requires withr::defer and adds a dependency
# # do.call(on.exit({
# # jaspBase::assignFunctionInPackage(tempRscales[["originalCreateRscales"]], "createRscales", "BayesFactor")
# #
# # s4env <- BayesFactor:::`.__T__compare:BayesFactor`
# # unlockBinding("BFlinearModel#missing#data.frame", s4env)
# # s4env$`BFlinearModel#missing#data.frame` <- newfun
# # lockBinding("BFlinearModel#missing#data.frame", s4env)
# # jaspBase::assignFunctionInPackage(s4env, ".__T__compare:BayesFactor", "BayesFactor")
# # }, add = TRUE, after = FALSE)
# }
# Dependencies ----
.BANOVAdataDependencies <- function() {
# seed is in here because basically everything depends on the seed
c("dependent", "randomFactors", "covariates", "fixedFactors", "betweenSubjectFactors", "repeatedMeasuresFactors", "repeatedMeasuresCells", "modelTerms",
"seed", "setSeed")
}
.BANOVAmodelSpaceDependencies <- function(modelPrior) {
c("modelPrior", "betaBinomialParameterA", "betaBinomialParameterB", "wilsonParameterLambda", "castilloParamerU", "enforcePrincipleOfMarginalityFixedEffects", "enforcePrincipleOfMarginalityRandomSlopes",
if (modelPrior == "custom") "customPriorSpecification"
)
}
.BANOVArscaleDependencies <- function(priorSpecificationMode) {
c("cauchyPriorScaleFixedEffects", "cauchyPriorScaleRandomEffects", "cauchyPriorScaleCovariates", "priorSpecificationMode",
if (priorSpecificationMode == "perTerm") "customPriorSpecification"
)
}
.BANOVAmodelPriorOptionsChanged <- function(state, options) {
!identical(state[["modelPriorOptions"]][.BANOVAmodelSpaceDependencies(options[["modelPrior"]])], options[.BANOVAmodelSpaceDependencies(options[["modelPrior"]])])
}
.BANOVAmarginalityOptionsUnchanged <- function(state, options) {
stateOpts <- state[["modelPriorOptions"]]
stateOpts[["enforcePrincipleOfMarginalityFixedEffects"]] == options[["enforcePrincipleOfMarginalityFixedEffects"]] &&
stateOpts[["enforcePrincipleOfMarginalityRandomSlopes"]] == options[["enforcePrincipleOfMarginalityRandomSlopes"]]
}
.BANOVAmodelBFTypeOrOrderChanged <- function(state, options) {
nms <- c("fixedFactors", "modelTerms", "randomFactors", "covariates", "seed", "setSeed")
nms <- intersect(names(options), nms) # excludes covariates for ANOVA
identical(state[nms], options[nms])
}
.BANOVAmodelTermsUnchanged <- function(state, options) {
identical(state[["modelTerms"]], options[["modelTerms"]])
}
# HF formulas ----
.BANOVAgetFormulaComponents <- function(x, what = c("components", "variables")) {
what <- match.arg(what)
if (what == "components") {
return(colnames(attr(terms(x), "factors")))
} else {
return(all.vars(x)[-1L])
}
}
.BANOVAgetModelTitlesWithAllTerms <- function(modelObjects, modelList, analysisType, hideNuisance) {
if (hideNuisance)
return(vapply(modelObjects, FUN = `[[`, FUN.VALUE = character(1L), "title"))
nullModelName <- gettext("Null model")
res <- gsub("(.*)~\\s+", "", vapply(modelList, function(x) if (is.null(x)) nullModelName else .BANOVAas.character.formula(x), character(1L)))
if (analysisType == "RM-ANOVA")
res <- gsub(.BANOVAsubjectName, "subject", res)
res[res == ""] <- nullModelName
return(jaspBase::gsubInteractionSymbol(res))
}
.BANOVAgenerateAllModelFormulas <- function(formula, nuisance = NULL, analysisType = "RM-ANOVA",
enforcePrincipleOfMarginalityFixedEffects = TRUE,
enforcePrincipleOfMarginalityRandomSlopes = FALSE,
rmFactors = NULL, legacy = FALSE
) {
# TODO: it might be nicer to represent the NULL model as
# y ~ 1, rather than NULL. Right now the null model is a formula
# when there are nuisance variables and otherwise NULL, which leads to a bunch
# of annoying if (is.null(model[i])) exceptions.
neverExclude <- paste0("^", nuisance, "$")
if (!legacy && analysisType == "RM-ANOVA" && is.null(rmFactors))
stop(".BANOVAgenerateAllModelFormulas called with invalid arguments: analysisType = \"RM-ANOVA\", legacy = FALSE, rmFactors = NULL", domain = NA)
modelSpace <- try(BayesFactor::enumerateGeneralModels(
formula,
whichModels = if (enforcePrincipleOfMarginalityFixedEffects) "withmain" else "all",
neverExclude = neverExclude)
)
nuisanceRandomSlopes <- NULL
if (isTryError(modelSpace))
.quitAnalysis(gettextf("The following error occured in BayesFactor::enumerateGeneralModels: %s",
.extractErrorMessage(modelSpace)))
if (analysisType == "RM-ANOVA" && !legacy) {
# adjust the nuisance to include the random slopes, only done here because BayesFactor::enumerateGeneralModels
# does not handle this properly
allPossibleSlopes <- labels(stats::terms(stats::as.formula(paste("y~", paste(rmFactors, collapse = "*")))))
# drop the most complex interaction
allPossibleSlopes <- allPossibleSlopes[-length(allPossibleSlopes)]
# if at least one interaction among the repeated measures is considered
if (length(allPossibleSlopes) > 0L) {
# add interaction with subject
nuisanceRandomSlopes <- paste0(allPossibleSlopes, ":", .BANOVAsubjectName)
if (enforcePrincipleOfMarginalityRandomSlopes) {
for (i in seq_along(modelSpace)) {
presentLabels <- labels(stats::terms(modelSpace[[i]]))
termsToAddChar <- intersect(allPossibleSlopes, presentLabels)
if (length(termsToAddChar) > 0L) {
termsToAddChar <- paste0(termsToAddChar, ":", .BANOVAsubjectName)
termsToAdd <- as.formula(paste("~ . +", paste0(termsToAddChar, collapse = " + ")))
modelSpace[[i]] <- update.formula(modelSpace[[i]], new = termsToAdd)
}
}
} else {
termsToAdd <- as.formula(paste("~ . +", paste0(nuisanceRandomSlopes, collapse = " + ")))
modelSpace <- lapply(modelSpace, update.formula, new = termsToAdd)
}
# for the reordering done below. termsToAdd always contains the most complex random effects
formula <- update.formula(formula, new = termsToAdd)
# update nuisance
nuisance <- c(nuisance, nuisanceRandomSlopes)
}
}
if (is.null(nuisance)) {
return(list(modelList = c(list(NULL), modelSpace), nuisance = nuisance, nuisanceRandomSlopes = nuisanceRandomSlopes))
} else {
# put the null-model first
i <- length(modelSpace)
modelSpace <- c(modelSpace[[i]], modelSpace[-i])
# BayesFactor has the nasty habit of changing the order of interactions whenever one of the components is
# a nuisance variable. Here we ensure the order matches that of the input formula (which matches the order a user entered the factors).
# see also https://github.com/jasp-stats/jasp-test-release/issues/1181
originalTerms <- terms(formula)
originalOrd <- attr(originalTerms, "order")
originalIdxNonInteraction <- which(originalOrd == 1L)
originalLabels <- attr(originalTerms, "term.labels")
originalLabelsPieces <- strsplit(originalLabels[originalOrd > 1L], ":")
originalLabelsPiecesSorted <- lapply(originalLabelsPieces, sort)
dependent <- all.vars(modelSpace[[1]])[1L]
for (i in seq_along(modelSpace)) {
term <- terms(modelSpace[[i]])
ord <- attr(term, "order") # interaction order
idxNonInteraction <- which(ord == 1L)
termLabels <- attr(term, "term.labels")
termLabels[idxNonInteraction] <- originalLabels[originalIdxNonInteraction][match(termLabels[idxNonInteraction], originalLabels[originalIdxNonInteraction])]
for (j in which(ord > 1L)) {
labelPieces <- strsplit(termLabels[j], ":")[[1L]]
sortedLabelPieces <- sort(labelPieces)
for (k in seq_along(originalLabelsPieces))
if (identical(sortedLabelPieces, originalLabelsPiecesSorted[[k]]))
termLabels[j] <- paste(originalLabelsPieces[[k]], collapse = ":")
}
modelSpace[[i]] <- as.formula(paste(dependent, "~", paste(termLabels, collapse = " + ")))
}
return(list(modelList = modelSpace, nuisance = nuisance, nuisanceRandomSlopes = nuisanceRandomSlopes))
}
}
.BANOVAcreateModelFormula <- function(dependent, modelTerms) {
rhs <- "" # right hand side of the formula
nuisance <- NULL
effects <- NULL
for (term in modelTerms) {
comp <- term[["components"]]
newValue <- paste(comp, collapse = ":")
rhs <- if (rhs != "") paste0(rhs, " + ", newValue) else newValue
if (!is.null(term[["isNuisance"]]) && term[["isNuisance"]]) {
nuisance <- c(nuisance, newValue)
} else {
effects <- c(effects, newValue)
}
}
model.formula <- formula(paste(dependent, "~", rhs))
# this would be cleaner ideal if BayesFactor::enumerateGeneralModels would handle the nuisance properly.
# if (isRMANOVA && !legacy) {
# randomSlopes <- paste0(rmFactors, ":", .BANOVAsubjectName)
# termsToAdd <- as.formula(paste("~ . +", paste0(randomSlopes, collapse = " + ")))
# model.formula <- update.formula(model.formula, termsToAdd)
# nuisance <- c(nuisance, randomSlopes)
# }
return(list(model.formula = model.formula, nuisance = nuisance, effects = effects))
}
.BANOVAgetModelFormulaFromBFobj <- function(BayesFactorObj, asCharacter = FALSE) {
out <- BayesFactorObj$bf@numerator[[1L]]@identifier$formula
if (asCharacter) {
return(out)
} else {
return(as.formula(out))
}
}
.BANOVAgetLevelInfo <- function(dataset, formula, dataTypes) {
levelNames <- .BANOVAmakeLabelList(formula, dataset, dataTypes)
levelCounts <- lengths(levelNames) # counts of the factors including interaction terms
return(list(levelCounts = levelCounts, levelNames = levelNames))
}
.BANOVAmakeLabelList <- function(formula, data, dataTypes, unreduce = TRUE, columnFilter = NULL) {
# from BayesFactor 0.9.12.4.2
# identical to BayesFactor:::makeLabelList, except that we return a named list
# with
# names: the variables that make up a particular (interaction) effect
# values: the parameter names as returned by BayesFactor::posterior(...)
terms = attr(terms(formula, data = data), "term.labels")
if (!is.null(columnFilter))
terms = terms[!BayesFactor:::filterVectorLogical(columnFilter, terms)]
if (unreduce)
dataTypes[dataTypes == "fixed"] = "random"
labelList = lapply(terms, function(term, data, dataTypes) {
effects = strsplit(term, ":", fixed = TRUE)[[1]]
my.names = BayesFactor:::design.names.intList(effects, data, dataTypes)
return(paste(term, "-", my.names, sep = ""))
}, data = data, dataTypes = dataTypes)
# this part is different from BayesFactor
names(labelList) <- terms
return(labelList)
}
.BANOVAgetDataTypes <- function(dataset, formula, whichRandom = NULL) {
# from BayesFactor:::lmBF
BayesFactor:::createDataTypes(
formula,
whichRandom = whichRandom,
data = dataset,
analysis = "lm"
)
}
.BANOVAgetRandomFactors <- function(options, analysisType) {
if (analysisType == "RM-ANOVA")
return(.BANOVAsubjectName)
else
return(unlist(options[["randomFactors"]]))
}
.BANOVAgetLevelsFromParamNames <- function(names) {
# NOTE: this works because base64 does not contain "-"
# split on first "-"; 2 implies output of length 2, i.e., only split once
out <- do.call(rbind, stringr::str_split(names, "-", 2L))
# continous variables have as level name the variable name
idx <- out[, 1L] == out[, 2L]
out[idx, 2L] <- ""
# change dots into spaces for aesthetic purposes
out[, 2L] <- gsub(".", " ", out[, 2L], fixed = TRUE)
colnames(out) <- c("parameter", "level")
return(out)
}
.BANOVAas.character.formula <- function(x, ...) {
# we could also extend the S3 function as.character
Reduce(paste, trimws(deparse(x)))
}
.BANOVAreorderFormulas <- function(x) {
# This function reorders the terms of a formula such that they are alphabetical
# e.g., a ~ c + b + c:b becomes a ~ b + b:c
# This is necessary because BayesFactor::enumerateGeneralModels always appends the nuisance terms
# and a ~ b + c != a ~ c + b
# so without this function the state does not get reused whenever a user modifies the nuisance terms
if (is.null(x))
return("NULL") # as a string so it can become one big character vector for match()
s <- strsplit(attr(stats::terms.formula(x), "term.labels"), ":")
for (i in which(lengths(s) > 1L))
s[[i]] <- paste0(sort(s[[i]]), collapse = ":")
return(paste(all.vars(x)[1L], "~", paste(sort(unlist(s)), collapse = " + ")))
}
.BANOVAreorderTerms <- function(x) {
# This function reorders the terms of a formula such that they are alphabetical
# it essentially does the same as .BANOVAreorderFormulas but expects x to be character
# x should be a character string of the form c("a" "a:b"), so not a ~ b + c:d.
# For example, output from `all.vars(formula)` or `.BANOVAgetFormulaComponents(formula)`.
s <- strsplit(x, ":")
for (i in which(lengths(s) > 1L))
s[[i]] <- paste0(sort(s[[i]]), collapse = ":")
return(unlist(s))
}
.BANOVAorderTermsByKnownOrder <- function(currentTerms, currentTermLabels, termLabelOrder) {
# This function reorders the terms of a formula such that they follow the order of termLabelOrder
# currentTerms: output of stats::terms(formula)
# currentTermLabels: output of labels(stats::formula(formula)), possibly after filtering out nuisance
# termLabelOrder: character vector of term labels in the desired order
currentFactors <- attr(currentTerms, "factors")
if (any(currentFactors > 1L)) { # TRUE implies may need to reorder some terms according to termLabelOrder
for (j in grep(":", currentTermLabels, fixed = TRUE)) {
tmp <- strsplit(currentTermLabels[j], ":", fixed = TRUE)[[1L]]
for (order in termLabelOrder) {
if (all(tmp %in% order) && length(tmp) == length(order)) {
currentTermLabels[j] <- paste(order, collapse = ":")
break
}
}
}
}
return(currentTermLabels)
}
# Single Model Inference (SMI) ----
.BANOVAsmi <- function(jaspResults, dataset, options, model) {
userWantsSMI <- any(unlist(options[c(
"singleModelPosteriorPlot", "singleModelQqPlot", "singleModelRsqPlot", "singleModelEstimates", "singleModelCriTable"
)]))
if (!userWantsSMI)
return()
if (!is.null(jaspResults[["containerSingleModel"]])) {
singleModelContainer <- jaspResults[["containerSingleModel"]]
} else {
singleModelContainer <- createJaspContainer(title = gettext("Single Model Inference"))
singleModelContainer$dependOn(c(
"singleModelTerms", "dependent", "samplingMethodMCMC", "samplesMCMC",
"repeatedMeasuresCells", "seed", "setSeed",
.BANOVArscaleDependencies(options[["priorSpecificationMode"]])
))
jaspResults[["containerSingleModel"]] <- singleModelContainer
singleModelContainer$position <- 7
}
singleModel <- jaspResults[["singleModelState"]]$object
if (!is.null(model[["models"]]) && is.null(singleModel) && length(options$singleModelTerms) > 0L && userWantsSMI) {
singleModel <- try(.BANOVAsmiSamplePosterior(dataset, options, model[["analysisType"]]))
if (isTryError(singleModel)) {
singleModelContainer$setError(gettextf("Error in single model inference:\n%s", singleModel))
singleModel <- NULL
} else {
singleModelState <- createJaspState(object = singleModel)
singleModelState$dependOn(optionsFromObject = singleModelContainer)
jaspResults[["singleModelState"]] <- singleModelState
}
}
.BANOVAsmiEstimates(singleModelContainer, options, singleModel)
.BANOVAsmiRsquared (singleModelContainer, options, singleModel)
.BANOVAsmiQqPlot (singleModelContainer, options, singleModel)
.BANOVAsmiRsqPlot (singleModelContainer, options, singleModel)
.BANOVAsmiPosteriorPlot(singleModelContainer, dataset, options, singleModel)
return()
}
.BANOVAsmiSamplePosterior <- function(dataset, options, analysisType) {
nIter <- if (options[["samplingMethodMCMC"]] == "auto") 1e3L else options[["samplesMCMC"]]
modelTerms <- options$singleModelTerms
dependent <- options$dependent
if (analysisType == "RM-ANOVA") {
modelTerms[[length(modelTerms) + 1L]] <- list(components = .BANOVAsubjectName, isNuisance = TRUE)
dependent <- .BANOVAdependentName
}
tempRScale <- .BANOVAgetRScale(options, analysisType)
rscaleFixed <- tempRScale[["rscaleFixed"]]
rscaleRandom <- tempRScale[["rscaleRandom"]]
rscaleCont <- tempRScale[["rscaleCont"]]
rscaleEffects <- tempRScale[["rscaleEffects"]]
formula <- .BANOVAcreateModelFormula(dependent, modelTerms)$model.formula
randomFactors <- .BANOVAgetRandomFactors(options, analysisType)
dataTypes <- .BANOVAgetDataTypes(dataset, formula, randomFactors)
levelInfo <- .BANOVAgetLevelInfo(dataset, formula, dataTypes)
# if all variables are continuous, then continuous variables keep their name. E.g., "contGamma" stays "contGamma".
# if any variables are factors, then continuous variables names are doubled. E.g., "contGamma" becomes "contGamma-contGamma".
allParamNames <- c(
"mu",
if (all(dataTypes == "continuous")) names(levelInfo$levelNames) else unlist(levelInfo$levelNames)
)
.setSeedJASP(options)
samples <- BayesFactor::lmBF(
formula = formula,
data = dataset,
whichRandom = unlist(randomFactors),
progress = TRUE,
posterior = TRUE,
rscaleFixed = rscaleFixed,
rscaleRandom = rscaleRandom,
rscaleCont = rscaleCont,
rscaleEffects = rscaleEffects,
iterations = nIter
)
types <- samples@model@dataTypes
# keep only relevant columns, drop sig2, g_xxx, ...
idx <- match(allParamNames, colnames(samples), nomatch = 0L)
samples <- samples[, idx, drop = FALSE]
# for some odd reason, Bayesfactor uses as column name contcor1-contcor1
# if there is a covariate AND fixed factors, but only contcor1 if all variables are continuous...
allContinuous <- all(types == "continuous")
if (allContinuous) {
cnms <- colnames(samples)[-1L] # omit the intercept (mu) which is not changed by Bayesfactor
colnames(samples)[-1L] <- paste0(cnms, "-", cnms)
}
# find out which parameters are random
isRandom <- logical(ncol(samples))
idx <- which(types == "random")
for (i in idx)
isRandom <- isRandom | startsWith(colnames(samples), names(types)[i])
means <- colMeans(samples)
sds <- matrixStats::colSds(samples)
names(means) <- names(sds) <- colnames(samples)
h <- (1 - options[["credibleInterval"]]) / 2
probs <- c(h, 1-h)
cri <- matrixStats::colQuantiles(samples, probs = probs)
densities <- .BANOVAfitDensity(samples, 2^9, FALSE)
tmp <- .BANOVAgetSMIResidRsq(samples, dataset, formula, nIter, options = options)
residmeans <- rowMeans(tmp$resid)
residquants <- matrixStats::rowQuantiles(tmp$resid, probs = c(h, 1-h))
# all information for q-q plot of residuals
residSumStats <- matrix(c(residmeans, residquants), nrow = length(residmeans), ncol = 3L,
dimnames = list(NULL, c("mean", "cri.2.5%", "cri.97.5%")))
# all information for r-squared density plot
rsqDens <- density(tmp$rsq, n = 2^11, from = 0, to = 1)
rsqCri <- quantile(tmp$rsq, probs = probs)
rsqMean <- mean(tmp$rsq)
return(list(
means = means, sds = sds, CRIs = cri, densities = densities$fit,
residSumStats = residSumStats, rsqDens = rsqDens, rsqCri = rsqCri, rsqMean = rsqMean,
allContinuous = allContinuous, isRandom = isRandom
))
}
.BANOVAsmiEstimates <- function(jaspContainer, options, model) {
if (!is.null(jaspContainer[["SMItablePosteriorEstimates"]]) || !options[["singleModelEstimates"]])
return()
estsTable <- createJaspTable(title = gettext("Single Model Posterior Summary"))
estsTable$position <- 1
estsTable$dependOn(c("singleModelEstimates", "credibleInterval"))
overTitle <- gettextf("%s%% Credible Interval", format(100 * options[["credibleInterval"]], digits = 3))
estsTable$addColumnInfo(name = "Variable", title = gettext("Variable"), type = "string")
estsTable$addColumnInfo(name = "Level", title = gettext("Level"), type = "string")
estsTable$addColumnInfo(name = "Mean", title = gettext("Mean"), type = "number")
estsTable$addColumnInfo(name = "SD", title = gettext("SD"), type = "number")
estsTable$addColumnInfo(name = "Lower", title = gettext("Lower"), type = "number", overtitle = overTitle)
estsTable$addColumnInfo(name = "Upper", title = gettext("Upper"), type = "number", overtitle = overTitle)
if (!(is.null(model) || estsTable$getError())) {
.BANOVAfillEstimatesTable(
jaspTable = estsTable,
mus = model$means,
sds = model$sds,
cri = model$CRIs,
hasNoLevels = model$allContinuous,
isRandom = model$isRandom
)
}
jaspContainer[["SMItablePosteriorEstimates"]] <- estsTable
return()
}
.BANOVAsmiPosteriorPlot <- function(jaspContainer, dataset, options, model) {
# meta wrapper for model averaged posterior plots, single model posterior plots, and Q-Q plots
if (!is.null(jaspContainer[["SMIposteriorPlot"]]) || !options$singleModelPosteriorPlot)
return()
posteriorPlotContainer <- createJaspContainer(title = gettext("Single Model Posterior Distributions"))
jaspContainer[["SMIposteriorPlot"]] <- posteriorPlotContainer
posteriorPlotContainer$position <- 2
posteriorPlotContainer$dependOn(c("singleModelPosteriorPlot", "singleModelGroupPosterior"))
if (is.null(model) || posteriorPlotContainer$getError()) {
posteriorPlotContainer[["dummyplot"]] <- createJaspPlot(title = gettext("Posterior distribution"), width = 400, height = 400,
plot = NULL)
} else {
.BANOVAfillPosteriorPlotContainer(
container = posteriorPlotContainer,
densities = model$densities[, -1L, ], # omit intercept
cris = model$CRIs[-1L, ], # omit intercept
isRandom = model$isRandom[-1L], # omit intercept
groupParameters = identical(options[["groupPosterior"]], "grouped")
)
}
return()
}
.BANOVAsmiQqPlot <- function(jaspContainer, options, model) {
if (!is.null(jaspContainer[["QQplot"]]) || !options$singleModelQqPlot)
return()
if (is.null(model) || jaspContainer$getError()) {
p <- NULL
} else {
p <- jaspGraphs::plotQQnorm(
residuals = model$residSumStats[,"mean"],
lower = model$residSumStats[,"cri.2.5%"],
upper = model$residSumStats[,"cri.97.5%"],
ablineColor = "darkred"
)
}
plot <- createJaspPlot(
title = gettext("Q-Q Plot"),
width = 400,
height = 400,
plot = p,
aspectRatio = 1
)
plot$dependOn("singleModelQqPlot")
plot$position <- 3
jaspContainer[["QQplot"]] <- plot
return()
}
.BANOVAsmiRsqPlot <- function(jaspContainer, options, model) {
if (!is.null(jaspContainer[["smirsqplot"]]) || !options$singleModelRsqPlot)
return()
if (is.null(model) || jaspContainer$getError()) {
p <- NULL
} else {
dd <- model$rsqDens
rsqCri <- model$rsqCri
df <- data.frame(x = dd$x, y = dd$y)
xName <- expression(R^2)
p <- jaspGraphs::PlotPriorAndPosterior(dfLines = df, xName = xName, CRI = rsqCri, drawCRItxt = FALSE)
}
plot <- createJaspPlot(
title = gettextf("Posterior R%s", "\u00B2"),
width = 400,
height = 400,
plot = p,
aspectRatio = 1
)
plot$dependOn("singleModelRsqPlot")
plot$position <- 4
jaspContainer[["smirsqplot"]] <- plot
return()
}
.BANOVAsmiRsquared <- function(jaspResults, options, model) {
if (!is.null(jaspResults[["tableSMICRI"]]) || !options[["singleModelCriTable"]])
return()
criTable <- createJaspTable(title = gettextf("Single Model R%s", "\u00B2"))
criTable$position <- 3.5
criTable$dependOn(c("singleModelCriTable", "credibleInterval"))
overTitle <- gettextf("%s%% Credible Interval", format(100 * options[["credibleInterval"]], digits = 3))
criTable$addColumnInfo(name = "rsq", title = "", type = "string")
criTable$addColumnInfo(name = "Mean", title = gettext("Mean"), type = "number")
criTable$addColumnInfo(name = "Lower", title = gettext("Lower"), type = "number", overtitle = overTitle)
criTable$addColumnInfo(name = "Upper", title = gettext("Upper"), type = "number", overtitle = overTitle)
if (!is.null(model)) {
cri <- model[["rsqCri"]]
df <- data.frame(
rsq = "R\u00B2",
Mean = model[["rsqMean"]],
Lower = cri[1L],
Upper = cri[2L],
row.names = NULL
)
criTable$setData(df)
} else {
criTable[["rsq"]] <- "R\u00B2"
}
jaspResults[["tableSMICRI"]] <- criTable
return()
}
# Citations ----
.BANOVAcitations <- c(
"MoreyEtal2015" = "Morey, R. D. & Rouder, J. N. (2015). BayesFactor (Version 0.9.10-2)[Computer software].",
"RouderEtal2012" = "Rouder, J. N., Morey, R. D., Speckman, P. L., Province, J. M., (2012) Default Bayes Factors for ANOVA Designs. Journal of Mathematical Psychology. 56. p. 356-374.",
"Jeffreys1938" = "Jeffreys, H. (1938). Significance tests when several degrees of freedom arise simultaneously. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 165, 161-198.",
"WestfallEtal1997" = "Westfall, P. H., Johnson, W. O., & Utts, J. M. (1997). A Bayesian perspective on the Bonferroni adjustment. Biometrika, 84, 419-427."
)
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.