Usage examples

  collapse = TRUE,
  comment = "#>"
knitr::opts_chunk$set(warning = FALSE, message = FALSE)

Minimal example

The main function is fit_cumhist() that takes a data frame with time-series as a first argument. In addition, you need to specify the name of the column that codes the perceptual state (state argument) and a column that holds either dominance phase duration (duration) or its onset (onset). The code below fits data using Gamma distribution (default family) for a single run of a single participant. By default, the function fits cumulative history time constant but uses default fixed mixed state value (mixed_state = 0.5), history mixing proportion (history_mix = 0.2), and initial history values (history_init = 0).

gamma_fit <- fit_cumhist(br_singleblock,

Alternatively, you specify onset of individual dominance phases that will be used to compute their duration.

gamma_fit <- fit_cumhist(br_singleblock,

Now you can look at the fitted value for history time constant using history_tau()


and main effect of history for both parameters of gamma distribution


The following model is fitted for the example above, see also companion vignette for details on cumulative history computation. $$Duration[i] \sim Gamma(shape[i], rate[i]) \ log(shape[i]) = \alpha^{shape} + \beta^{shape}H \cdot H{mix}[i] \ log(rate[i]) = \alpha^{rate} + \beta^{rate}H \cdot H{mix}[i] \ H_{mix}[i] = \text{cumulative_history}(\tau,\text{history_mixture}, \text{history_init})\ \alpha^{shape}, \alpha^{rate} \sim Normal(log(3), 5) \ \beta^{shape}_H, \beta^{rate}_H \sim Normal(0, 1) \ \tau \sim Normal(log(1), 0.75)$$


By default, fit_cumhist() function assumes that the time-series represent a single run with the consequence that history states are initialized only once at the very beginning. You can use run argument to pass the name of a column that specifies individual runs. In this case, history is initialize at the beginning of every run to avoid spill-over effects.

gamma_fit <- fit_cumhist(br_single_subject,

Experimental session

Experimental session specifies which time-series were measured together and is used to compute an average dominance phase duration that, in turn, is used when computing cumulative history: $\tau_H = \tau $, where $\tau$ is normalized time constant and $$ is the mean dominance phase duration. This can be used to account for changes in overall alternation rate between different sessions (days), as, for example, participants new to the stimuli tend to "speed up" over the course of days Suzuki & Grabowecky (2007). If you do not specify session parameter then a single mean dominance phase duration is computed for all runs of a single subject.

Random effect

The random_effect argument allows you to specify a name of the column that codes for a random effect, e.g., participant identity, bistable display (if different displays were used for a single participant), etc. If specified, it is used to fit a hierarchical model with random slopes for the history effect ($\beta_H$). Note that we if random independent intercepts are used as prior research suggest large differences in overall alternation rate between participants Brascamp et al. (2019).

Here, is the R code that specifies participants as random effect

gamma_fit <-  fit_cumhist(kde_two_observers,

And here is the corresponding model, specified for the shape parameter only as identical formulas are used for the rate parameter as well. Here, $R_i$ codes for a random effect level (participant identity) and a non-centered parametrization is used for the pooled random slopes.

$$Duration[i] \sim Gamma(shape[i], rate[i]) \ log(shape[i]) = \alpha[R_i] + \beta_H[R_i] \cdot H_{mix}[i] \ H_{mix}[i] = \text{cumulative_history}(\tau,\text{history_mixture}, \text{history_init})\ \alpha[R_i] \sim Normal(log(3), 5) \ \beta_H[R_i] = \beta^{pop}_H + \beta^{z}_H[R_i] \cdot \sigma^{pop}_H\ \beta^{pop}_H \sim Normal(0, 1) \ \beta^{z}_H[R_i] \sim Normal(0, 1) \ \sigma^{pop}_H \sim Exponential(1) \ \tau \sim Normal(log(1), 0.75)$$

Identical approach is take for $\tau$, if tau=' "1|random"' was specified (same holds for mixed_state=' "1|random"' and history_mix=' "1|random"' arguments, see below).

Fixed effects

fit_cumhist() functions allows you to specify multiple fixed effect terms as a vector of strings. The implementation is restricted to:

Although this limits usability of the fixed effects, these restrictions allowed for both a simpler model specification and a simpler underlying code. If you do require more complex models, please refer to companion vignette that provides several examples on writing model using Stan directly.

Once fitted, you can use fixef() function to extract a posterior distribution or its summary for each effect.

Cumulative history parameters

fit_cumhist() function takes four parameters for cumulative history computation (see also companion vignette):

Note that although history_init accepts only fixed values either a single value used for both states or a vector of two. In contrast, both fixed and fitted values can be used for the other three parameters. Here are possible function argument values

Once fitted, you can use history_tau(), history_mixed_state(), and history_mix() functions to obtain a posterior distribution or its summary for each parameter.

Distribution family

fit_cumhist currenly supports three distributions: 'gamma', 'lognormal', and 'normal'.


$$Duration[i] \sim Gamma(shape[i], rate[i])$$ For Gamma distribution independent linear models with a log link function are fitted for both shape and rate parameter. Priors for intercepts for both parameters are $\alpha ~ Normal(log(3), 5)$.


$$Duration[i] \sim LogNormal(\mu[i], \sigma)$$ The $\mu$ parameter is computed via a linear model with a log link function. Priors for the intercept are $\alpha ~ Normal(log(3), 5)$. Prior for $\sigma$ was $\sigma \sim Exponential(1)$.


$$Duration[i] \sim Normal(\mu[i], \sigma)$$ The $\mu$ parameter is computed via a linear model. Priors for the intercept are $\alpha ~ Normal(3, 5)$. Prior for $\sigma$ was $\sigma \sim Exponential(1)$.

Computing and using cumulative history

If you are interested in the cumulative history itself, you can extract from the fitted object via extract_history() function

H <- extract_history(gam_fit)

Alternatively, you can skip fittiong and compute history directly using predefined values via compute_history().

df <- compute_history(br_singleblock,

Try the bistablehistory package in your browser

Any scripts or data that you put into this service are public.

bistablehistory documentation built on Sept. 28, 2021, 5:10 p.m.