Alternative Covariance Structure of the Error Term

1. Getting Started

The data set (long format) consists of four columns.

Let resp be factor A (with levels $\alpha_i$ or $a_i$) and align be factor B (with levels $\beta_j$ or $b_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).

simData_eqVar <- function(seed=277, N=50, eqVar=T) {
  #' Input -
  #' seed:  random seed
  #' N:     the number of subjects (default is 50)
  #' eqVar: whether the variances are equal.
  #'        equal variance across the four conditions (a1b1, a2b1, a1b2, a2b2), but covariance differs.
  #'        
  #' Output - the data frame for a 2x2 repeated measures design
  set.seed(seed)
  if(eqVar) {
    x <- rnorm(N, 100, 6)     #set origin scores
    x1 <- rnorm(N, 0, 10)     #error variability
    y1 <- x+x1                #first set of base scores
    x2 <- rnorm(N, 0, 10)     #error variability
    y2 <- x+x2                #second set of base scores

    a1 <- rnorm(N, 0, 2.82)   #error variability
    a1b1 <- round(y1+a1)      # first condition

    a1 <- rnorm(N, 1.3, 2.82) #error variability and effect for factor A
    a1b2 <- round(y1+a1)      # third condition

    a2 <- rnorm(N, 15, 2.82)   #error variability and effect for factor B
    a2b1 <- round(y2+a2)       # second condition
    
    a2 <- rnorm(N, 16.3, 2.82) #error variability and effect for factors A and B
    a2b2 <- round(y2+a2)       # fourth condition
    
  } else { #eqVar==F
    x <- rnorm(N, 100, 4.24)  #set origin scores
    x1 <- rnorm(N, 0, 1.41)   #error variability
    a1b1 <- round(x+x1)       #first set of a1 scores
    x2 <- rnorm(N, 1.8, 1.41) #error variability
    a1b2 <- round(x+x2)       #second set of a1 scores
    
    #base scores for a2
    x1 <- rnorm(N, 0, 9.38)   #base error variability
    a2b1 <- x+x1
    a2b2 <- x+x1
    x2 <- rnorm(N, 16, 3.46)    #extra variability and effect
    a2b1 <- round(a2b1+x2)
    x2 <- rnorm(N, 17.8, 3.46)  #extra variability and effect
    a2b2 <- round(a2b2+x2)
  }
  
  data.frame("subj"=paste0("s",1:N),
             "RT"=c(a1b1,a2b1,a1b2,a2b2),
             "resp"=rep(rep(c("a1","a2"),2),each=N),
             "align"=rep(c("b1","b2"),each=N*2),
             stringsAsFactors=T)
}

df <- simData_eqVar()
dat <- matrix(df$RT,ncol=4,
              dimnames=list(paste0("s",1:length(unique(df$subj))),
                            c("onehand&A","twohands&A","onehand&M","twohands&M"))) #wide data format
cov(dat) %>% #sample covariance matrix
  kable(digits=3) %>% kable_styling(full_width=F)

	onehand&A	twohands&A	onehand&M	twohands&M
onehand&A	149.759	56.148	143.079	64.404
twohands&A	56.148	125.894	47.718	124.094
onehand&M	143.079	47.718	148.165	54.731
twohands&M	64.404	124.094	54.731	136.408

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 49  18872   385.1               
## 
## Error: subj:resp
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## resp       1   9800    9800   60.44 4.23e-10 ***
## Residuals 49   7945     162                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: subj:align
##           Df Sum Sq Mean Sq F value Pr(>F)  
## align      1  30.42  30.420   5.201  0.027 *
## Residuals 49 286.58   5.849                 
## ---
## 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   11.5  11.520   1.624  0.208
## Residuals  49  347.5   7.091

Method 1: Compare the model without each effect to the full model in (1).

\[\begin{equation} \tag{1} \mathcal{M}_{\mathrm{full}}:\ Y_{ijk}=\mu+s_k+\alpha_i+\beta_j+(\alpha\beta)_{ij}+\epsilon_{ijk},\quad \epsilon_{ijk}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,\ \sigma_\epsilon^2) \end{equation}\]

set.seed(277)
anovaBF(RT ~ resp * align + subj, data=df, whichRandom="subj", whichModels="top", progress=F)

## Bayes factor top-down analysis
## --------------
## When effect is omitted from resp + align + resp:align + subj , BF is...
## [1] Omit align:resp : 4.568858     ±5.56%
## [2] Omit align      : 6.038909     ±12.25%
## [3] Omit resp       : 8.616201e-24 ±4.05%
## 
## Against denominator:
##   RT ~ resp + align + resp:align + subj 
## ---
## Bayes factor type: BFlinearModel, JZS

Method 2: Include the random slopes in (2) and test each effect again.

\[\begin{equation} \tag{2} \mathcal{M}_{\mathrm{full}}^{'}:\ Y_{ijk}=\mu+s_k+\alpha_i+\beta_j+(\alpha\beta)_{ij}+(\alpha s)_{ik}+(\beta s)_{jk}+\epsilon_{ijk},\quad \epsilon_{ijk}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,\ \sigma_\epsilon^2) \end{equation}\]

Assume $\alpha_i$, $\beta_j$, and $(\alpha\beta)_{ij}$ are the fixed effects; $s_k$ are the subject-specific random effects.

The model that only omits AB in (2) versus the full (2).

randomEffects <- c("subj")

set.seed(277)
full <-   lmBF(RT~resp+align+subj+resp:align+resp:subj+align:subj, df, whichRandom=randomEffects, progress=F)

set.seed(277)
omitAB <- lmBF(RT~resp+align+subj+            resp:subj+align:subj, df, whichRandom=randomEffects, progress=F)
omitAB / full

## Bayes factor analysis
## --------------
## [1] resp + align + subj + resp:subj + align:subj : 1.907953 ±4.71%
## 
## Against denominator:
##   RT ~ resp + align + subj + resp:align + resp:subj + align:subj 
## ---
## Bayes factor type: BFlinearModel, JZS

The model that only omits A in (2) versus the full (2).

set.seed(277)
omitA <-  lmBF(RT~      align+subj+resp:align+resp:subj+align:subj, df, whichRandom=randomEffects, progress=F)
omitA / full

## Bayes factor analysis
## --------------
## [1] align + subj + resp:align + resp:subj + align:subj : 9.735705e-08 ±4.71%
## 
## Against denominator:
##   RT ~ resp + align + subj + resp:align + resp:subj + align:subj 
## ---
## Bayes factor type: BFlinearModel, JZS

The model that only omits B in (2) versus the full (2).

set.seed(277)
omitB <-  lmBF(RT~resp+      subj+resp:align+resp:subj+align:subj, df, whichRandom=randomEffects, progress=F)
omitB / full

## Bayes factor analysis
## --------------
## [1] resp + subj + resp:align + resp:subj + align:subj : 1.025312 ±4.75%
## 
## Against denominator:
##   RT ~ resp + align + subj + resp:align + resp:subj + align:subj 
## ---
## Bayes factor type: BFlinearModel, JZS

Method 3:

The model that excludes A, AB, and A:subj in (2) versus the full (2) that drops AB (i.e., an additive model).

set.seed(277)
notA <-   lmBF(RT~      align+subj+                     align:subj, df, whichRandom=randomEffects, progress=F)
notA / omitAB

## Bayes factor analysis
## --------------
## [1] align + subj + align:subj : 1.816505e-60 ±1.9%
## 
## Against denominator:
##   RT ~ resp + align + subj + resp:subj + align:subj 
## ---
## Bayes factor type: BFlinearModel, JZS

The model that excludes B, AB, and B:subj in (2) versus the full (2) that drops AB.

set.seed(277)
notB <-   lmBF(RT~resp+      subj+            resp:subj         , df, whichRandom=randomEffects, progress=F)
notB / omitAB

## Bayes factor analysis
## --------------
## [1] resp + subj + resp:subj : 53.61516 ±2.02%
## 
## Against denominator:
##   RT ~ resp + align + subj + resp:subj + align:subj 
## ---
## Bayes factor type: BFlinearModel, JZS

Method 4:

The model that excludes A and AB in (2) versus the full (2) that drops AB. Keep all random slopes.

set.seed(277)
exclA <-  lmBF(RT~      align+subj+            resp:subj+align:subj, df, whichRandom=randomEffects, progress=F)
exclA / omitAB

## Bayes factor analysis
## --------------
## [1] align + subj + resp:subj + align:subj : 9.89531e-08 ±2.03%
## 
## Against denominator:
##   RT ~ resp + align + subj + resp:subj + align:subj 
## ---
## Bayes factor type: BFlinearModel, JZS

The model that excludes B and AB in (2) versus the full (2) that drops AB. Keep all random slopes.

set.seed(277)
exclB <-  lmBF(RT~resp+      subj+            resp:subj+align:subj, df, whichRandom=randomEffects, progress=F)
exclB / omitAB

## Bayes factor analysis
## --------------
## [1] resp + subj + resp:subj + align:subj : 1.078329 ±2.05%
## 
## Against denominator:
##   RT ~ resp + align + subj + resp:subj + align:subj 
## ---
## Bayes factor type: BFlinearModel, JZS

Method 1’:

Assume $\alpha_i$, $\beta_j$, $(\alpha\beta)_{ij}$, and $s_k$ are the random effects.

randomEffects <- c("subj","resp","align")

set.seed(277)
full <-   lmBF(RT~resp+align+subj+resp:align                  , df, whichRandom=randomEffects, progress=F)

set.seed(277)
omitAB <- lmBF(RT~resp+align+subj                              , df, whichRandom=randomEffects, progress=F)
omitAB / full

## Bayes factor analysis
## --------------
## [1] resp + align + subj : 8.778604 ±7.96%
## 
## Against denominator:
##   RT ~ resp + align + subj + resp:align 
## ---
## Bayes factor type: BFlinearModel, JZS

set.seed(277)
omitA <-  lmBF(RT~      align+subj+resp:align                  , df, whichRandom=randomEffects, progress=F)
omitA / full

## Bayes factor analysis
## --------------
## [1] align + subj + resp:align : 0.4643873 ±5.98%
## 
## Against denominator:
##   RT ~ resp + align + subj + resp:align 
## ---
## Bayes factor type: BFlinearModel, JZS

set.seed(277)
omitB <-  lmBF(RT~resp+      subj+resp:align                  , df, whichRandom=randomEffects, progress=F)
omitB / full

## Bayes factor analysis
## --------------
## [1] resp + subj + resp:align : 2.445315 ±6.01%
## 
## Against denominator:
##   RT ~ resp + align + subj + resp:align 
## ---
## Bayes factor type: BFlinearModel, JZS

Method 2’:

The model that only omits AB in (2) versus the full (2). All are random effects.

set.seed(277)
full <-   lmBF(RT~resp+align+subj+resp:align+resp:subj+align:subj, df, whichRandom=randomEffects, progress=F)

set.seed(277)
omitAB <- lmBF(RT~resp+align+subj+            resp:subj+align:subj, df, whichRandom=randomEffects, progress=F)
omitAB / full

## Bayes factor analysis
## --------------
## [1] resp + align + subj + resp:subj + align:subj : 3.082228 ±18.61%
## 
## Against denominator:
##   RT ~ resp + align + subj + resp:align + resp:subj + align:subj 
## ---
## Bayes factor type: BFlinearModel, JZS

The model that only omits A in (2) versus the full (2). All are random effects.

set.seed(277)
omitA <-  lmBF(RT~      align+subj+resp:align+resp:subj+align:subj, df, whichRandom=randomEffects, progress=F)
omitA / full

## Bayes factor analysis
## --------------
## [1] align + subj + resp:align + resp:subj + align:subj : 0.09954129 ±59.89%
## 
## Against denominator:
##   RT ~ resp + align + subj + resp:align + resp:subj + align:subj 
## ---
## Bayes factor type: BFlinearModel, JZS

The model that only omits B in (2) versus the full (2). All are random effects.

set.seed(277)
omitB <-  lmBF(RT~resp+      subj+resp:align+resp:subj+align:subj, df, whichRandom=randomEffects, progress=F)
omitB / full

## Bayes factor analysis
## --------------
## [1] resp + subj + resp:align + resp:subj + align:subj : 2.272957 ±23.44%
## 
## Against denominator:
##   RT ~ resp + align + subj + resp:align + resp:subj + align:subj 
## ---
## Bayes factor type: BFlinearModel, JZS

Method 3’:

The model that excludes A, AB, and A:subj in (2) versus the full (2) that drops AB. All are random effects.

set.seed(277)
notA <-   lmBF(RT~      align+subj+                     align:subj, df, whichRandom=randomEffects, progress=F)
notA / omitAB

## Bayes factor analysis
## --------------
## [1] align + subj + align:subj : 3.784507e-61 ±4.34%
## 
## Against denominator:
##   RT ~ resp + align + subj + resp:subj + align:subj 
## ---
## Bayes factor type: BFlinearModel, JZS

The model that excludes B, AB, and B:subj in (2) versus the full (2) that drops AB. All are random effects.

set.seed(277)
notB <-   lmBF(RT~resp+      subj+            resp:subj         , df, whichRandom=randomEffects, progress=F)
notB / omitAB

## Bayes factor analysis
## --------------
## [1] resp + subj + resp:subj : 91.16869 ±5.84%
## 
## Against denominator:
##   RT ~ resp + align + subj + resp:subj + align:subj 
## ---
## Bayes factor type: BFlinearModel, JZS

Method 4’:

The model that excludes A and AB in (2) versus the full (2) that drops AB. All are random effects. Keep all random slopes.

set.seed(277)
exclA <-  lmBF(RT~      align+subj+            resp:subj+align:subj, df, whichRandom=randomEffects, progress=F)
exclA / omitAB

## Bayes factor analysis
## --------------
## [1] align + subj + resp:subj + align:subj : 1.293667e-07 ±26.36%
## 
## Against denominator:
##   RT ~ resp + align + subj + resp:subj + align:subj 
## ---
## Bayes factor type: BFlinearModel, JZS

The model that excludes B and AB in (2) versus the full (2) that drops AB. All are random effects. Keep all random slopes.

set.seed(277)
exclB <-  lmBF(RT~resp+      subj+            resp:subj+align:subj, df, whichRandom=randomEffects, progress=F)
exclB / omitAB

## Bayes factor analysis
## --------------
## [1] resp + subj + resp:subj + align:subj : 1.873388 ±5.54%
## 
## Against denominator:
##   RT ~ resp + align + subj + resp:subj + align:subj 
## ---
## Bayes factor type: BFlinearModel, JZS

2. New (Full) Model

Specify the full model in (3), which has no explicit terms for the subject-specific random effects at all. Compile the model in Stan. $\alpha_i$, $\beta_j$, and $(\alpha\beta)_{ij}$ are coded as random effects (in Sections 2, 3.1, 3.2, & 3.3) for $i,j\in\{1,2\}$.

\[\begin{equation} \tag{3} \mathcal{M}_{\mathrm{full}}^{*}:\ Y_{ijk}=\mu+\alpha_i+\beta_j+(\alpha\beta)_{ij}+\epsilon_{ijk}^{*} \end{equation}\]

\[\begin{equation} \pi(\mu)\propto 1 \end{equation}\]

\[\begin{equation} \alpha_i\mid g_A\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,g_A),\quad\beta_j\mid g_B\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,g_B),\quad(\alpha\beta)_{ij}\mid g_{AB}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,g_{AB}) \end{equation}\]

