The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

Fitting Bayesian Multilevel Single Case models using bmscstan

Michele Scandola

2022-09-04

The package bmscstan provides useful functions to fit Bayesian Multilevel Single Case models (BMSC) using as backend Stan (Carpenter et al. 2017).

This approach is based on the seminal approach of the Crawford’s tests (Crawford and Howell 1998; Crawford and Garthwaite 2005; Crawford et al. 2010), using a small control sample of individuals, to see whether the performance of the single case deviates from them. Unfortunately, Crawford’s tests are limited to a number of specific experimental designs that do not allow researchers to use complex experimental designs.

The BMSC approach is born mainly to deal with this problem: its purpose is, in fact, to allow the fitting of models with the same flexibility of a Multilevel Model, with single case and controls data.

The core function of the bmscstan package is BMSC, whose theoretical assumptions, and its validation are reported in (Scandola and Romano 2021).

The syntax used by the BMSC function is extremely similar to the syntax used in the lme4 package. However, the specification of random effects is limited, but it can cover the greatest part of cases (for further details, please see ?bmscstan::randomeffects).

Example on real data

In order to show an example on the use of the bmscstan package, the datasets in this package will be used.

In these datasets we have data coming from a Body Sidedness Effect paradigm (Ottoboni et al. 2005; Tessari et al. 2012), that is a Simon-like paradigm useful to measure body representation.

In this experimental paradigm, participants have to answer to a circle showed in the centre of the computer screen, superimposed to an irrelevant image of a left or right hand, or to a left or right foot.

The circle can be of two colors (e.g. red or blue), and participants have to press one button with the left when the circle is of a specific colour, and with the right hand when the circle is of the another colour.

When the irrelevant background image (foot or hand) is incongruent with the hand used to answer, the reaction times and frequency of errors are higher.

The two irrelevant backgrounds are administered in different experimental blocks.

This is considered an effect of the body representation.

In the package there are two datasets, one composed by 16 healthy control participants, and the other one by an individual affected by right unilateral brachial plexus lesion (however, s/he could independently press the keyboard buttons).

Explore the data

The datasets are called data.pt for the single case, and data.ctrl for the control group, and they can be loaded using data(BSE).

In these datasets there are the Reaction Times RT, a Body.District factor with levels FOOT and HAND, a Congruency factor (levels: Congruent, Incongruent), and a Side factor (levels: Left, Right). In the data.ctrl dataset there also is an ID factor, representing the different 16 control participants.

library(ggplot2)
library(bmscstan)

data(BSE)

str(data.pt)
#> 'data.frame':    467 obs. of  4 variables:
#>  $ RT           : int  562 424 411 491 439 593 504 483 467 413 ...
#>  $ Body.District: Factor w/ 2 levels "FOOT","HAND": 1 1 1 1 1 1 1 1 1 1 ...
#>  $ Congruency   : Factor w/ 2 levels "Congruent","Incongruent": 1 2 2 1 1 2 2 1 1 2 ...
#>  $ Side         : Factor w/ 2 levels "Left","Right": 1 2 1 2 1 1 2 1 2 2 ...

str(data.ctrl)
#> 'data.frame':    4049 obs. of  5 variables:
#>  $ RT           : int  785 641 938 841 486 425 408 394 611 387 ...
#>  $ Body.District: Factor w/ 2 levels "FOOT","HAND": 1 1 1 1 1 1 1 1 1 1 ...
#>  $ Congruency   : Factor w/ 2 levels "Congruent","Incongruent": 2 2 2 2 2 2 1 1 1 1 ...
#>  $ Side         : Factor w/ 2 levels "Left","Right": 1 1 1 1 2 1 1 1 2 2 ...
#>  $ ID           : Factor w/ 16 levels "HN_017","HN_019",..: 1 1 1 1 1 1 1 1 1 1 ...

ggplot(data.pt, aes(y = RT, x = Body.District:Side , fill = Congruency))+
  geom_boxplot()


ggplot(data.ctrl, aes(y = RT, x = Body.District:Side , fill = Congruency))+
  geom_boxplot()+
  facet_wrap( ~ ID , ncol = 4)

These data seem to have some outliers. Let see if they are normally distributed.

qqnorm(data.ctrl$RT, main = "Controls")
qqline(data.ctrl$RT)

qqnorm(data.pt$RT, main = "Single Case")
qqline(data.pt$RT)

They are not normally distributed. Outliers will be removed by using the boxplot.stats function.

out <- boxplot.stats( data.ctrl$RT )$out
data.ctrl <- droplevels( data.ctrl[ !data.ctrl$RT %in% out , ] )

out <- boxplot.stats( data.pt$RT )$out
data.pt <- droplevels( data.pt[ !data.pt$RT %in% out , ] )

qqnorm(data.ctrl$RT, main = "Controls")
qqline(data.ctrl$RT)

qqnorm(data.pt$RT, main = "Single Case")
qqline(data.pt$RT)

They are not perfect, but definitively better.

Deciding the contrasts and the random effects part

First of all, there is the necessity to think to our hypotheses, and setting the contrast matrices consequently.

In all cases, our factors have only two levels. Therefore, we set the factors with a Treatment Contrasts matrix, with baseline level for Side the Left level, for Congruency the Congruent level, and for Body.District the FOOT level.

In this way, each coefficient will represent the difference between the two levels.

contrasts( data.ctrl$Side )          <- contr.treatment( n = 2 )
contrasts( data.ctrl$Congruency )    <- contr.treatment( n = 2 )
contrasts( data.ctrl$Body.District ) <- contr.treatment( n = 2 )

contrasts( data.pt$Side )            <- contr.treatment( n = 2 )
contrasts( data.pt$Congruency )      <- contr.treatment( n = 2 )
contrasts( data.pt$Body.District )   <- contr.treatment( n = 2 )

The use of the BMSC function, for those who are used to lme4 or brms syntax should be straightforward.

In this case, we want to fit the following model:

RT ~ Body.District * Congruency * Side + (Congruency * Side | ID : Body.District)

Unfortunately, BMSC does not directly allow the syntax ID : Body.District in the specification of the random effects.

Therefore, it is necessary to create a new variable for ID : Body.District

data.ctrl$BD_ID <- interaction( data.ctrl$Body.District , data.ctrl$ID )

and the model would be:

RT ~ Body.District * Congruency * Side + (Congruency * Side | BD_ID)

