Distributions of Random Effects

1. Getting Started

Random Intercept-Only (RIO): JASP v 0.16.2 or below; anovaBF v 0.9.12-4.6 or below.

\[\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}\]

Maximal Set of Random Effects (MRE): Include the random slopes in (2); JASP v 0.16.3 or later; use lmBF.

\[\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}\]

Note: \((\alpha\beta s)_{ijk}\) is estimable without aggregation.

Consider a 2×2 repeated-measures design. In a fixed factor, two levels are usually represented as dummy variables (0 and 1) in frequentist approaches, while they are contrast-coded (\(-\sqrt{2}/2\) and \(\sqrt{2}/2\)) in Bayesian analyses. See the discussion here.

dist <- function(data, legacy=F, REML=T, nlme=T) {
  #' Input  -
  #' data:         two-way repeated-measures data whose columns are "DV", "ID", "A", and "B" in order
  #' legacy:       if FALSE (default), fit an MRE model; if TRUE, fit a RIO model
  #' REML:         if TRUE (default), restricted maximum likelihood; if FALSE, maximum likelihood estimate
  #' nlme:         if TRUE (default), call the nlme::lme; if FALSE, call the lme4::lmer function
  #' Output -      a list of the MCMC draws' means and random effects estimates
  if(REML) {method <- list("REML", T)} else {method <- list("ML", F)} #estimation method
  colnames(data)[1:4] <- c("DV", "ID", "A", "B") #rename variables
  levels(data$A) <- paste0("a",1:length(levels(data$A))) #rename levels
  levels(data$B) <- paste0("b",1:length(levels(data$B)))
  
  if(!legacy) { #MRE
    fit_BF <- BayesFactor::lmBF(DV ~ A * B + ID + A:ID + B:ID, data=data, 
                                whichRandom="ID", iterations=1e5, progress=F)
    
    if(nlme) { #nlme v 3.1-164
      fit_LME <- nlme::lme(fixed=DV ~ A * B, random=~A + B | ID, data=data, method=method[[1]],
                           control=nlme::lmeControl(opt="optim", optCtrl=list(maxit=1e5)))
    } else { #lme4 v 1.1-35.1
      fit_LME <- lme4::lmer(DV ~ A * B + (A + B | ID), data=data, REML=method[[2]],
                            control=lme4::lmerControl(optCtrl=list(maxfun=1e5)))
    }

  } else { #RIO
    fit_BF <- BayesFactor::lmBF(DV ~ A * B + ID, data=data, 
                                whichRandom="ID", iterations=1e5, progress=F)
    
    if(nlme) { #nlme
      fit_LME <- nlme::lme(fixed=DV ~ A * B, random=~1 | ID, data=data, method=method[[1]],
                           control=nlme::lmeControl(opt="optim", optCtrl=list(maxit=1e5)))
    } else { #lme4
      fit_LME <- lme4::lmer(DV ~ A * B + (1 | ID), data=data, REML=method[[2]],
                            control=lme4::lmerControl(optCtrl=list(maxfun=1e5)))
    }
  }
  
  if(nlme) {ranef <- nlme::ranef(fit_LME)} else {ranef <- lme4::ranef(fit_LME)$ID} #extract the random effects
  
  draws <- BayesFactor::posterior(fit_BF, iterations=1e5, progress=F) #sample from the posterior
  pars <- colnames(draws) #parameter names
  randomIntercept <- colMeans(draws[,grepl("^ID-", pars)]) #posterior means of the random intercepts
  if(!legacy) { #MRE
    randomSlope_a2 <- colMeans(draws[,grepl("^A:ID-a2", pars)]) #posterior means of the random slopes for a2
    randomSlope_b2 <- colMeans(draws[,grepl("^B:ID-b2", pars)]) #posterior means of the random slopes for b2
    MCMC <- list(randomIntercept, randomSlope_a2, randomSlope_b2)
  } else { #RIO
    MCMC <- list(randomIntercept)
  }
  
  list("MCMC"=MCMC, "ranef"=ranef)
}