\[\begin{equation} g_A\sim\text{Scale-inv-}\chi^2(1,h_A^2),\quad g_B\sim\text{Scale-inv-}\chi^2(1,h_B^2),\quad g_{AB}\sim\text{Scale-inv-}\chi^2(1,h_{AB}^2) \end{equation}\]

\[\begin{equation} \boldsymbol{\epsilon}_k=(\epsilon_{11k}^{*},\epsilon_{21k}^{*},\epsilon_{12k}^{*},\epsilon_{22k}^{*})^{\top}\overset{\text{i.i.d.}}{\sim}\boldsymbol{\mathcal{N}}_4(\mathbf{0},\mathbf{\Sigma}) \end{equation}\]

\[\begin{equation} \mathbf{\Sigma}\sim\mathcal{W}^{-1}(\mathbf{I}_{4},5) \end{equation}\]

Note: the Wishart prior is assumed for the precision matrix $\mathbf{\Omega}$ (which is the inverse of the covariance matrix) of the error term; thus, the inverse-Wishart prior is assumed for the covariance matrix $\mathbf{\Sigma}$ of the error term. ChatGPT (GPT-3.5) may confuse them (See pic.).

model.full <- "
data {
  int<lower=1> nS;        // number of subjects
  int<lower=2> nA;        // number of levels of the treatment effect A
  int<lower=2> nB;        // number of levels of the treatment effect B
  // vector[nA*nB] Y[nS];        // responses [removed features]
  array[nS] vector[nA*nB] Y;  // responses [Stan 2.26+ syntax for array declarations]
  real<lower=0> hA;       // (square root) scale parameter of the variance of the treatment effect A
  real<lower=0> hB;       // (square root) scale parameter of the variance of the treatment effect B
  real<lower=0> hAB;      // (square root) scale parameter of the variance of the interaction effect AB
  matrix[nA*nB,nA*nB] I;  // scale matrix of the inverse-Wishart prior
  real lwr;               // lower bound for the grand mean
  real upr;               // upper bound for the grand mean
}

parameters {
  real<lower=lwr, upper=upr> mu;     // grand mean
  cov_matrix[nA*nB] Sigma;           // error covariance matrix
  real<lower=0> gA;                  // variance of the treatment effect A
  real<lower=0> gB;                  // variance of the treatment effect B
  real<lower=0> gAB;                 // variance of the interaction effect AB
  vector[nA] tA;                     // treatment effect A
  vector[nB] tB;                     // treatment effect B
  matrix[nA,nB] tAB;                 // interaction effect AB
}

transformed parameters {
  matrix[nA,nB] muS;     // condition means
  for (i in 1:nA) {
    for (j in 1:nB) {
      muS[i,j] = mu + tA[i] + tB[j] + tAB[i,j];
    }
  }
}

model {
  target += uniform_lpdf(mu | lwr, upr);
  target += scaled_inv_chi_square_lpdf(gA | 1, hA);
  target += scaled_inv_chi_square_lpdf(gB | 1, hB);
  target += scaled_inv_chi_square_lpdf(gAB | 1, hAB);
  target += inv_wishart_lpdf(Sigma | nA*nB+1, I);
  target += normal_lpdf(tA | 0, sqrt(gA));
  target += normal_lpdf(tB | 0, sqrt(gB));
  target += normal_lpdf(to_vector(tAB) | 0, sqrt(gAB));
  for (k in 1:nS) {
    target += multi_normal_lpdf(Y[k] | to_vector(muS), Sigma);
  }
}
"

stanmodel.full <- stan_model(model_code=model.full)
datalist <- list("nA"=2, "nB"=2, "nS"=nrow(dat), "I"=diag(4), "Y"=dat,
                 "hA"=1, "hB"=1, "hAB"=1,
                 "lwr"=mean(dat)-6*sd(dat), "upr"=mean(dat)+6*sd(dat))
stanfit.full <- sampling(stanmodel.full, data=datalist,
                         iter=10000, warmup=1000, chains=2, seed=277, refresh=0)

rstan::summary(stanfit.full)[[1]] %>% kable(digits=4) %>% kable_classic(full_width=F)

	mean	se_mean	sd	2.5%	25%	50%	75%	97.5%	n_eff	Rhat
mu	108.3283	0.1899	10.1122	86.8097	104.3273	108.4582	112.5657	127.6423	2834.078	1.0008
Sigma[1,1]	149.1313	0.3163	30.6889	100.8404	127.3612	144.9396	166.4014	220.7443	9415.259	1.0000
Sigma[1,2]	55.6162	0.2166	21.2513	19.5680	40.7993	53.5925	68.1231	102.9179	9629.485	1.0002
Sigma[1,3]	142.5077	0.3099	29.9364	95.3400	121.3068	138.6213	159.5158	212.1285	9331.114	1.0001
Sigma[1,4]	63.8190	0.2299	22.4335	26.7094	48.2745	61.7318	76.7553	114.2998	9524.601	1.0002
Sigma[2,1]	55.6162	0.2166	21.2513	19.5680	40.7993	53.5925	68.1231	102.9179	9629.485	1.0002
Sigma[2,2]	125.5043	0.2538	25.9453	84.8078	107.0617	122.1224	140.1155	186.0828	10453.108	1.0005
Sigma[2,3]	47.1960	0.2099	20.7510	11.7366	32.8459	45.4795	59.4055	93.1816	9770.532	1.0001
Sigma[2,4]	123.6728	0.2603	26.2651	82.5447	105.0500	120.1788	138.4029	184.0509	10180.905	1.0005
Sigma[3,1]	142.5077	0.3099	29.9364	95.3400	121.3068	138.6213	159.5158	212.1285	9331.114	1.0001
Sigma[3,2]	47.1960	0.2099	20.7510	11.7366	32.8459	45.4795	59.4055	93.1816	9770.532	1.0001
Sigma[3,3]	147.7057	0.3112	30.3907	99.6435	126.1894	143.8310	164.7525	217.9685	9536.992	1.0001
Sigma[3,4]	54.1499	0.2219	21.7851	16.9937	39.0841	52.2865	66.9692	102.3805	9642.387	1.0001
Sigma[4,1]	63.8190	0.2299	22.4335	26.7094	48.2745	61.7318	76.7553	114.2998	9524.601	1.0002
Sigma[4,2]	123.6728	0.2603	26.2651	82.5447	105.0500	120.1788	138.4029	184.0509	10180.905	1.0005
Sigma[4,3]	54.1499	0.2219	21.7851	16.9937	39.0841	52.2865	66.9692	102.3805	9642.387	1.0001
Sigma[4,4]	135.9629	0.2757	28.0422	92.3012	116.0038	132.2575	151.8406	200.5357	10348.408	1.0006
gA	500.6777	229.2782	26146.9728	7.0395	29.3202	62.1959	151.7020	1681.8091	13005.208	1.0001
gB	10.9352	2.9407	252.1595	0.1675	0.4861	1.0266	2.5741	32.6467	7352.663	1.0004
gAB	3.7146	0.5179	19.0830	0.1501	0.3971	0.7711	1.7457	23.2088	1357.869	1.0036
tA[1]	-6.3516	0.1875	9.8978	-25.7639	-10.3671	-6.2426	-2.3725	14.5351	2787.298	1.0009
tA[2]	6.6626	0.1878	9.9239	-11.7855	2.4607	6.3275	10.4383	28.4324	2791.714	1.0007
tB[1]	-0.2333	0.0270	1.6801	-2.9889	-0.8067	-0.2210	0.3413	2.5456	3872.743	1.0005
tB[2]	0.2787	0.0292	1.6706	-2.3487	-0.3222	0.2194	0.8046	3.2893	3269.198	1.0006
tAB[1,1]	-0.1501	0.0311	1.5289	-3.3794	-0.5421	-0.0020	0.4899	1.9144	2421.410	1.0027
tAB[1,2]	-0.2886	0.0373	1.5864	-3.3707	-0.6854	-0.1576	0.3237	1.8718	1804.233	1.0030
tAB[2,1]	-0.0588	0.0445	1.7329	-2.2920	-0.7059	-0.2001	0.3307	3.2028	1515.291	1.0011
tAB[2,2]	0.5692	0.0461	1.8209	-1.5065	-0.1395	0.3703	0.9460	4.1456	1557.831	1.0007
muS[1,1]	101.5932	0.0132	1.7303	98.2023	100.4505	101.5675	102.7393	105.0462	17056.455	1.0000
muS[1,2]	101.9667	0.0133	1.7279	98.5905	100.8214	101.9421	103.1037	105.3999	16869.281	1.0000
muS[2,1]	114.6988	0.0116	1.5871	111.5583	113.6681	114.7062	115.7474	117.8258	18711.788	1.0002
muS[2,2]	115.8389	0.0118	1.6352	112.5685	114.7656	115.8483	116.9342	119.0514	19135.938	1.0002
lp__	-731.4280	0.0821	4.3594	-741.2512	-734.0340	-730.9670	-728.2978	-724.3035	2819.981	1.0022

Note that to compute the (log) marginal likelihood for a Stan model, we need to specify the model in a certain way. Instead of using “~” signs for specifying distributions, we need to directly use the (log) density functions. The reason for this is that when using the “~” sign, constant terms are dropped which are not needed for sampling from the posterior. However, for computing the marginal likelihood, these constants need to be retained. For instance, instead of writing y ~ normal(mu, sigma); we would need to write target += normal_lpdf(y | mu, sigma); (Gronau, 2021). See also the discussion regarding the implicit uniform prior for computing the log marginal likelihood.

The computation of marginal likelihoods based on bridge sampling requires a lot more posterior draws than usual. A good conservative rule of thump is perhaps 10-fold more draws (bridge_sampler.brmsfit {brms}).

