# Using brm to reproduce analysis of Liddell and Kruschke (2018)

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 )
``````
1 Like

Some corrections/clarifications. This is the actual formulation that I used for unequal variances.

``````fit_uneq <- brm(
formula = bf(resp ~ 1 + cond + (1|subID) + (1|item)) +
lf(disc ~ 0 + cond + (1|subID) + (1|item), cmc = FALSE),
data = df,
family = cumulative("probit"),
control = list(adapt_delta = 0.99, max_treedepth = 12),
iter = 2000, chains = 4, warmup = 1000,
cores = 4,
backend = "cmdstanr",
save_pars = save_pars(all = TRUE)
)
``````

The estimate that I get for the condition term is 0.65 [-0.05, 1.47], which is still far from 0 (with 3 divergent transitions out of 4000). Note that when I do the analysis with data containing a single item (in this case, the random effects do not come into place), things work quite well.

The brms model doesn’t match the data simulation.

This is how the simulation works: The latent variables are generated independently for the two conditions with

``````d = mvrnorm( n=nS , mu=mQ , Sigma=sigmaCov , empirical=TRUE )
``````

where `nS` is the number of subjects, `mQ` are the item means and `sigmaCov` is the item covariance matrix. `sigmaCov` is constructed so that the variances are unequal across conditions but the correlations are the same.

So in this simulation: (a) there are no item-level effects (as there are no correlations between subjects for the same item); (b) the covariances between items are unequal across conditions (since `Cov{item1, item2 | cond} = rho * sigma_cond^2)`; (c) the subject indexing starts at 1 for each condition.

Putting this together, I think the correct model specification is:

``````fit_uneq <- brm(
formula =
bf(resp ~ 1 + cond + (1 + cond || cond:subID)) +
lf(disc ~ 0 + cond, cmc = FALSE),
data = mutate(df, cond = factor(cond)), # condition is categorical, not numerical
family = cumulative(probit)
)
``````

There is no `(1 | item)` term because there are no item-level effects. The subject-level term is a bit complex. From right to left, `cond:subID` means that subject 1 in condition 1 is distinct from subject 1 in condition 2. (Or just assign unique subject IDs.) The `||` means there are no correlations between conditions at the subject level. And `1 + cond` means that each condition has its own std. dev parameter at the subject level. The last bit was hardest to figure out; in the simulation the correlations are the same but the covariances (between items for the same subject) are higher in the second condition.

Here is the summary of the fitted model:

``````#> Multilevel Hyperparameters:
#> ~cond:subID (Number of levels: 1000)
#>               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept)     2.02      0.11     1.80     2.24 1.00     1106     1741
#> sd(cond2)         7.69      0.58     6.60     8.89 1.00     1789     2828
#>
#> Regression Coefficients:
#>              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept[1]    -6.00      0.30    -6.60    -5.45 1.00     1785     2739
#> Intercept[2]    -3.77      0.18    -4.13    -3.43 1.00     1701     2282
#> Intercept[3]    -1.54      0.12    -1.77    -1.31 1.00     1461     2618
#> Intercept[4]     0.75      0.11     0.54     0.96 1.00     1253     2515
#> cond2           -0.06      0.42    -0.86     0.77 1.00     1285     2306
#> disc_cond2      -1.43      0.06    -1.56    -1.31 1.00     1919     2561
``````

Thanks a lot! I really appreciate your effort.

To demonstrate the authors’ point more clearly and showcase the value of probit models, I really need to play with a different data generation process that can simulate more realistic datasets, avoiding complex formulas.