============================================================================================================ PREVIOUS: https://rpubs.com/sherloconan/1033617
THIS: https://rpubs.com/sherloconan/1113820
MORE: https://github.com/zhengxiaoUVic/Bayesian
For testing the predictor of interest \(a\) in two-way repeated-measures designs,
the inclusion BFs of models containing the effect relative to models stripping of it
(always including random slopes [RS] and enforcing the principle of marginality for fixed effects) are
\[\begin{align*} \textit{BF}_{\text{incl}}\ \to\quad&\ [s+a+b+ab+sa+sb] \\ +&\ [s+a+b+sa+sb] \\ +&\ [s+a+sa+sb] \\ &\qquad\ / \tag{1} \\ &\ [s+b+sa+sb] \\ +&\ [s+sa+sb] \end{align*}\]
\[\begin{align*} \textit{BF}_{\text{incl}}^{\hspace{0.2em}\text{match}}\ \to\quad& \\ &\ [s+a+b+sa+sb] \\ +&\ [s+a+sa+sb] \\ &\qquad\ / \tag{2} \\ &\ [s+b+sa+sb] \\ +&\ [s+sa+sb] \end{align*}\]
Prior Model Probabilities, s.t. \(\sum_\xi p(\mathcal{M}_\xi)=1\).
Posterior Model Probabilities, s.t. \(\sum_\xi p(\mathcal{M}_\xi\mid\boldsymbol{y})=1\).
Prior Inclusion Probability, \(p(\mathcal{M}_\text{incl})=\sum_\eta p(\mathcal{M}_\eta)\)
Posterior Inclusion Probability, \[\begin{equation*} p(\mathcal{M}_\text{incl}\mid\boldsymbol{y})=\sum_\eta p(\mathcal{M}_\eta\mid\boldsymbol{y})=\frac{\sum_\eta p(\boldsymbol{y}\mid\mathcal{M}_\eta)\cdot p(\mathcal{M}_\eta)}{\sum_\xi p(\boldsymbol{y}\mid\mathcal{M}_\xi)\cdot p(\mathcal{M}_\xi)}=\frac{\sum_\eta\textit{BF}_{\eta\centerdot}^{\hspace{0.25em}\text{RS}}\cdot p(\mathcal{M}_\eta)}{\sum_\xi\textit{BF}_{\xi\centerdot}^{\hspace{0.25em}\text{RS}}\cdot p(\mathcal{M}_\xi)} \end{equation*}\]
\[\begin{equation*} \tag{3} \textit{BF}_{\hspace{0.1em}\text{incl}}=\frac{p(\mathcal{M}_\text{incl}\mid\boldsymbol{y})\ /\ \big(1-p(\mathcal{M}_\text{incl}\mid\boldsymbol{y})\big)}{p(\mathcal{M}_\text{incl})\ /\ \big(1-p(\mathcal{M}_\text{incl})\big)} \end{equation*}\]
\[\begin{equation*} \tag{4} \textit{BF}_{\hspace{0.1em}\text{incl}}^{\hspace{0.25em}\text{match}}=\frac{p(\mathcal{M}_{\text{incl}}^{\text{match}}\mid\boldsymbol{y})\ /\ p(\mathcal{M}_{\text{excl}}^{\text{match}}\mid\boldsymbol{y})}{p(\mathcal{M}_{\text{incl}}^{\text{match}})\ /\ p(\mathcal{M}_{\text{excl}}^{\text{match}})} \end{equation*}\]
inclusionBF <- function(data, legacy=F, match_models=F, prior_prob=rep(.2, 5), ...) {
#' Input -
#' data: two-way repeated-measures data whose column names must contain "DV", "ID", "A", and "B"
#' legacy: if FALSE (default), compare MRE models; if TRUE, compare RIO models
#' match_models: whether effects are across matched models (default is FALSE)
#' prior_prob: prior probabilities of [1] full, [2] additive, [3] A, [4] B,
#' and [5] null models (default is being equal)
#' Output -
#' a list of the JASP-replicating (v 0.16.3 or later) Bayes factors
#' from the "BayesFactor" dependency
if(round(sum(prior_prob), 3) != 1 || length(prior_prob) != 5) stop("Prior model probabilities must sum to 1!")
if(!legacy) {
BF_full <- BayesFactor::lmBF(DV ~ A * B + ID + A:ID + B:ID, data, whichRandom="ID", ...)
BF_additive <- BayesFactor::lmBF(DV ~ A + B + ID + A:ID + B:ID, data, whichRandom="ID", ...)
BF_A <- BayesFactor::lmBF(DV ~ A + ID + A:ID + B:ID, data, whichRandom="ID", ...)
BF_B <- BayesFactor::lmBF(DV ~ B + ID + A:ID + B:ID, data, whichRandom="ID", ...)
BF_0 <- BayesFactor::lmBF(DV ~ ID + A:ID + B:ID, data, whichRandom="ID", ...)
} else {
BF_full <- BayesFactor::lmBF(DV ~ A * B + ID, data, whichRandom="ID", ...)
BF_additive <- BayesFactor::lmBF(DV ~ A + B + ID, data, whichRandom="ID", ...)
BF_A <- BayesFactor::lmBF(DV ~ A + ID, data, whichRandom="ID", ...)
BF_B <- BayesFactor::lmBF(DV ~ B + ID, data, whichRandom="ID", ...)
BF_0 <- BayesFactor::lmBF(DV ~ ID, data, whichRandom="ID", ...)
}
BF <- data.frame(c(BF_full, BF_additive, BF_A, BF_B, BF_0) / BF_0) #compare to null
BF_weighted <- BF$bf * prior_prob
BF_avg <- sum(BF_weighted)
if(match_models) {
indices_incl <- list(c(2, 3), c(2, 4), 1) #indices of inclusion matched models
indices_excl <- list(c(4, 5), c(3, 5), 2) #indices of exclusion matched models
} else {
indices_incl <- list(c(1, 2, 3), c(1, 2, 4), 1) #indices of inclusion models
indices_excl <- list(c(4, 5), c(3, 5), 2:5) #indices of exclusion models
}
out <- data.frame("Effects" = c("A", "B", "AB"),
"Prior.Incl" = sapply(indices_incl, function(x) sum(prior_prob[x])), #prior inclusion prob
"Prior.Excl" = sapply(indices_excl, function(x) sum(prior_prob[x])), #prior exclusion prob
"Post.Incl" = sapply(indices_incl, function(x) sum(BF_weighted[x])) / BF_avg,
"Post.Excl" = sapply(indices_excl, function(x) sum(BF_weighted[x])) / BF_avg)
out$BF.Incl <- (out$Post.Incl * out$Prior.Excl) / (out$Post.Excl * out$Prior.Incl) #inclusion Bayes factors
OG <- data.frame("Models" = row.names(BF), "Prior.M" = prior_prob, "Post.M" = BF_weighted / BF_avg,
"BF10" = BF$bf, "error.per" = BF$error * 100) #JASP model comparison
OG$BF.M <- (OG$Post.M * (1 - prior_prob)) / ((1 - OG$Post.M) * prior_prob) #model Bayes factors
OG <- OG[order(OG$BF10, decreasing=T), c(1:3,6,4,5)]
row.names(OG) <- NULL
#RIO: random intercept-only; MRE: maximal set of random effects; BMA: Bayesian model averaging
if(!legacy) {list("MRE" = OG, "BMA" = out)} else {list("RIO" = OG, "BMA" = out)}
}
The data set (long format) consists of four columns.
Let resp be factor A (with levels \(\alpha_i\)) and align be factor B (with levels \(\beta_j\)) for \(i,j\in\{1,2\}\).
The Paradox: Factor B is statistically significant in the repeated-measures ANOVA but not supported in the anovaBF R function (version 0.9.12-4.6 or lower).
One can break the data into two parts, separated by the factor A, and conduct separate Bayesian tests for the factor of interest B in each condition of A. Under these circumstances, a very clear effect B emerges in both analyses. This discrepancy is not the product of a contrast between a liberal and a conservative test (Bub, Masson, & van Noordenne, 2021, p. 80).
df <- read.table("~/Documents/UVic/Project Meeting/paradox/UpDown3.txt", header=T, stringsAsFactors=T)
dat <- matrix(df$RT,ncol=4,byrow=T,
dimnames=list(paste0("s",1:31),
c("onehand&A","onehand&M","twohands&A","twohands&M"))) #wide data format
cov(dat[,c(1,3,2,4)]) %>% #sample covariance matrix
kable(digits=3) %>% kable_styling(full_width=F)| onehand&A | twohands&A | onehand&M | twohands&M | |
|---|---|---|---|---|
| onehand&A | 3937.495 | -207.186 | 4066.825 | -830.866 |
| twohands&A | -207.186 | 5202.157 | -565.871 | 5922.451 |
| onehand&M | 4066.825 | -565.871 | 4779.546 | -979.485 |
| twohands&M | -830.866 | 5922.451 | -979.485 | 7274.140 |
summary(aov(RT~(resp*align)+Error(subj/(resp*align)),df)) #repeated-measures ANOVA##
## Error: subj
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 30 270038 9001
##
## Error: subj:resp
## Df Sum Sq Mean Sq F value Pr(>F)
## resp 1 21479 21479 1.854 0.183
## Residuals 30 347540 11585
##
## Error: subj:align
## Df Sum Sq Mean Sq F value Pr(>F)
## align 1 11536 11536 28.22 9.64e-06 ***
## Residuals 30 12262 409
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: subj:resp:align
## Df Sum Sq Mean Sq F value Pr(>F)
## resp:align 1 1034 1033.6 5.203 0.0298 *
## Residuals 30 5960 198.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\[\begin{align} \tag{5} \text{Total}\ \mathit{SS}&=\mathit{SS}_\text{subjects}+\mathit{SS}_{\text{within-subjects}} \\ &=\mathit{SS}_\text{subjects}+[(\mathit{SS}_\text{A}+\mathit{SS}_\text{Error A})+(\mathit{SS}_\text{B}+\mathit{SS}_\text{Error B})+(\mathit{SS}_\text{AB}+\mathit{SS}_\text{Error AB})] \end{align}\]
\[\begin{align} \tag{6} \text{Total}\ \mathit{df}&=\mathit{df}_\text{subjects}+[(\mathit{df}_\text{A}+\mathit{df}_\text{Error A})+(\mathit{df}_\text{B}+\mathit{df}_\text{Error B})+(\mathit{df}_\text{AB}+\mathit{df}_\text{Error AB})]\\ &=(n-1)+[((a-1)+(a-1)(n-1))+((b-1)+(b-1)(n-1))+((a-1)(b-1)+(a-1)(b-1)(n-1))] \end{align}\]
options(digits=3)
colnames(df) <- c("DV", "ID", "A", "B")
set.seed(277)
inclusionBF(df, legacy=T, iterations=1e5, progress=F)## $RIO
## Models Prior.M Post.M BF.M BF10 error.per
## 1 A + ID 0.2 0.3615 2.264 2.100 1.08
## 2 A + B + ID 0.2 0.2646 1.439 1.537 2.28
## 3 ID 0.2 0.1721 0.832 1.000 0.00
## 4 B + ID 0.2 0.1265 0.579 0.735 8.80
## 5 A * B + ID 0.2 0.0753 0.326 0.438 1.65
##
## $BMA
## Effects Prior.Incl Prior.Excl Post.Incl Post.Excl BF.Incl
## 1 A 0.6 0.4 0.7014 0.299 1.566
## 2 B 0.6 0.4 0.4664 0.534 0.583
## 3 AB 0.2 0.8 0.0753 0.925 0.326inclusionBF(df, legacy=T, match_models=T, iterations=1e5, progress=F)## $RIO
## Models Prior.M Post.M BF.M BF10 error.per
## 1 A + ID 0.2 0.3673 2.322 2.089 1.272
## 2 A + B + ID 0.2 0.2648 1.441 1.507 0.804
## 3 ID 0.2 0.1758 0.853 1.000 0.000
## 4 B + ID 0.2 0.1192 0.541 0.678 0.662
## 5 A * B + ID 0.2 0.0729 0.314 0.415 0.885
##
## $BMA
## Effects Prior.Incl Prior.Excl Post.Incl Post.Excl BF.Incl
## 1 A 0.4 0.4 0.6321 0.295 2.143
## 2 B 0.4 0.4 0.3841 0.543 0.707
## 3 AB 0.2 0.2 0.0729 0.265 0.275
JASP v 0.16.3 or later. See van den Bergh, Wagenmakers, and Aust (2023).
Inclusion Bayes factors for matched models (“Baws factors”) have been featured since v 0.8.2.0.
set.seed(277)
inclusionBF(df, iterations=1e5, progress=F)## $MRE
## Models Prior.M Post.M BF.M BF10 error.per
## 1 A * B + ID + A:ID + B:ID 0.2 0.387661 2.532334 1810.06 5.59
## 2 B + ID + A:ID + B:ID 0.2 0.353610 2.188217 1651.07 1.16
## 3 A + B + ID + A:ID + B:ID 0.2 0.258405 1.393779 1206.54 27.48
## 4 ID + A:ID + B:ID 0.2 0.000214 0.000857 1.00 0.00
## 5 A + ID + A:ID + B:ID 0.2 0.000109 0.000437 0.51 3.28
##
## $BMA
## Effects Prior.Incl Prior.Excl Post.Incl Post.Excl BF.Incl
## 1 A 0.6 0.4 0.646 0.353824 1.22
## 2 B 0.6 0.4 1.000 0.000323 2060.85
## 3 AB 0.2 0.8 0.388 0.612339 2.53inclusionBF(df, match_models=T, iterations=1e5, progress=F)## $MRE
## Models Prior.M Post.M BF.M BF10 error.per
## 1 A * B + ID + A:ID + B:ID 0.2 0.420339 2.900581 1851.466 9.03
## 2 B + ID + A:ID + B:ID 0.2 0.370280 2.352031 1630.974 1.45
## 3 A + B + ID + A:ID + B:ID 0.2 0.209031 1.057087 920.718 5.33
## 4 ID + A:ID + B:ID 0.2 0.000227 0.000908 1.000 0.00
## 5 A + ID + A:ID + B:ID 0.2 0.000123 0.000494 0.543 6.67
##
## $BMA
## Effects Prior.Incl Prior.Excl Post.Incl Post.Excl BF.Incl
## 1 A 0.4 0.4 0.209 0.37051 0.565
## 2 B 0.4 0.4 0.579 0.00035 1653.259
## 3 AB 0.2 0.2 0.420 0.20903 2.011