Tips: Knit the R markdown in a fresh session without ‘rstan’ loaded. Otherwise, the compiler will generate the clang messages: ld: warning: -undefined dynamic_lookup may not work with chained fixups.

3.1 New (Omitted AB) Model

Specify the additive model in (4).

\[\begin{equation} \tag{4} \mathcal{M}_{\text{omit-}\textit{AB}}^{*}:\ Y_{ijk}=\mu+\alpha_i+\beta_j+\epsilon_{ijk}^{*} \end{equation}\]

model.omitAB <- "
data {
  int<lower=1> nS;        // number of subjects
  int<lower=2> nA;        // number of levels of the treatment effect A
  int<lower=2> nB;        // number of levels of the treatment effect B
  // vector[nA*nB] Y[nS];        // responses [removed features]
  array[nS] vector[nA*nB] Y;  // responses [Stan 2.26+ syntax for array declarations]
  real<lower=0> hA;       // (square root) scale parameter of the variance of the treatment effect A
  real<lower=0> hB;       // (square root) scale parameter of the variance of the treatment effect B
  matrix[nA*nB,nA*nB] I;  // scale matrix of the inverse-Wishart prior
  real lwr;               // lower bound for the grand mean
  real upr;               // upper bound for the grand mean
}

parameters {
  real<lower=lwr, upper=upr> mu;     // grand mean
  cov_matrix[nA*nB] Sigma;           // error covariance matrix
  real<lower=0> gA;                  // variance of the treatment effect A
  real<lower=0> gB;                  // variance of the treatment effect B
  vector[nA] tA;                     // treatment effect A
  vector[nB] tB;                     // treatment effect B
}

transformed parameters {
  matrix[nA,nB] muS;     // condition means
  for (i in 1:nA) {
    for (j in 1:nB) {
      muS[i,j] = mu + tA[i] + tB[j];
    }
  }
}

model {
  target += uniform_lpdf(mu | lwr, upr);
  target += scaled_inv_chi_square_lpdf(gA | 1, hA);
  target += scaled_inv_chi_square_lpdf(gB | 1, hB);
  target += inv_wishart_lpdf(Sigma | nA*nB+1, I);
  target += normal_lpdf(tA | 0, sqrt(gA));
  target += normal_lpdf(tB | 0, sqrt(gB));
  for (k in 1:nS) {
    target += multi_normal_lpdf(Y[k] | to_vector(muS), Sigma);
  }
}
"

stanmodel.omitAB <- stan_model(model_code=model.omitAB)
stanfit.omitAB <- sampling(stanmodel.omitAB, data=datalist[-8],
                           iter=10000, warmup=1000, chains=2, seed=277, refresh=0)

rstan::summary(stanfit.omitAB)[[1]] %>% kable(digits=4) %>% kable_classic(full_width=F)

	mean	se_mean	sd	2.5%	25%	50%	75%	97.5%	n_eff	Rhat
mu	107.8329	0.2106	10.2759	86.8231	103.8407	108.1250	112.1412	126.3842	2381.885	1.0001
Sigma[1,1]	149.8113	0.3621	30.8913	101.6782	127.9885	145.8478	167.0924	222.1465	7279.253	1.0000
Sigma[1,2]	56.0647	0.2464	21.6926	19.3233	41.2911	53.9543	68.3600	104.8603	7750.505	1.0000
Sigma[1,3]	143.1273	0.3548	30.1214	96.3507	121.9442	139.2978	159.9424	213.1272	7208.034	1.0001
Sigma[1,4]	64.2515	0.2624	22.9936	25.7425	48.4789	61.9661	77.3017	116.8069	7680.111	1.0000
Sigma[2,1]	56.0647	0.2464	21.6926	19.3233	41.2911	53.9543	68.3600	104.8603	7750.505	1.0000
Sigma[2,2]	126.0383	0.2724	25.9157	85.3441	107.7657	122.6632	140.8945	186.1463	9050.056	1.0001
Sigma[2,3]	47.5291	0.2396	21.0562	11.2537	33.2770	45.7950	59.5670	94.2544	7723.024	1.0000
Sigma[2,4]	124.3544	0.2794	26.2985	82.6833	105.8319	120.8993	139.1904	184.9151	8858.813	1.0001
Sigma[3,1]	143.1273	0.3548	30.1214	96.3507	121.9442	139.2978	159.9424	213.1272	7208.034	1.0001
Sigma[3,2]	47.5291	0.2396	21.0562	11.2537	33.2770	45.7950	59.5670	94.2544	7723.024	1.0000
Sigma[3,3]	148.2919	0.3542	30.5804	101.0046	126.6749	144.1858	165.3353	219.0978	7455.459	1.0001
Sigma[3,4]	54.3760	0.2532	22.2452	16.5826	39.3903	52.4066	67.2715	104.3041	7720.003	1.0000
Sigma[4,1]	64.2515	0.2624	22.9936	25.7425	48.4789	61.9661	77.3017	116.8069	7680.111	1.0000
Sigma[4,2]	124.3544	0.2794	26.2985	82.6833	105.8319	120.8993	139.1904	184.9151	8858.813	1.0001
Sigma[4,3]	54.3760	0.2532	22.2452	16.5826	39.3903	52.4066	67.2715	104.3041	7720.003	1.0000
Sigma[4,4]	136.9913	0.2950	28.1643	92.2610	117.2261	133.3968	152.6978	202.1085	9115.290	1.0001
gA	255.6183	15.7918	1138.4830	10.8598	31.8035	65.0156	155.1266	1625.1098	5197.422	1.0002
gB	5.1431	1.0174	90.9782	0.1677	0.4564	0.9093	2.2031	23.2145	7996.674	1.0001
tA[1]	-6.3150	0.2067	10.1276	-24.9989	-10.4472	-6.4370	-2.5031	14.0624	2399.634	1.0000
tA[2]	7.0364	0.2076	10.1410	-11.0479	2.8908	6.6956	10.7411	27.8081	2387.184	1.0001
tB[1]	-0.3261	0.0293	1.3815	-2.5923	-0.8221	-0.3364	0.1380	2.2232	2228.384	1.0002
tB[2]	0.3721	0.0294	1.3827	-1.8461	-0.1247	0.3470	0.8381	2.9787	2219.379	1.0002
muS[1,1]	101.1918	0.0129	1.7089	97.8218	100.0599	101.1936	102.3069	104.5363	17463.223	1.0001
muS[1,2]	101.8900	0.0132	1.7483	98.4611	100.7266	101.8878	103.0299	105.3136	17638.019	1.0001
muS[2,1]	114.5432	0.0119	1.6041	111.3864	113.4759	114.5270	115.6314	117.6676	18024.586	0.9999
muS[2,2]	115.2414	0.0117	1.5710	112.1427	114.1922	115.2397	116.3042	118.3176	18009.525	0.9999
lp__	-724.8943	0.0508	3.4035	-732.7079	-726.9609	-724.4740	-722.4391	-719.4517	4489.907	1.0002

3.2 New (Omitted A) Model

New Method 2:

Specify the model in (5.1), which violates the principle of marginality.

\[\begin{equation} \tag{5.1} \mathcal{M}_{\text{omit-}\textit{A}}^{*}:\ Y_{ijk}=\mu+\beta_j+(\alpha\beta)_{ij}+\epsilon_{ijk}^{*} \end{equation}\]

model.omitA <- "
data {
  int<lower=1> nS;        // number of subjects
  int<lower=2> nA;        // number of levels of the treatment effect A
  int<lower=2> nB;        // number of levels of the treatment effect B
  // vector[nA*nB] Y[nS];        // responses [removed features]
  array[nS] vector[nA*nB] Y;  // responses [Stan 2.26+ syntax for array declarations]
  real<lower=0> hB;       // (square root) scale parameter of the variance of the treatment effect B
  real<lower=0> hAB;      // (square root) scale parameter of the variance of the interaction effect AB
  matrix[nA*nB,nA*nB] I;  // scale matrix of the inverse-Wishart prior
  real lwr;               // lower bound for the grand mean
  real upr;               // upper bound for the grand mean
}

parameters {
  real<lower=lwr, upper=upr> mu;     // grand mean
  cov_matrix[nA*nB] Sigma;           // error covariance matrix
  real<lower=0> gB;                  // variance of the treatment effect B
  real<lower=0> gAB;                 // variance of the interaction effect AB
  vector[nB] tB;                     // treatment effect B
  matrix[nA,nB] tAB;                 // interaction effect AB
}

transformed parameters {
  matrix[nA,nB] muS;     // condition means
  for (i in 1:nA) {
    for (j in 1:nB) {
      muS[i,j] = mu + tB[j] + tAB[i,j];
    }
  }
}

model {
  target += uniform_lpdf(mu | lwr, upr);
  target += scaled_inv_chi_square_lpdf(gB | 1, hB);
  target += scaled_inv_chi_square_lpdf(gAB | 1, hAB);
  target += inv_wishart_lpdf(Sigma | nA*nB+1, I);
  target += normal_lpdf(tB | 0, sqrt(gB));
  target += normal_lpdf(to_vector(tAB) | 0, sqrt(gAB));
  for (k in 1:nS) {
    target += multi_normal_lpdf(Y[k] | to_vector(muS), Sigma);
  }
}
"

stanmodel.omitA <- stan_model(model_code=model.omitA)
stanfit.omitA <- sampling(stanmodel.omitA, data=datalist[-6],
                          iter=10000, warmup=1000, chains=2, seed=277, refresh=0)

rstan::summary(stanfit.omitA)[[1]] %>% kable(digits=4) %>% kable_classic(full_width=F)

	mean	se_mean	sd	2.5%	25%	50%	75%	97.5%	n_eff	Rhat
mu	108.8570	0.2074	7.2071	96.8994	105.6862	108.6326	111.6701	120.8202	1207.508	1.0003
Sigma[1,1]	150.1278	0.3350	31.4058	100.8737	127.4203	145.8401	168.2853	222.0518	8786.558	1.0000
Sigma[1,2]	55.9134	0.2483	21.7837	18.9731	40.8001	53.7640	68.7732	104.7440	7697.121	1.0003
Sigma[1,3]	143.4796	0.3249	30.6225	95.2906	121.4760	139.3979	161.2909	214.4826	8881.236	1.0000
Sigma[1,4]	64.1523	0.2657	23.0546	25.6775	48.0693	61.7420	77.8160	115.5844	7527.427	1.0003
Sigma[2,1]	55.9134	0.2483	21.7837	18.9731	40.8001	53.7640	68.7732	104.7440	7697.121	1.0003
Sigma[2,2]	126.2547	0.2771	26.5184	84.7503	107.3945	122.6952	141.2330	188.4768	9156.251	1.0001
Sigma[2,3]	47.4476	0.2365	21.1705	10.8630	32.8645	45.6839	60.0616	94.6711	8011.600	1.0002
Sigma[2,4]	124.4401	0.2867	26.9083	82.2036	105.1744	120.7883	139.8502	186.9690	8806.855	1.0001
Sigma[3,1]	143.4796	0.3249	30.6225	95.2906	121.4760	139.3979	161.2909	214.4826	8881.236	1.0000
Sigma[3,2]	47.4476	0.2365	21.1705	10.8630	32.8645	45.6839	60.0616	94.6711	8011.600	1.0002
Sigma[3,3]	148.6552	0.3237	31.0772	99.7615	126.1097	144.5331	166.7051	219.9229	9216.940	1.0000
Sigma[3,4]	54.4569	0.2529	22.3354	16.5605	39.1074	52.4123	67.6318	104.3596	7798.526	1.0002
Sigma[4,1]	64.1523	0.2657	23.0546	25.6775	48.0693	61.7420	77.8160	115.5844	7527.427	1.0003
Sigma[4,2]	124.4401	0.2867	26.9083	82.2036	105.1744	120.7883	139.8502	186.9690	8806.855	1.0001
Sigma[4,3]	54.4569	0.2529	22.3354	16.5605	39.1074	52.4123	67.6318	104.3596	7798.526	1.0002
Sigma[4,4]	136.8081	0.3056	28.8333	92.0410	116.0410	132.8955	153.3666	202.5049	8902.726	1.0000
gB	17.8383	3.3571	156.0513	0.1937	0.6547	1.6272	5.1314	103.4225	2160.708	1.0013
gAB	115.4627	16.8614	584.6195	13.1107	31.6390	53.2156	97.4966	471.3135	1202.157	1.0015
tB[1]	-0.1268	0.0839	3.2318	-5.7617	-0.8433	-0.0256	0.7768	4.8248	1485.123	1.0019
tB[2]	-0.0533	0.0854	3.1611	-5.2261	-0.8152	0.0177	0.8590	5.2014	1370.351	1.0023
tAB[1,1]	-6.7780	0.2024	6.6908	-18.1631	-9.3538	-6.4973	-3.7994	3.6868	1093.139	1.0012
tAB[1,2]	-6.5084	0.2025	6.7261	-18.0178	-9.0260	-6.1991	-3.5151	4.1199	1103.052	1.0014
tAB[2,1]	5.7758	0.2001	6.6597	-4.6387	3.0651	5.7678	8.6240	16.9349	1107.969	1.0011
tAB[2,2]	6.9719	0.1989	6.6753	-3.6385	4.2808	7.0438	9.8151	18.2590	1126.782	1.0013
muS[1,1]	101.9522	0.0130	1.7794	98.5025	100.7603	101.9368	103.1333	105.5200	18663.221	1.0000
muS[1,2]	102.2953	0.0130	1.7762	98.8627	101.1095	102.2596	103.4625	105.8949	18582.802	1.0000
muS[2,1]	114.5060	0.0117	1.6025	111.2631	113.4568	114.5469	115.5923	117.5385	18599.572	0.9999
muS[2,2]	115.7756	0.0120	1.6606	112.4457	114.6900	115.8031	116.8871	118.9614	19035.295	0.9999
lp__	-732.1785	0.0732	3.9217	-741.1771	-734.4846	-731.7256	-729.3525	-725.9118	2873.991	1.0002