For further details concerning the random effects available in bmscstan, please type ?bmscstan::randomeffect.

Fitting the BMSC model

At this point, fitting the model is easy, and it can be done with the use of a single function.

mdl <- BMSC(formula = RT ~ Body.District * Congruency * Side +
             (Congruency * Side | BD_ID),
             data_ctrl = data.ctrl,
             data_sc = data.pt,
             chains = 2,
             cores = 1,
             seed = 2020)
#> 
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.001497 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 14.97 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
#> Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
#> Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
#> Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
#> Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
#> Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
#> Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 190.083 seconds (Warm-up)
#> Chain 1:                178.987 seconds (Sampling)
#> Chain 1:                369.07 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 0.000939 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 9.39 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2: 
#> Chain 2: 
#> Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
#> Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
#> Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
#> Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
#> Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
#> Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
#> Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
#> Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
#> Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
#> Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
#> Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
#> Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
#> Chain 2: 
#> Chain 2:  Elapsed Time: 188.854 seconds (Warm-up)
#> Chain 2:                179.24 seconds (Sampling)
#> Chain 2:                368.094 seconds (Total)
#> Chain 2:

After fitting the model, we should check its quality by means of Posterior Predictive P-Values (Gelman 2013) with the bmscstan::pp_check function.

Thanks to this graphical function, we will see if the observed data and the data sampled from the posterior distributions of our BMSC model are similar.

If we observe strong deviations, it means that your BMSC model is not adequately describing your data. In this case, you might want to change the priors distribution (see the help page), change the random effects structure, or transform your dependent variable (using the logarithm or other math functions).

pp_check( mdl )

#> TableGrob (2 x 1) "arrange": 2 grobs
#>   z     cells    name           grob
#> 1 1 (1-1,1-1) arrange gtable[layout]
#> 2 2 (2-2,1-1) arrange gtable[layout]

In both the controls and the single case data, the Posterior Predictive P-Values check seems to adequately resemble the observed data.

A further control on our model is given by checking the Effective Sample Size (ESS) for each coefficient and the \(\hat{R}\) diagnostic index (Gelman and Rubin 1992).

The ESS is the “effective number of simulation draws” for any coefficient, namely the approximate number of independent draws, taking into account that the various simulations in a Monte Carlo Markov Chain (MCMC) are not independent each other. For further details, see an introductory book in Bayesian Statistics. A good ESS estimates should be \(ESS > 100\) or \(ESS > 10\%\) of the total draws (remembering that you should remove the burn-in simulations from the total iterations counting).

The \(\hat{R}\) is an index of the convergence of the MCMCs. In BMSC the default is 4. Usually, MCMCs are considered convergent when \(\hat{R} < 1.1\) (Stan default).

In order to check these values, the summary.BMSC function is needed (see next section).

The summary.BMSC output

The output of the brmscstan::summary.BMSC function is divided in four main parts:

  1. In the first part, the model and the selected priors are recalled.
  2. In the second part, the coefficients of the fixed effects for the control group are shown.
  3. In the third part, the coefficients of the fixed effects for the single case are shown.
  4. In the fourth and last part, the fixed effects coefficients for the difference between the single case and the control group are shown.
print( sum_mdl <- summary( mdl ) , digits = 3 )
#> 
#> Bayesian Multilevel Single Case model
#> 
#> RT ~ Body.District * Congruency * Side + (Congruency * Side | 
#>     BD_ID)
#> 
#> [1] "Priors for the regression coefficients: normal distribution; Dispersion parameter (scale or sigma): 10"
#> 
#> 
#>   Fixed Effects for the Control Group
#> 
#>                                    mean se_mean   sd   2.5%     25%    50%
#> (Intercept)                      202.25  0.1383 7.82 186.80 196.970 202.24
#> Body.District2                    24.21  0.1542 8.06   8.31  18.842  24.24
#> Congruency2                       16.55  0.1646 6.95   3.82  11.658  16.33
#> Side2                             19.90  0.1737 7.62   5.27  14.702  19.76
#> Body.District2:Congruency2        -9.28  0.0938 5.35 -19.65 -12.864  -9.20
#> Body.District2:Side2              -4.43  0.1034 5.31 -15.17  -7.930  -4.47
#> Congruency2:Side2                 -5.63  0.1440 7.49 -20.68 -10.574  -5.62
#> Body.District2:Congruency2:Side2   4.01  0.1325 6.99 -10.28  -0.645   4.01
#>                                      75%  97.5% n_eff Rhat     BF10
#> (Intercept)                      207.607 217.47  3195    1 4.21e+29
#> Body.District2                    29.657  40.32  2733    1     63.4
#> Congruency2                       21.151  30.82  1783    1     18.6
#> Side2                             24.979  35.16  1923    1       35
#> Body.District2:Congruency2        -5.658   1.02  3253    1     2.59
#> Body.District2:Side2              -0.831   6.03  2635    1     0.75
#> Congruency2:Side2                 -0.583   9.05  2708    1    0.975
#> Body.District2:Congruency2:Side2   8.703  17.68  2784    1    0.826
#> 
#>        mean se_mean    sd 2.5%  25%  50%  75% 97.5% n_eff Rhat
#> sigmaC 75.8  0.0108 0.874 74.1 75.2 75.8 76.4  77.6  6548    1
#> 
#> 
#>   Fixed Effects for the Patient
#> 
#>                                   mean se_mean    sd  2.5%   25%   50%    75%
#> (Intercept)                      394.3   0.103  6.52 381.3 390.0 394.5 398.77
#> Body.District2                    42.8   0.114  7.19  28.9  37.9  42.5  47.70
#> Congruency2                       46.6   0.116  7.31  31.9  41.7  46.7  51.44
#> Side2                             50.7   0.114  7.23  36.9  45.6  50.7  55.42
#> Body.District2:Congruency2       -27.1   0.137  8.68 -44.0 -33.0 -27.1 -21.52
#> Body.District2:Side2             -22.3   0.132  8.33 -38.4 -27.8 -22.2 -16.63
#> Congruency2:Side2                -12.6   0.139  8.80 -30.1 -18.8 -12.4  -6.65
#> Body.District2:Congruency2:Side2 -12.5   0.160 10.09 -32.4 -19.5 -12.4  -5.58
#>                                   97.5%     BF10
#> (Intercept)                      406.67 4.92e+56
#> Body.District2                    57.19    56187
#> Congruency2                       60.76  1103949
#> Side2                             65.29  5127730
#> Body.District2:Congruency2        -9.93     52.2
#> Body.District2:Side2              -6.30     31.6
#> Congruency2:Side2                  4.28     2.58
#> Body.District2:Congruency2:Side2   7.04        2
#> 
#> 
#>   Fixed Effects for the difference between the Patient and the Control Group
#> 
#>                                    mean se_mean   sd   2.5%   25%    50%    75%
#> (Intercept)                      192.09   0.136 7.77 176.93 186.9 192.07 197.44
#> Body.District2                    18.56   0.137 8.01   3.08  13.0  18.60  23.81
#> Congruency2                       30.03   0.134 7.55  15.06  25.0  30.10  34.98
#> Side2                             30.77   0.133 7.71  15.41  25.6  30.77  35.77
#> Body.District2:Congruency2       -17.87   0.117 7.90 -33.26 -23.1 -17.97 -12.49
#> Body.District2:Side2             -17.83   0.108 7.59 -32.41 -23.0 -17.86 -12.77
#> Congruency2:Side2                 -7.02   0.128 8.20 -23.03 -12.6  -7.03  -1.61
#> Body.District2:Congruency2:Side2 -16.47   0.125 8.45 -32.95 -22.2 -16.56 -10.90
#>                                    97.5% n_eff Rhat     BF10
#> (Intercept)                      207.115  3283    1 4.58e+27
#> Body.District2                    34.365  3397    1     13.4
#> Congruency2                       44.777  3190    1      347
#> Side2                             46.207  3370    1      849
#> Body.District2:Congruency2        -2.687  4550    1     11.7
#> Body.District2:Side2              -2.639  4974    1     9.96
#> Congruency2:Side2                  9.080  4117    1     1.17
#> Body.District2:Congruency2:Side2   0.374  4564    1      5.6
#> 
#>        mean se_mean   sd 2.5%  25%  50% 75% 97.5% n_eff Rhat
#> sigmaP 64.3  0.0381 2.63 59.3 62.5 64.2  66  69.7  4774    1

