============================================================================================================
Part 1: Using JAGS, https://rpubs.com/sherloconan/1015445 (unreliable estimates)

Part 2: Using Stan, https://rpubs.com/sherloconan/1024130

The paradox on real data, https://rpubs.com/sherloconan/1033617

The paradox on simulated data, https://rpubs.com/sherloconan/1036792

Three ways to standardize effects, https://rpubs.com/sherloconan/1071855

Three ways to standardize effects (continued), https://rpubs.com/sherloconan/1071863

GLS and BIC approximation, https://rpubs.com/sherloconan/1093319

MORE: https://osf.io/5dm28/

1. Getting Started

 

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

  • RT the dependent variable, which is the number of seconds that it took to complete a puzzle;

  • ID which is a participant identifier;

  • shape and color, which are two factors that describe the type of puzzle solved.

shape and color each have two levels, and each of 12 participants completed puzzles within combination of shape and color. The design is thus \(2\times2\) factorial within-subject (Morey, 2022). Let shape be factor A (with levels \(\alpha_i\)) and color be factor B (with levels \(\beta_j\)) for \(i,j\in\{1,2\}\).

data(puzzles)
dat <- matrix(puzzles$RT, ncol=4,
              dimnames=list(paste0("s",1:12),
                            c("round&mono","square&mono","round&color","square&color"))) #wide data format
cov(dat) %>% #sample covariance matrix
  kable(digits=3) %>% kable_styling(full_width=F)
round&mono square&mono round&color square&color
round&mono 4.727 3.273 4.545 4.545
square&mono 3.273 6.000 5.909 4.727
round&color 4.545 5.909 7.636 5.273
square&color 4.545 4.727 5.273 7.455

 

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 ~ shape * color + ID, data=puzzles, whichRandom="ID", whichModels="top", progress=F)
## Bayes factor top-down analysis
## --------------
## When effect is omitted from shape + color + shape:color + ID , BF is...
## [1] Omit color:shape : 2.780721  ±4.21%
## [2] Omit color       : 0.280169  ±11.35%
## [3] Omit shape       : 0.2494565 ±3.3%
## 
## Against denominator:
##   RT ~ shape + color + shape:color + ID 
## ---
## 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("ID")

set.seed(277)
full <-   lmBF(RT~shape+color+ID+shape:color+shape:ID+color:ID, puzzles, whichRandom=randomEffects, progress=F)

set.seed(277)
omitAB <- lmBF(RT~shape+color+ID+            shape:ID+color:ID, puzzles, whichRandom=randomEffects, progress=F)
omitAB / full
## Bayes factor analysis
## --------------
## [1] shape + color + ID + shape:ID + color:ID : 2.664999 ±3.24%
## 
## Against denominator:
##   RT ~ shape + color + ID + shape:color + shape:ID + color:ID 
## ---
## Bayes factor type: BFlinearModel, JZS

 

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

set.seed(277)
omitA <-  lmBF(RT~      color+ID+shape:color+shape:ID+color:ID, puzzles, whichRandom=randomEffects, progress=F)
omitA / full
## Bayes factor analysis
## --------------
## [1] color + ID + shape:color + shape:ID + color:ID : 0.5335185 ±3.21%
## 
## Against denominator:
##   RT ~ shape + color + ID + shape:color + shape:ID + color:ID 
## ---
## Bayes factor type: BFlinearModel, JZS

 

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