New Method 4:

Specify the model in (5.2), which follows the principle of marginality.

\[\begin{equation} \tag{5.2} \mathcal{M}_{\text{not-}\textit{A}}^{*}:\ Y_{ijk}=\mu+\beta_j+\epsilon_{ijk}^{*} \end{equation}\]

model.notA <- "
data {
  int<lower=1> nS;        // number of subjects
  int<lower=2> nA;        // number of levels of the treatment effect A
  int<lower=2> nB;        // number of levels of the treatment effect B
  // vector[nA*nB] Y[nS];        // responses [removed features]
  array[nS] vector[nA*nB] Y;  // responses [Stan 2.26+ syntax for array declarations]
  real<lower=0> hB;       // (square root) scale parameter of the variance of the treatment effect B
  matrix[nA*nB,nA*nB] I;  // scale matrix of the inverse-Wishart prior
  real lwr;               // lower bound for the grand mean
  real upr;               // upper bound for the grand mean
}

parameters {
  real<lower=lwr, upper=upr> mu;     // grand mean
  cov_matrix[nA*nB] Sigma;           // error covariance matrix
  real<lower=0> gB;                  // variance of the treatment effect B
  vector[nB] tB;                     // treatment effect B
}

transformed parameters {
  matrix[nA,nB] muS;     // condition means
  for (i in 1:nA) {
    for (j in 1:nB) {
      muS[i,j] = mu + tB[j];
    }
  }
}

model {
  target += uniform_lpdf(mu | lwr, upr);
  target += scaled_inv_chi_square_lpdf(gB | 1, hB);
  target += inv_wishart_lpdf(Sigma | nA*nB+1, I);
  target += normal_lpdf(tB | 0, sqrt(gB));
  for (k in 1:nS) {
    target += multi_normal_lpdf(Y[k] | to_vector(muS), Sigma);
  }
}
"

stanmodel.notA <- stan_model(model_code=model.notA)
stanfit.notA <- sampling(stanmodel.notA, data=datalist[-c(6,8)],
                         iter=10000, warmup=1000, chains=2, seed=277, refresh=0)

rstan::summary(stanfit.notA)[[1]] %>% kable(digits=4) %>% kable_classic(full_width=F)

	mean	se_mean	sd	2.5%	25%	50%	75%	97.5%	n_eff	Rhat
mu	109.1766	0.0372	2.4509	104.4597	107.6734	109.1837	110.6824	113.9660	4351.373	1.0001
Sigma[1,1]	202.9448	0.5593	51.7103	124.6038	166.1069	195.6027	231.0731	325.8826	8549.026	1.0001
Sigma[1,2]	13.7339	0.2521	26.9924	-37.3588	-3.9048	13.0558	30.4343	69.5340	11464.355	1.0003
Sigma[1,3]	201.8827	0.5714	52.7683	121.6861	164.1968	194.0698	230.6216	326.8394	8528.627	1.0000
Sigma[1,4]	20.2582	0.2605	28.0382	-32.4810	1.8351	19.5063	37.5653	78.4112	11587.453	1.0002
Sigma[2,1]	13.7339	0.2521	26.9924	-37.3588	-3.9048	13.0558	30.4343	69.5340	11464.355	1.0003
Sigma[2,2]	167.3041	0.4855	43.3241	102.8688	136.5036	160.7432	190.3168	270.3290	7964.395	1.0003
Sigma[2,3]	0.5119	0.2600	27.9348	-53.7753	-17.3895	0.2782	17.9663	56.8486	11547.594	1.0004
Sigma[2,4]	166.6138	0.4929	43.9135	101.3617	135.3652	159.7530	190.0752	271.3644	7938.383	1.0003
Sigma[3,1]	201.8827	0.5714	52.7683	121.6861	164.1968	194.0698	230.6216	326.8394	8528.627	1.0000
Sigma[3,2]	0.5119	0.2600	27.9348	-53.7753	-17.3895	0.2782	17.9663	56.8486	11547.594	1.0004
Sigma[3,3]	213.1832	0.5900	55.1531	128.6730	174.1119	204.7930	243.2051	343.1148	8739.596	1.0000
Sigma[3,4]	5.9115	0.2644	28.4823	-48.9107	-12.4665	5.3175	23.5872	63.7173	11600.262	1.0003
Sigma[4,1]	20.2582	0.2605	28.0382	-32.4810	1.8351	19.5063	37.5653	78.4112	11587.453	1.0002
Sigma[4,2]	166.6138	0.4929	43.9135	101.3617	135.3652	159.7530	190.0752	271.3644	7938.383	1.0003
Sigma[4,3]	5.9115	0.2644	28.4823	-48.9107	-12.4665	5.3175	23.5872	63.7173	11600.262	1.0003
Sigma[4,4]	180.0722	0.5079	45.9522	111.4526	147.4856	173.1799	204.8056	289.5344	8187.118	1.0004
gB	5.3845	0.4909	37.5494	0.1921	0.5677	1.1857	2.9158	30.6668	5851.447	1.0000
tB[1]	-0.5209	0.0303	1.4425	-3.3116	-1.0715	-0.5037	0.0200	2.1392	2261.479	0.9999
tB[2]	0.5472	0.0305	1.4406	-2.1274	-0.0181	0.5056	1.0794	3.2911	2236.316	1.0000
muS[1,1]	108.6557	0.0211	1.9999	104.7819	107.3134	108.6505	109.9881	112.6168	9020.619	1.0001
muS[1,2]	109.7238	0.0208	1.9970	105.8442	108.3967	109.7219	111.0472	113.6363	9260.724	1.0002
muS[2,1]	108.6557	0.0211	1.9999	104.7819	107.3134	108.6505	109.9881	112.6168	9020.619	1.0001
muS[2,2]	109.7238	0.0208	1.9970	105.8442	108.3967	109.7219	111.0472	113.6363	9260.724	1.0002
lp__	-735.4252	0.0389	2.9592	-742.1842	-737.2036	-735.0560	-733.2440	-730.7308	5772.335	1.0001

3.3 New (Omitted B) Model

New Method 2:

Specify the model in (6.1), which violates the principle of marginality.

\[\begin{equation} \tag{6.1} \mathcal{M}_{\text{omit-}\textit{B}}^{*}:\ Y_{ijk}=\mu+\alpha_i+(\alpha\beta)_{ij}+\epsilon_{ijk}^{*} \end{equation}\]

datalist[["Y"]] <- dat[,c(1,3,2,4)]
stanfit.omitB <- sampling(stanmodel.omitA, data=datalist[-6],
                          iter=10000, warmup=1000, chains=2, seed=277, refresh=0)

## Warning: There were 120 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.

## Warning: Examine the pairs() plot to diagnose sampling problems

rstan::summary(stanfit.omitB)[[1]] %>% kable(digits=4) %>% kable_classic(full_width=F)

	mean	se_mean	sd	2.5%	25%	50%	75%	97.5%	n_eff	Rhat
mu	108.4931	0.1239	9.1886	91.0056	104.4606	108.4408	112.4009	127.4397	5497.1836	1.0003
Sigma[1,1]	149.4138	0.3218	30.6246	100.9990	127.7086	145.4398	166.8429	219.3507	9058.6742	1.0000
Sigma[1,2]	142.6794	0.3132	29.8521	95.2112	121.6321	138.7485	159.6829	211.1870	9084.3159	1.0001
Sigma[1,3]	55.7526	0.2395	21.1661	19.8806	41.4522	53.6766	68.0854	103.2859	7811.3955	1.0002
Sigma[1,4]	63.9345	0.2521	22.3618	26.5407	48.5871	61.6417	76.8384	114.4616	7870.6033	1.0002
Sigma[2,1]	142.6794	0.3132	29.8521	95.2112	121.6321	138.7485	159.6829	211.1870	9084.3159	1.0001
Sigma[2,2]	147.6906	0.3132	30.2834	99.8715	126.3601	143.7970	164.8308	217.2387	9347.9158	1.0001
Sigma[2,3]	47.3628	0.2292	20.5667	11.0992	33.6723	45.5822	59.6569	92.5271	8049.4138	1.0000
Sigma[2,4]	54.3080	0.2427	21.6321	17.0919	39.6755	52.3743	66.9500	102.6401	7943.5161	1.0001
Sigma[3,1]	55.7526	0.2395	21.1661	19.8806	41.4522	53.6766	68.0854	103.2859	7811.3955	1.0002
Sigma[3,2]	47.3628	0.2292	20.5667	11.0992	33.6723	45.5822	59.6569	92.5271	8049.4138	1.0000
Sigma[3,3]	125.6863	0.2789	25.5545	85.1402	107.7046	122.5686	140.3558	184.6282	8393.6109	1.0002
Sigma[3,4]	123.7566	0.2769	25.7694	82.7389	105.5970	120.5939	138.0962	183.2224	8663.2001	1.0002
Sigma[4,1]	63.9345	0.2521	22.3618	26.5407	48.5871	61.6417	76.8384	114.4616	7870.6033	1.0002
Sigma[4,2]	54.3080	0.2427	21.6321	17.0919	39.6755	52.3743	66.9500	102.6401	7943.5161	1.0001
Sigma[4,3]	123.7566	0.2769	25.7694	82.7389	105.5970	120.5939	138.0962	183.2224	8663.2001	1.0002
Sigma[4,4]	135.9718	0.2906	27.4955	92.1040	116.6235	132.5999	151.4312	199.1441	8955.2541	1.0002
gB	235.5785	18.7341	1709.0713	2.6900	28.4211	59.8808	144.8551	1184.5939	8322.5604	1.0000
gAB	3.6821	0.9155	19.0295	0.1660	0.4287	0.7933	1.7061	34.8140	432.0878	1.0044
tB[1]	-6.4986	0.1283	9.1670	-25.5496	-10.3845	-6.2506	-2.2380	10.3776	5101.7428	1.0005
tB[2]	6.4109	0.1283	9.1851	-11.7681	2.2559	6.2558	10.3459	24.3199	5121.4081	1.0007
tAB[1,1]	-0.3999	0.0358	1.3093	-3.3970	-0.7646	-0.2376	0.2106	1.3804	1337.2952	1.0023
tAB[1,2]	-0.1378	0.1382	1.7916	-2.1725	-0.8631	-0.3915	0.1018	5.1176	167.9907	1.0146
tAB[2,1]	-0.1008	0.0354	1.2998	-3.0669	-0.4558	0.0426	0.4955	1.7228	1347.3872	1.0024
tAB[2,2]	0.9225	0.1485	1.8776	-0.9832	0.1365	0.5906	1.1478	6.3843	159.7590	1.0152
muS[1,1]	101.5945	0.0143	1.7408	98.1848	100.4349	101.5877	102.7521	105.0231	14872.2917	1.0003
muS[1,2]	114.7662	0.0134	1.5864	111.6020	113.7145	114.7738	115.8459	117.8205	13981.3329	1.0000
muS[2,1]	101.8937	0.0143	1.7376	98.5015	100.7311	101.8971	103.0536	105.2981	14760.7490	1.0004
muS[2,2]	115.8265	0.0128	1.6349	112.5595	114.7594	115.8373	116.9251	118.9876	16380.0605	0.9999
lp__	-727.8921	0.0631	3.6901	-736.1131	-730.2007	-727.5302	-725.1700	-721.7541	3416.9010	1.0008