In the second and fourth part of the output, we can observe a descriptive summary reporting the mean, the standard error, the standard deviation, the \(2.5\%\), \(25%\), \(50\%\), \(75\%\) and \(97.5\%\) of the posterior distributions of each coefficient. If we want the \(95\%\) Credible Interval, we can consider only the \(2.5\%\) and \(97.5\%\) extremes. Then, two diagnostic indexes are reported: the n_eff parameter, that is the ESS, and the Rhat (\(\hat{R}\)). Finally, the Savage-Dickey Bayes Factor is reported (BF10).

In the third part the diagnostic indexes are not reported because these coefficients are computed as marginal probabilities from the probabilities summarized in the second and fourth part.

Understanding the summary.BMSC output

Checking the diagnostic indexes

The first step should be controlling the diagnostic indexes.

In this model, all n_eff are greater than the \(10\%\) of the total iterations (default iterations: 4000, default warmup iterations: 2000, default chains: 4 = 800). Also, all \(\hat{R}s < 1.1\). Finally, we already saw that the Posterior Predictive P-values are showing that the model is representative of the data.

The Control Group results

Then, observing what the fixed effects of the Control group are showing is important before of seeing the differences with the single case.

In this analysis, there are 5 fixed effects which \(BF_{10}\) is greater than 3 (Raftery 1995).

tmp <- sum_mdl[[1]][sum_mdl[[1]]$BF10 > 3,c("BF10","mean","2.5%","97.5%")]

colnames(tmp) <- c("$BF_{10}$", "$\\mu$", "low $95\\%~CI$", "up $95\\%~CI$")

knitr::kable(
  tmp,
  digits = 3
)
\(BF_{10}\) \(\mu\) low \(95\%~CI\) up \(95\%~CI\)
(Intercept) 4.209343e+29 202.250 186.805 217.474
Body.District2 63.36141 24.205 8.306 40.317
Congruency2 18.62051 16.552 3.825 30.817
Side2 34.96145 19.904 5.274 35.161

We can have a general overview of the coefficients of the model with the plot.BMSC function.

plot( mdl , who = "control" )

The interaction between Body District and Congruency needs a further analysis to better understand the phenomenon. It comes useful the function pairwise.BMSC.

pp <- pairwise.BMSC(mdl = mdl , contrast = "Body.District2:Congruency2" ,
                    who = "control")

print( pp , digits = 3 )
#> 
#> Pairwise Bayesian Multilevel Single Case contrasts of coefficients divided by Body.District2:Congruency2 
#> 
#> 
#> 
#>   Marginal distributions
#> 
#>                  mean se_mean    sd 2.5% 25% 50% 75% 97.5% BF10 (not zero)
#> FOOT Incongruent  219   0.157  9.92  199 212 219 225   239        1.38e+25
#> FOOT Congruent    202   0.124  7.82  187 197 202 208   217        4.21e+29
#> HAND Incongruent  234   0.203 12.82  209 225 234 242   259        5.75e+18
#> HAND Congruent    226   0.166 10.50  206 219 227 233   247        7.27e+19
#> 
#> 
#> 
#>   Table of contrasts
#> 
#>                                       mean se_mean    sd   2.5%    25%    50%
#> FOOT Incongruent - FOOT Congruent    16.55   0.110  6.95   3.82  11.66  16.33
#> FOOT Incongruent - HAND Incongruent -14.93   0.147  9.31 -33.11 -21.16 -15.13
#> FOOT Incongruent - HAND Congruent    -7.65   0.163 10.30 -27.07 -14.76  -7.79
#> FOOT Congruent - HAND Incongruent   -31.48   0.178 11.26 -53.88 -38.89 -31.46
#> FOOT Congruent - HAND Congruent     -24.21   0.128  8.06 -40.32 -29.66 -24.24
#> HAND Incongruent - HAND Congruent     7.28   0.124  7.85  -6.72   1.65   6.84
#>                                         75% 97.5%  BF10
#> FOOT Incongruent - FOOT Congruent    21.151 30.82 18.19
#> FOOT Incongruent - HAND Incongruent  -8.713  3.77  3.27
#> FOOT Incongruent - HAND Congruent    -0.574 12.82  1.38
#> FOOT Congruent - HAND Incongruent   -23.938 -9.61 68.69
#> FOOT Congruent - HAND Congruent     -18.842 -8.31 56.38
#> HAND Incongruent - HAND Congruent    12.445 23.50  1.03