set.seed(277)
omitB <-  lmBF(RT~shape+      ID+shape:color+shape:ID+color:ID, puzzles, whichRandom=randomEffects, progress=F)
omitB / full
## Bayes factor analysis
## --------------
## [1] shape + ID + shape:color + shape:ID + color:ID : 0.4508712 ±3.2%
## 
## Against denominator:
##   RT ~ shape + color + ID + shape:color + shape:ID + color:ID 
## ---
## 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~      color+ID+                     color:ID, puzzles, whichRandom=randomEffects, progress=F)
notA / omitAB
## Bayes factor analysis
## --------------
## [1] color + ID + color:ID : 1.264948 ±1.75%
## 
## Against denominator:
##   RT ~ shape + color + ID + shape:ID + color:ID 
## ---
## 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~shape+      ID+            shape:ID         , puzzles, whichRandom=randomEffects, progress=F)
notB / omitAB
## Bayes factor analysis
## --------------
## [1] shape + ID + shape:ID : 2.993029 ±1.76%
## 
## Against denominator:
##   RT ~ shape + color + ID + shape:ID + color:ID 
## ---
## 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~      color+ID+            shape:ID+color:ID, puzzles, whichRandom=randomEffects, progress=F)
exclA / omitAB
## Bayes factor analysis
## --------------
## [1] color + ID + shape:ID + color:ID : 0.5450039 ±1.76%
## 
## Against denominator:
##   RT ~ shape + color + ID + shape:ID + color:ID 
## ---
## 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~shape+      ID+            shape:ID+color:ID, puzzles, whichRandom=randomEffects, progress=F)
exclB / omitAB
## Bayes factor analysis
## --------------
## [1] shape + ID + shape:ID + color:ID : 0.4599249 ±1.75%
## 
## Against denominator:
##   RT ~ shape + color + ID + shape:ID + color:ID 
## ---
## Bayes factor type: BFlinearModel, JZS

     

Method 1’:

  • Assume \(\alpha_i\), \(\beta_j\), \((\alpha\beta)_{ij}\), and \(s_k\) are the random effects.
randomEffects <- c("ID","shape","color")

set.seed(277)
full <-   lmBF(RT~shape+color+ID+shape:color                  , puzzles, whichRandom=randomEffects, progress=F)

set.seed(277)
omitAB <- lmBF(RT~shape+color+ID                              , puzzles, whichRandom=randomEffects, progress=F)
omitAB / full
## Bayes factor analysis
## --------------
## [1] shape + color + ID : 5.749677 ±7.08%
## 
## Against denominator:
##   RT ~ shape + color + ID + shape:color 
## ---
## Bayes factor type: BFlinearModel, JZS
set.seed(277)
omitA <-  lmBF(RT~      color+ID+shape:color                  , puzzles, whichRandom=randomEffects, progress=F)
omitA / full
## Bayes factor analysis
## --------------
## [1] color + ID + shape:color : 1.600381 ±5.1%
## 
## Against denominator:
##   RT ~ shape + color + ID + shape:color 
## ---
## Bayes factor type: BFlinearModel, JZS
set.seed(277)
omitB <-  lmBF(RT~shape+      ID+shape:color                  , puzzles, whichRandom=randomEffects, progress=F)
omitB / full
## Bayes factor analysis
## --------------
## [1] shape + ID + shape:color : 1.600381 ±5.1%
## 
## Against denominator:
##   RT ~ shape + color + ID + shape:color 
## ---
## 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~shape+color+ID+shape:color+shape:ID+color:ID, puzzles, whichRandom=randomEffects, progress=F)

set.seed(277)
omitAB <- lmBF(RT~shape+color+ID+            shape:ID+color:ID, puzzles, whichRandom=randomEffects, progress=F)
omitAB / full
## Bayes factor analysis
## --------------
## [1] shape + color + ID + shape:ID + color:ID : 6.419993 ±6.59%
## 
## Against denominator:
##   RT ~ shape + color + ID + shape:color + shape:ID + color:ID 
## ---
## 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~      color+ID+shape:color+shape:ID+color:ID, puzzles, whichRandom=randomEffects, progress=F)
omitA / full
## Bayes factor analysis
## --------------
## [1] color + ID + shape:color + shape:ID + color:ID : 1.845765 ±5.97%
## 
## Against denominator:
##   RT ~ shape + color + ID + shape:color + shape:ID + color:ID 
## ---
## 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~shape+      ID+shape:color+shape:ID+color:ID, puzzles, whichRandom=randomEffects, progress=F)
omitB / full
## Bayes factor analysis
## --------------
## [1] shape + ID + shape:color + shape:ID + color:ID : 1.839681 ±6.13%
## 
## Against denominator:
##   RT ~ shape + color + ID + shape:color + shape:ID + color:ID 
## ---
## 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~      color+ID+                     color:ID, puzzles, whichRandom=randomEffects, progress=F)
notA / omitAB
## Bayes factor analysis
## --------------
## [1] color + ID + color:ID : 1.245677 ±5.15%
## 
## Against denominator:
##   RT ~ shape + color + ID + shape:ID + color:ID 
## ---
## 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~shape+      ID+            shape:ID         , puzzles, whichRandom=randomEffects, progress=F)
notB / omitAB
## Bayes factor analysis
## --------------
## [1] shape + ID + shape:ID : 3.031622 ±5.16%
## 
## Against denominator:
##   RT ~ shape + color + ID + shape:ID + color:ID 
## ---
## 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~      color+ID+            shape:ID+color:ID, puzzles, whichRandom=randomEffects, progress=F)
exclA / omitAB
## Bayes factor analysis
## --------------
## [1] color + ID + shape:ID + color:ID : 0.6674978 ±5.45%
## 
## Against denominator:
##   RT ~ shape + color + ID + shape:ID + color:ID 
## ---
## 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~shape+      ID+            shape:ID+color:ID, puzzles, whichRandom=randomEffects, progress=F)
exclB / omitAB
## Bayes factor analysis
## --------------
## [1] shape + ID + shape:ID + color:ID : 0.5530826 ±5.43%
## 
## Against denominator:
##   RT ~ shape + color + ID + shape:ID + color:ID 
## ---
## Bayes factor type: BFlinearModel, JZS

   