New Method 4:

Specify the model in (6.2), which follows the principle of marginality.

\[\begin{equation} \tag{6.2} \mathcal{M}_{\text{not-}\textit{B}}^{*}:\ Y_{ijk}=\mu+\alpha_i+\epsilon_{ijk}^{*} \end{equation}\]

stanfit.notB <- sampling(stanmodel.notA, data=datalist[-c(6,8)],
                         iter=10000, warmup=1000, chains=2, seed=277, refresh=0)
rstan::summary(stanfit.notB)[[1]] %>% kable(digits=4) %>% kable_classic(full_width=F)

	mean	se_mean	sd	2.5%	25%	50%	75%	97.5%	n_eff	Rhat
mu	107.8775	0.2604	10.9702	84.9969	103.9839	108.1194	112.2339	127.4964	1774.731	1.0000
Sigma[1,1]	150.0577	0.3192	31.0598	101.5734	127.9586	145.9650	167.5284	223.7547	9466.631	0.9999
Sigma[1,2]	143.5478	0.3123	30.3542	96.0937	121.9501	139.4146	160.3996	215.6043	9444.570	0.9999
Sigma[1,3]	56.0581	0.2289	21.5328	19.4103	41.0595	54.2558	68.9823	104.1878	8851.935	1.0002
Sigma[1,4]	64.5230	0.2441	22.8993	26.1185	48.6196	62.3010	78.0660	115.8330	8800.354	1.0003
Sigma[2,1]	143.5478	0.3123	30.3542	96.0937	121.9501	139.4146	160.3996	215.6043	9444.570	0.9999
Sigma[2,2]	148.6907	0.3125	30.8699	100.6296	126.8574	144.4138	165.8906	221.8429	9759.844	0.9999
Sigma[2,3]	47.7204	0.2265	21.1169	11.2699	33.0858	45.9956	60.3098	94.2960	8690.748	1.0004
Sigma[2,4]	55.3286	0.2416	22.3933	17.1764	39.8021	53.4359	68.3684	105.6727	8594.376	1.0005
Sigma[3,1]	56.0581	0.2289	21.5328	19.4103	41.0595	54.2558	68.9823	104.1878	8851.935	1.0002
Sigma[3,2]	47.7204	0.2265	21.1169	11.2699	33.0858	45.9956	60.3098	94.2960	8690.748	1.0004
Sigma[3,3]	126.1989	0.2628	26.2839	85.2870	107.4949	122.5812	141.0154	187.0899	10004.555	1.0000
Sigma[3,4]	124.0858	0.2684	26.5588	82.7121	105.2919	120.4184	139.1337	187.2344	9788.922	1.0000
Sigma[4,1]	64.5230	0.2441	22.8993	26.1185	48.6196	62.3010	78.0660	115.8330	8800.354	1.0003
Sigma[4,2]	55.3286	0.2416	22.3933	17.1764	39.8021	53.4359	68.3684	105.6727	8594.376	1.0005
Sigma[4,3]	124.0858	0.2684	26.5588	82.7121	105.2919	120.4184	139.1337	187.2344	9788.922	1.0000
Sigma[4,4]	137.4135	0.2855	28.6196	92.3313	117.0271	133.4810	153.5063	204.3276	10047.633	1.0000
gB	510.2030	143.9676	16156.1558	11.9564	34.0326	69.5810	170.8378	1888.1053	12593.507	1.0000
tB[1]	-6.7055	0.2590	10.9156	-26.3295	-10.9011	-6.8143	-2.8560	15.8689	1775.570	1.0000
tB[2]	7.3136	0.2602	10.9255	-11.9903	3.0524	6.9680	11.0257	30.6909	1763.331	1.0000
muS[1,1]	101.1719	0.0131	1.7715	97.6936	99.9929	101.1809	102.3514	104.6189	18367.690	1.0000
muS[1,2]	115.1911	0.0123	1.6149	112.0308	114.1208	115.1852	116.2753	118.3385	17109.799	0.9999
muS[2,1]	101.1719	0.0131	1.7715	97.6936	99.9929	101.1809	102.3514	104.6189	18367.690	1.0000
muS[2,2]	115.1911	0.0123	1.6149	112.0308	114.1208	115.1852	116.2753	118.3385	17109.799	0.9999
lp__	-723.7619	0.0461	3.0351	-730.8352	-725.4967	-723.3758	-721.5812	-719.0658	4334.262	1.0001

4. Bayes Factor via bridge Sampling

Compute log marginal likelihoods for competing models via bridge sampling. Then, $\textit{BF}_{10}=\frac{p(\text{Data}\ \mid\ \mathcal{M}_1)}{p(\text{Data}\ \mid\ \mathcal{M}_0)}$.

set.seed(277)
M.full <- bridge_sampler(stanfit.full, silent=T)
summary(M.full)

## 
## Bridge sampling log marginal likelihood estimate 
## (method = "normal", repetitions = 1):
## 
##  -716.3043
## 
## Error Measures:
## 
##  Relative Mean-Squared Error: 0.0004457507
##  Coefficient of Variation: 0.02111281
##  Percentage Error: 2%
## 
## Note:
## All error measures are approximate.

set.seed(277)
M.omitAB <- bridge_sampler(stanfit.omitAB, silent=T)
summary(M.omitAB)

## 
## Bridge sampling log marginal likelihood estimate 
## (method = "normal", repetitions = 1):
## 
##  -715.5448
## 
## Error Measures:
## 
##  Relative Mean-Squared Error: 0.0002005177
##  Coefficient of Variation: 0.01416043
##  Percentage Error: 1%
## 
## Note:
## All error measures are approximate.

bf(M.omitAB, M.full) #New Method of testing AB

## Estimated Bayes factor in favor of M.omitAB over M.full: 2.13733

set.seed(277)
M.omitA <- bridge_sampler(stanfit.omitA, silent=T)
summary(M.omitA)

## 
## Bridge sampling log marginal likelihood estimate 
## (method = "normal", repetitions = 1):
## 
##  -719.7991
## 
## Error Measures:
## 
##  Relative Mean-Squared Error: 0.0002807846
##  Coefficient of Variation: 0.01675663
##  Percentage Error: 2%
## 
## Note:
## All error measures are approximate.

bf(M.omitA, M.full) #New Method 2 of testing A

## Estimated Bayes factor in favor of M.omitA over M.full: 0.03036

set.seed(277)
M.notA <- bridge_sampler(stanfit.notA, silent=T)
summary(M.notA)

## 
## Bridge sampling log marginal likelihood estimate 
## (method = "normal", repetitions = 1):
## 
##  -731.8172
## 
## Error Measures:
## 
##  Relative Mean-Squared Error: 6.676433e-05
##  Coefficient of Variation: 0.008170944
##  Percentage Error: 1%
## 
## Note:
## All error measures are approximate.

bf(M.notA, M.omitAB) #New Method 4 of testing A

## Estimated Bayes factor in favor of M.notA over M.omitAB: 0.00000

set.seed(277)
M.omitB <- bridge_sampler(stanfit.omitB, silent=T)
summary(M.omitB)

## 
## Bridge sampling log marginal likelihood estimate 
## (method = "normal", repetitions = 1):
## 
##  -715.7094
## 
## Error Measures:
## 
##  Relative Mean-Squared Error: 0.0002551457
##  Coefficient of Variation: 0.01597328
##  Percentage Error: 2%
## 
## Note:
## All error measures are approximate.

bf(M.omitB, M.full) #New Method 2 of testing B

## Estimated Bayes factor in favor of M.omitB over M.full: 1.81288

set.seed(277)
M.notB <- bridge_sampler(stanfit.notB, silent=T)
summary(M.notB)

## 
## Bridge sampling log marginal likelihood estimate 
## (method = "normal", repetitions = 1):
## 
##  -716.0738
## 
## Error Measures:
## 
##  Relative Mean-Squared Error: 0.0002693148
##  Coefficient of Variation: 0.01641081
##  Percentage Error: 2%
## 
## Note:
## All error measures are approximate.

bf(M.notB, M.omitAB) #New Method 4 of testing B

## Estimated Bayes factor in favor of M.notB over M.omitAB: 0.58918

5. Modeling Fixed Effects

Specify the full model in (3). Compile the model (which only works for $2\times2$) in Stan. $\alpha_i$, $\beta_j$, and $(\alpha\beta)_{ij}$ are now coded as fixed effects by projecting a set of $a$ main effects into $a-1$ parameters with the property that the marginal prior on all $a$ main effects is identical.

\[\begin{equation} \tag{3} \mathcal{M}_{\mathrm{full}}^{*}:\ Y_{ijk}=\mu+\alpha_i+\beta_j+(\alpha\beta)_{ij}+\epsilon_{ijk}^{*} \end{equation}\]

\[\begin{equation} \pi(\mu)\propto 1 \end{equation}\]

\[\begin{equation} \tag{$\star$} (\alpha_1^\star,\dotsb,\alpha_{a-1}^\star)=(\alpha_1,\dotsb,\alpha_a)\cdot\mathbf{Q},\quad\mathbf{I}_a-a^{-1}\mathbf{J}_a=\mathbf{Q}\cdot\mathbf{Q}^\top,\quad\dotsb \end{equation}\]

\[\begin{equation} \alpha_i^\star\mid g_A\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,g_A),\quad\beta_j^\star\mid g_B\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,g_B),\quad(\alpha\beta)_{ij}^\star\mid g_{AB}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,g_{AB}) \end{equation}\]

\[\begin{equation} g_A\sim\text{Scale-inv-}\chi^2(1,h_A^2),\quad g_B\sim\text{Scale-inv-}\chi^2(1,h_B^2),\quad g_{AB}\sim\text{Scale-inv-}\chi^2(1,h_{AB}^2) \end{equation}\]

\[\begin{equation} \boldsymbol{\epsilon}_k=(\epsilon_{11k}^{*},\epsilon_{21k}^{*},\epsilon_{12k}^{*},\epsilon_{22k}^{*})^{\top}\overset{\text{i.i.d.}}{\sim}\boldsymbol{\mathcal{N}}_4(\mathbf{0},\mathbf{\Sigma}) \end{equation}\]

\[\begin{equation} \mathbf{\Sigma}\sim\mathcal{W}^{-1}(\mathbf{I}_{4},5) \end{equation}\]

$\mathbf{Q}$ is an $a\times(a-1)$ matrix of the $a-1$ eigenvectors of unit length corresponding to the nonzero eigenvalues of the left side term. For example, the projected main effect $\alpha^\star=(\alpha_1-\alpha_2)/\sqrt{2}$ when $a=2$. In the other direction (given $\alpha_1+\alpha_2=0$), $\alpha_1=\alpha^\star/\sqrt{2}$ and $\alpha_2=-\alpha^\star/\sqrt{2}$. The same reduced parametrization applies to $\beta_j$.

As for the interaction effects, four variables are projected into one. $(\alpha\beta)^\star=((\alpha\beta)_{11}-(\alpha\beta)_{21}-(\alpha\beta)_{12}+(\alpha\beta)_{22})/2$. In the other direction (given $(\alpha\beta)_{11}+(\alpha\beta)_{12}=0$, $(\alpha\beta)_{21}+(\alpha\beta)_{22}=0$, and $(\alpha\beta)_{11}+(\alpha\beta)_{21}=0$), $(\alpha\beta)_{11}=(\alpha\beta)^\star/2$, $(\alpha\beta)_{21}=-(\alpha\beta)^\star/2$, $(\alpha\beta)_{12}=-(\alpha\beta)^\star/2$, and $(\alpha\beta)_{22}=(\alpha\beta)^\star/2$.

model.full.fixed <- "
data {
  int<lower=1> nS;        // number of subjects
  // vector[4] Y[nS];        // responses [removed features]
  array[nS] vector[4] Y;  // responses [Stan 2.26+ syntax for array declarations]
  real<lower=0> hA;       // (square root) scale parameter of the variance of the projected treatment effect A
  real<lower=0> hB;       // (square root) scale parameter of the variance of the projected treatment effect B
  real<lower=0> hAB;      // (square root) scale parameter of the variance of the projected interaction effect AB
  matrix[4,4] I;          // scale matrix of the inverse-Wishart prior
  real lwr;               // lower bound for the grand mean
  real upr;               // upper bound for the grand mean
}