viz <- function(df, legacy=F) {
  #' Output -      Plot the distributions of random effects
  fit <- dist(df, legacy=legacy)
  MCMC <- fit$MCMC
  ranef <- fit$ranef
  nS <- nrow(ranef) #number of subjects
  nBreaks <- round(nS/3) #number of cells for the histogram
  
  if(legacy) { #RIO
    par(mfrow=c(1,1))
    hist(ranef[,1], breaks=nBreaks, prob=T, col="white", yaxt="n", ylab="",
         xlab=paste0("Random Intercepts (",nS,")"), main="Linear Mixed Model")
    lines(density(ranef[,1]), col="red", lwd=2)

    
    hist(MCMC[[1]], breaks=nBreaks, prob=T, col="white", yaxt="n", ylab="",
         xlab=paste0("Random Intercepts (",nS,")"), main="JZS Priors")
    lines(density(MCMC[[1]]), col="red", lwd=2)
    
  } else { #MRE
    par(mfrow=c(1,3))
    hist(ranef[,1], breaks=nBreaks, prob=T, col="white", yaxt="n", ylab="",
         xlab=paste0("Random Intercepts (",nS,")"), main="LME")
    lines(density(ranef[,1]), col="red", lwd=2)

    hist(ranef[,"Aa2"], breaks=nBreaks, prob=T, col="white", yaxt="n", ylab="",
         xlab=paste0("Random Slopes for a2 (",nS,")"), main="LME")
    lines(density(ranef[,"Aa2"]), col="red", lwd=2)

    hist(ranef[,"Bb2"], breaks=nBreaks, prob=T, col="white", yaxt="n", ylab="",
         xlab=paste0("Random Slopes for b2 (",nS,")"), main="LME")
    lines(density(ranef[,"Bb2"]), col="red", lwd=2)

    hist(MCMC[[1]], breaks=nBreaks, prob=T, col="white", yaxt="n", ylab="",
         xlab=paste0("Random Intercepts (",nS,")"), main="JZS")
    lines(density(MCMC[[1]]), col="red", lwd=2)

    hist(MCMC[[2]], breaks=nBreaks, prob=T, col="white", yaxt="n", ylab="",
         xlab=paste0("Random Slopes for a2 (",nS,")"), main="JZS")
    lines(density(MCMC[[2]]), col="red", lwd=2)

    hist(MCMC[[3]], breaks=nBreaks, prob=T, col="white", yaxt="n", ylab="",
         xlab=paste0("Random Slopes for b2 (",nS,")"), main="JZS")
    lines(density(MCMC[[3]]), col="red", lwd=2)
  }
}

2.1. The Paradox on Real Data

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

Let resp be factor A (with levels \(\alpha_i\) or \(a_i\)) and align be factor B (with levels \(\beta_j\) or \(b_j\)) for \(i,j\in\{1,2\}\).

The Paradox: Factor B is statistically significant in the repeated-measures ANOVA but not supported in the anovaBF R function (version 0.9.12-4.6 or lower).

One can break the data into two parts, separated by the factor A, and conduct separate Bayesian tests for the factor of interest B in each condition of A. Under these circumstances, a very clear effect B emerges in both analyses. This discrepancy is not the product of a contrast between a liberal and a conservative test (Bub, Masson, & van Noordenne, 2021, p. 80).

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

	onehand&A	twohands&A	onehand&M	twohands&M
onehand&A	3937.495	-207.186	4066.825	-830.866
twohands&A	-207.186	5202.157	-565.871	5922.451
onehand&M	4066.825	-565.871	4779.546	-979.485
twohands&M	-830.866	5922.451	-979.485	7274.140

summary(aov(RT~(resp*align)+Error(subj/(resp*align)),df)) #repeated-measures ANOVA

