Theory Building in Psychology (Summer School 2026) • sdbuildR

These materials were developed for a stock-and-flow modelling workshop given at the Summer School Theory Building in Psychology (University of Amsterdam, July 2026; rendered materials can be accessed here). As the materials supported a lecture, they contain limited written explanations. For more detailed guidance on stock-and-flow modelling, see the Build vignette.

Package Setup

First install the package from GitHub:

if (!requireNamespace("remotes", quietly = TRUE)) install.packages("remotes")
remotes::install_github("kcevers/sdbuildR")

library(sdbuildR)

Conceptual Understanding of Stock-and-Flow Models

A Concrete Example: A Queue of People Waiting for Service

In iterative steps, we construct a simple model of a queue of people waiting for service. Either scroll through the steps below by unfolding the section, or use the widget.

Scroll through steps

Stock-and-flow diagram

Recap

Stock-and-flow models formalize the processes by which states change over time.
Stocks are accumulations of a quantity, and flows are the processes that increase or decrease those accumulations. The change in a stock is the net flow (inflow minus outflow) into that stock.
A model can have multiple stocks, and each stock can have multiple inflows and outflows.
Stocks are continuous quantities, meaning we can end up with fractional values (e.g., 3.5 people in a queue). This is a simplifying assumption.
When the outflow depends on the stock, the model contains a negative feedback loop: the more people in the queue, the greater the outflow, which in turn shrinks the queue. A positive feedback loop amplifies; a negative feedback loop dampens or balances.

Functional Forms

A functional form is the mathematical relationship specifying how one variable affects another. Above, we modified the service outflow from a constant to a linear function of the queue. This is visualized in the widget below:

eqn =

Applied Examples of Stock-and-Flow Models

The stock-and-flow framework can be applied to model a wide range of phenomena.

Stock-and-Flow Models to Build Psychological Theory

Similarly, we may use stock-and-flow models to formalize psychological theory. For example, we can model the dynamics of knowledge acquisition and forgetting:

Stock-and-flow models force us to specify how a state increases and decreases over time.

Building Stock-and-Flow Models in R

Load the stock-and-flow model of people waiting in a queue from the package’s model library:

sfm <- stockflow("queue")

sfm stands for stock-and-flow model.

Look at the model’s structure:

print(sfm)
#> 
#> ── Stock-and-Flow Model: Queue Model ───────────────────────────────────────────
#> 1 stock • 3 flows • 1 constant • 1 auxiliary
#> 
#> ── Stock-Flow Structure ──
#> people: + arrivals - leave - service
#> 
#> ── Other Variables ──
#> Constants: `service_rate`
#> Auxiliaries: `satisfaction`
#> 
#> ── Simulation Settings ──
#> Time: 0.0 to 10.0 hours (dt = 0.01) • euler • R
#> Simulation output: all variables

Or inspect it in data frame format:

as.data.frame(sfm, properties = "eqn")
#>       type         name                         eqn
#> 1    stock       people                           0
#> 2     flow     arrivals                           1
#> 3     flow        leave                   people^10
#> 4     flow      service       service_rate * people
#> 5 constant service_rate                         0.5
#> 6      aux satisfaction service / (service + leave)

Plot the stock-and-flow diagram:

plot(sfm, show_constants = TRUE)

Simulate the model and plot the resulting time series:

sfm |>
  simulate(save_length = 101) |>
  plot()

The pipe operator |> passes the result on its left into the function on its right. For example, x |> f(y) is interpreted as f(x, y).

The timeseries can also be inspected in data frame format:

sfm |>
  simulate() |>
  as.data.frame(direction = "wide") |>
  head()
#>   time     people    service        leave arrivals satisfaction
#> 1 0.00 0.00000000 0.00000000 0.000000e+00        1          NaN
#> 2 0.01 0.01000000 0.00500000 1.000000e-20        1            1
#> 3 0.02 0.01995000 0.00997500 9.986861e-18        1            1
#> 4 0.03 0.02985025 0.01492513 5.616681e-16        1            1
#> 5 0.04 0.03970100 0.01985050 9.727793e-15        1            1
#> 6 0.05 0.04950249 0.02475125 8.836307e-14        1            1