2. Data

 

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). The data set.

df <- read.table("~/Documents/UVic/Project Meeting/paradox/UpDown3.txt", header=T, stringsAsFactors=T)
dat2 <- matrix(df$RT,ncol=4,byrow=T,
               dimnames=list(paste0("s",1:31),
                             c("onehand&A","onehand&M","twohands&A","twohands&M"))) #wide data format
cov(dat2[,c(1,3,2,4)]) %>% #sample covariance matrix
  kable(digits=3) %>% kable_styling(full_width=F)
onehand&A twohands&A onehand&M twohands&M
onehand&A 3937.495 -207.186 4066.825 -830.866
twohands&A -207.186 5202.157 -565.871 5922.451
onehand&M 4066.825 -565.871 4779.546 -979.485
twohands&M -830.866 5922.451 -979.485 7274.140

 

Or, the simulation.

simData <- function(seed=277, N=50, eqVar=T, null=F, mult=1) {
  #' Input -
  #' seed:  random seed
  #' N:     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.
  #' null:  whether the effect B is not significant
  #' mult:  response multiplier (to check the sensitivity issue)
  #'        
  #' Output - long format 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, ifelse(null, 0, 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, ifelse(null, 15, 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, ifelse(null, 0, 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, ifelse(null, 16, 17.8), 3.46) #extra variability and effect  <==
    a2b2 <- round(a2b2+x2)
  }
  
  data.frame("RT"=c(a1b1,a2b1,a1b2,a2b2)*mult,
             "subj"=paste0("s",1:N),
             "resp"=rep(rep(c("a1","a2"),2),each=N),
             "align"=rep(c("b1","b2"),each=N*2),
             stringsAsFactors=T)
}


nS <- 50
dim <- list(paste0("s",1:nS),
            c("onehand&A","twohands&A","onehand&M","twohands&M"))

df_EqVar <- simData(N=nS)
dat_EqVar <- matrix(df_EqVar$RT, ncol=4, dimnames=dim) #wide data format
cov(dat_EqVar) %>% #Equal Variance
  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 
df_UnEqVar <- simData(seed=1, N=nS, eqVar=F)
dat_UnEqVar <- matrix(df_UnEqVar$RT, ncol=4, dimnames=dim)
cov(dat_UnEqVar) %>% #Unequal Variances
  kable(digits=3) %>% kable_styling(full_width=F)
onehand&A twohands&A onehand&M twohands&M
onehand&A 13.990 4.020 12.367 1.702
twohands&A 4.020 93.306 7.184 73.551
onehand&M 12.367 7.184 14.531 5.327
twohands&M 1.702 73.551 5.327 77.031 
df_EqBnull <- simData(seed=1, N=nS, null=T)
dat_EqBnull <- matrix(df_EqBnull$RT, ncol=4, dimnames=dim)
cov(dat_EqBnull) %>% #Equal Variance, Not Significant B
  kable(digits=3) %>% kable_styling(full_width=F)