The output of this function is divided in two parts:

It is also possible to plot the results of this function with the use of plot.pairwise.BMSC.

plot( pp )
#> [[1]]

#> 
#> [[2]]

Finally, it is possible to plot marginal posterior distributions for each effects with \(BF_{10} > 3\).

p1 <- pairwise.BMSC(mdl , contrast = "Body.District2" ,  who = "control" )

plot( p1 )[[1]] +
  ggtitle("Body District" , subtitle = "Marginal effects") 

plot( p1 )[[2]] +
  ggtitle("Body District" , subtitle = "Contrasts") 

p2 <- pairwise.BMSC(mdl , contrast = "Congruency2" ,  who = "control" )

plot( p2 )[[1]] +
  ggtitle("Congruency" , subtitle = "Marginal effects")

plot( p2 )[[2]] +
  ggtitle("Congruency" , subtitle = "Contrasts")

p3 <- pairwise.BMSC(mdl , contrast = "Side2" ,  who = "control" )

plot( p3 )[[1]] +
  ggtitle("Side" , subtitle = "Marginal effects")

plot( p3 )[[2]] +
  ggtitle("Side" , subtitle = "Contrasts")

The differences between the Control Group and the Single Case

Finally, the difference between the Control Group and the Single Case is of interest.

A general plot can be obtained in the following way, plotting both the Control Group and the Single Case:

plot( mdl ) +
  theme_bw( base_size = 18 )+
  theme( legend.position = "bottom",
         legend.direction = "horizontal")

or plotting only the difference

plot( mdl ,who = "delta" ) +
  theme_bw( base_size = 18 )

The relevant coefficients are:

tmp <- sum_mdl[[3]][sum_mdl[[3]]$BF10 > 3,c("BF10","mean","2.5%","97.5%")]

colnames(tmp) <- c("$BF_{10}$", "$\\mu$", "low $95\\%~CI$", "up $95\\%~CI$")

knitr::kable(
  tmp,
  digits = 3
)
\(BF_{10}\) \(\mu\) low \(95\%~CI\) up \(95\%~CI\)
(Intercept) 4.575384e+27 192.094 176.935 207.115
Body.District2 13.41511 18.559 3.085 34.365
Congruency2 346.7337 30.034 15.063 44.777
Side2 849.4772 30.769 15.410 46.207
Body.District2:Congruency2 11.69081 -17.865 -33.260 -2.687
Body.District2:Side2 9.958894 -17.834 -32.406 -2.639
Body.District2:Congruency2:Side2 5.596429 -16.468 -32.946 0.374

The Intercept coefficient is showing us that the single case is generally slower than the Control Sample (generally speaking, when you analyse healthy controls against a single case with a specific disease, the single case is slower).

All the main effects can be further analysed by simply looking at their estimates (knowing the contrast matrix and the direction of the estimate you can understand which level is greater than the other), or by means of the pairwise.BMSC function, if you also want marginal effects and automatic plots.

The interactions require the use of the pairwise.BMSC function.

The Body District : Congruency interaction:

p4 <- pairwise.BMSC(mdl , contrast = "Body.District2:Congruency2" ,
                    who = "delta")

print( p4 , digits = 3 )
#> 
#> Pairwise Bayesian Multilevel Single Case contrasts of coefficients divided by Body.District2:Congruency2 
#> 
#> 
#> 
#>   Marginal distributions
#> 
#>                  mean se_mean    sd 2.5% 25% 50% 75% 97.5% BF10 (not zero)
#> FOOT Congruent    192   0.123  7.77  177 187 192 197   207        4.58e+27
#> FOOT Incongruent  222   0.162 10.22  202 215 222 229   242        1.48e+23
#> HAND Incongruent  223   0.212 13.42  196 214 223 232   249        6.68e+19
#> HAND Congruent    211   0.170 10.74  190 203 211 218   231        3.18e+17
#> 
#> 
#> 
#>   Table of contrasts
#> 
#>                                        mean se_mean    sd   2.5%    25%     50%
#> FOOT Congruent - FOOT Incongruent   -30.034   0.119  7.55 -44.78 -34.98 -30.099
#> FOOT Congruent - HAND Incongruent   -30.728   0.189 11.93 -53.81 -38.76 -30.725
#> FOOT Congruent - HAND Congruent     -18.559   0.127  8.01 -34.36 -23.81 -18.599
#> FOOT Incongruent - HAND Incongruent  -0.694   0.163 10.30 -21.20  -7.62  -0.745
#> FOOT Incongruent - HAND Congruent    11.475   0.170 10.74  -9.17   4.27  11.558
#> HAND Incongruent - HAND Congruent    12.169   0.152  9.62  -6.86   5.85  12.116
#>                                        75%  97.5%   BF10
#> FOOT Congruent - FOOT Incongruent   -24.98 -15.06 273.48
#> FOOT Congruent - HAND Incongruent   -22.63  -6.68  29.91
#> FOOT Congruent - HAND Congruent     -13.04  -3.08  13.37
#> FOOT Incongruent - HAND Incongruent   6.32  19.99   1.06
#> FOOT Incongruent - HAND Congruent    18.59  32.80   1.87
#> HAND Incongruent - HAND Congruent    18.50  31.18   2.23

The pairwise.BMSC function shows that in all cases the marginal effects of the RTs where greater than zero, but the differences where present only in the comparison between FOOT Congruent and the other cases.

plot( p4 , type = "interval")
#> [[1]]
#> 
#> [[2]]

plot( p4 , type = "area")
#> [[1]]
#> 
#> [[2]]

plot( p4 , type = "hist")
#> [[1]]
#> 
#> [[2]]

In this case we can observe that the single case was more facilitated by the FOOT Congruent condition than the Control Group.

If the interpretation of the results is difficult, it can be useful look what happens in the Single Case marginal effects.

p5 <- pairwise.BMSC(mdl , contrast = "Body.District2:Congruency2" ,
                    who = "singlecase")

plot( p5 , type = "hist")[[1]]

The Body District : Congruency interaction:

p6 <- pairwise.BMSC(mdl , contrast = "Body.District2:Side2" , who = "delta")

