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$) 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
dat <- dat[,c(1,3,2,4)]
cov(dat) %>% #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

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 : 3.625492  ±5.35%
## [2] Omit align      : 1.621491  ±11.94%
## [3] Omit resp       : 0.4673661 ±3.91%
## 
## 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 : 0.4931613 ±12.44%
## 
## 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 : 2.099691 ±10.16%
## 
## 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 : 0.0005168662 ±12.96%
## 
## 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.977033e-29 ±7.64%
## 
## 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 : 2.583492e-05 ±9.52%
## 
## 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 : 2.148667 ±7.78%
## 
## 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 : 0.0007444065 ±13.86%
## 
## 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 : 9.131032 ±7.44%
## 
## 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 : 2.162451 ±5.24%
## 
## 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.343311 ±5.24%
## 
## 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 : 0.3859289 ±27.02%
## 
## 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.6553196 ±27.54%
## 
## 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 : 0.7108546 ±28.13%
## 
## 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 : 2.413584e-30 ±8.14%
## 
## 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 : 1.207287e-05 ±19.72%
## 
## 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.885159 ±23.56%
## 
## 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 : 0.0006992793 ±19.88%
## 
## 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)

## Warning: There were 11 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.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	324.4336	0.2250	14.3097	296.1189	316.0138	324.4232	332.8015	352.5539	4044.493	1.0000
Sigma[1,1]	3946.8278	10.6326	1037.3893	2403.2435	3211.5861	3781.7298	4509.1369	6405.6970	9519.377	1.0000
Sigma[1,2]	-222.6865	8.2844	840.9754	-1993.0571	-729.3603	-201.3563	295.9437	1411.3856	10305.040	1.0000
Sigma[1,3]	4085.6770	11.5472	1110.8081	2415.3310	3296.7067	3903.7430	4689.9835	6726.0893	9253.950	1.0000
Sigma[1,4]	-840.0288	10.0698	1002.7214	-3021.3697	-1428.1134	-783.9164	-198.5604	1044.5858	9915.694	0.9999
Sigma[2,1]	-222.6865	8.2844	840.9754	-1993.0571	-729.3603	-201.3563	295.9437	1411.3856	10305.040	1.0000
Sigma[2,2]	5206.9772	15.4870	1391.7493	3168.9505	4239.4247	4975.6319	5935.3169	8502.5306	8075.814	1.0000
Sigma[2,3]	-594.6007	9.1481	936.5349	-2594.9550	-1134.4925	-557.1269	-5.2831	1194.7914	10480.616	0.9999
Sigma[2,4]	5912.4850	17.7936	1605.3896	3545.5241	4795.1209	5647.8286	6740.8921	9756.4524	8140.191	1.0000
Sigma[3,1]	4085.6770	11.5472	1110.8081	2415.3310	3296.7067	3903.7430	4689.9835	6726.0893	9253.950	1.0000
Sigma[3,2]	-594.6007	9.1481	936.5349	-2594.9550	-1134.4925	-557.1269	-5.2831	1194.7914	10480.616	0.9999
Sigma[3,3]	4813.2716	13.0110	1269.2891	2916.3805	3911.4068	4604.7833	5516.8594	7838.9757	9517.024	1.0001
Sigma[3,4]	-1001.1696	10.8654	1108.6521	-3424.2648	-1640.7774	-929.1503	-290.0968	1056.6708	10411.087	0.9999
Sigma[4,1]	-840.0288	10.0698	1002.7214	-3021.3697	-1428.1134	-783.9164	-198.5604	1044.5858	9915.694	0.9999
Sigma[4,2]	5912.4850	17.7936	1605.3896	3545.5241	4795.1209	5647.8286	6740.8921	9756.4524	8140.191	1.0000
Sigma[4,3]	-1001.1696	10.8654	1108.6521	-3424.2648	-1640.7774	-929.1503	-290.0968	1056.6708	10411.087	0.9999
Sigma[4,4]	7246.8207	20.9749	1926.6029	4394.9331	5897.9099	6944.7017	8227.3139	11936.8666	8436.897	1.0001
gA	144.3016	30.2998	2524.1383	0.2078	0.7804	2.3383	11.7161	857.2128	6939.800	1.0002
gB	180.6154	25.9918	2415.9451	0.2266	1.1003	5.2840	44.4633	965.9891	8639.773	1.0001
gAB	247.3520	8.7315	723.7045	4.2849	48.2767	112.4143	245.8722	1278.4384	6869.740	1.0002
tA[1]	1.5605	0.1448	7.2165	-6.1980	-0.6640	0.2287	1.5450	22.4883	2484.375	1.0000
tA[2]	-1.7528	0.2060	8.1597	-23.1126	-1.5257	-0.2388	0.6768	5.9790	1569.241	1.0037
tB[1]	-2.6857	0.1513	8.0835	-22.4851	-3.8676	-0.6803	0.4404	7.6140	2854.469	1.0007
tB[2]	2.8015	0.1539	8.2533	-7.3139	-0.4299	0.7137	4.0134	23.4975	2875.682	1.0013
tAB[1,1]	0.9877	0.1476	9.6963	-17.9040	-4.1672	0.9331	5.6706	21.7818	4313.977	1.0000
tAB[1,2]	6.9848	0.1986	11.1898	-10.6758	-0.5469	5.4323	12.9915	33.1877	3173.391	1.0002
tAB[2,1]	-13.7812	0.1870	10.8985	-37.8987	-19.6848	-12.6279	-6.3167	2.9254	3397.954	1.0002
tAB[2,2]	6.0708	0.1385	10.2182	-15.3328	0.7192	5.8746	11.6741	26.0957	5440.724	1.0003
muS[1,1]	324.2961	0.1126	9.9800	305.0559	317.7516	324.2146	330.6261	344.7442	7853.230	1.0003
muS[1,2]	335.7804	0.1288	11.0776	314.5667	328.3932	335.5343	342.7547	358.9055	7402.156	1.0003
muS[2,1]	306.2139	0.1260	10.8868	283.4572	299.2285	306.5225	313.5826	326.5285	7459.652	1.0006
muS[2,2]	331.5532	0.1604	12.3024	305.1035	323.9522	332.1501	339.7565	354.1450	5885.441	1.0011
lp__	-760.3060	0.0959	4.6595	-770.3823	-763.2645	-759.9502	-756.9680	-752.3053	2360.456	1.0001

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)

