#' Visualize a prior on Cohen's d for a Bayesian t-test
#'
#' @param alternative A function object. The default is a Cauchy prior
#' with scaling parameter `sqrt(2) / 2` as is the default in package
#' `BayesFactor` (Morey & Rouder, 2015). This argument can also be a
#' scalar number, in which case it is assumed that the alternative
#' is a point hypothesis on Cohen's d (with Cohen's d = `prior`).
#' @param null Boolean. Plot a line on x = 0 representing the null
#' hypothesis?
#' @param from the left margin of the x-axis
#' @param to the right margin of the x-axis
#' @param xlab the label of the x-axis
#' @param frame.plot Should a frame be drawn?
#' @param col The color of the curve of the alternative hypothesis.
#' @param lwd The line width
#' @param ... additional parameters passed to function `curve`
#'
#' @examples
#' ## Standard cauchy prior:
#' visualize_prior()
#' ## Wide cauchy prior:
#' visualize_prior(function(x) dcauchy(x, scale = 1))
#' ## Normal prior with M = 0 and SD = 0.4
#' visualize_prior(function(x) dnorm(x, 0, 0.4))
#' ## Point prior on Cohen's d = 0.4
#' visualize_prior(0.4)
#'
#' @references
#'
#' Morey, R. D., & Rouder, J. N. (2015). BayesFactor: Computation of
#' bayes factors for common designs. Retrieved from
#' https://CRAN.R-project.org/package=BayesFactor
#'
#' Rouder, J. N., Speckman, P. L., Sun, D., Morey, R. D., & Iverson,
#' G. (2009). Bayesian t tests for accepting and rejecting the null
#' hypothesis. Psychonomic Bulletin & Review, 16(2), 225–237.
#'
#' @export
#'
visualize_prior <- function(alternative = function(x) dcauchy(x, scale = sqrt(2) / 2),
null = TRUE, from = -4, to = 4,
xlab = "Cohen's d", frame.plot = FALSE,
col = "red", lwd = 3, ...) {
if (class(alternative) != "function") {
plot(c(alternative, alternative), c(0, 1), type = "l",
las = 1, xlim = c(from, to), xlab = xlab,
frame.plot = frame.plot, yaxt = "n", col = "transparent",
ylab = "")
abline(v = alternative, col = col, lwd = lwd)
} else {
curve(alternative, las = 1, from = from, to = to, xlab = xlab,
frame.plot = frame.plot, yaxt = "n", col = col, ylab = "",
lwd = lwd, ...)
}
if (null) {
abline(v = 0, col = "#005A31", lwd = lwd)
legend("topleft", lwd = lwd, legend = c("Null", "Alternative"),
cex = 1.2, col = c("#005A31", "red"), lty = 1, bty = "n")
}
}
#' Visualize hypotheses' predictions for a Bayesian two-group t-test
#'
#' Visualizes the predictions of the null and an alternative hypothesis.
#' By specifying an observed effect (observed t-value), an illustration
#' of the ratio of the marginal likelihoods, i.e., an
#' illustration of the Bayes factor is displayed.
#'
#' @param alternative A function object. The default is a Cauchy prior
#' with scaling parameter `sqrt(2) / 2` as is the default in package
#' `BayesFactor` (Morey & Rouder, 2015). This argument can also be a
#' scalar number, in which case it is assumed that the alternative
#' is a point hypothesis on Cohen's d (with Cohen's d = `prior`).
#' If `alternative` is NULL, only the null hypothesis predictions
#' are drawn.
#' @param n1 The sample size in group 1
#' @param n2 The sample size in group 2
#' @param from the left margin of the x-axis
#' @param to the right margin of the x-axis
#' @param col Specify the coloring of the plot. col[1] specifies the
#' color for the null hypothesis; col[2] specifies the color for the
#' alternative hypothesis; col[3] specifies the color for the line
#' the illustrates the observed effect (only has an effect if an
#' observed effect size is specified via argument `observed_t`);
#' col[4] specifies the color of the points at the intersection of
#' the observed effect and the curves for the hypotheses (only has
#' an effect if an observed effect size is specified via argument
#' `observed_t`)
#' @param observed_t An observed t-value. Does not need to be specified
#' and defaults to NULL. If this argument is passed, the marginal
#' likelihoods of the null and the alternative hypothesis are drawn
#' and a Bayes factor is displayed.
#' @param frame.plot Should a frame be drawn?
#' @param xlab The label of the x-axis
#' @param BFx The x coordinate where the Bayes factor is displayed in
#' the plot. Only has an effect if an `observed_t` is passed.
#' @param BFy The y coordinate where the Bayes factor is displayed in
#' the plot. Only has an effect if an `observed_t` is passed.
#' @param BF10 Boolean; defaults to TRUE. Should a Bayes factor be
#' displayed quantifying the evidence in favor of the alternative
#' hypothesis? The BF01 is displayed if `BF10` is `FALSE`. Only has
#' an effect if an `observed_t` is passed.
#' @param legend_placement A keyword passed to `legend` to indicate
#' where the legend has to be placed. Defaults to "topleft". Only
#' has an effect if an `observed_t` is passed.
#' @param legend_content Specify the content of the legend.
#' legend_content[1] refers to the label for the null hypothesis;
#' legend_content[2] refers to the label for the alternative hypothesis;
#' legend_content[3] refers to the label for the observed effect;
#' (only has an effect if an
#' observed effect size is specified via argument `observed_t`);
#' @param ... additional parameters passed to `curve`
#'
#' @references
#'
#' Morey, R. D., & Rouder, J. N. (2015). BayesFactor: Computation of
#' Bayes factors for common designs. Retrieved from
#' https://CRAN.R-project.org/package=BayesFactor
#'
#' Rouder, J. N., Speckman, P. L., Sun, D., Morey, R. D., & Iverson,
#' G. (2009). Bayesian t tests for accepting and rejecting the null
#' hypothesis. Psychonomic Bulletin & Review, 16(2), 225–237.
#'
#' @examples
#' ## Using the default cauchy prior
#' visualize_predictions(n1 = 100, n2 = 100, observed_t = 3)
#' ## Using a wide Cauchy prior
#' visualize_predictions(function(x) dcauchy(x, scale = 1),
#' n1 = 100, n2 = 100, observed_t = 3)
#' ## Using a normal prior (M = 0, SD = 0.4)
#' visualize_predictions(function(x) dnorm(x, 0, 0.4),
#' n1 = 100, n2 = 100, observed_t = 3, lty = c(1, 1, 2))
#'
#' @export
#'
visualize_predictions <- function(alternative = function(x) dcauchy(x, scale = sqrt(2) / 2),
n1, n2, from = -6, to = 6, lwd = 3,
col = c("#005A31", "red", "black", "orange"),
observed_t = NULL, frame.plot = FALSE,
xlab = "t-value", BFx = from + 0.1,
BFy = 0.1, BF10 = TRUE,
legend_placement = "topleft",
legend_content = c("Null", "Alternative", "Observed"),
...) {
## 1. Draw null predictions
nullhyp <- function(x) dt(x, n1 + n2 - 2)
curve(nullhyp, from = from , to = to, col = col[1], lwd = lwd[1],
frame.plot = frame.plot, yaxt = "n", ylab = "", xlab = xlab, ...)
maxdensity <- nullhyp(0)
## Prepare the legend
options <- legend_content[1]
lty <- 1
col_legend <- col[1]
## 2. Draw alternative hypothesis predictions
if (!is.null(alternative)) {
predictions_alt <- sample_likelihoods(alternative, n1, n2, from, to)
points(predictions_alt$x, predictions_alt$y, type = "l", col = col[2],
lwd = lwd)
maxdensity <- max(maxdensity, predictions_alt$y)
## For legend:
options <- c(options, legend_content[2])
lty <- c(lty, 1)
col_legend <- c(col_legend, col[2])
}
## 3. Illustrate the observed effect and the Bayes factor
if (!is.null(observed_t)) {
null_likelihood <- nullhyp(observed_t)
lines(x = rep(observed_t, 2), y = c(0, maxdensity),
col = col[3], lwd = lwd, lty = 2)
points(observed_t, null_likelihood, pch = 19, col = col[4], cex = 2)
if (!is.null(alternative)) {
alt_likelihood <- marginal_likelihood(observed_t, alternative, n1, n2)
BF <- alt_likelihood / null_likelihood
displayBF <- "BF10 = "
if (BF10 == FALSE) {
BF <- 1 / BF
displayBF <- "BF01 = "
}
points(observed_t, alt_likelihood, pch = 19, col = col[4], cex = 2)
max_likelihood <- max(null_likelihood, alt_likelihood)
text(x = BFx, y = BFy, pos = 4,
labels = paste(displayBF, round(BF, 2)))
}
options <- c(options, legend_content[3])
lty <- c(lty, 2)
col_legend <- c(col_legend, col[3])
}
legend(legend_placement, lwd = 3, legend = options,
cex = 1.2, col = col_legend, lty = lty, bty = "n")
}
#' Generate marginal likelihoods in a given interval of t-values
#'
#' @param alternative A function object. The default is a Cauchy prior
#' with scaling parameter `sqrt(2) / 2` as is the default in package
#' `BayesFactor` (Morey & Rouder, 2015). This argument can also be a
#' scalar number, in which case it is assumed that the alternative
#' is a point hypothesis on Cohen's d (with Cohen's d = `prior`).
#' @param n1 The sample size in group 1
#' @param n2 The sample size in group 2
#' @param from The lowest t-value
#' @param to The largest t-value
#'
#' @return A data.frame of two columns: Column x contains the t-values,
#' column y contains the marginal likelihoods
#'
#' @examples
#' sample_likelihoods(n1 = 30, n2 = 30, from = -3, to = 3)
#' sample_likelihoods(function(x) dnorm(x, 0, 0.3), n1 = 30, n2 = 30)
#'
#' @export
#'
sample_likelihoods <- function(alternative = function(x) dcauchy(x, scale = sqrt(2) / 2),
n1, n2, from = -6, to = 6) {
observed_ts <- seq(from, to, .01) # hypothetical observed t-values
## Allow for a point null hypothesis
if (class(alternative) != "function") {
althyp <- function(x) dt(x, df = n1 + n2 - 2, ncp = ncp_from_d(alternative, n1, n2))
predictions <- althyp(observed_ts)
return(data.frame(x = observed_ts, y = predictions))
}
predictions <- vector(length = length(observed_ts)) # predict marginal likelihood of each t
for (i in 1:length(predictions)) {
predictions[i] <- marginal_likelihood(observed_ts[i], alternative, n1, n2)
}
return(data.frame(x = observed_ts, y = predictions))
}
#' Compute the marginal likelihood of an observed t-value
#'
#' This function can be used to compute a Bayes factor as the ratio
#' of two marginal likelihoods.
#'
#' @param observed_t The observed t-value.
#' @param prior A function object describing the a-priori plausibility
#' of all effect sizes. The default is a Cauchy prior as in package
#' BayesFactor (Morey & Rouder, 2015). Can also be a scalar, in
#' which case it is assumed that the alternative is a point
#' hypothesis on Cohen's d (with value `prior`).
#' @param n1 The sample size in group 1
#' @param n2 The sample size in group 2
#'
#' @return The marginal likelihood of the observed t-value given the
#' prior on Cohen's d.
#'
#' @examples
#' ## Compute Bayes factor for standard Cauchy prior:
#' n1 <- 100
#' n2 <- 100
#' sample1 <- rnorm(n1, 0.2, 1)
#' sample2 <- rnorm(n1, 0.2, 1)
#' tvalue <- t.test(sample1, sample2)$statistic
#' ml_alt <- marginal_likelihood(tvalue, n1 = n1, n2 = n2)
#' bf10 <- ml_alt / dt(tvalue, n1 + n2 - 2)
#' # Compare:
#' BayesFactor::ttestBF(sample1, sample2)
#'
#'
#' ## Normal prior - N(0, 0.3)
#' marginal_likelihood(3, function(x) dnorm(x, 0, 0.3), n1 = 30, n2 = 30)
#'
#' @export
#'
marginal_likelihood <- function(observed_t,
prior = function(x) dcauchy(x, scale = sqrt(2) / 2),
n1, n2) {
## allow for point alternative:
if (class(prior) != "function") {
althyp <- function(x) dt(x, df = n1 + n2 - 2, ncp = ncp_from_d(prior, n1, n2))
return(althyp(observed_t))
}
integrate(weighted_likelihood, lower = -Inf, upper = Inf,
observed_t = observed_t, prior = prior, n1 = n1, n2 = n2)$value
}
#' Compute the weighted likelihood of an observed t-value
#'
#' For an assumed true effect size (Cohen's d), compute the likelihood
#' of an observed t-value weighted by the a priori plausbility of this
#' effect size. This function is primarily used to compute the marginal
#' likelihood of the effect size given a continuous hypothesis (e.g.,
#' a Cauchy prior) by integrating in function `marginal_likelihood`
#' below.
#'
#' @param true_d A one-element vector, representing an assumed true
#' effect size d, conditioned on which the likelihood is computed.
#' @param observed_t The observed t-value.
#' @param prior A function object describing the a-priori plausibility
#' of all effect sizes. The default is a cauchy prior
#' as in package BayesFactor (Morey & Rouder, 2015)
#' @param n1 The sample size in group 1
#' @param n2 The sample size in group 2
#'
#' @return The weighted likelihood of the observed t-value.
#'
#' @export
#'
weighted_likelihood <- function(true_d, observed_t,
prior = function(x) dcauchy(x, scale = sqrt(2) / 2),
n1, n2) {
dfr <- n1 + n2 - 2 ## degrees of freedom for two-sample t-test
ncp = ncp_from_d(true_d, n1, n2)
return(dt(observed_t, dfr, ncp) * prior(true_d))
}
#' Compute non-centrality parameter for the non-central t distribution
#'
#' @param d Cohen's d
#' @param n1 Sample size in group 1
#' @param n2 Sample size in group 2
#'
#' @details See page 7 in Erdfelder, Faul, & Buchner (1996) for the
#' formula.
#'
#' @return The non-centrality parameter
#'
#' @references
#' Erdfelder, E., Faul, F., & Buchner, A. (1996). GPOWER: A general
#' power analysis program. Behavior research methods, instruments, &
#' computers, 28(1), 1-11.
#'
#' @examples
#' ncp_from_d(0.3, 50, 50)
#'
#' @export
#'
ncp_from_d <- function(d, n1, n2) {
ncp <- d * sqrt((n1 * n2^2 + n2 * n1^2) / ((n1 + n2)^2))
return(ncp)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.