parameters {
  real<lower=lwr, upper=upr> mu;     // grand mean
  cov_matrix[4] Sigma;               // error covariance matrix
  real<lower=0> gA;                  // variance of the projected treatment effect A
  real<lower=0> gB;                  // variance of the projected treatment effect B
  real<lower=0> gAB;                 // variance of the projected interaction effect AB
  real tAfix;                        // projected treatment effect A
  real tBfix;                        // projected treatment effect B
  real tABfix;                       // projected interaction effect AB
}

transformed parameters {
  matrix[2,2] muS;       // condition means
  for (i in 1:2) {
    for (j in 1:2) {
      muS[i,j] = mu + pow(-1,i+1) * tAfix / sqrt(2) 
                    + pow(-1,j+1) * tBfix / sqrt(2)
                    + pow(-1,i+j) * tABfix / 2;
    }
  }
}

model {
  target += uniform_lpdf(mu | lwr, upr);
  target += scaled_inv_chi_square_lpdf(gA | 1, hA);
  target += scaled_inv_chi_square_lpdf(gB | 1, hB);
  target += scaled_inv_chi_square_lpdf(gAB | 1, hAB);
  target += inv_wishart_lpdf(Sigma | 5, I);
  target += normal_lpdf(tAfix | 0, sqrt(gA));
  target += normal_lpdf(tBfix | 0, sqrt(gB));
  target += normal_lpdf(tABfix | 0, sqrt(gAB));
  for (k in 1:nS) {
    target += multi_normal_lpdf(Y[k] | to_vector(muS), Sigma);
  }
}
"
stanmodel.full.fixed <- stan_model(model_code=model.full.fixed)
datalist.fixed <- list("nS"=nrow(dat), "I"=diag(4), "Y"=dat,
                       "hA"=.5, "hB"=.5, "hAB"=.5,
                       "lwr"=mean(dat)-6*sd(dat),
                       "upr"=mean(dat)+6*sd(dat))

stanfit.full.fixed <- sampling(stanmodel.full.fixed, data=datalist.fixed,
                               iter=10000, warmup=1000, chains=2, seed=277, refresh=0)

set.seed(277)
M.full.fixed <- bridge_sampler(stanfit.full.fixed, silent=T)

Proposing a New (Omitted Fixed AB) Model

model.omitAB.fixed <- "
data {
  int<lower=1> nS;        // number of subjects
  // vector[4] Y[nS];        // responses [removed features]
  array[nS] vector[4] Y;  // responses [Stan 2.26+ syntax for array declarations]
  real<lower=0> hA;       // (square root) scale parameter of the variance of the projected treatment effect A
  real<lower=0> hB;       // (square root) scale parameter of the variance of the projected treatment effect B
  matrix[4,4] I;          // scale matrix of the inverse-Wishart prior
  real lwr;               // lower bound for the grand mean
  real upr;               // upper bound for the grand mean
}

parameters {
  real<lower=lwr, upper=upr> mu;     // grand mean
  cov_matrix[4] Sigma;               // error covariance matrix
  real<lower=0> gA;                  // variance of the projected treatment effect A
  real<lower=0> gB;                  // variance of the projected treatment effect B
  real tAfix;                        // projected treatment effect A
  real tBfix;                        // projected treatment effect B
}

transformed parameters {
  matrix[2,2] muS;       // condition means
  for (i in 1:2) {
    for (j in 1:2) {
      muS[i,j] = mu + pow(-1,i+1) * tAfix / sqrt(2) 
                    + pow(-1,j+1) * tBfix / sqrt(2);
    }
  }
}

model {
  target += uniform_lpdf(mu | lwr, upr);
  target += scaled_inv_chi_square_lpdf(gA | 1, hA);
  target += scaled_inv_chi_square_lpdf(gB | 1, hB);
  target += inv_wishart_lpdf(Sigma | 5, I);
  target += normal_lpdf(tAfix | 0, sqrt(gA));
  target += normal_lpdf(tBfix | 0, sqrt(gB));
  for (k in 1:nS) {
    target += multi_normal_lpdf(Y[k] | to_vector(muS), Sigma);
  }
}
"

stanmodel.omitAB.fixed <- stan_model(model_code=model.omitAB.fixed)
stanfit.omitAB.fixed <- sampling(stanmodel.omitAB.fixed, data=datalist.fixed[-6],
                                 iter=10000, warmup=1000, chains=2, seed=277, refresh=0)
set.seed(277)
M.omitAB.fixed <- bridge_sampler(stanfit.omitAB.fixed, silent=T)
bf(M.omitAB.fixed, M.full.fixed)

## Estimated Bayes factor in favor of M.omitAB.fixed over M.full.fixed: 1.39139

Proposing a New (Omitted Fixed A) Model

New Method 2, which violates the principle of marginality.

model.omitA.fixed <- "
data {
  int<lower=1> nS;        // number of subjects
  // vector[4] Y[nS];        // responses [removed features]
  array[nS] vector[4] Y;  // responses [Stan 2.26+ syntax for array declarations]
  real<lower=0> hB;       // (square root) scale parameter of the variance of the projected treatment effect B
  real<lower=0> hAB;      // (square root) scale parameter of the variance of the projected interaction effect AB
  matrix[4,4] I;          // scale matrix of the inverse-Wishart prior
  real lwr;               // lower bound for the grand mean
  real upr;               // upper bound for the grand mean
}

parameters {
  real<lower=lwr, upper=upr> mu;     // grand mean
  cov_matrix[4] Sigma;               // error covariance matrix
  real<lower=0> gB;                  // variance of the projected treatment effect B
  real<lower=0> gAB;                 // variance of the projected interaction effect AB
  real tBfix;                        // projected treatment effect B
  real tABfix;                       // projected interaction effect AB
}

transformed parameters {
  matrix[2,2] muS;       // condition means
  for (i in 1:2) {
    for (j in 1:2) {
      muS[i,j] = mu + pow(-1,j+1) * tBfix / sqrt(2)
                    + pow(-1,i+j) * tABfix / 2;
    }
  }
}

model {
  target += uniform_lpdf(mu | lwr, upr);
  target += scaled_inv_chi_square_lpdf(gB | 1, hB);
  target += scaled_inv_chi_square_lpdf(gAB | 1, hAB);
  target += inv_wishart_lpdf(Sigma | 5, I);
  target += normal_lpdf(tBfix | 0, sqrt(gB));
  target += normal_lpdf(tABfix | 0, sqrt(gAB));
  for (k in 1:nS) {
    target += multi_normal_lpdf(Y[k] | to_vector(muS), Sigma);
  }
}
"

stanmodel.omitA.fixed <- stan_model(model_code=model.omitA.fixed)
stanfit.omitA.fixed <- sampling(stanmodel.omitA.fixed, data=datalist.fixed[-4],
                                iter=10000, warmup=1000, chains=2, seed=277, refresh=0)
set.seed(277)
M.omitA.fixed <- bridge_sampler(stanfit.omitA.fixed, silent=T)
bf(M.omitA.fixed, M.full.fixed)

## Estimated Bayes factor in favor of M.omitA.fixed over M.full.fixed: 0.00000

New Method 4, which follows the principle of marginality.

model.notA.fixed <- "
data {
  int<lower=1> nS;        // number of subjects
  // vector[4] Y[nS];        // responses [removed features]
  array[nS] vector[4] Y;  // responses [Stan 2.26+ syntax for array declarations]
  real<lower=0> hB;       // (square root) scale parameter of the variance of the projected treatment effect B
  matrix[4,4] I;          // scale matrix of the inverse-Wishart prior
  real lwr;               // lower bound for the grand mean
  real upr;               // upper bound for the grand mean
}

parameters {
  real<lower=lwr, upper=upr> mu;     // grand mean
  cov_matrix[4] Sigma;               // error covariance matrix
  real<lower=0> gB;                  // variance of the projected treatment effect B
  real tBfix;                        // projected treatment effect B
}

transformed parameters {
  matrix[2,2] muS;       // condition means
  for (i in 1:2) {
    for (j in 1:2) {
      muS[i,j] = mu + pow(-1,j+1) * tBfix / sqrt(2);
    }
  }
}

model {
  target += uniform_lpdf(mu | lwr, upr);
  target += scaled_inv_chi_square_lpdf(gB | 1, hB);
  target += inv_wishart_lpdf(Sigma | 5, I);
  target += normal_lpdf(tBfix | 0, sqrt(gB));
  for (k in 1:nS) {
    target += multi_normal_lpdf(Y[k] | to_vector(muS), Sigma);
  }
}
"

stanmodel.notA.fixed <- stan_model(model_code=model.notA.fixed)
stanfit.notA.fixed <- sampling(stanmodel.notA.fixed, data=datalist.fixed[-c(4,6)],
                               iter=10000, warmup=1000, chains=2, seed=277, refresh=0)
set.seed(277)
M.notA.fixed <- bridge_sampler(stanfit.notA.fixed, silent=T)
bf(M.notA.fixed, M.omitAB.fixed)

## Estimated Bayes factor in favor of M.notA.fixed over M.omitAB.fixed: 0.00000

Proposing a New (Omitted Fixed B) Model

New Method 2, which violates the principle of marginality.

datalist.fixed[["Y"]] <- dat[,c(1,3,2,4)]
stanfit.omitB.fixed <- sampling(stanmodel.omitA.fixed, data=datalist.fixed[-4],
                                iter=10000, warmup=1000, chains=2, seed=277, refresh=0)
set.seed(277)
M.omitB.fixed <- bridge_sampler(stanfit.omitB.fixed, silent=T)
bf(M.omitB.fixed, M.full.fixed)

## Estimated Bayes factor in favor of M.omitB.fixed over M.full.fixed: 0.39077

New Method 4, which follows the principle of marginality.

stanfit.notB.fixed <- sampling(stanmodel.notA.fixed, data=datalist.fixed[-c(4,6)],
                               iter=10000, warmup=1000, chains=2, seed=277, refresh=0)
set.seed(277)
M.notB.fixed <- bridge_sampler(stanfit.notB.fixed, silent=T)
bf(M.notB.fixed, M.omitAB.fixed)

## Estimated Bayes factor in favor of M.notB.fixed over M.omitAB.fixed: 0.44566

summary(M.full.fixed)

## 
## Bridge sampling log marginal likelihood estimate 
## (method = "normal", repetitions = 1):
## 
##  -716.2848
## 
## Error Measures:
## 
##  Relative Mean-Squared Error: 3.847178e-05
##  Coefficient of Variation: 0.006202562
##  Percentage Error: 1%
## 
## Note:
## All error measures are approximate.

summary(M.omitAB.fixed)

## 
## Bridge sampling log marginal likelihood estimate 
## (method = "normal", repetitions = 1):
## 
##  -715.9545
## 
## Error Measures:
## 
##  Relative Mean-Squared Error: 2.816446e-05
##  Coefficient of Variation: 0.00530702
##  Percentage Error: 1%
## 
## Note:
## All error measures are approximate.

summary(M.omitA.fixed)

## 
## Bridge sampling log marginal likelihood estimate 
## (method = "normal", repetitions = 1):
## 
##  -732.3235
## 
## Error Measures:
## 
##  Relative Mean-Squared Error: 4.101512e-05
##  Coefficient of Variation: 0.006404305
##  Percentage Error: 1%
## 
## Note:
## All error measures are approximate.

summary(M.notA.fixed)

## 
## Bridge sampling log marginal likelihood estimate 
## (method = "normal", repetitions = 1):
## 
##  -731.7643
## 
## Error Measures:
## 
##  Relative Mean-Squared Error: 2.734897e-05
##  Coefficient of Variation: 0.005229624
##  Percentage Error: 1%
## 
## Note:
## All error measures are approximate.

summary(M.omitB.fixed)

## 
## Bridge sampling log marginal likelihood estimate 
## (method = "normal", repetitions = 1):
## 
##  -717.2244
## 
## Error Measures:
## 
##  Relative Mean-Squared Error: 2.965821e-05
##  Coefficient of Variation: 0.005445935
##  Percentage Error: 1%
## 
## Note:
## All error measures are approximate.

summary(M.notB.fixed)

## 
## Bridge sampling log marginal likelihood estimate 
## (method = "normal", repetitions = 1):
## 
##  -716.7626
## 
## Error Measures:
## 
##  Relative Mean-Squared Error: 1.960913e-05
##  Coefficient of Variation: 0.00442822
##  Percentage Error: 0%
## 
## Note:
## All error measures are approximate.

rstan::summary(stanfit.full.fixed)[[1]] %>% kable(digits=4) %>% kable_classic(full_width=F)

	mean	se_mean	sd	2.5%	25%	50%	75%	97.5%	n_eff	Rhat