print( p6 , digits = 3 )
#> 
#> Pairwise Bayesian Multilevel Single Case contrasts of coefficients divided by Body.District2:Side2 
#> 
#> 
#> 
#>   Marginal distributions
#> 
#>            mean se_mean    sd 2.5% 25% 50% 75% 97.5% BF10 (not zero)
#> FOOT Left   192   0.123  7.77  177 187 192 197   207        4.58e+27
#> FOOT Right  223   0.165 10.45  202 216 223 230   243        4.86e+20
#> HAND Left   211   0.170 10.74  190 203 211 218   231        3.18e+17
#> HAND Right  224   0.221 13.94  196 214 223 233   251        1.31e+16
#> 
#> 
#> 
#>   Table of contrasts
#> 
#>                            mean se_mean    sd   2.5%    25%     50%    75%
#> FOOT Left - FOOT Right  -30.769   0.122  7.71 -46.21 -35.77 -30.771 -25.57
#> FOOT Left - HAND Left   -18.559   0.127  8.01 -34.36 -23.81 -18.599 -13.04
#> FOOT Left - HAND Right  -31.494   0.192 12.17 -54.94 -39.68 -31.622 -23.22
#> FOOT Right - HAND Left   12.210   0.174 11.01  -9.23   4.96  12.320  19.59
#> FOOT Right - HAND Right  -0.725   0.163 10.34 -21.11  -7.61  -0.613   6.23
#> HAND Left - HAND Right  -12.935   0.155  9.82 -32.48 -19.51 -12.890  -6.32
#>                          97.5%    BF10
#> FOOT Left - FOOT Right  -15.41 1104.54
#> FOOT Left - HAND Left    -3.08   13.37
#> FOOT Left - HAND Right   -8.12   40.24
#> FOOT Right - HAND Left   33.42    2.06
#> FOOT Right - HAND Right  19.36    1.03
#> HAND Left - HAND Right    5.99    2.32

plot( p6 , type = "hist")[[1]] +
  theme_bw( base_size = 18)+
  theme( strip.text.y = element_text( angle = 0 ) )

In this case, we can see that the left - right difference in the single case is always present, with faster RTs in the left foot than in the other cases.

The Body District : Congruency : Side interaction:

p7 <- pairwise.BMSC(mdl ,
                    contrast = "Body.District2:Congruency2:Side2" ,
                    who = "delta")