## Warning: There were 93 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.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	321.9630	1.6726	53.6650	225.7598	302.8358	321.5437	340.2696	422.7417	1029.4572	1.0011
Sigma[1,1]	4000.7115	11.4673	1085.3345	2396.7869	3237.6184	3832.4329	4577.8796	6600.6531	8957.8630	0.9999
Sigma[1,2]	-280.0281	9.3890	872.6778	-2085.2815	-814.0803	-258.6365	256.7391	1452.6284	8639.1006	1.0006
Sigma[1,3]	4143.7640	12.4455	1162.8231	2417.4968	3319.5285	3960.9761	4750.3836	6923.9211	8729.6997	0.9999
Sigma[1,4]	-948.6435	11.2882	1064.0859	-3256.5590	-1577.9808	-884.6982	-262.9230	1037.7605	8885.9789	1.0006
Sigma[2,1]	-280.0281	9.3890	872.6778	-2085.2815	-814.0803	-258.6365	256.7391	1452.6284	8639.1006	1.0006
Sigma[2,2]	5349.5039	23.1298	1446.9863	3212.3466	4316.1792	5127.7342	6110.7483	8925.2425	3913.6660	1.0013
Sigma[2,3]	-669.3105	10.1342	968.8258	-2714.5372	-1253.9105	-618.0998	-51.3574	1147.3051	9139.2618	1.0005
Sigma[2,4]	6134.9982	20.4161	1681.2635	3644.9607	4932.0854	5881.7791	7029.9486	10154.2004	6781.5050	1.0012
Sigma[3,1]	4143.7640	12.4455	1162.8231	2417.4968	3319.5285	3960.9761	4750.3836	6923.9211	8729.6997	0.9999
Sigma[3,2]	-669.3105	10.1342	968.8258	-2714.5372	-1253.9105	-618.0998	-51.3574	1147.3051	9139.2618	1.0005
Sigma[3,3]	4887.6882	14.1069	1328.4039	2927.3795	3953.7853	4674.0083	5588.2380	8081.3969	8867.3612	0.9999
Sigma[3,4]	-1146.7672	12.0434	1169.2489	-3633.1074	-1847.3735	-1068.0803	-388.3398	1000.0870	9425.8018	1.0005
Sigma[4,1]	-948.6435	11.2882	1064.0859	-3256.5590	-1577.9808	-884.6982	-262.9230	1037.7605	8885.9789	1.0006
Sigma[4,2]	6134.9982	20.4161	1681.2635	3644.9607	4932.0854	5881.7791	7029.9486	10154.2004	6781.5050	1.0012
Sigma[4,3]	-1146.7672	12.0434	1169.2489	-3633.1074	-1847.3735	-1068.0803	-388.3398	1000.0870	9425.8018	1.0005
Sigma[4,4]	7602.2011	22.2909	2041.9404	4561.2474	6141.7739	7291.7959	8698.0289	12574.0767	8391.3008	1.0011
gA	6605.9694	1038.7580	87423.5133	28.7091	386.1198	918.7790	2427.1598	33877.4720	7083.1708	1.0000
gB	709.2544	142.9976	6683.0631	13.2667	50.9950	113.5367	289.8556	3857.4990	2184.2071	1.0009
tA[1]	25.1221	1.4972	50.1424	-61.1269	7.4978	23.5067	41.4310	117.3004	1121.5786	1.0009
tA[2]	-25.7164	1.4836	49.9872	-117.9584	-40.9870	-23.1011	-7.7088	58.4591	1135.1861	1.0009
tB[1]	-9.1585	0.6532	19.4286	-39.4020	-14.0434	-8.3438	-3.2154	18.9295	884.5694	1.0005
tB[2]	8.5777	0.6496	19.4087	-20.0624	3.3103	8.4851	14.2279	38.6041	892.7809	1.0007
muS[1,1]	337.9267	0.1143	11.8721	313.3298	330.3761	338.2547	345.8478	360.6333	10792.8311	1.0000
muS[1,2]	355.6629	0.1122	11.6758	331.6938	348.1865	356.0001	363.5103	377.4546	10823.7072	0.9999
muS[2,1]	287.0881	0.1934	13.1895	262.1802	278.4226	286.6652	295.1889	315.2235	4648.5934	1.0000
muS[2,2]	304.8244	0.1959	14.5144	277.3161	295.3039	304.3776	313.7803	335.8390	5487.9367	1.0001
lp__	-757.1612	0.0678	3.5119	-765.1268	-759.2708	-756.7720	-754.6013	-751.4855	2685.1562	1.0004

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)

## Warning: There were 6 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.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	324.3824	0.3457	15.7890	297.3118	316.8803	324.6615	332.6153	351.3327	2086.506	1.0007
Sigma[1,1]	3931.9297	10.5315	1044.5435	2383.0003	3190.4074	3770.4494	4499.6892	6430.2096	9837.217	1.0000
Sigma[1,2]	-215.7685	8.9497	834.3321	-1955.9093	-709.0112	-188.0173	304.6338	1413.9677	8690.791	1.0001
Sigma[1,3]	4070.5122	11.4928	1119.9809	2405.6058	3273.6939	3894.3006	4677.4852	6731.1156	9496.606	1.0000
Sigma[1,4]	-823.3815	10.6427	991.4596	-2964.9001	-1396.0038	-766.7649	-188.4179	1023.1853	8678.511	1.0000
Sigma[2,1]	-215.7685	8.9497	834.3321	-1955.9093	-709.0112	-188.0173	304.6338	1413.9677	8690.791	1.0001
Sigma[2,2]	5204.2346	14.8287	1379.8657	3162.5851	4230.3657	4985.1472	5943.3103	8489.0473	8658.949	1.0002
Sigma[2,3]	-584.9544	9.8763	928.1437	-2554.3408	-1128.3835	-544.8610	1.7090	1155.0207	8831.709	1.0000
Sigma[2,4]	5905.8469	17.2255	1593.3142	3558.2179	4773.4771	5638.7416	6745.5810	9718.3511	8555.730	1.0001
Sigma[3,1]	4070.5122	11.4928	1119.9809	2405.6058	3273.6939	3894.3006	4677.4852	6731.1156	9496.606	1.0000
Sigma[3,2]	-584.9544	9.8763	928.1437	-2554.3408	-1128.3835	-544.8610	1.7090	1155.0207	8831.709	1.0000
Sigma[3,3]	4796.0229	13.0686	1281.1690	2905.5066	3887.5331	4595.7549	5475.5201	7875.1682	9610.648	1.0000
Sigma[3,4]	-981.3468	11.6824	1097.9542	-3361.1640	-1613.5916	-923.5801	-276.1226	1025.8343	8832.864	1.0000
Sigma[4,1]	-823.3815	10.6427	991.4596	-2964.9001	-1396.0038	-766.7649	-188.4179	1023.1853	8678.511	1.0000
Sigma[4,2]	5905.8469	17.2255	1593.3142	3558.2179	4773.4771	5638.7416	6745.5810	9718.3511	8555.730	1.0001
Sigma[4,3]	-981.3468	11.6824	1097.9542	-3361.1640	-1613.5916	-923.5801	-276.1226	1025.8343	8832.864	1.0000
Sigma[4,4]	7234.1774	20.6011	1916.5162	4419.1734	5876.2079	6913.7259	8250.7080	11785.1656	8654.547	1.0001
gB	310.4821	89.1443	5896.9195	0.2204	1.0050	4.1493	37.6889	1143.8878	4375.855	1.0002
gAB	262.8415	13.1569	708.8701	9.6402	55.4732	121.1477	259.2642	1354.7475	2902.875	1.0002
tB[1]	-2.1178	0.2908	10.9264	-22.5944	-3.3541	-0.5877	0.4718	7.9273	1412.024	1.0002
tB[2]	2.9617	0.3255	11.4199	-7.2327	-0.4550	0.6099	3.4467	24.2577	1230.987	1.0014
tAB[1,1]	1.0447	0.1815	10.1959	-17.7941	-4.3771	1.0517	6.0243	21.5870	3156.928	1.0007
tAB[1,2]	7.2956	0.2188	11.5841	-10.7548	-0.1137	5.9655	13.2514	32.6087	2803.494	1.0000
tAB[2,1]	-14.7200	0.2026	11.4976	-39.8509	-20.5703	-13.5571	-7.3529	2.5396	3220.297	1.0009
tAB[2,2]	6.0115	0.1744	11.1620	-16.4973	0.6255	5.9226	12.0426	26.4140	4094.789	1.0000
muS[1,1]	323.3093	0.0761	9.4958	304.9241	317.0740	323.1174	329.4808	342.4858	15588.982	0.9999
muS[1,2]	334.6396	0.0830	10.3682	314.4526	327.7635	334.4353	341.4127	355.4694	15598.720	0.9999
muS[2,1]	307.5446	0.0914	10.2221	286.5492	300.9844	307.7238	314.3985	326.9614	12499.493	1.0001
muS[2,2]	333.3555	0.1011	11.2351	310.6320	326.2506	333.5962	340.7858	354.6897	12355.054	1.0001
lp__	-754.3054	0.0709	3.9784	-763.1486	-756.7420	-753.9276	-751.4661	-747.6115	3148.721	1.0008

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)

