I used the code of Liddell and Kruschke (2018) to reproduce their analysis for Fig.10-Left, where multiple Likert item responses are generated from distributions with equal means (but different standard deviations). The generated data have the following form:

resp item subID cond

4 1 1 1

3 2 1 1

4 3 1 1

…

- item = 1, 2, 3
- subjID = 1…500
- cond = 1, 2

According to the authors, the metric model applied to the average of the three items detects a difference d = 0.53 (Type I error). However, the authors also mention that their ordered-profit model finds no effect (the 95% CI includes 0 and d = 0). I’m trying to use a *brm* model to confirm this results but without success. Could someone help?

The model I fit assuming equal variances is the following:

```
fit_eq <- brm(
formula = resp ~ 1 + cond + (1|subID) + (1|item),
data = df,
family = cumulative("probit")
)
```

However, since variances between conditions are not equal, I also tried the following one:

```
fit_uneq <- brm(
formula = bf(resp ~ 1 + cond + (1|subID) + (1|item)) +
lf(disc ~ 0 + cond, cmc = FALSE),
data = df,
family = cumulative("probit")
)
```

Is my formulation correct?

Both models return a 95% CI (for the difference between conditions) that clearly does not contain zero. Consider that I have warnings about divergent transitions (even with adapt_delta=0.97) but I’m not an expert and don’t fully understand how to diagnose such problems.

Below, I attach the code used to generate the data (extracted from here where simIdx = "3QeqM”):

```
genGroupData <- function(
# specify means on underlying scales:
mQ = c( 1.0 , 1.5 , 2.5 , 3.0 ) , # length is number of questions
# specify standard deviations of responses on underlying scales:
sdQ = c( 1.5 , 2.0 , 2.5 , 3.0 ) , # length is number of questions
# specify correlations of responses on underlying scales
rQval = 0.9 ,
# specify thresholds:
thresh = matrix( c( 1.5 , 2.5 , 3.5 , 4.5 , 5.5 , 6.5 ,
1.5 , 2.5 , 3.5 , 4.5 , 5.5 , 6.5 ,
1.5 , 2.5 , 3.5 , 4.5 , 5.5 , 6.5 ,
1.5 , 2.5 , 3.5 , 4.5 , 5.5 , 6.5 ),
nrow=length(mQ) , byrow=TRUE ) ,
# specify number of subjects in the group:
nS = 100 ,
showSummary = FALSE
) { # generate data...
require(MASS)
nQ = length(mQ)
rQ = matrix( rQval , nrow=nQ , ncol=nQ )
diag(rQ) = 1.0
# construct covariance matrix:
sigmaCov = ( sdQ %o% sdQ ) * rQ
# construct latent data:
d = mvrnorm( n=nS , mu=mQ , Sigma=sigmaCov , empirical=TRUE )
# compute response data:
y = 0*d
for ( sIdx in 1:nrow(d) ) {
for ( qIdx in 1:ncol(d) ) {
y[sIdx,qIdx] = max( which( d[sIdx,qIdx] > c(-Inf,thresh[qIdx,]) ) )
}
}
# show data summary on console:
if ( showSummary ) {
apply( d , 2 , mean )
var( d )
# show response data on console:
apply( y , 2 , mean )
var( y )
cor( y )
}
return(y)
}
NperGroup = 500 # 500
rQvalGlobal = 0.8 # 0.8
# specify thresholds for all groups:
nQ=3 # number of questions (a.k.a. items)
thresh = matrix( rep( seq( 1.5 , 4.5 , 1 ) , nQ ), nrow=nQ , byrow=TRUE ) # thresh's of each Q
g1m = 4.2 # group 1 mean
g1sd = 1.0 # group 1 sd
g2m = 4.2
g2sd = 4.0
#fileNameRoot=paste0("OrdinalScaleGroupJags-Sim-",simIdx)
likertScaleDensMax = 0.7 # for use in histogram later
dataScaleDensMax = 0.7 # for use in histogram later
# Now generate the data
set.seed(47405)
y1 = genGroupData(
mQ = rep(0.0,nQ)+g1m , # mean of each Q
sdQ = rep(1.0,nQ)*g1sd , # sd of each Q
rQval = rQvalGlobal , # correlation of every pair of Q's
thresh = thresh , # thresh's of each Q
nS = NperGroup # number of subjects in this group
)
set.seed(47405)
y2 = genGroupData(
mQ = rep(0.0,nQ)+g2m , # mean of each Q
sdQ = rep(1.0,nQ)*g2sd , # sd of each Q
rQval = rQvalGlobal , # correlation of every pair of Q's
thresh = thresh , # thresh's of each Q
nS = NperGroup # number of subjects in this group
)
# Build a data frame in long format:
resp = item = subID = cond = rep( 0 , prod(dim(y1))+prod(dim(y2)) )
dfIdx = 1
for ( sIdx in 1:nrow(y1) ) { for ( qIdx in 1:ncol(y1) ) {
resp[dfIdx] = y1[sIdx,qIdx]
item[dfIdx] = qIdx
subID[dfIdx] = sIdx
cond[dfIdx] = 1
dfIdx = dfIdx+1
} }
for ( sIdx in 1:nrow(y2) ) { for ( qIdx in 1:ncol(y2) ) {
resp[dfIdx] = y2[sIdx,qIdx]
item[dfIdx] = qIdx
subID[dfIdx] = sIdx
cond[dfIdx] = 2
dfIdx = dfIdx+1
} }
df = data.frame( resp=resp , item=item , subID=subID , cond=cond )
```