## 
## Error: subj
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 30 270038    9001               
## 
## Error: subj:resp
##           Df Sum Sq Mean Sq F value Pr(>F)
## resp       1  21479   21479   1.854  0.183
## Residuals 30 347540   11585               
## 
## Error: subj:align
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## align      1  11536   11536   28.22 9.64e-06 ***
## Residuals 30  12262     409                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: subj:resp:align
##            Df Sum Sq Mean Sq F value Pr(>F)  
## resp:align  1   1034  1033.6   5.203 0.0298 *
## Residuals  30   5960   198.7                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

\[\begin{align} \tag{3} \text{Total}\ \mathit{SS}&=\mathit{SS}_\text{subjects}+\mathit{SS}_{\text{within-subjects}} \\ &=\mathit{SS}_\text{subjects}+[(\mathit{SS}_\text{A}+\mathit{SS}_\text{Error A})+(\mathit{SS}_\text{B}+\mathit{SS}_\text{Error B})+(\mathit{SS}_\text{AB}+\mathit{SS}_\text{Error AB})] \end{align}\]

\[\begin{align} \tag{4} \text{Total}\ \mathit{df}&=\mathit{df}_\text{subjects}+[(\mathit{df}_\text{A}+\mathit{df}_\text{Error A})+(\mathit{df}_\text{B}+\mathit{df}_\text{Error B})+(\mathit{df}_\text{AB}+\mathit{df}_\text{Error AB})]\\ &=(n-1)+[((a-1)+(a-1)(n-1))+((b-1)+(b-1)(n-1))+((a-1)(b-1)+(a-1)(b-1)(n-1))] \end{align}\]

2.2. Results: RIO

JASP v 0.16.2 or below.

set.seed(277)
viz(df, legacy=T)

2.3. Results: MRE

JASP v 0.16.3 or later. See van den Bergh, Wagenmakers, and Aust (2023).

set.seed(277)
viz(df)

3.1. The Paradox on Simulated Data

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 #number of subjects
df_EqVar <- simData(N=nS) #equal variance
dat_EqVar <- matrix(df_EqVar$RT, ncol=4, #wide data format
                    dimnames=list(paste0("s", 1:nS),
                                  c("onehand&A","twohands&A","onehand&M","twohands&M")))
cov(dat_EqVar) %>%
  kable(digits=3) %>% kable_styling(full_width=F)

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

summary(aov(RT~(resp*align)+Error(subj/(resp*align)),df_EqVar)) #repeated-measures ANOVA

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

\[\begin{align} \tag{3} \text{Total}\ \mathit{SS}&=\mathit{SS}_\text{subjects}+\mathit{SS}_{\text{within-subjects}} \\ &=\mathit{SS}_\text{subjects}+[(\mathit{SS}_\text{A}+\mathit{SS}_\text{Error A})+(\mathit{SS}_\text{B}+\mathit{SS}_\text{Error B})+(\mathit{SS}_\text{AB}+\mathit{SS}_\text{Error AB})] \end{align}\]

\[\begin{align} \tag{4} \text{Total}\ \mathit{df}&=\mathit{df}_\text{subjects}+[(\mathit{df}_\text{A}+\mathit{df}_\text{Error A})+(\mathit{df}_\text{B}+\mathit{df}_\text{Error B})+(\mathit{df}_\text{AB}+\mathit{df}_\text{Error AB})]\\ &=(n-1)+[((a-1)+(a-1)(n-1))+((b-1)+(b-1)(n-1))+((a-1)(b-1)+(a-1)(b-1)(n-1))] \end{align}\]

3.2. Results: RIO

JASP v 0.16.2 or below.

set.seed(277)
viz(df_EqVar, legacy=T)

3.3. Results: MRE

JASP v 0.16.3 or later. See van den Bergh, Wagenmakers, and Aust (2023).

set.seed(277)
viz(df_EqVar)

4.1. Real Data: puzzle

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\) or \(a_i\)) and color be factor B (with levels \(\beta_j\) or \(b_j\)) for \(i,j\in\{1,2\}\).

data(puzzles) #library(BayesFactor)
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

summary(aov(RT~shape*color+Error(ID/(shape*color)),puzzles)) #repeated-measures ANOVA