## Warning: There were 2 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.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	324.5258	1.1326	26.1520	286.0593	314.8697	324.9999	335.4637	364.3128	533.1287	1.0047
Sigma[1,1]	4161.9843	11.1959	1145.9406	2489.6258	3350.4254	3970.3259	4759.5468	6938.6587	10476.2292	1.0000
Sigma[1,2]	-351.0303	9.1061	890.5995	-2166.4306	-878.7208	-335.5261	189.4716	1407.1568	9565.3720	1.0000
Sigma[1,3]	4175.1842	11.5692	1171.0947	2451.9674	3341.1222	3982.2273	4790.8181	6988.1380	10246.6227	0.9999
Sigma[1,4]	-895.1376	10.5683	1046.0132	-3099.4280	-1502.9065	-845.1212	-245.0661	1089.4009	9796.2594	0.9999
Sigma[2,1]	-351.0303	9.1061	890.5995	-2166.4306	-878.7208	-335.5261	189.4716	1407.1568	9565.3720	1.0000
Sigma[2,2]	5428.3697	15.7987	1483.2217	3254.3780	4386.2162	5195.3503	6200.8016	8952.4151	8813.8871	1.0000
Sigma[2,3]	-582.1034	9.6574	965.7232	-2587.4815	-1152.5479	-558.2529	20.9815	1280.1912	9999.6995	1.0001
Sigma[2,4]	6049.8894	17.7054	1659.0400	3600.2437	4882.6608	5792.3560	6932.5972	10010.5246	8780.1923	1.0001
Sigma[3,1]	4175.1842	11.5692	1171.0947	2451.9674	3341.1222	3982.2273	4790.8181	6988.1380	10246.6227	0.9999
Sigma[3,2]	-582.1034	9.6574	965.7232	-2587.4815	-1152.5479	-558.2529	20.9815	1280.1912	9999.6995	1.0001
Sigma[3,3]	4830.9154	12.5040	1288.9451	2933.1820	3911.1136	4628.2080	5506.1159	7901.7761	10626.1028	0.9999
Sigma[3,4]	-927.9898	11.1100	1120.2456	-3274.3922	-1581.1361	-882.2728	-221.4778	1187.2046	10167.1099	1.0001
Sigma[4,1]	-895.1376	10.5683	1046.0132	-3099.4280	-1502.9065	-845.1212	-245.0661	1089.4009	9796.2594	0.9999
Sigma[4,2]	6049.8894	17.7054	1659.0400	3600.2437	4882.6608	5792.3560	6932.5972	10010.5246	8780.1923	1.0001
Sigma[4,3]	-927.9898	11.1100	1120.2456	-3274.3922	-1581.1361	-882.2728	-221.4778	1187.2046	10167.1099	1.0001
Sigma[4,4]	7328.0205	20.3749	1947.8506	4457.2911	5953.9365	7019.6965	8353.3850	12018.9037	9139.4579	1.0001
gB	1199.5421	452.3817	15037.5259	17.8189	71.9612	158.5955	395.4852	5006.9099	1104.9500	1.0022
tB[1]	-9.6570	1.1221	23.8785	-44.6830	-16.4110	-9.7694	-3.7578	23.6232	452.8328	1.0055
tB[2]	11.1035	1.1313	23.9773	-21.5031	3.9100	9.9394	16.6674	46.5721	449.2282	1.0059
muS[1,1]	314.8688	0.0821	10.9434	292.8773	307.7210	314.8517	322.0667	336.5457	17782.6318	0.9999
muS[1,2]	335.6294	0.0936	11.7163	312.1191	328.0705	335.6798	343.2441	358.5727	15674.5437	0.9999
muS[2,1]	314.8688	0.0821	10.9434	292.8773	307.7210	314.8517	322.0667	336.5457	17782.6318	0.9999
muS[2,2]	335.6294	0.0936	11.7163	312.1191	328.0705	335.6798	343.2441	358.5727	15674.5437	0.9999
lp__	-750.6344	0.0523	3.0627	-757.6283	-752.4650	-750.2407	-748.3980	-745.7810	3424.0600	1.0011

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)
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	324.7660	0.1682	13.1150	299.5784	316.7464	324.5216	332.5882	350.2166	6083.277	1.0005
Sigma[1,1]	3939.6211	10.5189	1051.7492	2397.9313	3201.9899	3770.7747	4483.2001	6436.1899	9997.260	1.0001
Sigma[1,2]	4080.9161	11.4201	1127.9760	2427.8173	3281.6941	3908.2569	4683.5924	6771.7922	9755.674	1.0001
Sigma[1,3]	-211.1970	8.0942	844.8737	-1954.8636	-719.1901	-190.5303	317.0762	1445.2736	10895.306	1.0001
Sigma[1,4]	-822.1191	9.6324	1005.9684	-2993.6337	-1408.4730	-773.6490	-179.5888	1030.5648	10906.822	1.0001
Sigma[2,1]	4080.9161	11.4201	1127.9760	2427.8173	3281.6941	3908.2569	4683.5924	6771.7922	9755.674	1.0001
Sigma[2,2]	4813.0508	12.8724	1291.1311	2921.1338	3902.8467	4608.6852	5494.5149	7865.7046	10060.589	1.0000
Sigma[2,3]	-586.8072	9.1976	941.9792	-2595.6812	-1140.9951	-551.0554	4.5691	1206.9915	10488.983	1.0001
Sigma[2,4]	-985.4674	10.9210	1119.7198	-3406.0039	-1630.2388	-925.0385	-266.0532	1076.3914	10512.249	1.0001
Sigma[3,1]	-211.1970	8.0942	844.8737	-1954.8636	-719.1901	-190.5303	317.0762	1445.2736	10895.306	1.0001
Sigma[3,2]	-586.8072	9.1976	941.9792	-2595.6812	-1140.9951	-551.0554	4.5691	1206.9915	10488.983	1.0001
Sigma[3,3]	5213.0528	15.2070	1372.1623	3179.4791	4236.8721	4996.6313	5939.1489	8550.3970	8141.904	1.0005
Sigma[3,4]	5918.1722	17.4774	1584.9204	3559.5527	4793.2915	5668.9480	6765.0565	9750.8178	8223.599	1.0005
Sigma[4,1]	-822.1191	9.6324	1005.9684	-2993.6337	-1408.4730	-773.6490	-179.5888	1030.5648	10906.822	1.0001
Sigma[4,2]	-985.4674	10.9210	1119.7198	-3406.0039	-1630.2388	-925.0385	-266.0532	1076.3914	10512.249	1.0001
Sigma[4,3]	5918.1722	17.4774	1584.9204	3559.5527	4793.2915	5668.9480	6765.0565	9750.8178	8223.599	1.0005
Sigma[4,4]	7250.7954	20.6727	1905.8023	4417.8522	5908.4479	6936.4508	8253.1417	11831.2176	8498.908	1.0004
gB	50.9211	9.1286	480.1744	0.1979	0.7420	2.1738	9.5130	323.8054	2766.892	1.0005
gAB	328.7440	11.7113	813.8048	35.5445	95.6726	173.9139	337.0359	1443.2714	4828.677	1.0014
tB[1]	0.6784	0.0980	5.6628	-6.0248	-0.7644	0.1403	1.2345	12.8908	3339.034	1.0006
tB[2]	-1.0475	0.1482	6.2996	-14.1900	-1.3281	-0.1698	0.6976	5.5380	1806.933	1.0009
tAB[1,1]	-0.5928	0.1376	10.9916	-21.1248	-6.6978	-1.1380	5.0848	24.1004	6377.400	1.0001
tAB[1,2]	-17.7514	0.1431	10.9731	-41.9668	-23.2939	-16.8631	-11.1595	1.1214	5881.841	1.0029
tAB[2,1]	10.4326	0.1432	11.6281	-9.7708	3.5072	9.4327	16.3896	37.2052	6591.423	1.0000
tAB[2,2]	7.5456	0.1451	11.5346	-17.4691	1.6719	8.2567	14.2887	28.4023	6322.038	1.0026
muS[1,1]	324.8516	0.0774	9.7636	305.7462	318.3236	324.8131	331.3332	344.1036	15930.671	1.0002
muS[1,2]	305.9671	0.0927	10.6987	284.4389	298.9882	306.0087	313.1088	326.8298	13320.496	1.0001
muS[2,1]	335.8770	0.0893	10.8466	314.9944	328.5818	335.7997	343.2356	357.0951	14762.359	1.0004
muS[2,2]	331.2641	0.1059	11.9754	307.1721	323.6051	331.3215	339.1692	354.7449	12798.620	1.0000
lp__	-755.2794	0.0743	4.1194	-764.5857	-757.7836	-754.8525	-752.2841	-748.5721	3073.896	1.0004

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)

