#
# Copyright (C) 2013-2020 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/>.
#
CorrelationBayesianInternal <- function(jaspResults, dataset=NULL, options, ...) {
ready <- length(options[["variables"]]) >= 2
pairs <- unique(options[["variablePairs"]])
options <- .parseAndStoreFormulaOptions(jaspResults, options, "priorWidth")
# 2. Data retrieval --------
#
if (ready && is.null(dataset))
dataset <- .corBayesReadData(dataset, options)
corModel <- .computeCorBayes(jaspResults = jaspResults, dataset = dataset, options = options, ready = ready)
# 3. Table: Get, compute and fill -----------
#
.createTableCorBayes(jaspResults = jaspResults, corModel = corModel,
dataset = dataset, options = options)
# 4. Matrix plot: Get, compute and draw ------
#
if (options[["matrixPlot"]])
.createMatrixPlotCorBayes(jaspResults = jaspResults, corModel = corModel,
dataset = dataset, options = options)
# 5. PairsPlot Container: (Optional) ------
#
if (options[["scatterPlot"]] || options[["priorPosteriorPlot"]] || options[["bfRobustnessPlot"]] ||
options[["bfSequentialPlot"]]) {
pairsPlotCollection <- .getPairPlotsContainerCorBayes(jaspResults, options)
if (!is.null(pairsPlotCollection)) {
pairsPlotCollection <- .initPlotContainerSubStructureCorBayes(pairsPlotCollection, options)
.fillPairsPlotsCorBayes(jaspResults, pairsPlotCollection, corModel, dataset, options)
}
}
}
.computeCorBayes <- function(jaspResults, dataset, options, ready = TRUE) {
if (!is.null(jaspResults[["corModel"]]))
return(jaspResults[["corModel"]]$object)
if (!ready)
return(NULL)
if (options[["naAction"]] == "listwise")
.hasErrors(dataset, type="observations", observations.amount='< 2', exitAnalysisIfErrors=TRUE)
result <- list()
pairs <- combn(options[["variables"]], 2, simplify=FALSE)
for (pair in pairs) {
var1 <- pair[1]
var2 <- pair[2]
# TODO(Alexander) check if they have the same name. Or just let it be, but highly inefficient though
#
pairName <- paste(sort(c(var1, var2)), collapse="-")
result[[pairName]] <- list()
for (method in c("pearson", "kendall")) {
errorMsg <- NULL
dataCheck <- .corBayesCheckPairsErrors(dataset, var1, var2)
if (!identical(dataCheck, FALSE)) {
errorMsg <- dataCheck[["message"]]
v1 <- NA
v2 <- NA
} else {
v1 <- dataset[[var1]]
v2 <- dataset[[var2]]
}
.setSeedJASP(options)
bfObject <- bstats::bcor.test("x"=v1, "y"=v2, "kappa"=options[["priorWidth"]],
"method"=method, "ciValue"=options[["ciLevel"]])
if (!is.null(errorMsg)) {
bfObject[["error"]] <- errorMsg
bfObject[["dataError"]] <- errorMsg
} else {
# TODO(Alexander): This can be removed once bstats is adapted
alternativeLocal <- switch(options[["alternative"]],
twoSided = "two.sided",
greater = "greater",
less = "less")
ci <- .computeCorCredibleInterval(bfObject, ciLevel=options[["ciLevel"]], method)
if (!is.null(ci[["error"]]) && isTRUE(bfObject[[alternativeLocal]][["tooPeaked"]]))
ci[["error"]] <- gettext("Posterior is too peaked")
bfObject <- modifyList(bfObject, ci)
}
result[[pairName]][[method]] <- bfObject
}
}
jaspResults[["corModel"]] <- createJaspState(result)
jaspResults[["corModel"]]$dependOn(c("naAction", "variables", "priorWidth", "ciLevel", "setSeed", "seed"))
return(result)
}
.corBayesCheckPairsErrors <- function(dataset, var1, var2) {
errors <- .hasErrors(dataset, message='short',
type=c('observations', 'variance', 'infinity', 'observationsPairwise'),
all.target=c(var1, var2), observations.amount='< 2',
observationsPairwise.amount=2)
return(errors)
}
.createTableCorBayes <- function(jaspResults, corModel, dataset, options) {
if (!is.null(jaspResults[["corBayesTable"]]))
return()
methodItems <- .getCorMethods(options)
corBayesTable <- createJaspTable(title=.getCorTableTitle(methodItems, bayes=TRUE))
corBayesTable$position <- 1
# TODO(Alexander): This can be removed once bstats is adapted
alternativeLocal <- switch(options[["alternative"]],
twoSided = "two.sided",
greater = "greater",
less = "less")
corBayesTable$dependOn(c("pearson", "kendall", "spearman", "alternative", "priorWidth", "variables",
"pairwiseDisplay","bayesFactorReport", "naAction",
"supportCorrelationFlagged", "ci", "ciLevel",
"sampleSize", "posteriorMedian", "bayesFactorType",
"setSeed", "seed"))
corBayesTable$showSpecifiedColumnsOnly <- TRUE
corBayesTable$position <- 1
corBayesTable$addCitation(.getCorCitations(methodItems, bayes=TRUE))
if (alternativeLocal=="greater")
corBayesTable$addFootnote(.getBfTableSidedFootnote(alternative="greater", analysis="correlation"))
if (alternativeLocal=="less")
corBayesTable$addFootnote(.getBfTableSidedFootnote(alternative="less", analysis="correlation"))
# Add legend footnote of BFs
if (options[["supportCorrelationFlagged"]]) {
tempNote <- .bfFlagTableFootnote(options)
corBayesTable$addFootnote(message=tempNote, symbol="*")
}
.addTableColumnMarkupCorBayes(corBayesTable, methodItems, options)
.fillTableCorBayes(corBayesTable, options, corModel)
jaspResults[["corBayesTable"]] <- corBayesTable
}
.addTableColumnMarkupCorBayes <- function(table, methodItems, options) {
if (options[["pairwiseDisplay"]])
.addPairwiseTableColumnMarkupCorBayes(table, methodItems, options)
else
.addMatrixTableColumnMarkupCorBayes(table, methodItems, options)
}
.addPairwiseTableColumnMarkupCorBayes <- function(table, methodItems, options) {
# Add the variables names paired
#
table$addColumnInfo(name="variable1", title="", type="string")
table$addColumnInfo(name="separator", title="", type="separator")
table$addColumnInfo(name="variable2", title="", type="string")
if (options[["sampleSize"]])
table$addColumnInfo(name="n", title="n", type="integer")
for (m in seq_along(methodItems)) {
methodName <- methodItems[m]
overTitle <- NULL
if (length(methodItems) > 1) {
# Overwrite the .cormethodNamesList with one that is broken up in overtitle and methodName
#
overTitle <- .corOverTitlesList[[methodName]]
.corMethodNamesList <- list(pearson=gettext("r"), spearman=gettext("rho"), kendall=gettext("tau B"))
}
# Add's "r", "rho", "tau B"
table$addColumnInfo(name=paste0(methodName, "stat"), title=.corMethodNamesList[[methodName]],
overtitle=overTitle, type="number")
if (options[["bayesFactorReport"]])
table$addColumnInfo(name=paste0(methodName, "bf"), title=.getBfTitle(options[["bayesFactorType"]], options[["alternative"]]), overtitle=overTitle, type="number")
# # TODO(ALEXANDER): Also report error %? Only useful for mcmc
# #
# if (options[["reportPercentageError"]])
# table$addColumnInfo(name=paste0(methodName, "bfUq"), title=bfTitle, overtitle=overTitle, type="number")
# TODO(Alexander): Do we want to show the posterior median?
#
# if (options[["posteriorMedian"]])
# table$addColumnInfo(name=paste0(methodName, "posteriorMedian"), title="Posterior Median", type="number")
#
if (options[["ci"]]) {
table$addColumnInfo(name=paste0(methodName, "lowerCi"), overtitle=overTitle, type="number",
title=gettextf("Lower %s%% CI", options[["ciLevel"]] * 100))
table$addColumnInfo(name=paste0(methodName, "upperCi"), overtitle=overTitle, type="number",
title=gettextf("Upper %s%% CI", options[["ciLevel"]] * 100))
}
}
# Calculate expected dimension
#
nVariables <- length(options[["variables"]])
if (nVariables <= 1) {
numberOfRows <- 1
} else {
numberOfRows <- choose(nVariables, 2)
}
if (is.null(methodItems) || length(methodItems)==0) {
numberOfColumns <- 3
} else {
numberOfColumns <- 3 + options[["sampleSize"]] + length(methodItems) * sum(options[["bayesFactorReport"]],
2*options[["ci"]])
}
table$setExpectedSize(numberOfRows, numberOfColumns)
}
.addMatrixTableColumnMarkupCorBayes <- function(table, methodItems, options) {
# Note(Alexander): Correlation Matrix
#
table$addColumnInfo(name="variable", title=gettext("Variable"), type="string", combine=TRUE)
table$addColumnInfo(name="itemColumn", title="", type="string")
nVariables <- length(options[["variables"]])
if (nVariables < 2) {
if (nVariables == 0)
title <- "..."
else
title <- options[["variables"]][1]
table$addColumnInfo(name="firstPlaceholder", type="number", title=title)
table$addColumnInfo(name="secondPlaceholder", type="number", title="...")
} else {
for (variable in options[["variables"]])
table$addColumnInfo(name=variable, type="number",title=variable)
}
if (length(methodItems)==0) {
nItems <- sum(options[["bayesFactorReport"]], 2*options[["ci"]], options[["sampleSize"]])
if (nItems <= 1) {
numberOfColumns <- nVariables
} else {
numberOfColumns <- 2 + nVariables
}
numberOfRows <- nVariables
} else {
itemsPerMethod <- sum(options[["bayesFactorReport"]], 2*options[["ci"]])
nVariablesEffective <- max(nVariables, 2)
if (itemsPerMethod==0) {
numberOfColumns <- 1 + nVariablesEffective
numberOfRows <- nVariablesEffective
} else {
numberOfColumns <- 2 + nVariablesEffective
numberOfRows <- nVariablesEffective * length(methodItems) *
(itemsPerMethod + 1) + options[["sampleSize"]] * nVariablesEffective
}
}
table$setExpectedSize(numberOfRows, numberOfColumns)
}
.fillTableCorBayes <- function(table, options, corModel) {
if (is.null(corModel)) {
if (options[["pairwiseDisplay"]])
.insertPairwiseDefaultEmptyTableCellsCorBayes(table, options)
else
.insertMatrixDefaultEmptyTableCellsCorBayes(table, options)
} else {
if (options[["pairwiseDisplay"]])
.fillPairwiseTableCorBayes(table, options, corModel)
else
.fillMatrixTableCorBayes(table, options, corModel)
}
}
.insertPairwiseDefaultEmptyTableCellsCorBayes <- function(table, options) {
nVariables <- length(options[["variables"]])
if (nVariables < 2) {
if (nVariables == 0)
var1 <- "..."
else
var1 <- options[["variables"]][1]
var2 <- "..."
table$addRows(list(variable1=var1, separator="\u2014", variable2=var2))
}
}
.insertMatrixDefaultEmptyTableCellsCorBayes <- function(table, options) {
variables <- unlist(options[["variables"]])
if (length(variables) == 1)
variables <- c(paste("1.", variables), "2. ...")
else
variables <- c("1. ...", "2. ...")
methodItems <- .getCorMethods(options)
itemNames <- .bSelectItems(options)
# TODO(Alexander): fix length is better
statColumnPerVar <- c()
for (methodItem in methodItems) {
statsToReport <- unlist(.bCorRowNames(options, itemNames, method=methodItem))
statsToReport <- statsToReport[statsToReport != gettext("n")] # we want to only report this once per variable. probably better solved in .bCorRowNames() tho
statColumnPerVar <- c(statColumnPerVar, statsToReport)
}
if (options[["sampleSize"]])
statColumnPerVar <- c(gettext("n"), statColumnPerVar)
emptyCells <- matrix(c("\u2014", "", ".", "\u2014"), nrow=2, byrow=TRUE)
for (i in seq_along(variables)) {
for (j in seq_along(statColumnPerVar)) {
table$addRows(list(.isNewGroup = (i == 2 && j == 1),
variable = variables[i],
itemColumn = statColumnPerVar[j],
firstPlaceholder = emptyCells[i, 1],
secondPlaceholder = emptyCells[i, 2]))
}
}
}
.fillPairwiseTableCorBayes <- function(table, options, corModel) {
pairs <- combn(options[["variables"]], 2, simplify=FALSE)
for (pair in pairs) {
var1 <- pair[1]
var2 <- pair[2]
pairName <- paste(sort(c(var1, var2)), collapse="-")
tempRow <- list("variable1"=var1, "separator"="-", "variable2"=var2)
itemNames <- .bSelectItems(options)
# TODO(Alexander): This can be removed once bstats is adapted
alternativeLocal <- switch(options[["alternative"]],
twoSided = "two.sided",
greater = "greater",
less = "less")
for (method in .getCorMethods(options)) {
bfObject <- corModel[[pairName]][[method]]
sidedObject <- bstats::getSidedObject(bfObject, alternative=alternativeLocal,
itemNames=itemNames)
rowObject <- sidedObject[itemNames]
sampleSize <- rowObject[["n"]]
reportBf <- rowObject[["bf"]]
if (options[["bayesFactorReport"]]) {
if (options[["bayesFactorType"]] == "BF01")
rowObject[["bf"]] <- 1/reportBf
else if (options[["bayesFactorType"]] == "LogBF10")
rowObject[["bf"]] <- log(reportBf)
}
names(rowObject) <- paste0(method, itemNames)
if (!is.null(sampleSize))
rowObject[["n"]] <- sampleSize
errorMessage <- bfObject[["error"]]
if (!is.null(errorMessage))
table$addFootnote(errorMessage, rowNames = pairName, colNames = paste0(method, "stat"))
else if (options[["bayesFactorReport"]] && options[["supportCorrelationFlagged"]] && !is.na(reportBf))
.addSupportedBfIndicators(table, reportBf, rowNames = pairName, colNames = paste0(method, "stat"))
tempRow <- modifyList(tempRow, rowObject)
}
table$addRows(tempRow, rowNames=pairName)
}
}
.fillMatrixTableCorBayes <- function(table, options, corModel) {
nVariables <- length(options[["variables"]])
for (i in 1:nVariables) {
emptyCellInfo <- vector("list", length=i)
var1 <- options[["variables"]][i]
var2Names <- options[["variables"]][1:i]
# Note(Alexander): Dash on the diagonal
#
emptyCellInfo[[i]] <- "\u2014"
itemNames <- .bSelectItems(options)
methodItems <- .getCorMethods(options)
for (m in seq_along(.getCorMethods(options))) {
# Row info collected here
#
methodName <- methodItems[m]
allItemInfo <- vector("list", length=length(itemNames))
names(allItemInfo) <- itemNames
for (item in itemNames)
allItemInfo[[item]] <- emptyCellInfo
# Note(Alexander): RETRIEVE info from state here
#
# TODO(Alexander): This can be removed once bstats is adapted
alternativeLocal <- switch(options[["alternative"]],
twoSided = "two.sided",
greater = "greater",
less = "less")
if (i > 1) {
for (j in 1:(i-1)) {
var2 <- var2Names[j]
pairName <- paste(sort(c(var1, var2)), collapse="-")
# Note(Alexander): Here use actual data to fill the table
#
bfObject <- corModel[[pairName]][[methodName]]
sidedObject <- bstats::getSidedObject(bfObject, alternative=alternativeLocal,
itemNames=itemNames)
reportBf <- sidedObject[["bf"]]
if (options[["bayesFactorType"]] == "BF01")
sidedObject[["bf"]] <- 1/reportBf
else if (options[["bayesFactorType"]] == "LogBF10")
sidedObject[["bf"]] <- log(reportBf)
# Here add to info to the collection
#
for (item in itemNames)
allItemInfo[[item]][j] <- sidedObject[[item]]
errorMessage <- bfObject[["error"]]
rowName <- paste0(var1, methodName, "stat")
if (!is.null(errorMessage))
table$addFootnote(errorMessage, rowNames = rowName, colNames = var2)
else if (options[["bayesFactorReport"]] && options[["supportCorrelationFlagged"]] && !is.na(reportBf))
.addSupportedBfIndicators(table, reportBf, rowNames = rowName, colNames = var2)
}
}
# Note(Alexander): Now fill each row. For the ith variable the structure is
# variable=var1,
# itemColumn,
# var2Name[1],
# var2Name[2],
# ...,
# var2Name[i-1],
# var2Name[i]
#
# To be used for the column "itemName"
#
if (m > 1 && options[["sampleSize"]]) {
allItemInfo[["n"]] <- NULL
# Note(Alexander): After the first time we added "n", remove it from the other
#
itemNames <- setdiff(itemNames, "n")
}
rowItemNames <- .bCorRowNames(options, itemNames, method = methodName)
for (k in seq_along(itemNames)) {
if (k == 1 && m == 1)
tempRow <- list(variable = paste0(i, ". ", var1), itemColumn = rowItemNames[[k]], .isNewGroup = TRUE) # First item gets the variable name
else
tempRow <- list(variable = "", itemColumn = rowItemNames[[k]], .isNewGroup = FALSE)
itemInfo <- as.list(allItemInfo[[k]])
names(itemInfo) <- var2Names
tempRow <- modifyList(tempRow, itemInfo)
table$addRows(tempRow, rowNames = paste0(var1, methodName, itemNames[k]))
} # End items loop
} # End loop over the methodItems
} # End var1 loop
}
.addSupportedBfIndicators <- function(table, bf, rowNames, colNames) {
symbol <- NULL
if (bf >= 100)
symbol <- "***"
else if (bf >= 30)
symbol <- "**"
else if (bf >= 10)
symbol <- "*"
if (!is.null(symbol))
table$addFootnote(symbol = symbol, rowNames=rowNames, colNames=colNames)
}
.createMatrixPlotCorBayes <- function(jaspResults, corModel, dataset, options) {
# a. Retrieve from state ----
#
if (!is.null(jaspResults[["matrixPlot"]]))
return()
nVariables <- length(options[["variables"]])
# b. Compute -----
#
methodItems <- .getCorMethods(options)
# c. Get ----
#
matrixPlot <- createJaspPlot(title=gettext("Bayesian Correlation Matrix Plot"))
matrixPlot$position <- 2
matrixDependencies <- c("variables", "matrixPlot", "matrixPlotDensity", "matrixPlotPosterior", "naAction", "setSeed", "seed", "pearson")
if (options[["matrixPlotPosterior"]])
matrixDependencies <- c(matrixDependencies, "spearman", "kendall", "alternative", "priorWidth")
matrixPlot$dependOn(matrixDependencies)
jaspResults[["matrixPlot"]] <- matrixPlot
if (!is.null(corModel)) {
# d. Draw -----
#
if (nVariables <= 2 && (options[["matrixPlotDensity"]] || options[["matrixPlotPosterior"]])) {
matrixPlot[["width"]] <- 580
matrixPlot[["height"]] <- 580
} else if (nVariables <= 2) {
matrixPlot[["width"]] <- 400
matrixPlot[["height"]] <- 400
} else {
matrixPlot[["width"]] <- 250 * nVariables + 20
matrixPlot[["height"]] <- 250 * nVariables + 20
}
if (length(methodItems) != 0 && options[["matrixPlot"]])
.drawMatrixPlotCorBayes(jaspResults, matrixPlot, corModel, dataset, options)
}
}
.drawMatrixPlotCorBayes <- function(jaspResults, matrixPlot, corModel, dataset, options) {
methodItems <- .getCorMethods(options)
vars <- options[["variables"]]
nVariables <- length(vars)
plotMat <- matrix(list(), nVariables, nVariables)
useRanks <- !options[["pearson"]]
for (row in seq_len(nVariables-1)) {
for (col in row:nVariables) {
var1 <- vars[row]
var2 <- vars[col]
pairName <- paste(sort(c(var1, var2)), collapse="-")
if (row < col) {
if (!options[["matrixPlotPosterior"]]) {
posteriorPlot <- list()
} else {
posteriorPlot <- .drawPosteriorPlotCorBayes(jaspResults, corModel, options, methodItems, purpose="matrix",
pairName)
}
plotMat[[col, row]] <- posteriorPlot
dataError <- corModel[[pairName]][[1]][["dataError"]]
if (!is.null(dataError)) {
# TODO(Alexander): This is probably not run here anyways, but try avoid displayerror
scatterPlot <- .displayError(errorMessage=dataError)
} else {
subData <- dataset[, c(var1, var2), drop=FALSE]
subData <- subData[complete.cases(subData), , drop=FALSE]
scatterPlot <- .bCorScatter(x=subData[[var2]], y=subData[[var1]],
useRanks = useRanks)
}
plotMat[[row, col]] <- scatterPlot
}
if (row == col) {
if (!options[["matrixPlotDensity"]]) {
densityPlot <- list()
} else {
densityPlot <-try(.bCorMarginalDistribution(variable = dataset[, vars[row], drop=TRUE],
varName = vars[row], options = options)
)
if (isTryError(densityPlot))
# TODO(Alexander): again with error case for matrix
densityPlot <- .displayError(errorMessage=.extractErrorMessage(densityPlot))
}
plotMat[[row, row]] <- densityPlot
}
}
}
if (options[["matrixPlotDensity"]]) {
densityPlot <- try(.bCorMarginalDistribution(variable = dataset[, vars[nVariables], drop=TRUE],
varName = vars[nVariables], options = options)
)
if (isTryError(densityPlot))
# TODO(Alexander): again with error case for matrix
densityPlot <- .displayError(errorMessage=.extractErrorMessage(densityPlot))
plotMat[[nVariables, nVariables]] <- densityPlot
}
# TODO(Alexander): Stolen from Descriptives.R
# slightly adjust the positions of the labels left and above the plots.
#
labelPos <- matrix(.5, 4, 2)
labelPos[1, 1] <- .55
labelPos[4, 2] <- .65
labels <- if (useRanks) .addRankToVariableName(vars) else vars
obj <- try(jaspGraphs::ggMatrixPlot(plotList = plotMat, leftLabels = labels, topLabels = labels,
scaleXYlabels = 0.9, labelPos=labelPos))
if (isTryError(obj)) {
matrixPlot$setError(.extractErrorMessage(obj))
} else {
matrixPlot$plotObject <- obj
}
return()
}
.drawPosteriorPlotCorBayes <- function(jaspResults, corModel, options, methodItems,
purpose=c("matrix", "pairs", "sumStat"), pairName = NULL) {
purpose <- match.arg(purpose)
plotResult <- ""
# TODO(Alexander): Check whether we can do .setErrorMessages() to plots in correlation matrix plot
if (purpose=="matrix") {
if (!options[["matrixPlotPosterior"]]) {
# Note(Alexander): This should never happen as, this shouldn't be called with matrixPlotPosterior <- FALSE
# TODO(Alexander) return empty list instead?
plotResult <- .displayError(errorMessage="")
return(plotResult)
}
} else if (purpose=="pairs" || purpose=="sumStat") {
if (!isTRUE(options[["priorPosteriorPlot"]])) {
# Note(Alexander): This should never happen as, this shouldn't be called with priorPosteriorPlot <- FALSE or NULL
# TODO(Alexander) return empty list instead?
plotResult <- .displayError(errorMessage="")
return(plotResult)
}
}
if (purpose == "sumStat") {
pairStats <- list()
pairStats[[options[["method"]]]] <- corModel
} else {
pairStats <- corModel[[pairName]]
}
if (purpose %in% c("matrix", "pairs"))
dataError <- pairStats[[1]][["error"]]
else if (purpose=="sumStat")
dataError <- pairStats[["error"]]
if (!is.null(dataError)) {
if (purpose=="matrix")
plotResult <- .displayError(errorMessage=dataError)
else
plotResult <- dataError
return(plotResult)
}
# a. No data errors, try to plot ----
#
postPlotValuesPerMethod <- .getPosteriorPlotValuesCorBayes(jaspResults, pairStats, options, pairName)
for (i in seq_along(postPlotValuesPerMethod)) {
postPlotValues <- postPlotValuesPerMethod[[i]]
error <- postPlotValues[["error"]]
if (postPlotValues[["tooPeaked"]]) {
error <- gettext("Posterior is too peaked")
break()
}
posteriorLine <- postPlotValues[["posteriorLine"]]
if (is.null(posteriorLine) || length(posteriorLine) == 1) {
# TODO(Alexander): This covers the case when posteriorLine <- "Could not compute posterior"
# RETURN FAILED PLOT
error <- gettext("Could not compute the posteriors")
break()
}
}
# b. Check for error ----
#
if (!is.null(error)) {
if (purpose=="matrix")
plotResult <- .displayError(errorMessage=error)
else
plotResult <- error
return(plotResult)
}
# c. No errors, combine plot info ------
#
xDomain <- postPlotValuesPerMethod[[1]][["xDomain"]] # these values are the same for all methods
domainLength <- length(xDomain)
dfPoints <- NULL
BF10 <- NULL
CRI <- NULL
medianPoint <- NULL
CRItxt <- NULL
if (purpose == "matrix") {
xName <- NULL
if (length(methodItems) == 1) {
# xName <- unlist(.corXNames[methodItems], recursive = FALSE)
#
dfLines <- data.frame(
x = xDomain,
y = postPlotValuesPerMethod[[1]][["posteriorLine"]]
)
} else {
# xName <- "Correlation Coefficient"
#
gLegend <- unlist(.corGLegendList[methodItems], use.names = FALSE, recursive = FALSE)
dfLines <- data.frame(
x = xDomain,
y = unlist(lapply(postPlotValuesPerMethod, function(x) return(x[["posteriorLine"]]))),
g = rep(gLegend, each=domainLength)
)
}
} else if (purpose %in% c("pairs", "sumStat")) {
if (purpose=="sumStat")
postPlotValues <- postPlotValuesPerMethod[[1]]
else if (purpose=="pairs")
postPlotValues <- postPlotValuesPerMethod[[options[["pairsMethod"]]]]
xName <- unlist(.corXNames[[methodItems]], recursive = FALSE)
gLegend <- c(gettext("Prior"), gettext("Posterior"))
dfLines <- data.frame(
x = xDomain,
y = c(postPlotValues[["priorLine"]], postPlotValues[["posteriorLine"]]),
g = rep(gLegend, each=domainLength)
)
# TODO(Alexander): This can be removed once bstats is adapted
alternativeLocal <- switch(options[["alternative"]],
twoSided = "two.sided",
greater = "greater",
less = "less")
if (isTRUE(options[["priorPosteriorPlotAdditionalEstimationInfo"]])) {
if (postPlotValues[["ciValue"]] != options[["ciLevel"]]) {
if (purpose=="pairs") {
methodName <- options[["pairsMethod"]]
sidedBfWithoutPost <- corModel[[pairName]][[methodName]][[alternativeLocal]]
} else if (purpose=="sumStat") {
methodName <- options[["method"]]
sidedBfWithoutPost <- corModel[[alternativeLocal]]
}
# Note(Alexander): Not stored
#
lowerCi <- sidedBfWithoutPost[["lowerCi"]]
upperCi <- sidedBfWithoutPost[["upperCi"]]
medianPoint <- sidedBfWithoutPost[["posteriorMedian"]]
} else {
lowerCi <- postPlotValues[["lowerCi"]]
upperCi <- postPlotValues[["upperCi"]]
medianPoint <- postPlotValues[["posteriorMedian"]]
}
ciLevel <- options[["ciLevel"]]
CRI <- c(lowerCi, upperCi)
CRItxt <- gettextf("%s%% CI:", ciLevel * 100)
}
if (isTRUE(options[["priorPosteriorPlotAdditionalTestingInfo"]])) {
dfPoints <- data.frame(
x = c(postPlotValues[["h0"]], postPlotValues[["h0"]]),
y = c(postPlotValues[["priorAtH0"]], postPlotValues[["posteriorAtH0"]]),
g = c(gettext("Prior"), gettext("Posterior"))
)
BF10 <- postPlotValues[["bf"]]
}
}
hypothesis <- switch(options[["alternative"]],
"twoSided"="equal",
"greater"="greater",
"less"="smaller"
)
plotResult <- try(jaspGraphs::PlotPriorAndPosterior(
dfLines, dfPoints, BF10, "CRI" = CRI, "CRItxt" = CRItxt, "median" = medianPoint, "xName" = xName,
"hypothesis" = hypothesis, bfType = "BF10")
)
return(plotResult)
}
.getPosteriorPlotValuesCorBayes <- function(jaspResults, bfObject, options, pair = NULL) {
indexName <- paste0("posteriorLine-", pair)
if (!is.null(jaspResults[[indexName]]))
return(jaspResults[[indexName]]$object)
results <- list()
# TODO(Alexander): This can be removed once bstats is adapted
alternativeLocal <- switch(options[["alternative"]],
twoSided = "two.sided",
greater = "greater",
less = "less")
for (method in names(bfObject)){
.setSeedJASP(options)
results[[method]] <- bstats::computeCorPosteriorLine(bfObject = bfObject[[method]],
alternative = alternativeLocal)
}
jaspResults[[indexName]] <- createJaspState(results)
if (is.null(pair)) # sumStats
qmlInputElements <- c("alternative", "priorWidth", "setSeed", "seed")
else # pairs
qmlInputElements <- c("naAction", "alternative", "priorWidth", "setSeed", "seed")
jaspResults[[indexName]]$dependOn(qmlInputElements)
return(results)
}
.getPairPlotsContainerCorBayes <- function(jaspResults, options) {
if (!.readyToCreatePairPlotsCorBayes(options))
return(NULL)
pairsPlotCollection <- jaspResults[["pairsPlotCollection"]]
if (is.null(pairsPlotCollection)) {
pairsPlotCollection <- createJaspContainer(title=gettext("Bayesian Correlation Pairwise Plots"))
pairsPlotCollection$dependOn("naAction")
pairsPlotCollection$position <- 3
jaspResults[["pairsPlotCollection"]] <- pairsPlotCollection
}
return(pairsPlotCollection)
}
.readyToCreatePairPlotsCorBayes <- function(options) {
pairsIndeces <- .getPairsIndeces(options)
if (length(pairsIndeces) == 0)
return(FALSE)
plotItems <- .getCorPlotItems(options, bayes=TRUE, sumStat=FALSE)
if (length(plotItems) == 0)
return(FALSE)
return(TRUE)
}
.initPlotContainerSubStructureCorBayes <- function(pairsPlotCollection, options) {
plotItems <- .getCorPlotItems(options, bayes=TRUE, sumStat=FALSE)
bfPlotPriorPosteriorDependencies <- c("pairsMethod", "priorWidth", "alternative", "setSeed", "seed")
bfPlotDependencies <- c(bfPlotPriorPosteriorDependencies, "bayesFactorType")
if (options[["priorPosteriorPlot"]] && options[["priorPosteriorPlotAdditionalEstimationInfo"]])
bfPlotPriorPosteriorDependencies <- c(bfPlotPriorPosteriorDependencies, "ciLevel")
plotItemDependencies <- list(
"scatterPlot"=c("scatterPlot", "scatterPlotAdditionalInfo", "pairsMethod"),
"priorPosteriorPlot"=c("priorPosteriorPlot", bfPlotPriorPosteriorDependencies,
"priorPosteriorPlotAdditionalTestingInfo", "priorPosteriorPlotAdditionalEstimationInfo"),
"bfRobustnessPlot"=c("bfRobustnessPlot", "bfRobustnessPlotAdditionalInfo", bfPlotDependencies),
"bfSequentialPlot"=c("bfSequentialPlot", "bfSequentialPlotAdditionalInfo", bfPlotDependencies)
)
pairs <- options[["variablePairs"]]
for (i in seq_along(pairs)) {
# Loop over plotItems
pair <- pairs[[i]]
var1 <- pair[1]
var2 <- pair[2]
# TODO(Alexander): Think about when var1==var2
#
pairName <- paste(sort(c(var1, var2)), collapse="-")
pairContainerTitle <- paste(var1, var2, sep=" - ")
pairContainer <- pairsPlotCollection[[pairName]]
if (is.null(pairContainer)) {
pairContainer <- createJaspContainer(title = pairContainerTitle)
pairContainer$dependOn(optionContainsValue = list(pairs = unname(pair)))
pairsPlotCollection[[pairName]] <- pairContainer
}
for (i in seq_along(plotItems)) {
item <- plotItems[i]
itemTitle <- .bfPlotTitles(item)
itemPlot <- pairContainer[[item]]
if (is.null(itemPlot)) {
itemPlot <- createJaspPlot(title=itemTitle, width=530, height=400)
itemPlot$dependOn(options = plotItemDependencies[[item]])
itemPlot$position <- switch(item,
"scatterPlot"=1,
"priorPosteriorPlot"=2,
"bfRobustnessPlot"=3,
"bfSequentialPlot"=4)
pairContainer[[item]] <- itemPlot
}
}
}
return(pairsPlotCollection)
}
.fillPairsPlotsCorBayes <- function(jaspResults, pairsPlotCollection, corModel, dataset, options) {
plotItems <- .getCorPlotItems(options, bayes=TRUE, sumStat=FALSE)
pairs <- options[["variablePairs"]]
# TODO(Alexander): Get indeces for useable pairs instead
pairsIndeces <- .getPairsIndeces(options)
# TODO(Alexander): This can be removed once bstats is adapted
alternativeLocal <- switch(options[["alternative"]],
twoSided = "two.sided",
greater = "greater",
less = "less")
thisMethod <- options[["pairsMethod"]]
for (i in pairsIndeces) {
currentPair <- pairs[[i]]
var1 <- currentPair[1]
var2 <- currentPair[2]
pairName <- paste(sort(c(var1, var2)), collapse="-")
pairContainer <- pairsPlotCollection[[pairName]]
pairStats <- corModel[[pairName]]
dataError <- pairStats[[1]][["dataError"]]
if (!is.null(dataError)) {
pairContainer$setError(dataError)
next
}
for (i in seq_along(plotItems)) {
item <- plotItems[i]
jaspPlotResult <- pairContainer[[item]]
plotResult <- jaspPlotResult$plotObject
if (is.null(plotResult)) {
jaspPlotResult$status <- "running"
plot <- NULL
if (item == "scatterPlot") {
subData <- dataset[, c(var1, var2), drop=FALSE]
subData <- subData[complete.cases(subData), , drop=FALSE]
plot <- try(.bCorScatter(x=subData[[var1]], y=subData[[var2]], xName=var1, yName=var2,
useRanks = options[["pairsMethod"]] != "pearson"))
} else if (item == "priorPosteriorPlot") {
plot <- .drawPosteriorPlotCorBayes(jaspResults, corModel, options, thisMethod, purpose="pairs", "pairName"=pairName)
} else if (item == "bfRobustnessPlot") {
plot <- .drawBfRobustnessPlotCorBayes(pairStats[[thisMethod]], options, thisMethod)
} else if (item == "bfSequentialPlot") {
plot <- .drawBfSequentialPlotCorBayes(dataset[[var1]], dataset[[var2]], pairStats[[thisMethod]],
options)
}
.checkAndSetPlotCorBayes(plot, jaspPlotResult)
}
}
}
return()
}
.drawBfRobustnessPlotCorBayes <- function(bfObject, options, method) {
.setSeedJASP(options)
robustnessValuesPerHypothesis <- bstats::computeCorRobustnessLine(bfObject)
# TODO(Alexander): This can be removed once bstats is adapted
alternativeLocal <- switch(options[["alternative"]],
twoSided = "two.sided",
greater = "greater",
less = "less")
robustnessValues <- bstats::getSidedObject(robustnessValuesPerHypothesis, alternative=alternativeLocal,
itemNames=c("kappaDomain", "kappa"))
robustnessLine <- robustnessValues[["robustnessLine"]]
if (is.character(robustnessLine)) {
plotResult <- robustnessLine
} else {
xLine <- robustnessValues[["kappaDomain"]]
yLine <- robustnessLine
xPoint <- NULL
yPoint <- NULL
# TODO(Alexander): This can be removed once bstats is adapted
alternativeLocal <- switch(options[["alternative"]],
twoSided = "two.sided",
greater = "greater",
less = "less")
if (isTRUE(options[["bfRobustnessPlotAdditionalInfo"]])) {
maxBf <- robustnessValues[["robustnessMaxBf"]]
kappaOfMaxBf <- robustnessValues[["robustnessKappaOfMaxBf"]]
userBf <- bfObject[[alternativeLocal]][["bf"]]
userKappa <- bfObject[["kappa"]]
xPoint <- c(kappaOfMaxBf, userKappa)
yPoint <- c(maxBf, userBf)
}
plotResult <- .plotRobustnessCor(xLine=xLine, yLine=yLine, xPoint=xPoint, yPoint=yPoint,
bfType=options[["bayesFactorType"]], alternative=alternativeLocal,
nDigits=4)
}
return(plotResult)
}
.plotRobustnessCor <- function(xLine, yLine, xPoint, yPoint, bfType, alternative, nDigits=4,
pointColors=c("red", "grey", "black", "white"), pointFill=c("grey", "black", "white")) {
# TODO(Alexander): Perhaps add an argument
# category=c("max", "user")
#
bfPlotType <- switch(bfType,
"BF10"="BF10",
"LogBF10"="BF10",
"BF01"="BF01"
)
logYLine <- try(log(yLine))
if (isTryError(logYLine))
return(logYLine)
if (bfType=="BF01")
logYLine <- -logYLine
dfLines <- data.frame(
x = xLine,
y = logYLine
)
hypothesisjaspGraphsName <- switch(alternative,
"twoSided"="equal",
"greater"="greater",
"less"="smaller")
if (is.null(xPoint)) {
dfPoints <- NULL
pointColors <- NULL
pointFill <- NULL
} else {
nPoints <- length(xPoint)
pointColors <- pointColors[1:nPoints]
pointFill <- pointFill[1:nPoints]
bfLegendLabel <- jaspGraphs::getBFSubscripts(bfPlotType, hypothesis=hypothesisjaspGraphsName)[1]
legendText1 <- vector("character", length(xPoint))
legendText2 <- vector("character", length(xPoint))
for (i in seq_along(xPoint)) {
if (i==1) {
legendText1[i] <- gettextf("max %s", bfLegendLabel)
legendText1[i] <- gsub(pattern = "\\s+", "~", legendText1[i])
legendText2[i] <- gettextf("%1$s at kappa==%2$s", format(yPoint[i], digits=nDigits), format(xPoint[i], digits=nDigits))
legendText2[i] <- gsub(pattern = "\\s+", "~", legendText2[i])
} else if (i==2) {
legendText1[i] <- gettext("user prior")
legendText1[i] <- paste0("\"", legendText1[i], "\"")
legendText2[i] <- gettextf("%1$s at kappa==%2$s", format(yPoint[i], digits=nDigits), format(xPoint[i], digits=nDigits))
legendText2[i] <- gsub(pattern = "\\s+", "~", legendText2[i])
}
# TODO(Alexander): Do something here, with user-prior, wide etc etc. Also loop over colours
}
legendText1 <- paste0("paste(", legendText1, ", ':')")
# TODO(Alexander): Try and think of how to default
logYPoint <- log(.recodeBFtype(yPoint, newBFtype = bfPlotType, oldBFtype = "BF10"))
dfPoints <- data.frame(
x = xPoint,
y = logYPoint,
g = legendText1,
label1 = jaspGraphs::parseThis(legendText1),
label2 = jaspGraphs::parseThis(legendText2),
stringsAsFactors = FALSE
)
}
plotResult <- try(jaspGraphs::PlotRobustnessSequential(
dfLines = dfLines,
xName = bquote(paste(.(gettext("Stretched beta prior width")), ~kappa)),
dfPoints = dfPoints,
bfType = bfPlotType,
pointColors = pointColors,
pointFill = pointFill,
hypothesis = hypothesisjaspGraphsName,
plotLineOrPoint = "line"
))
return(plotResult)
}
.drawBfSequentialPlotCorBayes <- function(v1, v2, bfObject, options) {
.setSeedJASP(options)
# TODO(Alexander): This can be removed once bstats is adapted
alternativeLocal <- switch(options[["alternative"]],
twoSided = "two.sided",
greater = "greater",
less = "less")
sequentialValuesPerHyp <- bstats::computeCorSequentialLine(x=v1, y=v2, bfObject=bfObject)
sequentialValues <- bstats::getSidedObject(sequentialValuesPerHyp, alternative=alternativeLocal,
itemNames="nDomain")
sequentialLine <- sequentialValues[["sequentialLine"]]
if (is.character(sequentialLine)) {
return(sequentialLine)
} else {
y <- log(sequentialLine)
if (options[["bayesFactorType"]]=="BF01")
y <- -y
dfLines <- data.frame(
x = sequentialValues[["nDomain"]],
y = y
)
BF <- NULL
hypothesis <- switch(options[["alternative"]],
"twoSided"="equal",
"greater"="greater",
"less"="smaller"
)
bfType <- switch(options[["bayesFactorType"]],
"LogBF10"="BF10",
"BF10"="BF10",
"BF01"="BF01"
)
# TODO(Alexander): This can be removed once bstats is adapted
alternativeLocal <- switch(options[["alternative"]],
twoSided = "two.sided",
greater = "greater",
less = "less")
if (options[["bfSequentialPlotAdditionalInfo"]]) {
BF <- bfObject[[alternativeLocal]][["bf"]]
BF <- .recodeBFtype(BF, newBFtype = bfType, oldBFtype = "BF10")
}
plotResult <- try(jaspGraphs::PlotRobustnessSequential(
dfLines = dfLines,
xName = gettext("n"),
BF = BF,
bfType = bfType,
hypothesis = hypothesis
))
}
return(plotResult)
}
.checkAndSetPlotCorBayes <- function(triedPlot, jaspPlotResult) {
if (isTryError(triedPlot)) {
jaspPlotResult$setError(.extractErrorMessage(triedPlot))
} else if (is.character(triedPlot)) {
jaspPlotResult$setError(triedPlot)
} else {
jaspPlotResult$plotObject <- triedPlot
}
}
.computeCorCredibleInterval <- function(bfObject, ciLevel, method) {
if (method == "pearson") {
return(bstats::computePearsonCredibleInterval("betaA"=bfObject[["betaA"]], "betaB"=bfObject[["betaB"]],
"ciValue"=ciLevel))
} else if (method == "kendall") {
return(bstats::computeKendallCredibleInterval("n"=bfObject[["two.sided"]][["n"]], "tauObs"=bfObject[["two.sided"]][["stat"]],
"kappa"=bfObject[["kappa"]], "var"=1, "ciValue"=ciLevel,
"h0"=bfObject[["h0"]]))
} else if (method == "spearman") {
# TODO(Johnny):
print("NO THIS IS NOT IT, STILL NEED TO MAKE THIS, THIS IS JUST A PLACEHOLDER")
return(bstats::computePearsonCredibleInterval("betaA"=bfObject[["betaA"]], "betaB"=bfObject[["betaB"]],
"ciValue"=ciLevel))
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.