A Simplified Model of Burnout

We now use stock-and-flow modeling to formalize a simplified model of burnout. We first draw the target phenomenon:

Create an empty stock-and-flow model:

sfm <- stockflow()
print(sfm)
#> 
#> ── Stock-and-Flow Model ────────────────────────────────────────────────────────
#> ℹ Empty model without any variables.
#> 
#> ── Simulation Settings ──
#> Time: 0 to 100 seconds (dt = 0.01) • euler • R
#> Simulation output: stocks only

Change name of the model:

sfm <- meta(sfm, name = "Burnout")

Defining the Time Horizon and Time Unit

Change simulation settings:

sfm <- sim_settings(sfm,

  # Run simulation for 6 months (~ 180 days)
  stop = round(365 / 2), time_unit = "days",

  # Return all variables in output (not just stocks)
  only_stocks = FALSE
)

print(sfm)
#> 
#> ── Stock-and-Flow Model: Burnout ───────────────────────────────────────────────
#> ℹ Empty model without any variables.
#> 
#> ── Simulation Settings ──
#> Time: 0.0 to 182.0 days (dt = 0.01) • euler • R
#> Simulation output: all variables

Adding a Stock

sfm <- stock(sfm, name = engagement, eqn = .3)

Plot stock-and-flow diagram:

sfm |> plot()

Simulate and visualise timeseries:

sfm |>
  simulate() |>
  plot(animation = "time")

Adding an Outflow

sfm <- constant(sfm, decay_rate, eqn = .05) |>
  flow(decay, eqn = decay_rate, from = engagement)

Plot stock-and-flow diagram:

sfm |> plot()

Simulate and visualise timeseries:

sfm |>
  simulate() |>
  plot(animation = "time")

Stock-Dependent Outflow

sfm <- update(sfm, decay, eqn = decay_rate * engagement)

Plot stock-and-flow diagram:

sfm |> plot()

Simulate and visualise timeseries:

sfm |>
  simulate() |>
  plot(animation = "time")

Adding an Inflow

sfm <- constant(sfm, enjoyment, eqn = .3) |>
  flow(motivation, eqn = enjoyment, to = engagement)

Plot stock-and-flow diagram:

sfm |> plot()

Simulate and visualise timeseries:

sfm |>
  simulate() |>
  plot(animation = "time")

The stock settles at a steady state where inflow == outflow.

Dynamic Inflow Rate

What if work enjoyment is not static, but erodes over time? Put differently, what if work enjoyment is not a constant, but a stock?

sfm <- change_type(sfm, enjoyment, new_type = "stock")
sfm |> plot()

Adding an Outflow from Enjoyment

sfm <- sfm |>
  flow(overcommitment,
    eqn = enjoyment * new_projects,
    from = enjoyment
  ) |>
  aux(new_projects, eqn = .1 * engagement)

Plot stock-and-flow diagram:

sfm |> plot()

Simulate and visualise timeseries:

sfm |>
  simulate() |>
  plot(animation = "time")

Recap

Iteratively, we have built a stock-and-flow model with:

A stock for engagement
A constant outflow, updated to a stock-dependent outflow
A constant inflow, updated to an inflow with a dynamic rate

In other words, the inflow rate changed from an exogenous to an endogenous variable; from a constant to a stock that erodes over time.

The table below provides an overview of each model revision and the behaviour it produces.

Connecting Equations to Model Behaviour
Panel	Stocks	Constants	Recovery eqn (inflow)	Depletion eqn (outflow)	Interpretation	Behaviour
A	engagement				No process of change	Static
B	engagement	`decay_rate`		`decay_rate`	Engagement decreases at a constant rate	Linear decrease
C	engagement	`decay_rate`		`decay_rate * engagement`	Engagement decreases at a rate proportional to its current value	Exponential decrease towards zero
D	engagement	`decay_rate`, `enjoyment`	`enjoyment`	`decay_rate * engagement`	Engagement changes at a rate equal to a constant minus a rate proportional to its current value	Stability when recovery and depletion are equal
E	engagement, enjoyment	`decay_rate`	`enjoyment`	`decay_rate * engagement`	Engagement recovers at a rate which itself changes over time	Rise and collapse