## Warning: There were 2629 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

## Warning: The largest R-hat is 1.23, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#r-hat

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess

## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

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	316.3952	1.4242	53.3077	217.1720	298.6348	319.5790	335.7031	406.6356	1401.0926	1.0092
Sigma[1,1]	4617.3468	285.1040	1337.5583	2539.7789	3566.2861	4410.2878	5665.9734	7468.8339	22.0099	1.0980
Sigma[1,2]	4481.7774	204.0236	1292.7376	2467.8560	3492.7071	4321.6307	5466.9623	7502.1175	40.1476	1.0606
Sigma[1,3]	-390.5109	11.5171	883.8528	-2317.1135	-814.2082	-382.1428	74.6686	1379.6070	5889.4337	1.0005
Sigma[1,4]	-1414.3518	49.2912	1257.9844	-4154.3592	-2093.4336	-1380.3207	-702.7670	1196.0824	651.3461	1.0051
Sigma[2,1]	4481.7774	204.0236	1292.7376	2467.8560	3492.7071	4321.6307	5466.9623	7502.1175	40.1476	1.0606
Sigma[2,2]	5071.8828	104.2199	1364.6441	2972.9479	4089.8832	5006.5826	5782.8540	8449.1464	171.4501	1.0234
Sigma[2,3]	-602.8546	35.4297	919.1777	-2503.2490	-1023.1755	-717.1860	-57.7689	1274.6877	673.0765	1.0084
Sigma[2,4]	-1082.6429	100.8475	1297.6628	-3688.1991	-1803.3440	-1169.6535	-345.5975	1734.6368	165.5747	1.0165
Sigma[3,1]	-390.5109	11.5171	883.8528	-2317.1135	-814.2082	-382.1428	74.6686	1379.6070	5889.4337	1.0005
Sigma[3,2]	-602.8546	35.4297	919.1777	-2503.2490	-1023.1755	-717.1860	-57.7689	1274.6877	673.0765	1.0084
Sigma[3,3]	5227.6740	313.7778	1594.8365	3310.9095	3889.3156	4924.2135	6069.8157	9156.0216	25.8338	1.0810
Sigma[3,4]	6058.9199	568.2586	2214.5434	3251.4712	4305.2532	5821.7666	7284.3374	11270.2321	15.1872	1.1370
Sigma[4,1]	-1414.3518	49.2912	1257.9844	-4154.3592	-2093.4336	-1380.3207	-702.7670	1196.0824	651.3461	1.0051
Sigma[4,2]	-1082.6429	100.8475	1297.6628	-3688.1991	-1803.3440	-1169.6535	-345.5975	1734.6368	165.5747	1.0165
Sigma[4,3]	6058.9199	568.2586	2214.5434	3251.4712	4305.2532	5821.7666	7284.3374	11270.2321	15.1872	1.1370
Sigma[4,4]	8118.5459	815.7402	3033.9184	3718.4887	5834.1200	7846.4494	9819.2814	15064.6601	13.8326	1.1524
gB	7661.3251	1614.9561	107912.2055	0.7256	183.7051	970.8173	2858.2445	38917.4504	4464.9817	1.0009
tB[1]	26.1621	1.5250	54.1001	-54.4973	0.2016	20.4786	44.7831	132.4219	1258.5498	1.0253
tB[2]	-25.2760	1.4986	53.8116	-129.9902	-44.5907	-19.8871	-0.9118	58.6662	1289.3003	1.0137
muS[1,1]	342.5573	2.0132	16.9598	309.3053	334.0330	342.2541	353.7713	374.2416	70.9704	1.0437
muS[1,2]	291.1192	8.4855	24.0343	252.7249	273.8613	286.2527	305.9240	337.8376	8.0225	1.2952
muS[2,1]	342.5573	2.0132	16.9598	309.3053	334.0330	342.2541	353.7713	374.2416	70.9704	1.0437
muS[2,2]	291.1192	8.4855	24.0343	252.7249	273.8613	286.2527	305.9240	337.8376	8.0225	1.2952
lp__	-757.0845	0.1863	2.9500	-763.8858	-758.7800	-756.4767	-755.0283	-752.3513	250.6964	1.0208

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):
## 
##  -721.8924
## 
## Error Measures:
## 
##  Relative Mean-Squared Error: 0.001107519
##  Coefficient of Variation: 0.03327942
##  Percentage Error: 3%
## 
## 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):
## 
##  -726.2256
## 
## Error Measures:
## 
##  Relative Mean-Squared Error: 0.0009735747
##  Coefficient of Variation: 0.03120216
##  Percentage Error: 3%
## 
## 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: 0.01313

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

