contr.bayes.R
#' Orthonormal Contrast Matrices for Bayesian Estimation
#'
#' Returns a design or model matrix of orthonormal contrasts such that the
#' marginal prior on all effects is identical. Implementation from Singmann
#' \& Gronau's \href{https://github.com/bayesstuff/bfrms/}{\code{bfrms}},
#' following the description in Rouder, Morey, Speckman, \& Province (2012, p. 363).
#'
#' Though using this factor coding scheme might obscure the interpretation of
#' parameters, it is essential for correct estimation of Bayes factors for
#' contrasts and multi-level order restrictions. See info on specifying correct
#' priors for factors with more than 2 levels in
#' \href{https://easystats.github.io/bayestestR/articles/bayes_factors.html}{the Bayes factors vignette}.
#'
#' @param n a vector of levels for a factor, or the number of levels.
#' @param contrasts logical indicating whether contrasts should be computed.
#'
#' @references Rouder, J. N., Morey, R. D., Speckman, P. L., \& Province, J. M.
#' (2012). Default Bayes factors for ANOVA designs. *Journal of Mathematical
#' Psychology*, 56(5), 356-374. https://doi.org/10.1016/j.jmp.2012.08.001
#'
#' @return A \code{matrix} with n rows and k columns, with k=n-1 if contrasts is
#' \code{TRUE} and k=n if contrasts is \code{FALSE}.
#'
#' @examples
#' \dontrun{
#' contr.bayes(2) # Q_2 in Rouder et al. (2012, p. 363)
#' # [,1]
#' # [1,] -0.7071068
#' # [2,] 0.7071068
#'
#' contr.bayes(5) # equivalent to Q_5 in Rouder et al. (2012, p. 363)
#' # [,1] [,2] [,3] [,4]
#' # [1,] 0.0000000 0.8944272 0.0000000 0.0000000
#' # [2,] 0.0000000 -0.2236068 -0.5000000 0.7071068
#' # [3,] 0.7071068 -0.2236068 -0.1666667 -0.4714045
#' # [4,] -0.7071068 -0.2236068 -0.1666667 -0.4714045
#' # [5,] 0.0000000 -0.2236068 0.8333333 0.2357023
#'
#' ## check decomposition
#' Q3 <- contr.bayes(3)
#' Q3 %*% t(Q3)
#' # [,1] [,2] [,3]
#' # [1,] 0.6666667 -0.3333333 -0.3333333
#' # [2,] -0.3333333 0.6666667 -0.3333333
#' # [3,] -0.3333333 -0.3333333 0.6666667
#' ## 2/3 on diagonal and -1/3 on off-diagonal elements
#' }
#' @export
contr.bayes <- function(n, contrasts = TRUE) {
# validate n
if (length(n) <= 1L) {
if (is.numeric(n) && length(n) == 1L && n > 1L) {
TRUE # all good
} else {
stop("not enough degrees of freedom to define contrasts")
}
} else {
n <- length(n)
}
# make factor coding
cont <- diag(n)
if (contrasts) {
a <- n
I_a <- diag(a)
J_a <- matrix(1, nrow = a, ncol = a)
Sigma_a <- I_a - J_a / a
cont <- eigen(Sigma_a)$vectors[, seq_len(a - 1), drop = FALSE]
}
cont
}