onehand&A twohands&A onehand&M twohands&M
onehand&A 114.792 22.940 110.880 30.898
twohands&A 22.940 116.725 27.049 108.957
onehand&M 110.880 27.049 126.092 33.063
twohands&M 30.898 108.957 33.063 118.719 
df_UnEqBnull <- simData(seed=1, N=nS, eqVar=F, null=T)
dat_UnEqBnull <- matrix(df_UnEqBnull$RT, ncol=4, dimnames=dim)
cov(dat_UnEqBnull) %>% #Unequal Variances, Not Significant B
  kable(digits=3) %>% kable_styling(full_width=F)
onehand&A twohands&A onehand&M twohands&M
onehand&A 13.990 4.020 12.331 2.108
twohands&A 4.020 93.306 6.388 73.000
onehand&M 12.331 6.388 13.824 4.813
twohands&M 2.108 73.000 4.813 76.498

   

3. Generalized Least Squares

 

Within-Subject Model ?= Between-Subjects Model + Alternative Covariance Structure of the Error Term

Specify the full model in (3), which has no terms for the subject-specific random effects at all but an unstructured error covariance matrix (or correlation matrix).

\[\begin{align} \tag{3} \mathcal{M}_{\mathrm{full}}^{*}:\ Y_{ijk}=&\mu+\alpha_i+\beta_j+(\alpha\beta)_{ij}+\epsilon_{ijk}^{*} \\ &\mathrm{for}\ i=1,\dotsb,a,\ j=1,\dotsb,b,\mathrm{ and}\ k=1,\dotsb,n. \end{align}\]

\[\begin{equation} \tag{4} \mathbf{V}_k\equiv\text{Cov}(\epsilon_{ijk}^*)=\sigma^2 \begin{pmatrix} 1 & \rho_{1,2} & \dotsb & \rho_{1,ab}\\ \rho_{2,1} & 1 & \dotsb & \rho_{2,ab}\\ \vdots & \vdots & \rho_{p,q} & \vdots\\ \rho_{ab,1} & \rho_{ab,2} & \dotsb & 1 \end{pmatrix}_{ab\times ab} \end{equation}\]

\[\begin{equation} \tag{5} \mathit{BF}_{01}\approx\text{exp}\left\{\frac{\text{BIC}(\mathcal{M}_1)-\text{BIC}(\mathcal{M}_0)}{2}\right\} \end{equation}\]

 

The correlation between any two observations within the same subject \(k\) can be different in (4). \(1+ab(ab-1)/2=7\) parameters in a \(2\times2\) unstructured matrix.