## 
## Bridge sampling log marginal likelihood estimate 
## (method = "normal", repetitions = 1):
## 
##  -721.852
## 
## Error Measures:
## 
##  Relative Mean-Squared Error: 0.0005551869
##  Coefficient of Variation: 0.0235624
##  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: 1.04124

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

## 
## Bridge sampling log marginal likelihood estimate 
## (method = "normal", repetitions = 1):
## 
##  -729.1096
## 
## Error Measures:
## 
##  Relative Mean-Squared Error: 0.0001976231
##  Coefficient of Variation: 0.01405785
##  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.05591

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

## 
## Bridge sampling log marginal likelihood estimate 
## (method = "normal", repetitions = 1):
## 
##  -722.179
## 
## Error Measures:
## 
##  Relative Mean-Squared Error: 0.0003785044
##  Coefficient of Variation: 0.01945519
##  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: 0.75083

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

## 
## Bridge sampling log marginal likelihood estimate 
## (method = "normal", repetitions = 1):
## 
##  -732.1003
## 
## Error Measures:
## 
##  Relative Mean-Squared Error: 0.002616586
##  Coefficient of Variation: 0.05115258
##  Percentage Error: 5%
## 
## 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.00281

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)

## Warning: There were 495 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

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)

## Warning: There were 1447 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

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess

## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

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: 0.05786

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.92596

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)

## Warning: There were 2 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

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.09754

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)

## Warning: There were 5374 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

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess

## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

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.00122

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)

## Warning: There were 1613 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

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess

## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

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.00489

summary(M.full.fixed)

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

summary(M.omitAB.fixed)

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

summary(M.omitA.fixed)

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

summary(M.notA.fixed)

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

summary(M.omitB.fixed)

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

summary(M.notB.fixed)