mu	108.4500	0.0106	1.3753	105.7504	107.5202	108.4503	109.3696	111.1492	16703.699	0.9999
Sigma[1,1]	149.6417	0.3424	30.9193	100.6656	127.8076	145.6060	166.9457	222.0967	8153.252	0.9999
Sigma[1,2]	55.7927	0.2291	20.8781	19.6157	41.4648	53.9795	67.9252	102.3033	8304.950	1.0000
Sigma[1,3]	142.9304	0.3343	30.1637	95.4873	121.5778	138.8859	159.7269	213.3493	8140.802	0.9999
Sigma[1,4]	63.9465	0.2457	22.1101	26.6414	48.7323	61.9316	76.8350	113.4527	8100.989	1.0000
Sigma[2,1]	55.7927	0.2291	20.8781	19.6157	41.4648	53.9795	67.9252	102.3033	8304.950	1.0000
Sigma[2,2]	125.6237	0.2648	25.8152	84.7811	107.2380	122.4140	140.3503	185.6564	9501.961	1.0001
Sigma[2,3]	47.4106	0.2227	20.3659	11.4477	33.5148	45.7831	59.5226	92.0157	8359.522	1.0000
Sigma[2,4]	123.7917	0.2718	26.1422	82.1853	105.2211	120.5960	138.7046	184.6552	9253.267	1.0001
Sigma[3,1]	142.9304	0.3343	30.1637	95.4873	121.5778	138.8859	159.7269	213.3493	8140.802	0.9999
Sigma[3,2]	47.4106	0.2227	20.3659	11.4477	33.5148	45.7831	59.5226	92.0157	8359.522	1.0000
Sigma[3,3]	148.0119	0.3333	30.6336	99.7337	126.4437	143.9567	165.1908	219.6003	8446.755	0.9999
Sigma[3,4]	54.3358	0.2386	21.4609	17.3155	39.6560	52.4491	66.9400	101.9434	8089.407	1.0000
Sigma[4,1]	63.9465	0.2457	22.1101	26.6414	48.7323	61.9316	76.8350	113.4527	8100.989	1.0000
Sigma[4,2]	123.7917	0.2718	26.1422	82.1853	105.2211	120.5960	138.7046	184.6552	9253.267	1.0001
Sigma[4,3]	54.3358	0.2386	21.4609	17.3155	39.6560	52.4491	66.9400	101.9434	8089.407	1.0000
Sigma[4,4]	136.1448	0.2884	27.9422	91.5868	116.3961	132.5957	152.2300	201.2187	9386.177	1.0001
gA	479.9789	114.9119	11933.0217	10.9713	32.1386	65.0760	159.2141	1852.4952	10783.780	1.0002
gB	2.8725	1.0642	106.4528	0.0525	0.1619	0.3379	0.8418	9.4987	10006.079	1.0001
gAB	3.7879	2.2380	218.9112	0.0436	0.1293	0.2825	0.7315	8.6365	9568.294	1.0001
tAfix	-9.6045	0.0103	1.3002	-12.1554	-10.4825	-9.6071	-8.7389	-7.0542	15803.542	1.0000
tBfix	-0.4534	0.0018	0.2344	-0.9217	-0.6099	-0.4473	-0.2942	-0.0060	16255.118	0.9999
tABfix	0.2984	0.0025	0.3161	-0.2883	0.0841	0.2813	0.5006	0.9610	15447.910	0.9999
muS[1,1]	101.4872	0.0137	1.7330	98.1026	100.3066	101.5033	102.6414	104.8681	15907.654	0.9999
muS[1,2]	101.8301	0.0139	1.7367	98.3817	100.6722	101.8443	102.9856	105.2263	15689.540	0.9999
muS[2,1]	114.7716	0.0122	1.5845	111.6871	113.7128	114.7637	115.8285	117.8923	16734.723	1.0000
muS[2,2]	115.7113	0.0126	1.6207	112.5536	114.6304	115.7161	116.7912	118.8942	16599.431	1.0000
lp__	-722.3186	0.0373	3.1092	-729.3018	-724.1950	-722.0265	-720.0570	-717.2286	6930.221	1.0002

rstan::summary(stanfit.omitAB.fixed)[[1]] %>% kable(digits=4) %>% kable_classic(full_width=F)

	mean	se_mean	sd	2.5%	25%	50%	75%	97.5%	n_eff	Rhat
mu	108.2034	0.0101	1.3680	105.5225	107.2872	108.1987	109.1194	110.9136	18483.768	0.9999
Sigma[1,1]	149.3278	0.3397	30.6205	101.0113	127.7770	145.4553	166.4474	219.9658	8124.245	1.0000
Sigma[1,2]	55.8521	0.2316	21.1796	19.5212	41.3489	54.0481	68.3045	102.6094	8360.823	0.9999
Sigma[1,3]	142.7174	0.3299	29.8249	95.8596	121.7705	138.9482	159.4438	210.7569	8171.997	1.0000
Sigma[1,4]	64.0664	0.2477	22.4336	25.7776	48.6127	61.9570	77.1627	113.8462	8204.833	1.0000
Sigma[2,1]	55.8521	0.2316	21.1796	19.5212	41.3489	54.0481	68.3045	102.6094	8360.823	0.9999
Sigma[2,2]	125.9291	0.2676	25.9227	85.0598	107.5279	122.7561	140.6685	186.6609	9381.855	1.0000
Sigma[2,3]	47.3826	0.2242	20.6246	10.9912	33.4308	45.8687	59.5840	92.3077	8461.474	1.0000
Sigma[2,4]	124.1650	0.2748	26.3453	82.5764	105.2804	121.0102	139.2252	185.2692	9189.126	1.0000
Sigma[3,1]	142.7174	0.3299	29.8249	95.8596	121.7705	138.9482	159.4438	210.7569	8171.997	1.0000
Sigma[3,2]	47.3826	0.2242	20.6246	10.9912	33.4308	45.8687	59.5840	92.3077	8461.474	1.0000
Sigma[3,3]	147.8517	0.3274	30.2468	100.1754	126.4719	144.0688	164.8239	216.7622	8536.839	1.0000
Sigma[3,4]	54.3021	0.2394	21.7824	16.6852	39.4988	52.5718	67.2718	102.4822	8276.600	1.0000
Sigma[4,1]	64.0664	0.2477	22.4336	25.7776	48.6127	61.9570	77.1627	113.8462	8204.833	1.0000
Sigma[4,2]	124.1650	0.2748	26.3453	82.5764	105.2804	121.0102	139.2252	185.2692	9189.126	1.0000
Sigma[4,3]	54.3021	0.2394	21.7824	16.6852	39.4988	52.5718	67.2718	102.4822	8276.600	1.0000
Sigma[4,4]	136.7981	0.2916	28.2975	91.9256	116.6263	133.2146	153.0314	202.0326	9418.830	1.0000
gA	380.5931	35.8196	3769.7393	10.9690	31.7685	64.9912	157.1937	1897.2298	11075.973	0.9999
gB	2.3427	0.4118	37.0579	0.0512	0.1582	0.3290	0.8225	9.6028	8097.761	1.0002
tAfix	-9.4976	0.0097	1.3015	-12.0464	-10.3598	-9.4989	-8.6390	-6.9205	17991.038	1.0000
tBfix	-0.4368	0.0018	0.2369	-0.9171	-0.5939	-0.4293	-0.2739	0.0084	16801.477	1.0000
muS[1,1]	101.1787	0.0130	1.7245	97.8377	100.0150	101.1642	102.3014	104.6383	17609.723	0.9999
muS[1,2]	101.7965	0.0135	1.7559	98.4007	100.6295	101.7964	102.9419	105.3162	16888.342	0.9999
muS[2,1]	114.6104	0.0111	1.5832	111.4707	113.5401	114.6101	115.6864	117.7037	20193.701	0.9999
muS[2,2]	115.2281	0.0108	1.5563	112.1093	114.1746	115.2394	116.2899	118.2657	20787.529	0.9999
lp__	-720.1962	0.0365	2.9105	-726.7733	-721.9048	-719.8441	-718.0792	-715.5474	6374.630	1.0000

rstan::summary(stanfit.omitA.fixed)[[1]] %>% kable(digits=4) %>% kable_classic(full_width=F)

	mean	se_mean	sd	2.5%	25%	50%	75%	97.5%	n_eff	Rhat
mu	109.2812	0.0202	2.0177	105.2964	107.9333	109.2772	110.6292	113.2303	9981.438	1.0002
Sigma[1,1]	206.8205	0.5602	53.5877	127.0696	168.7277	198.7171	235.6029	335.0461	9148.912	1.0002
Sigma[1,2]	13.2893	0.2454	27.1565	-37.9389	-4.5028	12.5418	30.1255	69.3432	12249.144	1.0000
Sigma[1,3]	203.8634	0.5640	53.7166	123.0849	165.3771	195.9240	233.1711	332.0930	9070.951	1.0002
Sigma[1,4]	19.8837	0.2561	28.3691	-33.3651	1.1684	18.8418	37.2557	78.7190	12268.354	1.0000
Sigma[2,1]	13.2893	0.2454	27.1565	-37.9389	-4.5028	12.5418	30.1255	69.3432	12249.144	1.0000
Sigma[2,2]	165.4435	0.4428	42.5919	102.0329	135.2259	158.9126	188.6020	266.5112	9252.257	1.0006
Sigma[2,3]	1.5115	0.2484	27.7195	-52.7484	-16.2112	1.1915	18.9076	57.9360	12449.528	1.0000
Sigma[2,4]	164.9472	0.4583	43.5195	100.3895	134.1947	158.2628	188.0103	267.9925	9015.616	1.0006
Sigma[3,1]	203.8634	0.5640	53.7166	123.0849	165.3771	195.9240	233.1711	332.0930	9070.951	1.0002
Sigma[3,2]	1.5115	0.2484	27.7195	-52.7484	-16.2112	1.1915	18.9076	57.9360	12449.528	1.0000
Sigma[3,3]	213.1574	0.5746	55.2226	129.6416	173.2614	205.4179	243.7432	344.4928	9236.704	1.0002
Sigma[3,4]	6.8476	0.2578	28.7027	-48.8217	-11.4704	6.2071	24.8465	65.3632	12394.579	0.9999
Sigma[4,1]	19.8837	0.2561	28.3691	-33.3651	1.1684	18.8418	37.2557	78.7190	12268.354	1.0000
Sigma[4,2]	164.9472	0.4583	43.5195	100.3895	134.1947	158.2628	188.0103	267.9925	9015.616	1.0006
Sigma[4,3]	6.8476	0.2578	28.7027	-48.8217	-11.4704	6.2071	24.8465	65.3632	12394.579	0.9999
Sigma[4,4]	178.7501	0.4833	46.0559	110.8607	146.2923	171.8026	202.9378	288.7311	9079.446	1.0005
gB	10.5917	7.6382	1017.2566	0.0643	0.2268	0.5108	1.3386	14.7055	17736.893	1.0000
gAB	1.8966	0.1868	18.3836	0.0420	0.1252	0.2672	0.6694	8.8325	9683.103	1.0001
tBfix	-0.6666	0.0027	0.3429	-1.3706	-0.8947	-0.6581	-0.4307	-0.0253	15879.464	0.9999
tABfix	0.1012	0.0034	0.3914	-0.6676	-0.1455	0.0876	0.3350	0.9198	13188.439	1.0000
muS[1,1]	108.8604	0.0208	2.0711	104.7561	107.4771	108.8426	110.2579	112.9635	9931.071	1.0003
muS[1,2]	109.7020	0.0199	2.0082	105.7324	108.3702	109.7076	111.0463	113.5909	10164.060	1.0002
muS[2,1]	108.7592	0.0199	2.0103	104.8021	107.4231	108.7512	110.1047	112.6969	10203.812	1.0002
muS[2,2]	109.8032	0.0208	2.0760	105.6973	108.4167	109.8012	111.2072	113.8527	9955.099	1.0002
lp__	-736.4026	0.0358	2.9328	-743.0624	-738.1936	-736.0600	-734.2516	-731.6994	6709.660	1.0002

rstan::summary(stanfit.notA.fixed)[[1]] %>% kable(digits=4) %>% kable_classic(full_width=F)

	mean	se_mean	sd	2.5%	25%	50%	75%	97.5%	n_eff	Rhat