## 
## Error: ID
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 11  226.5   20.59               
## 
## Error: ID:shape
##           Df Sum Sq Mean Sq F value Pr(>F)  
## shape      1   12.0  12.000   7.543  0.019 *
## Residuals 11   17.5   1.591                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: ID:color
##           Df Sum Sq Mean Sq F value  Pr(>F)   
## color      1   12.0  12.000    13.9 0.00334 **
## Residuals 11    9.5   0.864                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: ID:shape:color
##             Df Sum Sq Mean Sq F value Pr(>F)
## shape:color  1    0.0   0.000       0      1
## Residuals   11   30.5   2.773

\[\begin{align} \tag{3} \text{Total}\ \mathit{SS}&=\mathit{SS}_\text{subjects}+\mathit{SS}_{\text{within-subjects}} \\ &=\mathit{SS}_\text{subjects}+[(\mathit{SS}_\text{A}+\mathit{SS}_\text{Error A})+(\mathit{SS}_\text{B}+\mathit{SS}_\text{Error B})+(\mathit{SS}_\text{AB}+\mathit{SS}_\text{Error AB})] \end{align}\]

\[\begin{align} \tag{4} \text{Total}\ \mathit{df}&=\mathit{df}_\text{subjects}+[(\mathit{df}_\text{A}+\mathit{df}_\text{Error A})+(\mathit{df}_\text{B}+\mathit{df}_\text{Error B})+(\mathit{df}_\text{AB}+\mathit{df}_\text{Error AB})]\\ &=(n-1)+[((a-1)+(a-1)(n-1))+((b-1)+(b-1)(n-1))+((a-1)(b-1)+(a-1)(b-1)(n-1))] \end{align}\]

4.2. Results: RIO

JASP v 0.16.2 or below.

set.seed(277)
viz(puzzles, legacy=T)

4.3. Results: MRE

JASP v 0.16.3 or later. See van den Bergh, Wagenmakers, and Aust (2023).

set.seed(277)
viz(puzzles)

5.1. Simulated Unaggregated Data

For \(i=1,\dotsb,N;\quad j,k\in\{1,2\};\quad r=1,\dotsb,n\), \(\qquad\qquad\)(note that the indicators differ from Eq. 1 & 2)

\[\begin{align*} \text{Level 1:}\hspace{3em} Y_{ijkr}&\sim\mathcal{N}(M_{ijk},\ 1)\tag{5} \\ M_{ijk}&=\mu_i+b_{i1}\cdot C_1+b_{i2}\cdot C_2\tag{6} \\ \text{Level 2:}\hspace{3.95em}\mu_i&\sim\mathcal{N}(0,\ \sigma^2)\tag{7} \\ b_{i1}&\sim\mathcal{N}(B_1,\ \sigma^2)\tag{8} \\ b_{i2}&\sim\mathcal{N}(B_2,\ \sigma^2)\tag{9} \end{align*}\]

\(\mu_i,\ b_{i1},\ b_{i2}:\hspace{0.5em}\) subject-level effects

\(B_1, B_2:\hspace{2.1em}\) population-level effects

\(\sigma:\hspace{4.3em}\) between-subjects standard deviation of effects

\(C_1, C_2:\hspace{2.2em}\) contrast-coded independent variables and interaction \(C_1C_2\)

\(\overline{Y}_{ijk\cdot}:\hspace{2.8em}\) aggregated data at the level of the design cell

Citation: Oberauer, K. (2022). The importance of random slopes in mixed models for Bayesian hypothesis testing. Psychological Science, 33, 648-665. https://doi.org/10.1177/09567976211046884

N <- 20; n <- 30; C1 <- C2 <- c(-.5, .5); B1 <- 0; B2 <- .5; sd <- .5

