Skip to contents

Your first simulation

To perform your first simulation you will need:

  1. A data.frame specifiying the experiment design, and
  2. A list with the parameters for the model you’ll be using.

The design data.frame

We will use a data.frame to specify our experimental design. Let’s specify a blocking design.

library(calmr)

my_blocking <- data.frame(
  Group = c("Exp", "Control"),
  Phase1 = c("10A>(US)", "10C>(US)"),
  R1 = c(FALSE, FALSE),
  Phase2 = c("10AB>(US)", "10AB>(US)"),
  R2 = c(FALSE, FALSE),
  Test = c("1A#/1B#", "1A#/1B#"),
  R3 = c(FALSE, FALSE)
)
my_blocking
#>     Group   Phase1    R1    Phase2    R2    Test    R3
#> 1     Exp 10A>(US) FALSE 10AB>(US) FALSE 1A#/1B# FALSE
#> 2 Control 10C>(US) FALSE 10AB>(US) FALSE 1A#/1B# FALSE

A few rules about the design data.frame:

  1. Each row represents a group.
  2. Its first column contains the group labels.
  3. The remaining columns are organized in pairs (trials in a phase, and whether to randomize them)

The trials in each phase column are specified using a very rigid notation. A few observations about it:

  1. Trials are preceded by a number. That number represents the number of times that trial will be given in each phase. “10A>(US)” means that the “A>(US)” trial will be given 10 times.
  2. The presence and absence of the unconditioned stimulus are not denoted with the traditional “+” and “-” symbols. Instead, here we use parenthesis to denote “complex” stimuli. These can be thought of as an element with a complex name (i.e., with more than one character). As such, “(US)” specifies a single element to represent the US.
  3. In the same vein, multiple characters with no parentheses denote individual elements. For example, “AB” implies the presence of two stimuli, A and B.
  4. The “>” character is used as a separator of the “expectations” and “correction” steps within the trial. “10A>(US)” means that the model generates an expectation with A only, but learns from the co-occurrence of A and the US.
  5. The “/” character is used as a trial separator (it does not imply randomization by itself). Thus, “1A/1B” specifies that a single “A” trial and a single “B” trial will be given during that phase. Recall that randomization of trials within a phase is specified by the column after it (above, R1, R2, and R3).
  6. The “#” character is used to denote probe trials. In contrast to real life, probe trials here entail no update of the model’s associations. As such, probe trials can be used to track the development of key associations, with no repercussion to what the model learns on normal training trials.

If you want to check whether your trial string will work with the simulator, you can use phase_parser. The function returns a list with a lot of information, so let’s print only some of the fields.

# not specifying a number of AB trials. Bad practice!
trial <- phase_parser("AB/10AC")
trial$general_info[c("trial_names", "trial_repeats")]
#> $trial_names
#> [1] "AB" "AC"
#> 
#> $trial_repeats
#> [1]  1 10
# considering a configural cue for elements AB
trial <- phase_parser("10AB(AB)(US)")
trial$general_info[c("func2nomi")]
#> $func2nomi
#>    A    B   AB   US 
#>  "A"  "B" "AB" "US"

The parameters list

Now we need to pick a model and its parameters.

To get the models currently supported in calmr, you can call supported_models().

supported_models()
#> [1] "HDI2020" "HD2022"  "RW1972"  "MAC1975" "PKH1982" "SM2007"  "RAND"   
#> [8] "ANCCR"

After picking a model, you can get default parameters for your design with get_parameters.

my_pars <- get_parameters(my_blocking, model = "RW1972")
# Increasing the beta parameter for US presentations
my_pars$betas_on["US"] <- .6
my_pars
#> $alphas
#>   A   B   C  US 
#> 0.4 0.4 0.4 0.4 
#> 
#> $betas_on
#>   A   B   C  US 
#> 0.4 0.4 0.4 0.6 
#> 
#> $betas_off
#>   A   B   C  US 
#> 0.4 0.4 0.4 0.4 
#> 
#> $lambdas
#>  A  B  C US 
#>  1  1  1  1

The (optional) additional options

If you opt to, you can also pass a list with simulation options. Simulation options are entirely optional, and if you don’t pass them you will just get a warning.

We can get a default one with get_exp_opts, and modify them accordingly.

my_opts <- get_exp_opts()
# Increase number of simulations to average out randomization effects
my_opts$iterations <- 5
my_opts
#> $iterations
#> [1] 5
#> 
#> $miniblocks
#> [1] TRUE

Simulating

With all of the above, we can run our simulation using the run_experiment function.