print( p7 , digits = 3 )
#> 
#> Pairwise Bayesian Multilevel Single Case contrasts of coefficients divided by Body.District2:Congruency2:Side2 
#> 
#> 
#> 
#>   Marginal distributions
#> 
#>                        mean se_mean    sd 2.5% 25% 50% 75% 97.5%
#> FOOT Congruent Left     192   0.123  7.77  177 187 192 197   207
#> FOOT Incongruent Right  246   0.215 13.60  219 237 246 255   272
#> FOOT Incongruent Left   222   0.162 10.22  202 215 222 229   242
#> FOOT Congruent Right    223   0.165 10.45  202 216 223 230   243
#> HAND Incongruent Left   223   0.212 13.42  196 214 223 232   249
#> HAND Congruent Right    224   0.221 13.94  196 214 223 233   251
#> HAND Congruent Left     211   0.170 10.74  190 203 211 218   231
#> HAND Incongruent Right  212   0.273 17.28  178 200 212 224   246
#>                        BF10 (not zero)
#> FOOT Congruent Left           4.58e+27
#> FOOT Incongruent Right        2.42e+19
#> FOOT Incongruent Left         1.48e+23
#> FOOT Congruent Right          4.86e+20
#> HAND Incongruent Left         6.68e+19
#> HAND Congruent Right          1.31e+16
#> HAND Congruent Left           3.18e+17
#> HAND Incongruent Right        3.14e+13
#> 
#> 
#> 
#>   Table of contrasts
#> 
#>                                                     mean se_mean    sd   2.5%
#> FOOT Congruent Left - FOOT Incongruent Right    -53.7876   0.193 12.21 -77.68
#> FOOT Congruent Left - FOOT Incongruent Left     -30.0338   0.119  7.55 -44.78
#> FOOT Congruent Left - FOOT Congruent Right      -30.7692   0.122  7.71 -46.21
#> FOOT Congruent Left - HAND Incongruent Left     -30.7277   0.189 11.93 -53.81
#> FOOT Congruent Left - HAND Congruent Right      -31.4940   0.192 12.17 -54.94
#> FOOT Congruent Left - HAND Congruent Left       -18.5590   0.127  8.01 -34.36
#> FOOT Congruent Left - HAND Incongruent Right    -20.1792   0.255 16.12 -52.47
#> FOOT Incongruent Right - FOOT Incongruent Left   23.7537   0.164 10.35   2.84
#> FOOT Incongruent Right - FOOT Congruent Right    23.0183   0.159 10.07   3.18
#> FOOT Incongruent Right - HAND Incongruent Left   23.0599   0.225 14.22  -4.97
#> FOOT Incongruent Right - HAND Congruent Right    22.2936   0.220 13.92  -4.07
#> FOOT Incongruent Right - HAND Congruent Left     35.2286   0.223 14.13   7.47
#> FOOT Incongruent Right - HAND Incongruent Right  33.6083   0.203 12.86   8.86
#> FOOT Incongruent Left - FOOT Congruent Right     -0.7354   0.165 10.44 -21.52
#> FOOT Incongruent Left - HAND Incongruent Left    -0.6939   0.163 10.30 -21.20
#> FOOT Incongruent Left - HAND Congruent Right     -1.4602   0.218 13.76 -28.01
#> FOOT Incongruent Left - HAND Congruent Left      11.4748   0.170 10.74  -9.17
#> FOOT Incongruent Left - HAND Incongruent Right    9.8546   0.243 15.38 -20.54
#> FOOT Congruent Right - HAND Incongruent Left      0.0415   0.215 13.60 -25.85
#> FOOT Congruent Right - HAND Congruent Right      -0.7248   0.163 10.34 -21.11
#> FOOT Congruent Right - HAND Congruent Left       12.2102   0.174 11.01  -9.23
#> FOOT Congruent Right - HAND Incongruent Right    10.5900   0.239 15.13 -18.95
#> HAND Incongruent Left - HAND Congruent Right     -0.7663   0.208 13.16 -26.76
#> HAND Incongruent Left - HAND Congruent Left      12.1687   0.152  9.62  -6.86
#> HAND Incongruent Left - HAND Incongruent Right   10.5485   0.202 12.78 -15.08
#> HAND Congruent Right - HAND Congruent Left       12.9350   0.155  9.82  -5.99
#> HAND Congruent Right - HAND Incongruent Right    11.3148   0.196 12.38 -13.27
#> HAND Congruent Left - HAND Incongruent Right     -1.6202   0.232 14.69 -31.20
#>                                                     25%     50%    75%  97.5%
#> FOOT Congruent Left - FOOT Incongruent Right    -61.927 -53.675 -45.63 -29.69
#> FOOT Congruent Left - FOOT Incongruent Left     -34.984 -30.099 -24.98 -15.06
#> FOOT Congruent Left - FOOT Congruent Right      -35.772 -30.771 -25.57 -15.41
#> FOOT Congruent Left - HAND Incongruent Left     -38.764 -30.725 -22.63  -6.68
#> FOOT Congruent Left - HAND Congruent Right      -39.684 -31.622 -23.22  -8.12
#> FOOT Congruent Left - HAND Congruent Left       -23.813 -18.599 -13.04  -3.08
#> FOOT Congruent Left - HAND Incongruent Right    -30.906 -20.355  -9.40  11.92
#> FOOT Incongruent Right - FOOT Incongruent Left   16.797  24.051  30.59  44.02
#> FOOT Incongruent Right - FOOT Congruent Right    16.473  23.039  29.78  42.58
#> FOOT Incongruent Right - HAND Incongruent Left   13.422  23.041  32.91  51.05
#> FOOT Incongruent Right - HAND Congruent Right    12.954  21.825  31.63  50.18
#> FOOT Incongruent Right - HAND Congruent Left     25.952  34.950  44.59  63.34
#> FOOT Incongruent Right - HAND Incongruent Right  24.735  33.725  42.30  58.79
#> FOOT Incongruent Left - FOOT Congruent Right     -7.756  -0.800   6.22  19.60
#> FOOT Incongruent Left - HAND Incongruent Left    -7.622  -0.745   6.32  19.99
#> FOOT Incongruent Left - HAND Congruent Right    -10.869  -1.647   8.03  25.65
#> FOOT Incongruent Left - HAND Congruent Left       4.266  11.558  18.59  32.80
#> FOOT Incongruent Left - HAND Incongruent Right   -0.571   9.806  20.22  39.64
#> FOOT Congruent Right - HAND Incongruent Left     -9.159   0.189   8.94  27.30
#> FOOT Congruent Right - HAND Congruent Right      -7.613  -0.613   6.23  19.36
#> FOOT Congruent Right - HAND Congruent Left        4.961  12.320  19.59  33.42
#> FOOT Congruent Right - HAND Incongruent Right     0.229  10.446  20.46  40.37
#> HAND Incongruent Left - HAND Congruent Right     -9.407  -0.805   7.73  25.79
#> HAND Incongruent Left - HAND Congruent Left       5.852  12.116  18.50  31.18
#> HAND Incongruent Left - HAND Incongruent Right    2.022  10.846  19.09  35.23
#> HAND Congruent Right - HAND Congruent Left        6.319  12.890  19.51  32.48
#> HAND Congruent Right - HAND Incongruent Right     2.858  11.075  19.69  35.80
#> HAND Congruent Left - HAND Incongruent Right    -11.635  -1.425   8.42  26.58
#>                                                    BF10
#> FOOT Congruent Left - FOOT Incongruent Right    3225.05
#> FOOT Congruent Left - FOOT Incongruent Left      273.48
#> FOOT Congruent Left - FOOT Congruent Right      1104.54
#> FOOT Congruent Left - HAND Incongruent Left       29.91
#> FOOT Congruent Left - HAND Congruent Right        40.24
#> FOOT Congruent Left - HAND Congruent Left         13.37
#> FOOT Congruent Left - HAND Incongruent Right       3.49
#> FOOT Incongruent Right - FOOT Incongruent Left    15.09
#> FOOT Incongruent Right - FOOT Congruent Right     13.84
#> FOOT Incongruent Right - HAND Incongruent Left     5.49
#> FOOT Incongruent Right - HAND Congruent Right      4.96
#> FOOT Incongruent Right - HAND Congruent Left      33.14
#> FOOT Incongruent Right - HAND Incongruent Right   48.31
#> FOOT Incongruent Left - FOOT Congruent Right       1.02
#> FOOT Incongruent Left - HAND Incongruent Left      1.06
#> FOOT Incongruent Left - HAND Congruent Right       1.41
#> FOOT Incongruent Left - HAND Congruent Left        1.87
#> FOOT Incongruent Left - HAND Incongruent Right     1.87
#> FOOT Congruent Right - HAND Incongruent Left       1.30
#> FOOT Congruent Right - HAND Congruent Right        1.03
#> FOOT Congruent Right - HAND Congruent Left         2.06
#> FOOT Congruent Right - HAND Incongruent Right      1.87
#> HAND Incongruent Left - HAND Congruent Right       1.31
#> HAND Incongruent Left - HAND Congruent Left        2.23
#> HAND Incongruent Left - HAND Incongruent Right     1.83
#> HAND Congruent Right - HAND Congruent Left         2.35
#> HAND Congruent Right - HAND Incongruent Right      1.87
#> HAND Congruent Left - HAND Incongruent Right       1.49

plot( p7 , type = "hist")[[1]] +
  theme_bw( base_size = 18)+
  theme( strip.text.y = element_text( angle = 0 ) )

Here we can see that the effect was pushed by the facilitation that the single case had in the Left Congruent Foot condition compared to the Control Group.

Efficient approximate leave-one-out cross-validation for fitted BMSC models

The bmscstan package has wrapper functions to interface with the loo package, to diagnostic and compare BMSC models.

Leaving-One-Out scores, diagnostics and comparisons are separately computed for the Control group and the Single Case data.

In order to see the Leaving-One-Out and the Pareto smoothed importance sampling (PSIS), it is possible to use the function loo.BMSC:

