#Project: lmhyp - Informed hypothesis testing for regression
#Script purpose: Function to examine complex hypothesis for lm objects with a minimal prior and BF as output
#Code: Anton Ohlsson Collentine
#*************************************
#Hypothesis testing function----
#*************************************
#'Testing Informed Hypotheses
#'
#'Test competing hypotheses about coefficients in an \code{\link{lm}}-object.
#'
#'This function is based on a method by Mulder (2014), a modification of the
#'fractional Bayes factor approach. In essence, it uses a number of observations
#'equal to the number of predictors in the \code{lm} model to construct a
#'minimally informative prior, and the remainder of the observations are then
#'used to test the hypotheses.
#'
#'Hypotheses are specified using the variable names from the \code{lm} object.
#'If a hypothesis involves multiple variables it is usually preferable to
#'standardize relevant variables before fitting the model with \code{lm} to
#'facilitate interpretation. This is done simply by substracting the mean of a
#'variable from each observation and dividing by the standard deviation. A
#'simple option for achieving this is to use the \code{\link{scale}} function.
#'
#'Multiple hypotheses can be specified at the same time by separating them with
#'a semicolon. It is advisable to only specify competing hypotheses in this way,
#'that is, hypotheses regarding the same variables, e.g., \dQuote{X1 > 0; X1 <
#'0; X1 = 0}. If specifying multiple hypotheses and comparing against a value it
#'is currently only possible to compare against the same value, i.e., \dQuote{X1
#'= 0; X1 = 2} is not functional input. This is because the prior is centered
#'around the input value (or zero if no input value), which is not possible in
#'the case of several comparison values.
#'
#'Parentheses can be used to compare multiple variables with the same variable
#'or value. For example, \dQuote{(X1, X2) > X3} is read as \dQuote{X1 > X3 and
#'X2 > X3}. Each variable should only be specified once in a single hypothesis.
#'
#'An alternative to specifying hypotheses is to input the string
#'\dQuote{exploratory}. This will compare the hypotheses \dQuote{X < 0; X = 0; X
#'> 0} for all independent variables in the regression object, including the
#'intercept.
#'
#'For each specified hypothesis the posterior probability is output. If the
#'hypotheses are not exhaustive (i.e., do not cover the entire parameter space)
#'this includes the posterior probability of the complement to the input
#'hypotheses. The complement is the hypothesis that neither of the input
#'hypotheses is true. For example, inputting \dQuote{X1 > 0; X1 < 0} gives
#'posterior probabilities for only for these hypotheses, whereas inputting
#'\dQuote{(X1, X2) > 0} gives posterior probabilities for \dQuote{(X1, X2) > 0}
#'and \dQuote{not (X1, X2) > 0}.
#'
#'If not using the \dQuote{exploratory} option, it is possible to specify prior
#'probabilities for the input hypotheses. By default these are equal (priorprob
#'= 1). Prior probabilites can both be input as probabilites, e.g., c(0.2, 0.3,
#'0.5) or relative weights of each hypothesis, e.g, c(2, 3, 5). If the input
#'probabilites do not sum to 1 they will simply be normalized. Prior
#'probabilities must be specified for all hypotheses, including the complement
#'if one exists.
#'
#'By saving the test as an object it is also possible to access the
#'\code{BF_matrix} which compares the hypotheses directly against each other
#'(see examples). This matrix divides the row hypothesis by each column
#'hypothesis and, assuming equal prior probabilities, can be interpreted as
#'\dQuote{given the data, [row hypothesis] is [value] times as likely as [column
#'hypothesis]}.
#'
#'Supplementary output intended to provide a deeper understanding of the
#'underlying method and primary output can be printed when the test has been
#'saved as an object. Calling \code{BF_computation} illustrates the computation
#'of the Bayes factor of each hypothesis against the unconstrained hypothesis.
#'In the output \code{c(E)} is the prior density, \code{c(I|E)} the prior
#'probability, \code{c} the product of these two, and columns prefixed with "f"
#'the equivalent for the posterior. \code{B(t,u)} is the Bayes factor of
#'hypothesis t against the unconstrained (u) and \code{PP(t)} is the
#'posterior probability of t. Finally, calling \code{BFu_CI} provides 90\%
#'credibility intervals for those few cases where the Bayes factor was
#'calculated numerically.
#'
#'@section References: Mulder, J. (2014). Prior adjusted default Bayes factors
#' for testing (in) equality constrained hypotheses. Computational Statistics &
#' Data Analysis, 71, 448-463.
#'
#'@examples
#'###Standardize variables and fit the linear model
#'dt <- as.data.frame(scale(mtcars[, c(1, 3:4, 6)]))
#'fit <- lm(mpg ~ disp + hp + wt, data = dt)
#'
#'###Exploratory analysis
#'test_hyp(fit, "exploratory")
#'
#'###Define hypotheses based on theory and test them
#'hyp <- "(wt, hp) > disp > 0; (wt, hp) > disp = 0"
#'res <- test_hyp(fit, hyp)
#'res
#'
#'###Bayes factor comparison of hypotheses
#'res$BF_matrix
#'
#'@param object A regression model object fit using the \code{lm} function.
#'@param hyp A string specifying hypotheses to be tested using the variable
#' names of the \code{lm} object, or the string \dQuote{exploratory}, see
#' details.
#'@param priorprob Vector of prior probabilities for the input hypotheses, by
#' default equal.
#'@param mcrep Integer specifying the number of iterations if no analytical
#' solutions is possible. This is rare and only the case if the rank of the
#' constraint matrix is less than its number of rows.
#'
#'@export test_hyp
test_hyp <- function(object, hyp, priorprob = 1, mcrep = 1e6){
#1) initial setup and checks of input----
varnames <- variable.names(object)
if(is.null(varnames)) stop("Please input proper linear model object")
varnames <- gsub("(\\(Intercept\\))", "Intercept", varnames)
betahat <- object$coefficients
k <- length(varnames)
n <- length(object$fitted.values)
b <- (k + 1) / n
if(hyp == "exploratory"){
if(!isTRUE(priorprob == 1)) stop('"exploratory" option requires priorprobs = 1 (default) ')
hyp <- vector("list", length(varnames))
BF_matrices <- vector("list", length(varnames))
list_post_prob <- vector("list", length(varnames))
for(IV in seq_along(varnames)){
hyp[[IV]] <- paste(varnames[IV], "< 0;", varnames[IV], "= 0;", varnames[IV], "> 0", sep = " ")
}
}
for(loop in seq_along(hyp)){ #If 'exploratory' option loop over all variables, else specified hyp
hyp2 <- gsub("[ \n]", "", hyp[[loop]])
hyp2 <- gsub("(\\(Intercept\\))", "Intercept", hyp2)
if(!grepl("^[0-9a-zA-Z><=,;().-]+$", hyp2)) stop("Impermissable characters in hypotheses.")
if(grepl("[><=]{2,}", hyp2)) stop("Do not use combined comparison signs e.g., '>=' or '=='")
step1 <- unlist(strsplit(hyp2, split = "[<>=,;()]"))
input_vars <- step1[grep("[a-zA-Z]+", step1)]
if(!all(input_vars %in% varnames)) stop("Hypothesis variable(s) not in object, check spelling")
hyp_out <- if(length(hyp) > 1) c("X < 0", "X = 0", "X > 0") else unlist(strsplit(hyp2, split = ";"))
for(no in seq_along(hyp_out)){names(hyp_out)[no] <- paste0("H", no)}
if(!isTRUE(priorprob == 1) && length(priorprob) < length(hyp_out)) stop("Use default equal priorprob or specify priorprob for all hypotheses")
hyps <- unlist(strsplit(hyp2, split = ";"))
BFu <- out_c_E <- out_f_E <- out_c_i_e <- out_f_i_e <- out_ci_lb <- out_ci_ub <- rep(NA, length = length(hyps))
BFip_posterior <- if(any(!grepl("=", hyps))) {rep(NA, sum(!grepl("=", hyps)))} else{NULL}
if(!is.null(BFip_posterior)) {
R_i_all <- r_i_all <- BFi_draws_post <- BFi_draws_prior <- vector("list", length = length(BFip_posterior))
ineq_marker <- 0
BFip_prior <- BFip_posterior
}
for(h in seq_along(hyps)){
hyp2 <- hyps[[h]]
step2 <- unlist(strsplit(hyp2, split = "[<>=,()]"))
hyp_vars <- step2[grep("[a-zA-Z]+", step2)]
if(any(duplicated(hyp_vars))) stop("Variables should occur only once in a hypothesis. Check semicolons.")
#2)hyp-to-matrices----
framer <- function(x){
pos_comparisons <- unlist(gregexpr("[<>=]", x))
leftside <- rep(NA, length(pos_comparisons) + 1)
rightside <- rep(NA, length(pos_comparisons) + 1)
pos1 <- c(-1, pos_comparisons)
pos2 <- c(pos_comparisons, nchar(x) + 1)
for(i in seq_along(pos1)){
leftside[i] <- substring(x, pos1[i] + 1, pos1[i+1] - 1)
rightside[i] <- substring(x, pos2[i] + 1, pos2[i+1] - 1)
}
leftside <- leftside[-length(leftside)]
rightside <- rightside[-length(rightside)]
comparisons <- substring(x, pos_comparisons, pos_comparisons)
data.frame(left = leftside, comp = comparisons, right = rightside, stringsAsFactors = FALSE)
}
framed <- framer(hyp2)
if(any(grepl(",", framed$left)) || any(grepl(",", framed$right))){
if(nrow(framed) > 1){
for(r in 1:(nrow(framed)-1)){
if(all.equal(framed$right[r], framed$left[r+1])){
if(substring(framed$right[r], 1, 1) == "(") {
framed$right[r] <- sub("),.+", ")", framed$right[r])
framed$left[r+1] <- sub(".+),", "", framed$left[r +1])
} else{
framed$right[r] <- sub(",.+", "", framed$right[r])
framed$left[r+1] <- sub("[^,]+,", "", framed$left[r+1])
}
}
}
}
commas_left <- framed$left[grep(",", framed$left)]
commas_right <- framed$right[grep(",", framed$right)]
if(isTRUE(any(!grepl("\\(.+)", commas_left))) || isTRUE(any(!grepl("\\(.+)", commas_right))) ||
isTRUE(any(grepl(").+", commas_left))) || isTRUE(any(grepl(").+", commas_right))) ||
isTRUE(any(grepl(".+\\(", commas_left))) || isTRUE(any(grepl(".+\\(", commas_right)))) {
stop("Incorrect hypothesis syntax or extra character, check specification")
}
framed$left <- gsub("[()]", "", framed$left)
framed$right <- gsub("[()]", "", framed$right)
commas <- unique(c(grep(",", framed$left), grep(",", framed$right)))
if(length(commas) > 0){
multiples <- vector("list", length = length(commas))
for(r in seq_along(commas)){
several <- framed[commas,][r, ]
if(several$comp == "="){
several <- c(several$left, several$right)
separate <- unlist(strsplit(several, split = ","))
if(any(grepl("^$", several))) stop("Misplaced comma in hypothesis")
converted_equality <- paste(separate, collapse = "=")
multiples[[r]] <- framer(converted_equality)
} else{
leftvars <- unlist(strsplit(several$left, split = ","))
rightvars <- unlist(strsplit(several$right, split = ","))
if(any(grepl("^$", leftvars)) || any(grepl("^$", rightvars))) stop("Misplaced comma in hypothesis")
left <- rep(leftvars, length.out = length(rightvars)*length(leftvars))
right <- rep(rightvars, each = length(leftvars))
comp <- rep(several$comp, length(left))
multiples[[r]] <- data.frame(left = left, comp = comp, right = right, stringsAsFactors = FALSE)
}
}
framed <- framed[-commas,]
multiples <- do.call(rbind, multiples)
framed <- rbind(multiples, framed)
}
}
equality <- framed[framed$comp == "=",]
inequality <- framed[!framed$comp == "=",]
#****Equality part string-to-matrix
if(nrow(equality) == 0) {
R_e <- r_e <- NULL
} else{
outcomes <- suppressWarnings(apply(equality[, -2], 2, as.numeric))
outcomes <- matrix(outcomes, ncol = 2, byrow = FALSE)
if(any(rowSums(is.na(outcomes)) == 0)) stop("Value compared with value rather than variable, e.g., '2 = 2', check hypotheses")
rows <- which(rowSums(is.na(outcomes)) < 2)
specified <- t(outcomes[rows,])
specified <- specified[!is.na(specified)]
r_e <- ifelse(rowSums(is.na(outcomes)) == 2, 0, specified)
r_e <- matrix(r_e)
var_locations <- apply(equality[, -2], 2, function(x) ifelse(x %in% varnames, match(x, varnames), 0))
var_locations <- matrix(var_locations, ncol = 2)
R_e <- matrix(rep(0, nrow(equality)*length(varnames)), ncol = length(varnames))
for(i in seq_along(r_e)){
if(!all(var_locations[i, ] > 0)){
R_e[i, var_locations[i,]] <- 1
} else{
R_e[i, var_locations[i,]] <- c(1, -1)
}
}
}
#****Inequality part string-to-matrix
if(nrow(inequality) == 0) {
R_i <- r_i <- NULL
} else{
outcomes <- suppressWarnings(apply(inequality[, -2], 2, as.numeric))
outcomes <- matrix(outcomes, ncol = 2, byrow = FALSE)
if(any(rowSums(is.na(outcomes)) == 0)) stop("Value compared with value rather than variable, e.g., '2 > 2', check hypotheses")
cols <- which(rowSums(is.na(outcomes)) < 2)
specified <- t(outcomes[cols,])
specified <- specified[!is.na(specified)]
r_i <- ifelse(rowSums(is.na(outcomes)) == 2, 0, specified)
r_i <- matrix(r_i)
leq <- which(inequality$comp == "<")
var_locations <- apply(inequality[, -2], 2, function(x) ifelse(x %in% varnames, match(x, varnames), 0))
var_locations <- matrix(var_locations, ncol = 2)
R_i <- matrix(rep(0, nrow(inequality)*length(varnames)), ncol = length(varnames))
for(i in seq_along(r_i)){
if(!all(var_locations[i, ] > 0)){
if(var_locations[i, 1] == 0){
if(i %in% leq){
value <- 1
} else{
r_i[i] <- r_i[i]*-1
value <- -1
}
} else{
if(i %in% leq){
r_i[i] <- r_i[i]*-1
value <- -1
} else{
value <- 1
}
}
R_i[i, var_locations[i,]] <- value
} else{
value <- if(i %in% leq) c(-1, 1) else c(1, -1)
R_i[i, var_locations[i,]] <- value
}
}
}
#3)check comparisons----------------
if(is.null(R_i)){
comparisons <- "only equality"
} else if(is.null(R_e)){
comparisons <- "only inequality"
} else{
comparisons <- "both comparisons"
}
#set prior mean
R_ei <- rbind(R_e,R_i)
r_ei <- rbind(r_e,r_i)
Rr_ei <- cbind(R_ei,r_ei)
beta_zero <- MASS::ginv(R_ei)%*%r_ei
if(nrow(Rr_ei) > 1){
rref_ei <- pracma::rref(Rr_ei)
nonzero <- rref_ei[,k+1]!=0
if(max(nonzero)>0){
row1 <- max(which(nonzero==T))
if(sum(abs(rref_ei[row1,1:k]))==0){
stop("Default prior mean cannot be constructed from constraints.")
}
}
}
if(comparisons == "only equality"){
#**3.1)only-equality----
delta <- R_e %*% betahat
delta_zero <- R_e %*% beta_zero
#Scale matrix components
X <- model.matrix(object)
RX <- R_e %*% solve((t(X) %*% X)) %*% t(R_e)
s2 <- sum((model.frame(object)[[1]] - X %*% betahat)^2)
#Scale matrix for posterior t-distribution
scale_post <- matrix(s2 * RX / (n - k), ncol = nrow(R_e))
#Scale matrix for prior t-distribution
scale_prior <- matrix(s2 * RX / (n*b - k), nrow(R_e))
#Hypothesis test
log_BF <- mvtnorm::dmvt(x = t(r_e), delta = delta, sigma = scale_post, df = n - k, log = TRUE) -
mvtnorm::dmvt(x = t(r_e), delta = delta_zero, sigma = scale_prior, df = n*b - k, log = TRUE)
f_E <- mvtnorm::dmvt(x = t(r_e), delta = delta, sigma = scale_post, df = n - k, log = FALSE)
c_E <- mvtnorm::dmvt(x = t(r_e), delta = delta_zero, sigma = scale_prior, df = n*b - k, log = FALSE)
c_i_e <- f_i_e <- NA
BF <- exp(log_BF)
ci_lb <- ci_ub <- NA
} else if(comparisons == "only inequality"){
#**3.2)only-inequality----
ineq_marker <- ineq_marker + 1
r_i <- as.vector(r_i)
if(Matrix::rankMatrix(R_i)[[1]] == nrow(R_i)){
delta <- as.vector(R_i %*% betahat)
delta_zero <- as.vector(R_i %*% beta_zero)
#Scale matrix components
X <- model.matrix(object)
RX <- R_i %*% solve(t(X) %*% X) %*% t(R_i)
s2 <- sum((model.frame(object)[[1]] - X %*% betahat)^2)
#Scale matrix for posterior t-distribution
scale_post <- s2 * RX / (n - k)
#Scale matrix for prior t-distribution
scale_prior <- s2 * RX / (n*b - k)
if(nrow(scale_post) == 1){ #If univariate
prior_prob <- pt((r_i - delta_zero) / sqrt(scale_prior), df = n*b - k, lower.tail = FALSE)[1]
posterior_prob <- pt((r_i - delta) / sqrt(scale_post), df = n - k, lower.tail = FALSE)[1]
BF <- posterior_prob / prior_prob #prior
} else { #if multivariate
prior_prob <- mvtnorm::pmvt(lower = r_i, upper = Inf, delta = delta_zero, sigma = scale_prior, df = n*b - k, type = "shifted")[1]
posterior_prob <- mvtnorm::pmvt(lower = r_i, upper = Inf, delta = delta, sigma = scale_post, df = n - k, type = "shifted")[1]
BF <- posterior_prob / prior_prob
}
ci_lb <- ci_ub <- NA
ci_draws_post <- ci_draws_prior <- NULL
} else{ #No transformation is possible.
if(!is.numeric(mcrep) || !mcrep %% 1 == 0) stop("Input for mcrep should be an integer")
#Scale matrix for posterior t-distribution.
scale_post <- vcov(object)
#Scale matrix for prior t-distribution
scale_prior <- vcov(object) * (n - k) / (n*b - k)
draws_post <- mvtnorm::rmvt(n = mcrep, delta = betahat, sigma = scale_post, df = n - k)
draws_prior <- mvtnorm::rmvt(n = mcrep, delta = beta_zero, sigma = scale_prior, df = n*b - k)
prior_prob <- mean(apply(draws_prior%*%t(R_i) > rep(1, mcrep)%*%t(r_i), 1, prod))
posterior_prob <- mean(apply(draws_post%*%t(R_i) > rep(1, mcrep)%*%t(r_i), 1, prod))
BF <- posterior_prob / prior_prob
#Credibility interval
x_post <- sum(apply(draws_post%*%t(R_i) > rep(1, mcrep)%*%t(r_i), 1, prod))
x_prior <- sum(apply(draws_prior%*%t(R_i) > rep(1, mcrep)%*%t(r_i), 1, prod))
ci_draws_post <- rbeta(1e4, x_post, 1 + mcrep - x_post)
ci_draws_prior <- rbeta(1e4, x_prior, 1 + mcrep - x_prior)
BF_draws <- ci_draws_post / ci_draws_prior
ci_lb <- quantile(sort(BF_draws), 0.05)
ci_ub <- quantile(sort(BF_draws), 0.95)
}
c_i_e <- prior_prob
f_i_e <- posterior_prob
f_E <- c_E <- NA
BFip_prior[ineq_marker] <- prior_prob
BFip_posterior[ineq_marker] <- posterior_prob
R_i_all[[ineq_marker]] <- R_i
r_i_all[[ineq_marker]] <- matrix(r_i)
BFi_draws_post[ineq_marker] <- list(ci_draws_post)
BFi_draws_prior[ineq_marker] <- list(ci_draws_prior)
} else{
#**3.3)both-comparisons----
#****Equality
q_e <- nrow(R_e)
#Scale matrix for posterior t-distribution
scale_post <- vcov(object)
#Scale matrix for prior t-distribution
scale_prior <- vcov(object) * (n - k) / (n*b - k)
#a)Transformation matrix
D <- diag(k) - t(R_e) %*% solve(R_e %*% t(R_e)) %*% R_e
D2 <- unique(round(D, 5))
D2 <- D2[as.logical(rowSums(D2 != 0)),]
Tm <- rbind(R_e, D2)
#b)
w_post <- Tm %*% betahat
w_prior <- Tm %*% beta_zero
K_post <- Tm %*% scale_post %*% t(Tm)
K_prior <- Tm %*% scale_prior %*% t(Tm)
#Equality BF
log_BF <- mvtnorm::dmvt(x = t(r_e), delta = w_post[1:q_e], sigma = matrix(K_post[1:q_e, 1:q_e], ncol = q_e), df = n - k, log = TRUE) -
mvtnorm::dmvt(x = t(r_e), delta = w_prior[1:q_e], sigma = matrix(K_prior[1:q_e, 1:q_e], ncol = q_e), df = n*b - k, log = TRUE)
f_E <- mvtnorm::dmvt(x = t(r_e), delta = w_post[1:q_e], sigma = matrix(K_post[1:q_e, 1:q_e], ncol = q_e), df = n - k, log = FALSE)
c_E <- mvtnorm::dmvt(x = t(r_e), delta = w_prior[1:q_e], sigma = matrix(K_prior[1:q_e, 1:q_e], ncol = q_e), df = n*b - k, log = FALSE)
BFe <- exp(log_BF)
#****Inequality
R_iv <- R_i %*% MASS::ginv(D2)
r_iv <- r_i - R_i %*% MASS::ginv(R_e) %*% r_e
##Partitioning
#Posterior part
w_1_post <- w_post[1:q_e]
w_2_post <- w_post[(q_e + 1):k]
K_11_post <- K_post[1:q_e, 1:q_e]
K_12_post <- K_post[1:q_e, (q_e + 1):k]
K_21_post <- K_post[(q_e + 1):k, 1:q_e]
K_22_post <- K_post[(q_e + 1):k, (q_e + 1):k]
#prior part
w_1_prior <- w_prior[1:q_e]
w_2_prior <- w_prior[(q_e + 1):k]
K_11_prior <- K_prior[1:q_e, 1:q_e]
K_12_prior <- K_prior[1:q_e, (q_e + 1):k]
K_21_prior <- K_prior[(q_e + 1):k, 1:q_e]
K_22_prior <- K_prior[(q_e + 1):k, (q_e + 1):k]
#Conditional mean vectors and scale matrices
#posterior
w_2g1_post <- w_2_post + K_21_post %*% solve(K_11_post) %*% matrix(r_e - w_1_post)
K_2g1_post <- as.vector((n - k + (t(matrix(r_e - w_1_post)) %*% solve(K_11_post) %*% matrix(r_e - w_1_post))) /
(n - k + q_e)) * (K_22_post - K_21_post %*% solve(K_11_post) %*% t(K_21_post))
#prior
w_2g1_prior <- w_2_prior + K_21_prior %*% solve(K_11_prior) %*% matrix(r_e - w_1_prior)
K_2g1_prior <- as.vector((n*b - k + (t(matrix(r_e - w_1_prior)) %*% solve(K_11_prior) %*% matrix(r_e - w_1_prior))) /
(n*b - k + q_e)) * (K_22_prior - K_21_prior %*% solve(K_11_prior) %*% t(K_21_prior))
if(Matrix::rankMatrix(R_iv)[[1]] == nrow(R_iv)){
r_iv <- as.vector(r_iv)
delta_post <- as.vector(R_iv %*% w_2g1_post)
delta_prior <- as.vector(R_iv %*% w_2g1_prior)
scale_post_trans <- R_iv %*% K_2g1_post %*% t(R_iv)
scale_prior_trans <- R_iv %*% K_2g1_prior %*% t(R_iv)
if(nrow(scale_post_trans) == 1){ #If univariate
posterior_prob <- pt((r_iv - delta_post) / sqrt(scale_post_trans), df = n - k + q_e, lower.tail = FALSE)[1]
prior_prob <- pt((r_iv - delta_prior) / sqrt(scale_prior_trans), df = n*b - k + q_e, lower.tail = FALSE)[1]
BFi <- posterior_prob / prior_prob
} else { #if multivariate
posterior_prob <- mvtnorm::pmvt(lower = r_iv, upper = Inf, delta = delta_post, sigma = scale_post_trans, df = n - k + q_e, type = "shifted")[1]
prior_prob <- mvtnorm::pmvt(lower = r_iv, upper = Inf, delta = delta_prior, sigma = scale_prior_trans, df = n*b - k + q_e, type = "shifted")[1]
BFi <- posterior_prob / prior_prob
}
ci_lb <- ci_ub <- NA
} else{
if(!is.numeric(mcrep) || !mcrep %% 1 == 0) stop("Input for mcrep should be an integer")
#Draw from prior and posterior
draws_post <- mvtnorm::rmvt(n = mcrep, delta = w_2g1_post, sigma = K_2g1_post, df = n - k + q_e)
draws_prior <- mvtnorm::rmvt(n = mcrep, delta = w_2g1_prior, sigma = K_2g1_prior, df = n*b - k + q_e)
posterior_prob <- mean(apply(draws_post%*%t(R_iv) > rep(1,mcrep)%*%t(r_iv),1,prod))
prior_prob <- mean(apply(draws_prior%*%t(R_iv) > rep(1,mcrep)%*%t(r_iv),1,prod))
BFi <- posterior_prob / prior_prob
#Credibility interval
x_post <- sum(apply(draws_post%*%t(R_iv) > rep(1,mcrep)%*%t(r_iv),1,prod))
x_prior <- sum(apply(draws_prior%*%t(R_iv) > rep(1,mcrep)%*%t(r_iv),1,prod))
ci_draws_post <- rbeta(1e4, x_post, 1 + mcrep - x_post)
ci_draws_prior <- rbeta(1e4, x_prior, 1 + mcrep - x_prior)
BF_draws <- ci_draws_post / ci_draws_prior
ci_lb <- quantile(sort(BF_draws), 0.05)
ci_ub <- quantile(sort(BF_draws), 0.95)
}
c_i_e <- prior_prob
f_i_e <- posterior_prob
#Total BF
BF <- BFe * BFi
ci_lb <- ci_lb*BFe
ci_ub <- ci_ub*BFe
} #end 'both comparisons' option
out_c_E[h] <- c_E
out_f_E[h] <- f_E
out_c_i_e[h] <- c_i_e
out_f_i_e[h] <- f_i_e
out_ci_lb[h] <- ci_lb
out_ci_ub[h] <- ci_ub
BFu[h] <- BF #hypothesis vs. unconstrained
names(BFu)[[h]] <- paste0("H", h)
} #end loop over all hypotheses
if(!is.null(BFip_posterior)){
if(length(BFip_posterior) == 1){
BFc <- (1 - BFip_posterior) / (1 - BFip_prior)
comp_fie <- (1 - BFip_posterior)
comp_cie <- (1 - BFip_prior)
comp_fe <- comp_ce <- NA
#Credibility interval
if(length(Filter(Negate(is.null), BFi_draws_post)) == 0){
BFc_ci_lb <- BFc_ci_ub <- NA
}else{
BFc_draws_post <- 1 - unlist(BFi_draws_post)
BFc_draws_prior <- 1 - unlist(BFi_draws_prior)
BFc_draws <- BFc_draws_post / BFc_draws_prior
BFc_ci_lb <- quantile(sort(BFc_draws), 0.05)
BFc_ci_ub <- quantile(sort(BFc_draws), 0.95)
}
} else{
R_i_overlap <- do.call(rbind, R_i_all)
r_i_overlap <- do.call(rbind, r_i_all)
ineq_draws_prior <- mvtnorm::rmvt(n = 1e4, delta = beta_zero, sigma = vcov(object) * (n - k) / (n*b - k), df = (n*b - k))
exhaustive <- mean(rowSums(ineq_draws_prior%*%t(R_i_overlap) > rep(1, 1e4)%*%t(r_i_overlap)) > 0)
if(exhaustive == 1){
BFc <- comp_cie <- comp_fie <- comp_fe <- comp_ce <- BFc_ci_lb <- BFc_ci_ub <- NULL
} else{
overlap <- mean(apply(ineq_draws_prior%*%t(R_i_overlap) > rep(1, 1e4)%*%t(r_i_overlap), 1, prod))
if(overlap == 0){
BFc <- (1 - sum(BFip_posterior)) / (1 - sum(BFip_prior))
comp_fie <- (1 - sum(BFip_posterior))
comp_cie <- (1 - sum(BFip_prior))
#Credibility interval
if(length(Filter(Negate(is.null), BFi_draws_post)) == 0){
BFc_ci_lb <- BFc_ci_ub <- NA
}else{
analytic_BFip_posterior <- BFip_posterior[which(unlist(lapply(BFi_draws_post, is.null)) == 1)]
num_draws_BFip_posterior <- Filter(Negate(is.null),BFi_draws_post)
BFc_draws_post <- 1 - sum(analytic_BFip_posterior) - Reduce(`-`, num_draws_BFip_posterior)
analytic_BFip_prior <- BFip_prior[which(unlist(lapply(BFi_draws_prior, is.null)) == 1)]
num_draws_BFip_prior <- Filter(Negate(is.null),BFi_draws_prior)
BFc_draws_prior <- 1 - sum(analytic_BFip_prior) - Reduce(`-`, num_draws_BFip_prior)
BFc_draws <- BFc_draws_post / BFc_draws_prior
BFc_ci_lb <- quantile(sort(BFc_draws), 0.05)
BFc_ci_ub <- quantile(sort(BFc_draws), 0.95)
}
} else{
ineq_draws_posterior <- mvtnorm::rmvt(n = 1e4, delta = betahat, sigma = vcov(object), df = n - k)
constraints_prior <- Map(function(Ri, ri){apply(ineq_draws_prior%*%t(Ri) > rep(1,1e4)%*%t(ri), 1, prod)}, R_i_all, r_i_all)
constraints_posterior <- Map(function(Ri, ri){apply(ineq_draws_posterior%*%t(Ri) > rep(1,1e4)%*%t(ri), 1, prod)}, R_i_all, r_i_all)
prob_prior <- mean(Reduce(`+`, constraints_prior) > 0)
prob_posterior <- mean(Reduce(`+`, constraints_posterior) > 0)
BFc <- (1 - prob_posterior) / (1 - prob_prior)
comp_fie <- (1 - prob_posterior)
comp_cie <- (1 - prob_prior)
#Credibility interval
x_post <- 1e4 - sum(Reduce(`+`, constraints_posterior) > 0)
x_prior <- 1e4 - sum(Reduce(`+`, constraints_prior) > 0)
ci_draws_post <- rbeta(1e4, x_post, 1 + 1e4 - x_post)
ci_draws_prior <- rbeta(1e4, x_prior, 1 + 1e4 - x_prior)
BFc_draws <- ci_draws_post / ci_draws_prior
BFc_ci_lb <- quantile(sort(BFc_draws), 0.05)
BFc_ci_ub <- quantile(sort(BFc_draws), 0.95)
}
comp_fe <- comp_ce <- NA
}
}
} else{
BFc <- 1
comp_ce <- 1
comp_fe <- 1
comp_fie <- comp_cie <- BFc_ci_lb <- BFc_ci_ub <- NA
}
if(!is.null(BFc)){names(BFc) <- "Hc"}
BFu <- c(BFu, BFc)
if(!isTRUE(priorprob == 1) && !is.null(BFc) && length(BFu) > length(priorprob)) stop("Hypotheses have complement, add priorprob for it or use equal priorprobs.")
if(length(priorprob) > length(BFu)) stop("Too many priorprobs specified")
if(!sum(priorprob) == 1){
priorprob <- priorprob / sum(priorprob)
warning(paste(c("priorprobs did not sum to 1 and have been normalized. Used priorprobs: ",
round(priorprob, 4)), collapse = " "))
}
out_hyp_prob <- BFu*priorprob / sum(BFu*priorprob)
BF_matrix <- matrix(rep(BFu, length(BFu)), ncol = length(BFu), byrow = TRUE)
BF_matrix <- t(BF_matrix / BFu)
colnames(BF_matrix) <- rownames(BF_matrix) <- names(BFu)
BF_matrix[is.nan(BF_matrix)] <- 0
diag(BF_matrix) <- 1
if(length(hyp) > 1){
BF_matrices[[loop]] <- round(BF_matrix, digits = 3)
list_post_prob[[loop]] <- out_hyp_prob
}
} #end exploratory loop
if(length(hyp) > 1){
matrix_post_prob <- matrix(do.call(rbind, list_post_prob), ncol = 3)
colnames(matrix_post_prob) <- c("H1", "H2", "H3")
varnames <- gsub("Intercept", "(Intercept)", varnames)
names(BF_matrices) <- rownames(matrix_post_prob) <- varnames
out <- list(BF_matrix = BF_matrices, post_prob = matrix_post_prob,
hypotheses = hyp_out, BF_computation = "Not available when option 'exploratory' chosen.",
BFu_CI = "Not available when option 'exploratory' chosen.")
} else{
out_c_i_e <- c(out_c_i_e, comp_cie)
out_f_i_e <- c(out_f_i_e, comp_fie)
out_c_E <- c(out_c_E, comp_ce)
out_f_E <- c(out_f_E, comp_fe)
out_ci_lb <- c(out_ci_lb, BFc_ci_lb)
out_ci_ub <- c(out_ci_ub, BFc_ci_ub)
BF_computation <- as.matrix(data.frame(out_c_E, out_c_i_e, c = out_c_E*out_c_i_e, out_f_E, out_f_i_e,
f = out_f_E*out_f_i_e, BFu, PP_t = out_hyp_prob))
colnames(BF_computation) <- c("c(E)", "C(I|E)", "c", "f(E)", "f(I|E)", "f", "B(t,u)", "PP(t)")
BFu_ci <- as.matrix(data.frame(BFu, out_ci_lb, out_ci_ub))
colnames(BFu_ci) <- c("B(t,u)", "lb. (5%)", "ub. (95%)")
out <- list(BF_matrix = round(BF_matrix, digits = 3),post_prob = out_hyp_prob,
hypotheses = hyp_out, BF_computation = BF_computation, BFu_CI = BFu_ci)
}
class(out) <- "hyp"
out
}
#End----
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.