set.seed(277)
s <- rnorm(N, 0, sd)   #random intercepts
b1 <- rnorm(N, B1, sd) #random slopes of factor A
b2 <- rnorm(N, B2, sd) #random slopes of factor B
M <- s + cbind(b1, b2) %*% rbind(rep(C1, 2), rep(C2, each=2)) #cell means
df <- data.frame("DV"=rnorm(N*4*n, c(M), 1),
                 "ID"=rep(paste0("s",1:N), 4*n),
                 "A"=rep(rep(c("a1","a2"), each=N), 2*n),
                 "B"=rep(rep(c("b1","b2"), each=N*2), n),
                 "Trial"=rep(paste0("r",1:n), each=N*4),
                 stringsAsFactors=T)
df_aggreg <- aggregate(DV ~ A + B + ID, data=df, FUN=mean) #aggregate across trials
df_aggreg <- df_aggreg[,c(4,3,1,2)]

dat <- matrix(df_aggreg$DV, ncol=4, byrow=T,
              dimnames=list(unique(df_aggreg$ID),
                            c("a1b1","a2b1","a1b2","a2b2"))) #wide data format
cov(dat) %>% #sample covariance matrix of the aggregated data
  kable(digits=3) %>% kable_styling(full_width=F)

	a1b1	a2b1	a1b2	a2b2
a1b1	0.588	0.264	0.395	0.127
a2b1	0.264	0.325	0.211	0.235
a1b2	0.395	0.211	0.520	0.323
a2b2	0.127	0.235	0.323	0.434

summary(aov(DV~A*B+Error(ID/(A*B)),df)) #repeated-measures ANOVA

## 
## Error: ID
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 19  709.3   37.33               
## 
## Error: ID:A
##           Df Sum Sq Mean Sq F value Pr(>F)
## A          1   5.93   5.930   0.619  0.441
## Residuals 19 181.94   9.576               
## 
## Error: ID:B
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## B          1  174.7   174.7   21.04 0.000201 ***
## Residuals 19  157.7     8.3                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: ID:A:B
##           Df Sum Sq Mean Sq F value Pr(>F)
## A:B        1  0.263  0.2633   0.327  0.574
## Residuals 19 15.287  0.8046               
## 
## Error: Within
##             Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 2320   2190  0.9441

\[\begin{align} \tag{3} \text{Total}\ \mathit{SS}&=\mathit{SS}_\text{subjects}+\mathit{SS}_{\text{within-subjects}} \\ &=\mathit{SS}_\text{subjects}+[(\mathit{SS}_\text{A}+\mathit{SS}_\text{Error A})+(\mathit{SS}_\text{B}+\mathit{SS}_\text{Error B})+(\mathit{SS}_\text{AB}+\mathit{SS}_\text{Error AB})] \end{align}\]

\[\begin{align} \tag{4} \text{Total}\ \mathit{df}&=\mathit{df}_\text{subjects}+[(\mathit{df}_\text{A}+\mathit{df}_\text{Error A})+(\mathit{df}_\text{B}+\mathit{df}_\text{Error B})+(\mathit{df}_\text{AB}+\mathit{df}_\text{Error AB})]\\ &=(n-1)+[((a-1)+(a-1)(n-1))+((b-1)+(b-1)(n-1))+((a-1)(b-1)+(a-1)(b-1)(n-1))] \end{align}\]

5.2. Results: RIO

JASP v 0.16.2 or below.

set.seed(277)
viz(df, legacy=T)

5.3. Results: MRE

JASP v 0.16.3 or later. See van den Bergh, Wagenmakers, and Aust (2023).

set.seed(277)
viz(df)

Distributions of Random Effects

Random Intercept-Only (RIO) and Maximal Random Effects (MRE) Models

Zhengxiao Wei, Farouk S. Nathoo, Michael E. J. Masson

2023-11-20

1. Getting Started

2.1. The Paradox on Real Data

2.2. Results: RIO

2.3. Results: MRE

3.1. The Paradox on Simulated Data

3.2. Results: RIO

3.3. Results: MRE

4.1. Real Data: puzzle

4.2. Results: RIO

4.3. Results: MRE

5.1. Simulated Unaggregated Data

5.2. Results: RIO

5.3. Results: MRE