print( loo1 <- BMSC_loo( mdl ) )
#> 
#> Leave-One-Out Cross-Validation using PSIS-LOO for the single case
#> 
#> 
#> Computed from 4000 by 449 log-likelihood matrix
#> 
#>          Estimate   SE
#> elpd_loo  -2509.0 15.2
#> p_loo         6.0  0.4
#> looic      5018.0 30.4
#> ------
#> Monte Carlo SE of elpd_loo is 0.0.
#> 
#> All Pareto k estimates are good (k < 0.5).
#> See help('pareto-k-diagnostic') for details.
#> 
#> Leave-One-Out Cross-Validation using PSIS-LOO for the control group
#> 
#> 
#> Computed from 4000 by 3933 log-likelihood matrix
#> 
#>          Estimate   SE
#> elpd_loo -22637.5 49.8
#> p_loo        70.3  1.8
#> looic     45275.0 99.7
#> ------
#> Monte Carlo SE of elpd_loo is 0.1.
#> 
#> All Pareto k estimates are good (k < 0.5).
#> See help('pareto-k-diagnostic') for details.

plot( loo1 )

Model comparison can be done by means of the BMSC_loo_compare function:

mdl.null <- BMSC(formula = RT ~ 1 +
             (Congruency * Side | BD_ID),
             data_ctrl = data.ctrl,
             data_sc = data.pt,
             cores = 1,
             chains = 2,
             seed = 2021)
#> 
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.001076 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 10.76 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
#> Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
#> Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
#> Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
#> Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
#> Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
#> Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 152.371 seconds (Warm-up)
#> Chain 1:                163.833 seconds (Sampling)
#> Chain 1:                316.204 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 0.000684 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 6.84 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2: 
#> Chain 2: 
#> Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
#> Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
#> Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
#> Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
#> Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
#> Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
#> Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
#> Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
#> Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
#> Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
#> Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
#> Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
#> Chain 2: 
#> Chain 2:  Elapsed Time: 176.451 seconds (Warm-up)
#> Chain 2:                168.809 seconds (Sampling)
#> Chain 2:                345.26 seconds (Total)
#> Chain 2:

print( loo2 <- BMSC_loo( mdl.null ) )
#> 
#> Leave-One-Out Cross-Validation using PSIS-LOO for the single case
#> 
#> 
#> Computed from 4000 by 449 log-likelihood matrix
#> 
#>          Estimate   SE
#> elpd_loo  -2484.7 16.4
#> p_loo         1.9  0.2
#> looic      4969.4 32.7
#> ------
#> Monte Carlo SE of elpd_loo is 0.0.
#> 
#> All Pareto k estimates are good (k < 0.5).
#> See help('pareto-k-diagnostic') for details.
#> 
#> Leave-One-Out Cross-Validation using PSIS-LOO for the control group
#> 
#> 
#> Computed from 4000 by 3933 log-likelihood matrix
#> 
#>          Estimate   SE
#> elpd_loo -22636.0 49.7
#> p_loo        57.0  1.5
#> looic     45272.1 99.4
#> ------
#> Monte Carlo SE of elpd_loo is 0.1.
#> 
#> All Pareto k estimates are good (k < 0.5).
#> See help('pareto-k-diagnostic') for details.

plot( loo2 )


BMSC_loo_compare( list( loo1, loo2 ) )
#> 
#> Leave-One-Out Cross-Validation model comparison for the single case
#> 
#>        elpd_diff se_diff
#> model2   0.0       0.0  
#> model1 -24.3       7.9  
#> 
#> Leave-One-Out Cross-Validation model comparison for the control group
#> 
#>        elpd_diff se_diff
#> model2   0.0       0.0  
#> model1 -24.3       7.9

Further details on LOO, PSIS and their use can be found in the loo package and in Vehtari, Gelman, and Gabry (2017) and Vehtari et al. (2015).

A binomial case

In this section, a brief example on how to use the package for binomial data.

We start simulating the data.

######################################
# simulation of controls' group data
######################################

# Number of levels for each condition and trials
NCond  <- 2
Ntrials <- 20
NSubjs  <- 40

betas <- c( 0.5 , 0 )

data.sim <- expand.grid(
  trial      = 1:Ntrials,
  ID         = factor(1:NSubjs),
  Cond      = factor(1:NCond)
)

### d.v. generation
y <- rep( times = nrow(data.sim) , NA )

# cheap simulation of individual random intercepts
set.seed(1)
rsubj <- rnorm(NSubjs , sd = 0.1)

for( i in 1:length( levels( data.sim$ID ) ) ){
  
  sel <- which( data.sim$ID == as.character(i) )
  
  mm  <- model.matrix(~ 1 + Cond , data = data.sim[ sel , ] )
  
  set.seed(1 + i)
  y[sel] <- mm %*% as.matrix(betas + rsubj[i]) +
    rnorm( n = Ntrials * NCond )
  
}

data.sim$y <- y
data.sim$bin <- sapply(
  LaplacesDemon::invlogit(data.sim$y),
  function(x) rbinom( 1, 1, x)
  )

data.sim.bin <- aggregate( bin ~ Cond * ID, data = data.sim, FUN = sum)
data.sim.bin$n <- aggregate( bin ~ Cond * ID,
                             data = data.sim, FUN = length)$bin

######################################
# simulation of patient data
######################################

betas.pt <- c( 0 , 2  )

data.pt <- expand.grid(
  trial      = 1:Ntrials,
  Cond      = factor(1:NCond)
)

### d.v. generation
mm  <- model.matrix(~ 1 + Cond , data = data.pt )

set.seed(5)
data.pt$y <- (mm %*% as.matrix(betas.pt + betas) +
  rnorm( n = Ntrials * NCond ))[,1]
data.pt$bin <- sapply(
  LaplacesDemon::invlogit(data.pt$y),
  function(x) rbinom( 1, 1, x)
  )

data.pt.bin <- aggregate( bin ~ Cond, data = data.pt, FUN = sum)
data.pt.bin$n <- aggregate( bin ~ Cond,
                             data = data.pt, FUN = length)$bin


plot(x = data.sim.bin$Cond, y = data.sim.bin$bin, ylim = c(0,20))
points(x = data.pt.bin$Cond, y = data.pt.bin$bin, col = "red")

The boxplot represents the control participants, the red dot the single case.

Now, we can specify the model:

cbind(bin, n) ~ Cond

The right-hand side of the formula follows the usual lmer- and brms-like syntax. In the left-hand side of the formula, brms and lme4 have divergent notations.

In future, the bmscstan package will be able to use both notations, for the moment it is necessary the lme4 notation cbind(bin, n) where:

mdlBin <- BMSC(formula = cbind(bin, n) ~ 1 + Cond,
            data_ctrl = data.sim.bin, data_sc = data.pt.bin, seed = 2022,
            chains = 2,
            family = "binomial", cores = 1)