my_experiment <- run_experiment(
  my_blocking,
  model = "RW1972",
  parameters = my_pars,
  options = my_opts,
)
slotNames(my_experiment)
#> [1] "arguments" "design"    "results"
# access arguments (see how many iterations we ran?)
my_experiment@arguments
#> # A tibble: 10 × 6
#>    model  iteration group   experience    mapping           parameters      
#>    <chr>      <int> <chr>   <list>        <list>            <list>          
#>  1 RW1972         1 Control <df [22 × 7]> <named list [12]> <named list [4]>
#>  2 RW1972         1 Exp     <df [22 × 7]> <named list [12]> <named list [4]>
#>  3 RW1972         2 Control <df [22 × 7]> <named list [12]> <named list [4]>
#>  4 RW1972         2 Exp     <df [22 × 7]> <named list [12]> <named list [4]>
#>  5 RW1972         3 Control <df [22 × 7]> <named list [12]> <named list [4]>
#>  6 RW1972         3 Exp     <df [22 × 7]> <named list [12]> <named list [4]>
#>  7 RW1972         4 Control <df [22 × 7]> <named list [12]> <named list [4]>
#>  8 RW1972         4 Exp     <df [22 × 7]> <named list [12]> <named list [4]>
#>  9 RW1972         5 Control <df [22 × 7]> <named list [12]> <named list [4]>
#> 10 RW1972         5 Exp     <df [22 × 7]> <named list [12]> <named list [4]>
# access results
results(my_experiment)
#> $rs
#> # A tibble: 704 × 9
#>    model  group trial trial_type phase  s1    s2    block_size value
#>    <chr>  <chr> <int> <chr>      <chr>  <chr> <chr>      <dbl> <dbl>
#>  1 RW1972 Exp       1 A>(US)     Phase1 A     A              1     0
#>  2 RW1972 Exp       2 A>(US)     Phase1 A     A              1     0
#>  3 RW1972 Exp       3 A>(US)     Phase1 A     A              1     0
#>  4 RW1972 Exp       4 A>(US)     Phase1 A     A              1     0
#>  5 RW1972 Exp       5 A>(US)     Phase1 A     A              1     0
#>  6 RW1972 Exp       6 A>(US)     Phase1 A     A              1     0
#>  7 RW1972 Exp       7 A>(US)     Phase1 A     A              1     0
#>  8 RW1972 Exp       8 A>(US)     Phase1 A     A              1     0
#>  9 RW1972 Exp       9 A>(US)     Phase1 A     A              1     0
#> 10 RW1972 Exp      10 A>(US)     Phase1 A     A              1     0
#> # ℹ 694 more rows
#> 
#> $vs
#> # A tibble: 704 × 9
#>    model  group trial trial_type phase  s1    s2    block_size value
#>    <chr>  <chr> <int> <chr>      <chr>  <chr> <chr>      <dbl> <dbl>
#>  1 RW1972 Exp       1 A>(US)     Phase1 A     A              1     0
#>  2 RW1972 Exp       2 A>(US)     Phase1 A     A              1     0
#>  3 RW1972 Exp       3 A>(US)     Phase1 A     A              1     0
#>  4 RW1972 Exp       4 A>(US)     Phase1 A     A              1     0
#>  5 RW1972 Exp       5 A>(US)     Phase1 A     A              1     0
#>  6 RW1972 Exp       6 A>(US)     Phase1 A     A              1     0
#>  7 RW1972 Exp       7 A>(US)     Phase1 A     A              1     0
#>  8 RW1972 Exp       8 A>(US)     Phase1 A     A              1     0
#>  9 RW1972 Exp       9 A>(US)     Phase1 A     A              1     0
#> 10 RW1972 Exp      10 A>(US)     Phase1 A     A              1     0
#> # ℹ 694 more rows

The run_experiment returns a lot of information. The str function will be your best friend, and there is plenty of documentation to go about.

Of course, if you are an advanced R user you will be able to dig into the data straight away.

However, the package also includes some methods to get a quick look at the results.

Plotting

Let’s use plot method to create some plots. Each model supports different types of plots according to the results they can produce (e.g., associations, responses, saliences, etc.)

# get all the plots for the experiment
plots <- plot(my_experiment)
names(plots)
#> [1] "Exp - Response Strength (RW1972)"       
#> [2] "Control - Response Strength (RW1972)"   
#> [3] "Exp - Association Strength (RW1972)"    
#> [4] "Control - Association Strength (RW1972)"
# or get a specific type of plot
specific_plot <- plot(my_experiment, type = "vs")
names(specific_plot)
#> [1] "Exp - Association Strength (RW1972)"    
#> [2] "Control - Association Strength (RW1972)"
# show which plots are supported by the model we are using
supported_plots("RW1972")
#> [1] "rs" "vs"

In this case, the RW model supports both associations and responses.

Stimulus associations

The columns below are the phases of the design and the rows denote the source of the association. The colors within each panel determine the target of the association.

plot(my_experiment, type = "vs")
#> $`Exp - Association Strength (RW1972)`

#> 
#> $`Control - Association Strength (RW1972)`

Responding

Fairly similar to the above, but this time responding is a function of the stimuli presented on a trial.

plot(my_experiment, type = "rs")
#> $`Exp - Response Strength (RW1972)`

#> 
#> $`Control - Response Strength (RW1972)`

Graphing

You can also take a look at the state of the model’s associations at any point during training, using the graph method on your experiment.

my_graph_opts <- get_graph_opts("small")
graph(my_experiment, t = 20, graph_opts = my_graph_opts)
#> $RW1972
#> $RW1972$`Exp - Associations (RW1972)`

#> 
#> $RW1972$`Control - Associations (RW1972)`

Final thoughts

The calmr package was designed to simulate quickly: write your design and get a glance at model predictions. That said, the package also has some features for advanced R users, so make sure to check more advanced vignettes when you are ready.