## 
## Bridge sampling log marginal likelihood estimate 
## (method = "normal", repetitions = 1):
## 
##  -732.7993
## 
## Error Measures:
## 
##  Relative Mean-Squared Error: 0.0001645306
##  Coefficient of Variation: 0.01282695
##  Percentage Error: 1%
## 
## 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	325.0409	0.0943	8.7048	307.6073	319.3771	325.0223	330.7564	342.0622	8516.170	1.0001
Sigma[1,1]	3989.8251	13.0054	1065.4184	2426.2960	3236.1909	3826.2366	4536.5310	6527.7185	6711.095	1.0000
Sigma[1,2]	-251.0728	10.7584	859.7961	-2000.1905	-768.8178	-240.1650	269.6804	1496.0763	6387.007	1.0002
Sigma[1,3]	4130.7483	14.2844	1141.2795	2452.8679	3323.6345	3942.2047	4725.8493	6850.5547	6383.557	1.0000
Sigma[1,4]	-870.7396	12.9612	1025.5778	-3070.2902	-1456.9906	-821.8356	-228.1072	1075.4097	6261.064	1.0002
Sigma[2,1]	-251.0728	10.7584	859.7961	-2000.1905	-768.8178	-240.1650	269.6804	1496.0763	6387.007	1.0002
Sigma[2,2]	5279.4880	18.9865	1415.2435	3199.6229	4284.1982	5054.7844	6019.1882	8632.6203	5556.113	0.9999
Sigma[2,3]	-620.9291	12.1432	954.7157	-2632.8236	-1191.5139	-589.3017	-30.7792	1236.0249	6181.288	1.0002
Sigma[2,4]	5993.7271	21.9836	1631.8976	3584.0242	4845.5326	5738.9666	6878.7373	9839.4048	5510.455	0.9999
Sigma[3,1]	4130.7483	14.2844	1141.2795	2452.8679	3323.6345	3942.2047	4725.8493	6850.5547	6383.557	1.0000
Sigma[3,2]	-620.9291	12.1432	954.7157	-2632.8236	-1191.5139	-589.3017	-30.7792	1236.0249	6181.288	1.0002
Sigma[3,3]	4861.5779	16.2652	1306.0781	2947.0762	3933.7925	4644.1210	5546.3410	8015.1226	6447.956	1.0000
Sigma[3,4]	-1028.9214	14.5566	1131.5076	-3458.7491	-1686.1945	-977.5701	-327.7517	1102.2540	6042.212	1.0002
Sigma[4,1]	-870.7396	12.9612	1025.5778	-3070.2902	-1456.9906	-821.8356	-228.1072	1075.4097	6261.064	1.0002
Sigma[4,2]	5993.7271	21.9836	1631.8976	3584.0242	4845.5326	5738.9666	6878.7373	9839.4048	5510.455	0.9999
Sigma[4,3]	-1028.9214	14.5566	1131.5076	-3458.7491	-1686.1945	-977.5701	-327.7517	1102.2540	6042.212	1.0002
Sigma[4,4]	7340.5307	25.9748	1959.2659	4448.9856	5949.8954	7030.7553	8384.6186	11914.0370	5689.618	0.9999
gA	642.2077	296.4809	39181.2588	0.0538	0.2242	0.8171	7.0023	1800.4606	17464.790	1.0000
gB	931.3240	186.3985	17923.7125	17.3691	59.7568	127.3723	318.8842	3485.4309	9246.363	1.0001
gAB	441.1502	125.3022	10162.6613	0.3839	14.2784	33.9417	89.5614	1123.6558	6578.056	1.0001
tAfix	3.6565	0.2546	9.7505	-2.3490	-0.2575	0.2422	1.5869	37.7440	1466.635	1.0007
tBfix	-13.3700	0.0343	2.8149	-18.6657	-15.2388	-13.4207	-11.6194	-7.7141	6753.302	1.0009
tABfix	6.8695	0.0527	2.5568	0.2973	5.6368	7.2217	8.5807	11.0220	2355.603	1.0002
muS[1,1]	321.6072	0.1507	10.0338	303.4033	314.9556	321.0816	327.4295	343.8049	4431.082	1.0002
muS[1,2]	333.6457	0.1910	11.2431	313.4138	326.4189	332.7329	339.9385	359.9406	3465.830	1.0001
muS[2,1]	309.5666	0.2030	11.0669	282.9468	303.9031	310.5352	316.6321	328.4362	2971.781	1.0003
muS[2,2]	335.3442	0.2645	12.8866	301.4868	329.4884	336.9432	343.5095	355.9797	2373.328	1.0006
lp__	-748.4755	0.0701	3.6816	-756.4762	-750.8672	-748.0682	-745.7889	-742.4584	2758.818	1.0007

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	322.5986	1.0252	10.3658	303.4100	315.9074	321.9248	328.5464	347.0021	102.2373	1.0203
Sigma[1,1]	3985.5698	54.8138	1129.5918	2357.6824	3172.6788	3808.2835	4617.2207	6657.9535	424.6809	1.0052
Sigma[1,2]	-262.3969	30.6517	873.2037	-2082.5770	-753.9007	-261.2353	249.9594	1528.1480	811.5630	1.0013
Sigma[1,3]	4111.8334	51.8824	1194.2385	2420.9943	3240.9626	3932.4067	4759.7697	6949.9480	529.8374	1.0044
Sigma[1,4]	-913.5457	35.3186	1065.5882	-3231.7044	-1498.9452	-871.3663	-269.7705	1213.4150	910.2728	1.0012
Sigma[2,1]	-262.3969	30.6517	873.2037	-2082.5770	-753.9007	-261.2353	249.9594	1528.1480	811.5630	1.0013
Sigma[2,2]	5352.3131	29.4044	1446.9778	3228.7195	4339.8756	5107.2768	6140.0138	8849.1022	2421.5827	1.0021
Sigma[2,3]	-629.8416	32.5089	963.8731	-2672.0006	-1185.8781	-572.5744	-59.0591	1312.0077	879.0953	1.0008
Sigma[2,4]	6112.0822	34.3800	1690.6873	3654.9727	4920.0544	5843.7307	7004.7917	10204.5929	2418.3311	1.0019
Sigma[3,1]	4111.8334	51.8824	1194.2385	2420.9943	3240.9626	3932.4067	4759.7697	6949.9480	529.8374	1.0044
Sigma[3,2]	-629.8416	32.5089	963.8731	-2672.0006	-1185.8781	-572.5744	-59.0591	1312.0077	879.0953	1.0008
Sigma[3,3]	4845.8145	48.8723	1348.8712	2932.6659	3864.4623	4622.8820	5556.4286	8138.0391	761.7566	1.0033
Sigma[3,4]	-1088.1710	38.4475	1167.7528	-3643.4604	-1748.0517	-996.5325	-353.1476	1219.2247	922.4977	1.0009
Sigma[4,1]	-913.5457	35.3186	1065.5882	-3231.7044	-1498.9452	-871.3663	-269.7705	1213.4150	910.2728	1.0012
Sigma[4,2]	6112.0822	34.3800	1690.6873	3654.9727	4920.0544	5843.7307	7004.7917	10204.5929	2418.3311	1.0019
Sigma[4,3]	-1088.1710	38.4475	1167.7528	-3643.4604	-1748.0517	-996.5325	-353.1476	1219.2247	922.4977	1.0009
Sigma[4,4]	7545.5058	44.4709	2061.2685	4594.5637	6066.4607	7201.3460	8632.2816	12644.0114	2148.4142	1.0021
gA	6085.8401	794.5050	77807.8203	0.1225	211.2850	704.5368	1979.8959	26071.5742	9590.7633	1.0001
gB	757.6710	103.2837	11643.9397	12.6291	53.6309	125.2945	299.7210	3486.0410	12709.7337	1.0002
tAfix	30.9337	1.4504	16.8453	-0.3671	22.4550	34.4139	42.7997	57.2131	134.8994	1.0164
tBfix	-12.7808	0.0987	3.2633	-19.0295	-14.8740	-12.9132	-10.7563	-6.0253	1092.8016	1.0009
muS[1,1]	335.4346	0.5171	13.5910	306.7859	327.6432	336.4280	344.4412	360.2628	690.7889	1.0016
muS[1,2]	353.5095	0.4197	12.9174	326.8762	345.3074	354.6972	362.1084	377.1793	947.2268	1.0014
muS[2,1]	291.6878	1.9897	17.4575	263.0198	279.8089	289.0912	301.1281	334.4886	76.9845	1.0276
muS[2,2]	309.7626	2.0875	19.0338	278.0734	296.8591	307.0101	320.3238	356.7624	83.1385	1.0255
lp__	-749.7676	0.0483	3.0042	-756.5925	-751.5019	-749.4328	-747.5768	-744.9891	3873.7220	1.0002

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	325.4501	0.0637	8.5968	308.6939	319.7579	325.4468	331.2576	342.1753	18227.787	1.0001
Sigma[1,1]	4001.7726	10.8034	1095.6232	2419.3559	3236.1399	3811.2908	4573.2325	6679.4333	10285.017	0.9999
Sigma[1,2]	-285.8025	8.2525	851.0897	-2024.3592	-785.0434	-266.8482	245.9667	1373.4153	10635.968	1.0002
Sigma[1,3]	4144.0480	11.7465	1171.8122	2449.7183	3332.3901	3940.4028	4744.5238	6995.2628	9951.762	1.0000
Sigma[1,4]	-908.3642	9.9423	1019.2002	-3095.1022	-1499.4126	-851.7564	-255.9829	970.2324	10508.626	1.0001
Sigma[2,1]	-285.8025	8.2525	851.0897	-2024.3592	-785.0434	-266.8482	245.9667	1373.4153	10635.968	1.0002
Sigma[2,2]	5298.3400	13.9827	1409.3279	3208.2579	4289.5471	5081.7240	6054.7848	8651.7869	10158.851	1.0001
Sigma[2,3]	-663.7015	8.9805	949.5893	-2673.0387	-1217.3138	-615.2984	-59.6837	1099.4664	11180.679	1.0002
Sigma[2,4]	6012.1812	16.1765	1627.2581	3601.9910	4851.6551	5754.0668	6857.9436	9933.5390	10119.165	1.0001
Sigma[3,1]	4144.0480	11.7465	1171.8122	2449.7183	3332.3901	3940.4028	4744.5238	6995.2628	9951.762	1.0000
Sigma[3,2]	-663.7015	8.9805	949.5893	-2673.0387	-1217.3138	-615.2984	-59.6837	1099.4664	11180.679	1.0002
Sigma[3,3]	4877.7332	13.3022	1334.9546	2936.2819	3952.2076	4647.3267	5558.1520	8082.6517	10071.246	1.0000
Sigma[3,4]	-1074.9097	10.8422	1129.5288	-3500.8609	-1721.3350	-1001.8207	-342.9146	981.0976	10853.220	1.0002
Sigma[4,1]	-908.3642	9.9423	1019.2002	-3095.1022	-1499.4126	-851.7564	-255.9829	970.2324	10508.626	1.0001
Sigma[4,2]	6012.1812	16.1765	1627.2581	3601.9910	4851.6551	5754.0668	6857.9436	9933.5390	10119.165	1.0001
Sigma[4,3]	-1074.9097	10.8422	1129.5288	-3500.8609	-1721.3350	-1001.8207	-342.9146	981.0976	10853.220	1.0002
Sigma[4,4]	7357.6809	19.2211	1953.0367	4453.9474	5966.0288	7049.7167	8372.7524	12053.6537	10324.423	1.0001
gB	873.8697	104.3223	11437.3486	16.9231	60.8990	128.8510	323.3922	3997.5864	12019.777	1.0002
gAB	248.0621	34.1201	3191.1164	3.8347	18.0793	40.5261	104.1968	1244.6046	8747.105	1.0002
tBfix	-13.4842	0.0233	2.7718	-18.7446	-15.3022	-13.5702	-11.7362	-7.7896	14091.586	0.9999
tABfix	7.4953	0.0178	2.0093	3.2332	6.2837	7.5762	8.8265	11.1851	12779.082	1.0000
muS[1,1]	319.6630	0.0624	8.5735	302.6988	314.0109	319.6895	325.3723	336.4884	18863.183	1.0000
muS[1,2]	331.2372	0.0712	9.2251	313.1340	325.0859	331.2637	337.4708	349.1236	16788.451	1.0001
muS[2,1]	312.1677	0.0618	8.5448	295.3037	306.5010	312.1572	317.8298	329.0652	19112.779	1.0001
muS[2,2]	338.7325	0.0698	9.1325	320.7056	332.6699	338.7550	344.8902	356.2031	17126.452	1.0001
lp__	-744.7839	0.0397	3.0918	-751.8748	-746.6019	-744.4201	-742.5398	-739.8686	6066.855	1.0000

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	325.3195	0.0982	11.2487	303.2104	317.8866	325.3239	332.8368	347.3174	13122.144	1.0000
Sigma[1,1]	4150.6761	11.9831	1145.1797	2492.1713	3345.2789	3960.9748	4735.3388	6944.6088	9132.924	1.0001
Sigma[1,2]	-336.2336	8.6226	899.1307	-2201.4646	-867.6634	-322.1488	213.2139	1443.4081	10873.376	1.0000
Sigma[1,3]	4164.2302	12.4849	1165.5645	2467.0805	3340.4680	3976.7387	4774.7561	7008.3626	8715.762	1.0001
Sigma[1,4]	-878.3080	10.1879	1059.0216	-3119.9292	-1486.6742	-830.4490	-218.4696	1109.5698	10805.426	1.0000
Sigma[2,1]	-336.2336	8.6226	899.1307	-2201.4646	-867.6634	-322.1488	213.2139	1443.4081	10873.376	1.0000
Sigma[2,2]	5449.7289	15.3020	1514.6740	3244.7072	4386.3002	5197.7965	6241.2029	9086.6969	9798.102	1.0000
Sigma[2,3]	-572.6696	9.3904	979.1795	-2641.5522	-1143.5021	-545.7366	25.0886	1345.3892	10873.069	1.0000
Sigma[2,4]	6069.8362	17.1970	1696.7788	3588.3836	4878.1414	5803.8765	6937.7639	10140.1449	9735.223	1.0000
Sigma[3,1]	4164.2302	12.4849	1165.5645	2467.0805	3340.4680	3976.7387	4774.7561	7008.3626	8715.762	1.0001
Sigma[3,2]	-572.6696	9.3904	979.1795	-2641.5522	-1143.5021	-545.7366	25.0886	1345.3892	10873.069	1.0000
Sigma[3,3]	4822.9891	13.4766	1284.2664	2934.9814	3922.4609	4611.2467	5507.8513	7935.2146	9081.382	1.0001
Sigma[3,4]	-915.7335	10.9607	1139.3032	-3345.3560	-1580.8908	-861.4454	-215.1876	1261.8467	10804.412	1.0000
Sigma[4,1]	-878.3080	10.1879	1059.0216	-3119.9292	-1486.6742	-830.4490	-218.4696	1109.5698	10805.426	1.0000
Sigma[4,2]	6069.8362	17.1970	1696.7788	3588.3836	4878.1414	5803.8765	6937.7639	10140.1449	9735.223	1.0000
Sigma[4,3]	-915.7335	10.9607	1139.3032	-3345.3560	-1580.8908	-861.4454	-215.1876	1261.8467	10804.412	1.0000
Sigma[4,4]	7349.0582	20.0523	1994.6354	4419.3964	5939.0007	7035.5050	8369.8624	12186.7824	9894.653	1.0000
gB	884.6462	86.7082	10062.1146	15.9160	68.8230	154.8318	395.6831	4231.2290	13466.596	0.9999
tBfix	-14.6190	0.0441	3.7443	-21.4655	-17.0949	-14.7839	-12.3657	-6.7330	7202.615	1.0001
muS[1,1]	314.9823	0.0933	11.1305	292.9974	307.5778	314.9502	322.3411	337.0867	14246.800	1.0000
muS[1,2]	335.6567	0.1116	11.9665	311.5453	327.8864	335.8109	343.6339	358.6823	11491.299	0.9999
muS[2,1]	314.9823	0.0933	11.1305	292.9974	307.5778	314.9502	322.3411	337.0867	14246.800	1.0000
muS[2,2]	335.6567	0.1116	11.9665	311.5453	327.8864	335.8109	343.6339	358.6823	11491.299	0.9999
lp__	-747.3300	0.0363	2.8301	-753.8211	-748.9649	-746.9727	-745.2807	-742.8811	6092.863	0.9999

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	317.5686	0.4291	11.7877	293.6524	309.3724	317.9777	325.5003	340.1666	754.6592	1.0040
Sigma[1,1]	4071.3331	42.2074	1134.1441	2477.6468	3270.7948	3856.9446	4623.0630	6872.4153	722.0357	1.0043
Sigma[1,2]	4252.3669	51.4022	1269.7522	2499.6347	3356.0283	4012.6659	4945.6698	7298.9072	610.2043	1.0057
Sigma[1,3]	-232.1627	33.9914	885.1166	-1994.7756	-743.1689	-232.1149	291.0636	1615.2296	678.0505	1.0032
Sigma[1,4]	-840.2921	39.4178	1113.6007	-3258.7594	-1410.8059	-780.3715	-172.7226	1204.4681	798.1312	1.0020
Sigma[2,1]	4252.3669	51.4022	1269.7522	2499.6347	3356.0283	4012.6659	4945.6698	7298.9072	610.2043	1.0057
Sigma[2,2]	5383.8291	71.0469	1600.5828	3154.2732	4262.5218	5079.4480	6256.1750	9427.3143	507.5346	1.0070
Sigma[2,3]	-748.0475	42.3970	1023.4782	-2812.1780	-1405.0201	-732.3035	-94.9608	1249.8996	582.7575	1.0044
Sigma[2,4]	-811.1272	49.6522	1247.6187	-3242.9597	-1613.4410	-842.7540	-31.0966	1646.0297	631.3744	1.0041
Sigma[3,1]	-232.1627	33.9914	885.1166	-1994.7756	-743.1689	-232.1149	291.0636	1615.2296	678.0505	1.0032
Sigma[3,2]	-748.0475	42.3970	1023.4782	-2812.1780	-1405.0201	-732.3035	-94.9608	1249.8996	582.7575	1.0044
Sigma[3,3]	5292.4137	49.9698	1417.7650	3170.2382	4363.8539	4983.7887	5992.5434	8846.0621	804.9937	1.0012
Sigma[3,4]	5871.2879	59.5037	1688.8689	3417.2210	4779.8294	5541.6843	6688.7313	10128.4131	805.5711	1.0007
Sigma[4,1]	-840.2921	39.4178	1113.6007	-3258.7594	-1410.8059	-780.3715	-172.7226	1204.4681	798.1312	1.0020
Sigma[4,2]	-811.1272	49.6522	1247.6187	-3242.9597	-1613.4410	-842.7540	-31.0966	1646.0297	631.3744	1.0041
Sigma[4,3]	5871.2879	59.5037	1688.8689	3417.2210	4779.8294	5541.6843	6688.7313	10128.4131	805.5711	1.0007
Sigma[4,4]	7439.1740	78.7928	2218.8358	4265.0367	5942.5622	7027.5670	8459.4345	13070.7635	793.0074	1.0005
gB	3952.0131	1675.9374	220892.8692	0.0996	0.3994	2.1846	309.1141	10388.2820	17371.9188	1.0000
gAB	171.0462	24.5957	1573.8658	0.0849	2.0785	23.1468	77.0935	847.8156	4094.6396	1.0003
tBfix	12.3164	0.9275	20.1224	-2.2135	-0.1529	0.6775	22.7302	61.5014	470.6595	1.0002
tABfix	5.8659	0.2343	4.1899	-0.4783	1.2744	6.5374	9.2861	12.7937	319.8005	1.0045
muS[1,1]	329.2105	0.6269	15.9704	301.1049	318.3516	327.6583	338.7635	364.1807	648.8914	1.0024
muS[1,2]	305.9266	0.8320	18.8400	263.1943	295.1347	309.2297	319.8716	335.0360	512.7769	1.0015
muS[2,1]	323.3446	0.7823	18.0550	293.3341	311.0516	321.2591	333.4898	363.7076	532.6528	1.0034
muS[2,2]	311.7925	0.9361	21.1468	263.3307	299.5463	315.9160	326.6412	342.6491	510.2976	1.0005
lp__	-753.9804	0.1409	3.3039	-761.2751	-756.0343	-753.7629	-751.6479	-748.3595	549.8921	1.0021

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	314.6838	0.4117	11.9781	292.0943	306.5684	314.5587	322.3841	338.4631	846.3784	1.0003
Sigma[1,1]	4247.7137	32.2665	1255.5813	2507.7557	3354.1546	4034.7307	4919.5039	7254.7410	1514.2099	1.0005
Sigma[1,2]	4201.2519	35.5229	1232.5858	2439.0395	3295.7715	4001.3388	4868.9792	7104.6811	1203.9730	1.0006
Sigma[1,3]	-367.6653	25.1199	973.2926	-2460.0089	-906.1508	-331.0974	235.5555	1550.9430	1501.2449	1.0004
Sigma[1,4]	-1263.6444	47.7385	1369.6868	-4317.0307	-2027.9230	-1165.6059	-379.6766	1115.5405	823.1984	1.0015
Sigma[2,1]	4201.2519	35.5229	1232.5858	2439.0395	3295.7715	4001.3388	4868.9792	7104.6811	1203.9730	1.0006
Sigma[2,2]	4903.3284	42.3126	1347.4932	2929.4143	3942.5715	4682.4914	5652.1045	8185.8157	1014.1777	1.0007
Sigma[2,3]	-564.1260	27.2782	1005.2603	-2674.3130	-1141.0839	-526.1943	11.9252	1426.9880	1358.0814	1.0015
Sigma[2,4]	-911.3036	46.1168	1400.6785	-3852.4087	-1700.7884	-868.9347	-75.2091	1669.3404	922.4831	1.0019
Sigma[3,1]	-367.6653	25.1199	973.2926	-2460.0089	-906.1508	-331.0974	235.5555	1550.9430	1501.2449	1.0004
Sigma[3,2]	-564.1260	27.2782	1005.2603	-2674.3130	-1141.0839	-526.1943	11.9252	1426.9880	1358.0814	1.0015
Sigma[3,3]	5582.3701	46.7303	1610.2904	3304.5423	4459.1897	5278.8312	6425.8257	9352.5803	1187.4377	1.0007
Sigma[3,4]	6551.8647	63.6010	2033.8074	3694.3992	5029.5935	6175.3565	7622.7108	11501.3899	1022.5665	1.0005
Sigma[4,1]	-1263.6444	47.7385	1369.6868	-4317.0307	-2027.9230	-1165.6059	-379.6766	1115.5405	823.1984	1.0015
Sigma[4,2]	-911.3036	46.1168	1400.6785	-3852.4087	-1700.7884	-868.9347	-75.2091	1669.3404	922.4831	1.0019
Sigma[4,3]	6551.8647	63.6010	2033.8074	3694.3992	5029.5935	6175.3565	7622.7108	11501.3899	1022.5665	1.0005
Sigma[4,4]	8763.0621	101.5860	2726.2480	4965.0293	6700.2147	8311.4559	10215.2448	15396.7646	720.2162	1.0006
gB	7765.9087	1122.7821	136183.5041	0.2772	314.0626	1081.2776	3019.9945	35422.3873	14711.5380	1.0001
tBfix	38.6226	1.5654	21.2133	-0.3804	27.7347	43.2327	53.6707	71.3655	183.6377	1.0079
muS[1,1]	341.9941	1.1816	17.9768	302.4034	331.1278	343.6966	354.3856	373.5866	231.4753	1.0068
muS[1,2]	287.3735	1.1598	20.3418	251.8476	273.0330	285.2679	300.4780	332.7824	307.6253	1.0035
muS[2,1]	341.9941	1.1816	17.9768	302.4034	331.1278	343.6966	354.3856	373.5866	231.4753	1.0068
muS[2,2]	287.3735	1.1598	20.3418	251.8476	273.0330	285.2679	300.4780	332.7824	307.6253	1.0035
lp__	-752.8758	0.0629	2.7816	-759.2151	-754.4954	-752.4964	-750.9074	-748.5301	1956.9290	1.0009

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).