# R-code for example analyses of Applications 1 to 8 (Section 4) in BFpack: # Flexible Bayes Factor Testing of Scientific Theories in R by Joris Mulder, # Donald R. Williams, Xin Gu, Andrew Tomarken, Florian Boing-Messing, Anton # Olsson-Collentine, Marlyne Bosman-Meijerink, Janosch Menke, Robbie van Aert, # Jean-Paul Fox, Yes Rossen, Eric-Jan Wagenmakers, Herbert Hoijtink, and Caspar # van Lissa. Journal of Statistical Software. # Application 1: two-sample t test on the birthwt data library("MASS") library("bain") library("BFpack") smoke0 <- subset(MASS::birthwt, smoke == 0, select = bwt) smoke1 <- subset(MASS::birthwt, smoke == 1, select = bwt) ttest1 <- bain::t_test(smoke0, smoke1, alternative = "greater", var.equal = TRUE) print(ttest1) constraints1 <- "difference = 0; difference > 0" BF1 <- BF(ttest1, hypothesis = constraints1, complement = FALSE) summary(BF1) # Bayes factor of H2 against H1 in the exploratory test BF1$BFtu_exploratory[3]/BF1$BFtu_exploratory[1] # Application 2: anova model on the tvprices data aov2 <- aov(price ~ anchor * motivation, data = tvprices) summary(aov2) # Perform exploratory Bayes factor tests for ANOVA design constraints2 <- "anchorrounded > 0 & motivationlow < 0 & anchorrounded:motivationlow < 0; anchorrounded = 0 & motivationlow = 0 & anchorrounded:motivationlow = 0" set.seed(1234) BF2 <- BF(aov2, hypothesis = constraints2) summary(BF2) # Perform analysis via a regression analysis lm2 <- lm(price ~ anchor * motivation, data = tvprices) summary(lm2) # Perform exploratory Bayes factor tests for ANOVA design BF2 <- BF(lm2) print(BF2) # Application 3: Group variance testing on the accuracy data bartlett3 <- bartlett_test(x = attention$accuracy, g = attention$group) print(bartlett3) # Specify informative hypotheses on the group variances based on substantive # expectations get_estimates(bartlett3) constraints3 <- c("Controls = TS < ADHD; Controls < TS = ADHD; Controls = TS = ADHD") set.seed(358) # Perform the Bayesian hypothesis tests BF3 <- BF(bartlett3, hypothesis = constraints3) summary(BF3) # Application 4: logistic regression model on the wilson data glm4 <- glm(sent ~ ztrust + zfWHR + zAfro + glasses + attract + maturity + tattoos, family = binomial(), data = wilson) # summary of classical tests summary(glm4) summary(glm4)$coefficients # Perform the Bayesian hypothesis tests set.seed(123) constraints4 <- "ztrust > (zfWHR, zAfro) > 0; ztrust > zfWHR = zAfro = 0" BF4 <- BF(glm4, hypothesis = constraints4) summary(BF4) # exploratory Bayesian test using output from classical test ct <- lmtest::coeftest(glm4) BF(ct[, 1], Sigma = diag(ct[, 2]^2), n = attr(ct, "nobs")) # Application 5: Multivariate regression model on the fmri data mlm5a <- lm(cbind(Superficial, Middle, Deep) ~ Face + Vehicle, data = fmri) summary(mlm5a) # Specify informative hypotheses on the coefficients across predictor and # dependent variables based on substantive expectations constraints5a <- "Face_on_Deep = Face_on_Superficial = Face_on_Middle < 0 < Vehicle_on_Deep = Vehicle_on_Superficial = Vehicle_on_Middle; Face_on_Deep < Face_on_Superficial = Face_on_Middle < 0 < Vehicle_on_Deep = Vehicle_on_Superficial = Vehicle_on_Middle" # Perform the Bayesian hypothesis tests set.seed(123) BF5a <- BF(mlm5a, hypothesis = constraints5a) summary(BF5a) # Application 5b (with random missing observations) Specify constrained # hypotheses constraints5b <- "Face_on_Deep = Face_on_Superficial = Face_on_Middle < 0; Face_on_Deep < Face_on_Superficial = Face_on_Middle < 0" # Fit model on complete data mlm5b <- lm(cbind(Superficial, Middle, Deep) ~ Face + Vehicle, data = fmri) # Compute Bayes factors and posterior probabilities based on complete data BF5b <- BF(mlm5b, hypothesis = constraints5b) print(BF5b) # Generate random missings fmri_missing <- fmri set.seed(1234) for (i in 1:10) { fmri_missing[sample(1:nrow(fmri), 1), sample(1:ncol(fmri), 1)] <- NA } # List-wise delete rows that contain at least one missing observation fmri_listdel <- fmri_missing[!is.na(apply(fmri_missing, 1, sum)), ] mlm5b_listdel <- lm(cbind(Superficial, Middle, Deep) ~ Face + Vehicle, data = fmri_listdel) # Execute Bayes factor test based on list-wise deleted data BF5b_listdel <- BF(mlm5b_listdel, hypothesis = constraints5b) print(BF5b_listdel) # Create 500 complete imputed data sets using the mice package M <- 500 library("mice") mice_fmri <- mice::mice(data = fmri_missing, m = M, meth = c("norm", "norm", "norm", "norm", "norm"), diagnostics = F, printFlag = F) # Execute tests based on complete imputed data. Extract the posterior and # prior quantities of the extended Savage-Dickey density ratio from the # Specification matrix for the imputed data. relmeas_all <- matrix(unlist(lapply(1:M, function(m) { mlm5b_m <- lm(cbind(Superficial, Middle, Deep) ~ Face + Vehicle, data = mice::complete(mice_fmri, m)) BF5b_m <- BF(mlm5b_m, hypothesis = constraints5b) # Return four quantities of SD-ratio per constrained hypothesis c(BF5b_m$BFtable_confirmatory[, 1:4]) })), ncol = M) # Compute Monte Carlo estimate of the four quantities for the constrained # hypotheses relmeas <- matrix(apply(relmeas_all, 1, mean), nrow = 3) # Give appropriate names row.names(relmeas) <- c("H1", "H2", "H3") colnames(relmeas) <- c("complex=", "complex>", "fit=", "fit>") # Compute Bayes factors using Equation 4 BF_tu_confirmatory <- relmeas[, 3] * relmeas[, 4]/(relmeas[, 1] * relmeas[, 2]) # Compute Posterior probabilities using Equation 3. PHP <- BF_tu_confirmatory/sum(BF_tu_confirmatory) print(PHP) # Application 6: Correlation analysis on the memory data memoryHC <- subset(memory, Group == "HC")[, -7] memorySZ <- subset(memory, Group == "SZ")[, -7] Hmisc::rcorr(as.matrix(memoryHC)) Hmisc::rcorr(as.matrix(memorySZ)) set.seed(123) # Execute unconstrained Bayesian correlation analysis using cor_test cor6 <- cor_test(memoryHC, memorySZ) print(cor6) # Perform Bayes factor test of the order hypothesis against its complement constraints6 <- "Del_with_Im_in_g1 > Del_with_Im_in_g2 & Del_with_Wmn_in_g1 > Del_with_Wmn_in_g2 & Del_with_Cat_in_g1 > Del_with_Cat_in_g2 & Del_with_Fas_in_g1 > Del_with_Fas_in_g2 & Del_with_Rat_in_g1 > Del_with_Rat_in_g2 & Im_with_Wmn_in_g1 > Im_with_Wmn_in_g2 & Im_with_Cat_in_g1 > Im_with_Cat_in_g2 & Im_with_Fas_in_g1 > Im_with_Fas_in_g2 & Im_with_Rat_in_g1 > Im_with_Rat_in_g2 & Wmn_with_Cat_in_g1 > Wmn_with_Cat_in_g2 & Wmn_with_Fas_in_g1 > Wmn_with_Fas_in_g2 & Wmn_with_Rat_in_g1 > Wmn_with_Rat_in_g2 & Cat_with_Fas_in_g1 > Cat_with_Fas_in_g2 & Cat_with_Rat_in_g1 > Cat_with_Rat_in_g2 & Fas_with_Rat_in_g1 > Fas_with_Rat_in_g2" BF6 <- BF(cor6, hypothesis = constraints6) summary(BF6) # Application 7. Testing intraclass correlations on the timssICC data timssICC_subset <- subset(timssICC, groupNL11 == 1 | groupHR11 == 1 | groupDE11 == 1 | groupDK11 == 1) # Fit a random intercept model with country specific random effect variances # across schools lmer7 <- lme4::lmer(math ~ -1 + gender + weight + lln + groupNL11 + (0 + groupNL11 | schoolID) + groupHR11 + (0 + groupHR11 | schoolID) + groupDE11 + (0 + groupDE11 | schoolID) + groupDK11 + (0 + groupDK11 | schoolID), data = timssICC_subset) print(lmer7) # Perform Bayes factor test on intraclass correlations across countries in 2011 set.seed(123) constraints7 <- "groupNL11 0; culture > location > 0; location > culture > 0" # create object for the unconstrained estimates estimates8 <- remdyad8$coef # create object for the error covariance matrix of the estimates covmatrix8 <- remdyad8$cov # create object for sample size (number of events in this case) samplesize8 <- remdyad8$m # compute Bayes factors using BF.default BF8 <- BF(estimates8, Sigma = covmatrix8, n = samplesize8, hypothesis = constraints8) # print output summary(BF8)