Dependence on Initial Condition

If you did not complete all steps above, load the model from the model library:

sfm <- stockflow("burnout")

Rather than a fixed value, initialize the engagement stock with a random value between 0 and 4:

sfm <- update(sfm, engagement, eqn = runif(1, min = 0, max = 4))

Equations are stored as expressions. Note the initial value of engagement:

as.data.frame(sfm, properties = "eqn")
#>       type           name                        eqn
#> 1    stock     engagement runif(1, min = 0, max = 4)
#> 2    stock      enjoyment                        0.3
#> 3     flow          decay    decay_rate * engagement
#> 4     flow     motivation                  enjoyment
#> 5     flow overcommitment   enjoyment * new_projects
#> 6 constant     decay_rate                       0.05
#> 7      aux   new_projects           0.1 * engagement

This means that each simulation will have a different initial value for engagement:

sfm |>
  simulate() |>
  plot()

sfm |>
  simulate() |>
  plot()

Systematically Investigating Dependence on Initial Conditions

Update simulation settings to prepare for ensemble simulations:

sfm <- sim_settings(sfm,
  # Save only 50 time points
  save_length = 50,

  # Save only stocks (not flows or auxiliaries)
  only_stocks = TRUE,

  # Save individual simulation runs (not just summary statistics)
  save_sims = TRUE
)

Run 50 simulations from different initial conditions:

sims <- ensemble(sfm, n = 50)
#> Starting ensemble simulation in "R" with 50 simulations.
#> ✔ Ensemble simulation completed in 11.495 seconds.

Summary statistics:

plot(sims)

Individual simulation runs:

plot(sims, which = "sims")

Parameter Sensitivity

Ensemble simulations can also be used to explore the sensitivity of the model to changes in parameter values.

# Reload model from library to reset any changes made above
sfm <- stockflow("burnout") |>
  update(engagement, eqn = runif(1, min = 0, max = 4)) |>
  sim_settings(save_length = 50, only_stocks = TRUE, save_sims = TRUE)

# Specify a range of values of decay_rate to vary
conditions <- list(decay_rate = c(0.005, 0.01, 0.05, 0.1, 0.2))

# Run 50 simulations for each condition
sims <- ensemble(sfm, n = 50, conditions = conditions)
#> Starting ensemble simulation in "R" with 250 simulations in total.
#> ℹ 5 conditions x 50 simulations per condition.
#> ✔ Ensemble simulation completed in 56.3373 seconds.

plot(sims, which = "sims", condition_display = "slider")

Assignment: Expressing Theoretical Ideas with Equations

Set-up

Let’s return to the simpler model with a static recovery rate:

sfm <- sfm |>
  # Remove the outflow and auxiliary
  discard(c(overcommitment, new_projects)) |>
  # Change recovery_rate back to a constant
  change_type(enjoyment, new_type = "constant")

sfm |> plot(show_constants = TRUE)

This model has a constant inflow. The functional form of the motivation inflow can be changed to explore different theoretical assumptions about how motivation depends on engagement.

Functional Forms

eqn =

Questions

Explore the functional forms above:

What theoretical ideas does each functional form express?
How do they differ in terms of the dynamics they produce?
Is the behaviour dependent on the initial condition of energy?

Parallel Simulations

If you have a multi-core machine, you can run simulations in parallel to speed up simulations:

# Set up parallel execution (3 cores less than the total number of available cores)
n_workers <- max(c(1, parallel::detectCores() - 3))
future::plan(future::multisession, workers = n_workers)

If configured, ensemble() will automatically run simulations in parallel. Otherwise, it will run sequentially. Run 100 simulations:

sims <- ensemble(sfm, n = 100, save_sims = TRUE)
#> Starting ensemble simulation in "R" with 100 simulations.
#> ✔ Ensemble simulation completed in 20.5051 seconds.

Individual simulation runs:

plot(sims, which = "sims")

If parallel execution was configured above, restore sequential execution:

future::plan(future::sequential)