mu	109.1869	0.0211	1.9947	105.2541	107.8860	109.2010	110.4879	113.1266	8954.234	0.9999
Sigma[1,1]	205.0389	0.5896	52.7557	125.7556	167.7452	196.9731	233.2771	330.1519	8006.596	0.9999
Sigma[1,2]	14.1051	0.2735	27.1296	-38.0364	-3.4987	13.3848	30.4902	70.5311	9837.846	1.0000
Sigma[1,3]	203.0054	0.5987	53.3643	122.7407	165.1427	195.0499	231.7477	329.0090	7945.485	0.9999
Sigma[1,4]	19.6064	0.2859	28.5254	-34.8860	1.0022	18.8995	36.9981	79.3413	9955.801	1.0001
Sigma[2,1]	14.1051	0.2735	27.1296	-38.0364	-3.4987	13.3848	30.4902	70.5311	9837.846	1.0000
Sigma[2,2]	166.2117	0.4738	43.3740	101.3953	135.4195	159.4702	189.4047	269.4888	8381.524	1.0000
Sigma[2,3]	1.7179	0.2768	27.8622	-52.8962	-15.8849	1.6493	18.9921	57.4098	10134.479	1.0000
Sigma[2,4]	166.3413	0.4888	44.2438	100.5908	134.7890	159.3242	190.3828	271.6590	8192.500	1.0000
Sigma[3,1]	203.0054	0.5987	53.3643	122.7407	165.1427	195.0499	231.7477	329.0090	7945.485	0.9999
Sigma[3,2]	1.7179	0.2768	27.8622	-52.8962	-15.8849	1.6493	18.9921	57.4098	10134.479	1.0000
Sigma[3,3]	213.1800	0.6136	55.2651	130.0646	173.9752	205.1101	243.3013	342.3510	8111.239	0.9999
Sigma[3,4]	6.0261	0.2849	28.7641	-50.3868	-12.2555	5.7162	23.7424	64.2874	10195.732	1.0000
Sigma[4,1]	19.6064	0.2859	28.5254	-34.8860	1.0022	18.8995	36.9981	79.3413	9955.801	1.0001
Sigma[4,2]	166.3413	0.4888	44.2438	100.5908	134.7890	159.3242	190.3828	271.6590	8192.500	1.0000
Sigma[4,3]	6.0261	0.2849	28.7641	-50.3868	-12.2555	5.7162	23.7424	64.2874	10195.732	1.0000
Sigma[4,4]	180.7359	0.5123	46.4958	111.3787	147.6344	173.6313	206.3348	290.1020	8237.290	1.0000
gB	3.7604	0.4610	47.8701	0.0627	0.2268	0.5094	1.3533	16.6678	10783.500	1.0000
tBfix	-0.6627	0.0029	0.3411	-1.3526	-0.8929	-0.6581	-0.4233	-0.0227	14117.614	1.0000
muS[1,1]	108.7184	0.0213	2.0102	104.8064	107.3927	108.7266	110.0532	112.6832	8942.503	0.9999
muS[1,2]	109.6555	0.0211	2.0083	105.6942	108.3445	109.6616	110.9802	113.6323	9078.139	1.0000
muS[2,1]	108.7184	0.0213	2.0102	104.8064	107.3927	108.7266	110.0532	112.6832	8942.503	0.9999
muS[2,2]	109.6555	0.0211	2.0083	105.6942	108.3445	109.6616	110.9802	113.6323	9078.139	1.0000
lp__	-733.8407	0.0320	2.6883	-739.9221	-735.4314	-733.5051	-731.8705	-729.6045	7045.181	1.0002

rstan::summary(stanfit.omitB.fixed)[[1]] %>% kable(digits=4) %>% kable_classic(full_width=F)

	mean	se_mean	sd	2.5%	25%	50%	75%	97.5%	n_eff	Rhat
mu	108.3936	0.0104	1.4220	105.5942	107.4398	108.3922	109.3502	111.1632	18540.992	0.9999
Sigma[1,1]	149.9667	0.3166	31.1374	101.1537	128.1426	145.8097	167.6594	221.2492	9670.116	1.0008
Sigma[1,2]	143.2740	0.3079	30.3621	95.2452	121.7684	139.2396	160.6858	214.2947	9724.946	1.0009
Sigma[1,3]	56.4575	0.2199	21.7401	19.4710	41.6477	54.3653	68.7734	105.2259	9777.997	1.0003
Sigma[1,4]	64.6337	0.2343	23.0621	26.3942	48.8567	62.3735	77.6042	117.0655	9685.418	1.0003
Sigma[2,1]	143.2740	0.3079	30.3621	95.2452	121.7684	139.2396	160.6858	214.2947	9724.946	1.0009
Sigma[2,2]	148.5536	0.3098	30.8295	99.9957	126.9253	144.4918	165.8372	220.5612	9905.949	1.0008
Sigma[2,3]	47.8269	0.2155	21.2229	11.1803	33.5028	45.9945	60.0168	94.6999	9701.955	1.0002
Sigma[2,4]	55.2347	0.2284	22.4648	16.4483	39.9060	53.1881	68.1937	105.1951	9675.311	1.0003
Sigma[3,1]	56.4575	0.2199	21.7401	19.4710	41.6477	54.3653	68.7734	105.2259	9777.997	1.0003
Sigma[3,2]	47.8269	0.2155	21.2229	11.1803	33.5028	45.9945	60.0168	94.6999	9701.955	1.0002
Sigma[3,3]	126.0460	0.2667	26.1627	84.9244	107.6596	122.3031	140.9861	186.8371	9621.762	1.0001
Sigma[3,4]	123.9896	0.2738	26.5074	82.2950	105.2946	120.4058	138.8651	185.9100	9374.271	1.0002
Sigma[4,1]	64.6337	0.2343	23.0621	26.3942	48.8567	62.3735	77.6042	117.0655	9685.418	1.0003
Sigma[4,2]	55.2347	0.2284	22.4648	16.4483	39.9060	53.1881	68.1937	105.1951	9675.311	1.0003
Sigma[4,3]	123.9896	0.2738	26.5074	82.2950	105.2946	120.4058	138.8651	185.9100	9374.271	1.0002
Sigma[4,4]	136.8840	0.2910	28.5357	92.4470	117.0291	133.0129	152.6571	203.3381	9614.644	1.0002
gB	572.7021	91.1082	10114.5705	12.2132	35.2657	71.8266	173.1695	1928.4532	12324.815	1.0001
gAB	2.6284	0.9737	92.1880	0.0432	0.1280	0.2682	0.6730	7.8169	8964.073	1.0002
tBfix	-9.9995	0.0096	1.3261	-12.5786	-10.8867	-10.0132	-9.1279	-7.3529	18887.626	0.9999
tABfix	0.2459	0.0026	0.3222	-0.3571	0.0259	0.2344	0.4479	0.9312	15514.600	1.0001
muS[1,1]	101.4458	0.0131	1.7831	97.9364	100.2501	101.4375	102.6271	104.9933	18488.980	0.9999
muS[1,2]	115.3415	0.0116	1.6182	112.1943	114.2630	115.3457	116.4392	118.5122	19414.203	0.9999
muS[2,1]	101.2000	0.0127	1.7549	97.7422	100.0182	101.1977	102.3586	104.6807	18959.774	0.9999
muS[2,2]	115.5873	0.0124	1.6829	112.2914	114.4628	115.6001	116.7083	118.8627	18428.226	0.9999
lp__	-723.0302	0.0360	2.8852	-729.6975	-724.7125	-722.6927	-720.9455	-718.4576	6440.065	1.0002

rstan::summary(stanfit.notB.fixed)[[1]] %>% kable(digits=4) %>% kable_classic(full_width=F)

	mean	se_mean	sd	2.5%	25%	50%	75%	97.5%	n_eff	Rhat
mu	108.1757	0.0104	1.4234	105.3644	107.2349	108.1692	109.1261	111.0067	18903.788	0.9999
Sigma[1,1]	150.2811	0.2984	31.2444	101.0405	127.9688	146.3854	167.7333	222.9841	10965.463	1.0000
Sigma[1,2]	143.7526	0.2904	30.4787	95.5531	121.9417	139.7821	160.9329	214.2889	11013.879	1.0000
Sigma[1,3]	55.9999	0.2185	21.6316	18.7928	41.1729	54.0416	68.9414	104.4579	9799.572	1.0002
Sigma[1,4]	64.4692	0.2345	22.9206	25.6599	48.6183	62.0877	78.0937	115.6079	9550.569	1.0001
Sigma[2,1]	143.7526	0.2904	30.4787	95.5531	121.9417	139.7821	160.9329	214.2889	11013.879	1.0000
Sigma[2,2]	148.8621	0.2914	30.9285	99.9288	126.8695	144.8913	166.4787	220.6707	11264.042	1.0000
Sigma[2,3]	47.6817	0.2112	21.0780	10.4789	33.5161	45.9191	60.2244	94.0055	9956.480	1.0004
Sigma[2,4]	55.3055	0.2247	22.3070	17.0219	40.0397	53.2866	68.5307	104.5872	9852.530	1.0002
Sigma[3,1]	55.9999	0.2185	21.6316	18.7928	41.1729	54.0416	68.9414	104.4579	9799.572	1.0002
Sigma[3,2]	47.6817	0.2112	21.0780	10.4789	33.5161	45.9191	60.2244	94.0055	9956.480	1.0004
Sigma[3,3]	126.0777	0.2692	25.9874	84.7512	107.4989	122.7973	141.0239	185.8633	9316.220	1.0002
Sigma[3,4]	123.9887	0.2697	26.3362	81.7652	105.3325	120.6364	139.0998	184.0340	9535.556	1.0002
Sigma[4,1]	64.4692	0.2345	22.9206	25.6599	48.6183	62.0877	78.0937	115.6079	9550.569	1.0001
Sigma[4,2]	55.3055	0.2247	22.3070	17.0219	40.0397	53.2866	68.5307	104.5872	9852.530	1.0002
Sigma[4,3]	123.9887	0.2697	26.3362	81.7652	105.3325	120.6364	139.0998	184.0340	9535.556	1.0002
Sigma[4,4]	137.3578	0.2852	28.4536	92.1258	117.3185	133.6024	153.4221	202.5832	9955.870	1.0002
gB	400.1721	42.2534	3922.1923	12.3564	33.9245	68.2649	164.8698	1850.4451	8616.552	1.0000
tBfix	-9.9307	0.0098	1.3211	-12.4802	-10.8132	-9.9394	-9.0506	-7.3324	18167.201	0.9999
muS[1,1]	101.1536	0.0134	1.7778	97.7042	99.9546	101.1533	102.3406	104.6391	17638.122	0.9999
muS[1,2]	115.1977	0.0114	1.6239	111.9925	114.1161	115.1860	116.2745	118.4135	20224.673	0.9999
muS[2,1]	101.1536	0.0134	1.7778	97.7042	99.9546	101.1533	102.3406	104.6391	17638.122	0.9999
muS[2,2]	115.1977	0.0114	1.6239	111.9925	114.1161	115.1860	116.2745	118.4135	20224.673	0.9999
lp__	-720.7700	0.0321	2.7003	-726.9239	-722.3951	-720.4247	-718.8059	-716.5229	7068.818	1.0000

6. Summary

\[\begin{align*} \mathcal{M}_{\mathrm{full}}:&\ Y_{ijk}=\mu+s_k+\alpha_i+\beta_j+(\alpha\beta)_{ij}+\epsilon_{ijk}\tag{1}\\ \mathcal{M}_{\mathrm{full}}^{'}:&\ Y_{ijk}=\mu+s_k+\alpha_i+\beta_j+(\alpha\beta)_{ij}+(\alpha s)_{ik}+(\beta s)_{jk}+\epsilon_{ijk}\tag{2}\\ \mathcal{M}_{\mathrm{full}}^{*}:&\ Y_{ijk}=\mu+\alpha_i+\beta_j+(\alpha\beta)_{ij}+\epsilon_{ijk}^{*}\tag{3} \end{align*}\]

Classical: Repeated-measures ANOVA
JASP: The reduced model excludes each effect term and the interaction in (2). Version 0.16.3 or later seems to match Method 4.
Method 1: Compare the model without each effect to the full model in (1). BayesFactor version 0.9.12-4.6.
Method 2: Include the random slopes in (2) and test each effect again. The reduced model only omits each effect term in (2).
Method 3: The reduced model excludes each effect term, the interaction, and the random slope for such an effect in (2).
Method 4: The reduced model excludes each effect term and the interaction in (2).
New Method 2: The full model in (3) has no explicit terms for the subject-specific random effects. The reduced model only omits each effect term in (3).
New Method 4: The reduced model excludes each effect term and the interaction in (3).