Nothing
R/BF.mvt_test.R
In BFpack: Flexible Bayes Factor Testing of Scientific Expectations
Defines functions
BF.mvt_test
get_estimates.mvt_test
mvt_test
Documented in
mvt_test
#' @title Multivariate Student t test
#' @description First step to performs a Bayesian multivariate one sample Student t test using the
#' (adjusted) fractional Bayes factor using the \code{BF()} function.
#'
#'@details \code{X} must be a data matrix and \code{null}
#'must be a vector of the assumed null values of the variables.
#'
#'@param X a data matrix with the variables in the columns.
#'@param Y an optional data matrix with the variables in the columns.
#'@param null a vector of the null values of the variables.
#'@param paired a logical indicating whether you want a multivariate paired t-test.
#'@param ... further arguments to be passed to or from methods.
#'
#'@return An object that can be applied to the \code{BF()}.
#'
#'@references Mulder, J. and Gu, X. (2023). Bayesian Testing of Scientific
#'Expectations under Multivariate Normal Linear Models. Multivariate Behavioral
#'Research, 57, 767-783. DOI: 10.1080/00273171.2021.1904809.
#'
#'@examples
#'
#'mvt_fmri <- mvt_test(fmri[,1:2],null = c(0,0))
#'BF(mvt_fmri)
#'
#'# the same test can be executed via the lm() function
#'intercept <- rep(1,nrow(fmri))
#'lm1 <- lm(cbind(Face,Vehicle) ~ -1 + intercept, data=fmri)
#'BF(lm1,hypothesis="intercept_on_Face=intercept_on_Vehicle=0")
#'
#' @rdname mvt_test
#' @export
mvt_test <- function(X, Y, null = NULL, paired = FALSE, ...){
if(missing(Y)){
X <- as.data.frame(as.matrix(X))
p <- ncol(X)
n <- nrow(X)
if(is.null(null)){
null <- rep(0,p)
}
if(length(null) != p){
stop("'null' must be a vector of length equal to the number of variables.")
}
if(p == 1){
out <- t_test(x=X,mu=null)
}else{
intercept <- rep(1,n)
varnames <- colnames(X)
formu <- as.formula(paste0("cbind(",paste0(varnames,collapse = ","),") ~ -1 + intercept"))
out <- lm(formu,data=X)
out$null <- null
out$numpop <- 1
out$paired <- paired
class(out) <- "mvt_test"
#
}
}else{
if(paired==TRUE){
p <- ncol(X)
n <- nrow(X)
if(is.null(null)){
null <- rep(0,p)
}
if(length(null) != p){
stop("'null' must be a vector of length equal to the number of variables.")
}
if(dim(X)[1]!=dim(Y)[1] || dim(X)[2]!=dim(Y)[2]){
stop("X and Y must have same dimension in case of paired test")
}
if(p == 1){
out <- t_test(x=X-Y,mu=null,paired=paired)
}else{
intercept <- rep(1,n)
varnames <- colnames(X)
formu <- as.formula(paste0("cbind(",paste0(varnames,collapse = ","),") ~ -1 + intercept"))
df.XY <- as.data.frame(X-Y)
colnames(df.XY) <- varnames
out <- lm(formu,data=df.XY)
out$null <- null
out$numpop <- 1
out$paired <- paired
class(out) <- "mvt_test"
#
}
}else{ # independent samples multivariate t test
# equal covariances assumed
p <- ncol(X)
n1 <- nrow(X)
n2 <- nrow(Y)
N <- n1 + n2
if(is.null(null)){
null <- rep(0,p)
}
if(length(null) != p){
stop("'null' must be a vector of length equal to the number of variables.")
}
if(dim(X)[2]!=dim(Y)[2]){
stop("X and Y must have same number of columns")
}
if(p == 1){
out <- t_test(x=X,y=Y,mu=null,paired=paired)
}else{
intercept <- rep(1,N)
difference <- c(rep(1,n1),rep(0,n2))
varnames <- colnames(X)
formu <- as.formula(paste0("cbind(",paste0(varnames,collapse = ","),") ~ -1 + intercept + difference"))
df.XY <- as.data.frame(rbind(X,Y))
colnames(df.XY) <- varnames
out <- lm(formu,data=df.XY)
out$null <- null
out$numpop <- 2
out$paired <- paired
class(out) <- "mvt_test"
#
}
}
}
return(out)
}
#' @method get_estimates mvt_test
#' @export
get_estimates.mvt_test <- function(x, ...){
class(x) <- "lm"
x_est <- get_estimates(x)
if(x$numpop==1){
if(x$paired==FALSE){
names(x_est$estimate) <- colnames(x_est$Sigma[[1]]) <-
row.names(x_est$Sigma[[1]]) <- gsub("intercept_on_","",names(x_est$estimate))
}else{
names(x_est$estimate) <- colnames(x_est$Sigma[[1]]) <-
row.names(x_est$Sigma[[1]]) <- gsub("intercept_on_","difference_",names(x_est$estimate))
}
}else{
P <- ncol(x$coefficients)
select_diff <- 2*(1:P)
names(x_est$estimate) <- colnames(x_est$Sigma[[1]]) <-
row.names(x_est$Sigma[[1]]) <- gsub("difference_on_","difference_",names(x_est$estimate))
x_est$estimate <- x_est$estimate[select_diff]
x_est$Sigma[[1]] <- x_est$Sigma[[1]][select_diff,select_diff]
}
x_est
}
#' @method BF mvt_test
#' @export
BF.mvt_test <- function(x,
hypothesis = NULL,
prior.hyp.explo = NULL,
prior.hyp.conf = NULL,
prior.hyp = NULL,
complement = TRUE,
log = FALSE,
cov.prob = .95,
BF.type = NULL,
...) {
if(!(cov.prob>0 & cov.prob<1)){
stop("The argument 'cov.prob' is a coverage probability for the interval estimates that should lie between 0 and 1. The default is 0.95.")
}
if(x$numpop==1 & x$paired==FALSE){
parameters <- "means"
}else{
parameters <- "differences"
}
if(x$numpop == 1){
P <- length(x$coefficients)
names1 <- colnames(x$coefficients)
names2 <- paste0("intercept_on_",names1)
}else{
P <- ncol(x$coefficients)
names1 <- colnames(x$coefficients)
names2 <- paste0("difference_on_",names1)
}
#exploratory test of joint equality to null
hypothesis.explo <- paste0(unlist(lapply(1:P,function(p){
paste0(names2[p],"=",x$null[p])
})),collapse = " & ")
x1 <- x
class(x1) <- "lm"
BF.explo <- BF(x1,
hypothesis=hypothesis.explo,
prior.hyp.conf=prior.hyp.explo,
log=log,
cov.prob=cov.prob,
BF.type=BF.type)
BF.explo$BFtu_confirmatory <- t(as.matrix(BF.explo$BFtu_confirmatory))
BF.explo$PHP_confirmatory <- t(as.matrix(BF.explo$PHP_confirmatory))
row.names(BF.explo$BFtu_confirmatory) <- row.names(BF.explo$PHP_confirmatory) <- parameters
colnames(BF.explo$BFtu_confirmatory) <- c("BF0u","BFuu")
colnames(BF.explo$PHP_confirmatory) <- c("Pr(=null)","Pr(not null)")
if(!is.null(hypothesis)){
variable.names <- colnames(x$coefficients)
if(x$numpop==1 | x$paired==TRUE){
add.name <- "intercept_on_"
rem.name <- ""
if(x$numpop==1 & x$paired==TRUE){
rem.name <- "difference_"
}
}else{
add.name <- "difference_on_"
rem.name <- "difference_"
}
hypothesis.updated <- hypothesis
for(varname in variable.names){
hypothesis.updated <- gsub(paste0(rem.name,varname),paste0(add.name,varname),hypothesis.updated)
}
}else{
hypothesis.updated <- NULL
}
BF.conf <- BF(x1,
hypothesis=hypothesis.updated,
prior.hyp.conf=prior.hyp.conf,
log=log,
complement=complement,
BF.type=BF.type)
if(!is.null(hypothesis)){
BF.conf$hypotheses <- gsub(add.name,rem.name,BF.conf$hypotheses)
}
BF.explo_estimates <- BF.explo$estimates
if(x$numpop == 1){
row.names(BF.explo_estimates) <- names(get_estimates(x)$estimate)
}else{
BF.explo_estimates <- BF.explo_estimates[(1:P)*2,]
row.names(BF.explo_estimates) <- names(get_estimates(x)$estimate)
}
BFlm_out <- list(
BFtu_exploratory=BF.explo$BFtu_confirmatory,
PHP_exploratory=BF.explo$PHP_confirmatory,
BFtu_confirmatory=BF.conf$BFtu_confirmatory,
PHP_confirmatory=BF.conf$PHP_confirmatory,
BFmatrix_confirmatory=BF.conf$BFmatrix_confirmatory,
BFtable_confirmatory=BF.conf$BFtable_confirmatory,
prior.hyp.explo=BF.explo$prior.hyp.conf,
prior.hyp.conf=BF.conf$prior.hyp.conf,
hypotheses=BF.conf$hypotheses,
estimates=BF.explo_estimates,
model=x1,
bayesfactor=BF.conf$bayesfactor,
parameter=parameters,
log = BF.conf$log,
fraction_number_groupIDs = BF.conf$fraction_number_groupIDs,
fraction_groupID_observations = BF.conf$fraction_groupID_observations,
fraction_groupID = BF.conf$fraction_groupID,
call=match.call())
class(BFlm_out) <- "BF"
return(BFlm_out)
}
Try the BFpack package in your browser
Any scripts or data that you put into this service are public.
BFpack documentation built on April 12, 2025, 1:40 a.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 BF.mvt_test get_estimates.mvt_test mvt_test
Documented in mvt_test
#' @title Multivariate Student t test
#' @description First step to performs a Bayesian multivariate one sample Student t test using the
#' (adjusted) fractional Bayes factor using the \code{BF()} function.
#'
#'@details \code{X} must be a data matrix and \code{null}
#'must be a vector of the assumed null values of the variables.
#'
#'@param X a data matrix with the variables in the columns.
#'@param Y an optional data matrix with the variables in the columns.
#'@param null a vector of the null values of the variables.
#'@param paired a logical indicating whether you want a multivariate paired t-test.
#'@param ... further arguments to be passed to or from methods.
#'
#'@return An object that can be applied to the \code{BF()}.
#'
#'@references Mulder, J. and Gu, X. (2023). Bayesian Testing of Scientific
#'Expectations under Multivariate Normal Linear Models. Multivariate Behavioral
#'Research, 57, 767-783. DOI: 10.1080/00273171.2021.1904809.
#'
#'@examples
#'
#'mvt_fmri <- mvt_test(fmri[,1:2],null = c(0,0))
#'BF(mvt_fmri)
#'
#'# the same test can be executed via the lm() function
#'intercept <- rep(1,nrow(fmri))
#'lm1 <- lm(cbind(Face,Vehicle) ~ -1 + intercept, data=fmri)
#'BF(lm1,hypothesis="intercept_on_Face=intercept_on_Vehicle=0")
#'
#' @rdname mvt_test
#' @export
mvt_test <- function(X, Y, null = NULL, paired = FALSE, ...){
if(missing(Y)){
X <- as.data.frame(as.matrix(X))
p <- ncol(X)
n <- nrow(X)
if(is.null(null)){
null <- rep(0,p)
}
if(length(null) != p){
stop("'null' must be a vector of length equal to the number of variables.")
}
if(p == 1){
out <- t_test(x=X,mu=null)
}else{
intercept <- rep(1,n)
varnames <- colnames(X)
formu <- as.formula(paste0("cbind(",paste0(varnames,collapse = ","),") ~ -1 + intercept"))
out <- lm(formu,data=X)
out$null <- null
out$numpop <- 1
out$paired <- paired
class(out) <- "mvt_test"
#
}
}else{
if(paired==TRUE){
p <- ncol(X)
n <- nrow(X)
if(is.null(null)){
null <- rep(0,p)
}
if(length(null) != p){
stop("'null' must be a vector of length equal to the number of variables.")
}
if(dim(X)[1]!=dim(Y)[1] || dim(X)[2]!=dim(Y)[2]){
stop("X and Y must have same dimension in case of paired test")
}
if(p == 1){
out <- t_test(x=X-Y,mu=null,paired=paired)
}else{
intercept <- rep(1,n)
varnames <- colnames(X)
formu <- as.formula(paste0("cbind(",paste0(varnames,collapse = ","),") ~ -1 + intercept"))
df.XY <- as.data.frame(X-Y)
colnames(df.XY) <- varnames
out <- lm(formu,data=df.XY)
out$null <- null
out$numpop <- 1
out$paired <- paired
class(out) <- "mvt_test"
#
}
}else{ # independent samples multivariate t test
# equal covariances assumed
p <- ncol(X)
n1 <- nrow(X)
n2 <- nrow(Y)
N <- n1 + n2
if(is.null(null)){
null <- rep(0,p)
}
if(length(null) != p){
stop("'null' must be a vector of length equal to the number of variables.")
}
if(dim(X)[2]!=dim(Y)[2]){
stop("X and Y must have same number of columns")
}
if(p == 1){
out <- t_test(x=X,y=Y,mu=null,paired=paired)
}else{
intercept <- rep(1,N)
difference <- c(rep(1,n1),rep(0,n2))
varnames <- colnames(X)
formu <- as.formula(paste0("cbind(",paste0(varnames,collapse = ","),") ~ -1 + intercept + difference"))
df.XY <- as.data.frame(rbind(X,Y))
colnames(df.XY) <- varnames
out <- lm(formu,data=df.XY)
out$null <- null
out$numpop <- 2
out$paired <- paired
class(out) <- "mvt_test"
#
}
}
}
return(out)
}
#' @method get_estimates mvt_test
#' @export
get_estimates.mvt_test <- function(x, ...){
class(x) <- "lm"
x_est <- get_estimates(x)
if(x$numpop==1){
if(x$paired==FALSE){
names(x_est$estimate) <- colnames(x_est$Sigma[[1]]) <-
row.names(x_est$Sigma[[1]]) <- gsub("intercept_on_","",names(x_est$estimate))
}else{
names(x_est$estimate) <- colnames(x_est$Sigma[[1]]) <-
row.names(x_est$Sigma[[1]]) <- gsub("intercept_on_","difference_",names(x_est$estimate))
}
}else{
P <- ncol(x$coefficients)
select_diff <- 2*(1:P)
names(x_est$estimate) <- colnames(x_est$Sigma[[1]]) <-
row.names(x_est$Sigma[[1]]) <- gsub("difference_on_","difference_",names(x_est$estimate))
x_est$estimate <- x_est$estimate[select_diff]
x_est$Sigma[[1]] <- x_est$Sigma[[1]][select_diff,select_diff]
}
x_est
}
#' @method BF mvt_test
#' @export
BF.mvt_test <- function(x,
hypothesis = NULL,
prior.hyp.explo = NULL,
prior.hyp.conf = NULL,
prior.hyp = NULL,
complement = TRUE,
log = FALSE,
cov.prob = .95,
BF.type = NULL,
...) {
if(!(cov.prob>0 & cov.prob<1)){
stop("The argument 'cov.prob' is a coverage probability for the interval estimates that should lie between 0 and 1. The default is 0.95.")
}
if(x$numpop==1 & x$paired==FALSE){
parameters <- "means"
}else{
parameters <- "differences"
}
if(x$numpop == 1){
P <- length(x$coefficients)
names1 <- colnames(x$coefficients)
names2 <- paste0("intercept_on_",names1)
}else{
P <- ncol(x$coefficients)
names1 <- colnames(x$coefficients)
names2 <- paste0("difference_on_",names1)
}
#exploratory test of joint equality to null
hypothesis.explo <- paste0(unlist(lapply(1:P,function(p){
paste0(names2[p],"=",x$null[p])
})),collapse = " & ")
x1 <- x
class(x1) <- "lm"
BF.explo <- BF(x1,
hypothesis=hypothesis.explo,
prior.hyp.conf=prior.hyp.explo,
log=log,
cov.prob=cov.prob,
BF.type=BF.type)
BF.explo$BFtu_confirmatory <- t(as.matrix(BF.explo$BFtu_confirmatory))
BF.explo$PHP_confirmatory <- t(as.matrix(BF.explo$PHP_confirmatory))
row.names(BF.explo$BFtu_confirmatory) <- row.names(BF.explo$PHP_confirmatory) <- parameters
colnames(BF.explo$BFtu_confirmatory) <- c("BF0u","BFuu")
colnames(BF.explo$PHP_confirmatory) <- c("Pr(=null)","Pr(not null)")
if(!is.null(hypothesis)){
variable.names <- colnames(x$coefficients)
if(x$numpop==1 | x$paired==TRUE){
add.name <- "intercept_on_"
rem.name <- ""
if(x$numpop==1 & x$paired==TRUE){
rem.name <- "difference_"
}
}else{
add.name <- "difference_on_"
rem.name <- "difference_"
}
hypothesis.updated <- hypothesis
for(varname in variable.names){
hypothesis.updated <- gsub(paste0(rem.name,varname),paste0(add.name,varname),hypothesis.updated)
}
}else{
hypothesis.updated <- NULL
}
BF.conf <- BF(x1,
hypothesis=hypothesis.updated,
prior.hyp.conf=prior.hyp.conf,
log=log,
complement=complement,
BF.type=BF.type)
if(!is.null(hypothesis)){
BF.conf$hypotheses <- gsub(add.name,rem.name,BF.conf$hypotheses)
}
BF.explo_estimates <- BF.explo$estimates
if(x$numpop == 1){
row.names(BF.explo_estimates) <- names(get_estimates(x)$estimate)
}else{
BF.explo_estimates <- BF.explo_estimates[(1:P)*2,]
row.names(BF.explo_estimates) <- names(get_estimates(x)$estimate)
}
BFlm_out <- list(
BFtu_exploratory=BF.explo$BFtu_confirmatory,
PHP_exploratory=BF.explo$PHP_confirmatory,
BFtu_confirmatory=BF.conf$BFtu_confirmatory,
PHP_confirmatory=BF.conf$PHP_confirmatory,
BFmatrix_confirmatory=BF.conf$BFmatrix_confirmatory,
BFtable_confirmatory=BF.conf$BFtable_confirmatory,
prior.hyp.explo=BF.explo$prior.hyp.conf,
prior.hyp.conf=BF.conf$prior.hyp.conf,
hypotheses=BF.conf$hypotheses,
estimates=BF.explo_estimates,
model=x1,
bayesfactor=BF.conf$bayesfactor,
parameter=parameters,
log = BF.conf$log,
fraction_number_groupIDs = BF.conf$fraction_number_groupIDs,
fraction_groupID_observations = BF.conf$fraction_groupID_observations,
fraction_groupID = BF.conf$fraction_groupID,
call=match.call())
class(BFlm_out) <- "BF"
return(BFlm_out)
}
Try the BFpack package in your browser
Any scripts or data that you put into this service are public.
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.