#> 
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 6.7e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.67 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
#> Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
#> Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
#> Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
#> Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
#> Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
#> Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 1.03 seconds (Warm-up)
#> Chain 1:                0.822 seconds (Sampling)
#> Chain 1:                1.852 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 3.9e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.39 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2: 
#> Chain 2: 
#> Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
#> Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
#> Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
#> Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
#> Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
#> Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
#> Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
#> Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
#> Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
#> Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
#> Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
#> Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
#> Chain 2: 
#> Chain 2:  Elapsed Time: 1.075 seconds (Warm-up)
#> Chain 2:                0.783 seconds (Sampling)
#> Chain 2:                1.858 seconds (Total)
#> Chain 2:

print( summary( mdlBin ) , digits = 3 )
#> 
#> Bayesian Multilevel Single Case model
#> 
#> cbind(bin, n) ~ 1 + Cond
#> 
#> [1] "Priors for the regression coefficients: normal distribution; Dispersion parameter (scale or sigma): 10"
#> 
#> 
#>   Fixed Effects for the Control Group
#> 
#>                mean se_mean     sd   2.5%     25%    50%    75% 97.5% n_eff
#> (Intercept) 0.41315 0.00176 0.0709  0.278  0.3649 0.4106 0.4602 0.556  1617
#> Cond2       0.00683 0.00249 0.0993 -0.191 -0.0622 0.0115 0.0729 0.201  1594
#>             Rhat  BF10
#> (Intercept)    1 11519
#> Cond2          1 0.011
#> 
#> NULL
#> 
#> 
#>   Fixed Effects for the Patient
#> 
#>               mean se_mean    sd   2.5%     25%   50%    75% 97.5%   BF10
#> (Intercept)  0.224 0.00729 0.461 -0.653 -0.0885 0.214  0.524  1.15 0.0484
#> Cond2       10.291 0.09191 5.813  2.762  5.8346 8.987 13.624 24.76    555
#> 
#> 
#>   Fixed Effects for the difference between the Patient and the Control Group
#> 
#>               mean se_mean    sd  2.5%    25%    50%    75%  97.5% n_eff Rhat
#> (Intercept) -0.189  0.0101 0.467 -1.08 -0.504 -0.196  0.119  0.771  2146    1
#> Cond2       10.284  0.1692 5.811  2.69  5.845  8.968 13.619 24.747  1180    1
#>               BF10
#> (Intercept) 0.0514
#> Cond2          699
#> 
#> NULL

Conclusions

In this vignette we have seen how to use the package bmscstan and its functions to analyse and make sense of Single Case data.

The output of the main functions is rich of information, and the Bayesian Inference can be done by taking into account the Savage-Dickey \(BF_{10}\), or the \(95\%\) CI (see Kruschke 2014 for further details).

In this vignette there is almost no discussion concerning how to test the Single Case fixed effects (third part of the main output), but it was used to better understand what happens in the differences between the single case and the control group.

However, if your hypotheses focus on the behaviour of the patient, and not only on the differences between single case and the control group, it will be important to analyse in detail also that part.

References

Carpenter, Bob, Andrew Gelman, Matthew D. Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. 2017. “<i>Stan</i> : A Probabilistic Programming Language.” Journal of Statistical Software 76 (1). https://doi.org/10.18637/jss.v076.i01.
Crawford, John R., and Paul H. Garthwaite. 2005. “Evaluation of criteria for classical dissociations in single-case studies by Monte Carlo simulation.” Neuropsychology 19 (5): 664–78. https://doi.org/10.1037/0894-4105.19.5.664.
Crawford, John R., Paul H Garthwaite, Sara Porter, and Sara Porter Point. 2010. “Point and interval estimates of effect sizes for the case-controls design in neuropsychology : Rationale , methods , implementations , and proposed reporting standards Point and interval estimates of effect sizes for the case-controls design in neuropsych” 3294. https://doi.org/10.1080/02643294.2010.513967.
Crawford, John R., and David C Howell. 1998. “Comparing an Individual ’ s Test Score Against Norms Derived from Small Samples.” The Clinical Neuropsychologist 12 (4): 482–86. https://doi.org/10.1076/clin.12.4.482.7241.
Gelman, Andrew. 2013. “Understanding posterior p-values.” Electronic Journal of Statistics, 1–8. https://doi.org/10.1214/13-EJS854.
Gelman, Andrew, and Donald B Rubin. 1992. “Inference from Iterative Simulation Using Multiple Sequences.” Statistical Science 7 (4): 457–72.
Kruschke, John K. 2014. Doing Bayesian data analysis: A tutorial with R, JAGS, and Stan, second edition. 2nd ed. Elsevier Inc. https://doi.org/10.1016/B978-0-12-405888-0.09999-2.
Ottoboni, Giovanni, Alessia Tessari, Roberto Cubelli, and Carlo Umiltà. 2005. “Is handedness recognition automatic? A study using a Simon-like paradigm.” Journal of Experimental Psychology: Human Perception and Performance 31 (4): 778–89. https://doi.org/10.1037/0096-1523.31.4.778.
Raftery, Adrian E. 1995. “Bayesian Model Selection in Social Research.” In Sociological Methodology, edited by P V Marsden, 25:111–63. 1995. Cambridge: Blackwells.
Scandola, Michele, and Daniele Romano. 2021. “Bayesian multilevel single case models using ‘Stan’. A new tool to study single cases in neuropsychology.” Neuropsychologia 156 (June): 107834. https://doi.org/10.1016/j.neuropsychologia.2021.107834.
Tessari, Alessia, Giovanni Ottoboni, Giulia Baroni, Ed Symes, and Roberto Nicoletti. 2012. “Is access to the body structural description sensitive to a body part’s significance for action and cognition? A study of the sidedness effect using feet.” Experimental Brain Research 218 (4): 515–25. https://doi.org/10.1007/s00221-012-3045-4.
Vehtari, Aki, Andrew Gelman, and Jonah Gabry. 2017. “Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC.” Statistics and Computing 27 (5): 1413–32. https://doi.org/10.1007/s11222-016-9696-4.
Vehtari, Aki, Daniel Simpson, Andrew Gelman, Yuling Yao, and Jonah Gabry. 2015. “Pareto Smoothed Importance Sampling,” no. July (July). https://arxiv.org/abs/1507.02646.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.