get_BIC_approx <- function(df, changeNames=T, meth="ML") {
  #' package -    nlme, lme4
  #' Input -
  #' df:          long format data frame for a 2x2 repeated-measures design
  #' changeNames: columns must be dependent variable,
  #'              subject identifier, Factor A, and Factor B
  #' meth:        estimation method, either "ML" or "REML"
  #' 
  #' Output -     BIC approximations of Bayes factors
  colnames(df) <- c("DV", "ID", "A", "B")
  
  #Testing the interaction effect AB
  gls_full <- nlme::gls(DV ~ A * B, data=df, method=meth, #Generalized Least Squares
                        correlation=nlme::corSymm(form=~1 | ID))
  gls_omitAB <- nlme::gls(DV ~ A + B, data=df, method=meth,
                          correlation=nlme::corSymm(form=~1 | ID))
  
  # lme_full <- nlme::lme(fixed=DV ~ A * B, random=~A + B | ID, #Linear Mixed-Effects
  #                       data=df, method=meth,
  #                       control=nlme::lmeControl(opt="optim", optCtrl=list(maxit=1e5)))
  # lme_omitAB <- nlme::lme(fixed=DV ~ A + B, random=~A + B | ID,
  #                         data=df, method=meth,
  #                         control=nlme::lmeControl(opt="optim", optCtrl=list(maxit=1e5)))
  lme_full <- lme4::lmer(DV ~ A * B + (A + B | ID), data=df, REML=(meth=="REML"), #Linear Mixed-Effects
                         control=lme4::lmerControl(optCtrl=list(maxfun=1e5)))   #BEWARE: boundary (singular) fit
  lme_omitAB <- lme4::lmer(DV ~ A + B + (A + B | ID), data=df, REML=(meth=="REML"),
                           control=lme4::lmerControl(optCtrl=list(maxfun=1e5)))
  
  #Testing the main effect A
  gls_omitA <- nlme::gls(DV ~ B, data=df, method=meth,
                         correlation=nlme::corSymm(form=~1 | ID))
  # lme_omitA <- nlme::lme(fixed=DV ~ B, random=~A + B | ID, #BEWARE: convergence error code = 1
  #                        data=df, method=meth,
  #                        control=nlme::lmeControl(opt="optim", optCtrl=list(maxit=1e5)))
  lme_omitA <- lme4::lmer(DV ~ B + (A + B | ID), data=df, REML=(meth=="REML"),
                          control=lme4::lmerControl(optCtrl=list(maxfun=1e5)))
  
  #Testing the main effect B
  gls_omitB <- nlme::gls(DV ~ A, data=df, method=meth,
                         correlation=nlme::corSymm(form=~1 | ID))
  # lme_omitB <- nlme::lme(fixed=DV ~ A, random=~A + B | ID,
  #                        data=df, method=meth,
  #                        control=nlme::lmeControl(opt="optim", optCtrl=list(maxit=1e5)))
  lme_omitB <- lme4::lmer(DV ~ A + (A + B | ID), data=df, REML=(meth=="REML"),
                          control=lme4::lmerControl(optCtrl=list(maxfun=1e5)))
  
  data.frame("BF_AB_GLS"=exp((BIC(gls_full)-BIC(gls_omitAB))/2),
             "BF_A_GLS"=exp((BIC(gls_omitAB)-BIC(gls_omitA))/2),
             "BF_B_GLS"=exp((BIC(gls_omitAB)-BIC(gls_omitB))/2),
             "BF_AB_LME"=exp((BIC(lme_full)-BIC(lme_omitAB))/2),
             "BF_A_LME"=exp((BIC(lme_omitAB)-BIC(lme_omitA))/2),
             "BF_B_LME"=exp((BIC(lme_omitAB)-BIC(lme_omitB))/2))
}
get_BIC_approx(puzzles) #example data
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
##   BF_AB_GLS  BF_A_GLS BF_B_GLS BF_AB_LME  BF_A_LME BF_B_LME
## 1  6.928203 0.4170107  0.03802  6.928203 0.3245586 0.331885

 

get_BIC_approx(df) #the paradox on real data
##   BF_AB_GLS BF_A_GLS    BF_B_GLS BF_AB_LME BF_A_LME     BF_B_LME
## 1 0.5335279 0.037708 0.006152088 0.9336214  4.39568 0.0003829248

 

get_BIC_approx(df_EqVar) #the paradox on simulated data, equal variance
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
##   BF_AB_GLS     BF_A_GLS BF_B_GLS BF_AB_LME     BF_A_LME BF_B_LME
## 1  6.119624 4.365432e-08 1.638499  5.675178 2.666837e-08 1.462907

 

get_BIC_approx(df_UnEqVar) #the paradox on simulated data, unequal variances
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues
##   BF_AB_GLS    BF_A_GLS     BF_B_GLS BF_AB_LME     BF_A_LME     BF_B_LME
## 1  5.628207 1.63522e-14 6.974966e-05  6.041731 8.318143e-16 0.0006722155

 

get_BIC_approx(df_EqBnull) #the paradox on simulated data, equal variance and null B
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
##   BF_AB_GLS     BF_A_GLS BF_B_GLS BF_AB_LME    BF_A_LME BF_B_LME
## 1   13.7268 8.602489e-07 13.44601  13.66466 1.01005e-06 13.50502

 

get_BIC_approx(df_UnEqBnull) #the paradox on simulated data, unequal variances and null B
## boundary (singular) fit: see help('isSingular')
##   BF_AB_GLS     BF_A_GLS BF_B_GLS BF_AB_LME     BF_A_LME BF_B_LME
## 1  5.870316 1.591414e-14  4.73571  6.299052 8.206495e-16 14.14213