Stock-and-flow models represent systems as states (stocks) that accumulate over time with processes (flows) that change these variables. In this vignette, we will demonstrate how to create stock-and-flow models from scratch using sdbuildR. It covers the basics of stock-and-flow modelling in the context of psychology with an example of burnout. Note that this vignette serves as online supplemental material A accompanying the paper Formalizing Psychological Theory with sdbuildR: A Stock-and-Flow Modelling Tutorial in R by Evers et al. (under review). To reproduce the figures in the paper, please see the bottom of the corresponding .Rmd file.
Stock-and-Flow Models
Stock-and-flow models conceptualize systems in terms of quantities that accumulate (i.e., stocks) and the processes (i.e., flows) that change them over time. Stocks are like the amount of water in a bathtub, whereas flows are like the water flowing in and out of the tub. Inflows – water from the tap – raise the stock, while outflows – water through the drain – lower it. Without an outflow, the water remains in the bathtub; without an inflow, the bathtub stays empty. Stocks can have multiple inflows and outflows, but each flow can only originate in one stock and terminate in one stock. Flow rates are expressed in units over time (e.g., litre per minute). The net rate of change in the water level is the difference between the inflows and outflows. The current water level (i.e., the value of the stock) is the result of how much water has flowed in and out of the bath tub. In this way, a stock functions as a memory of past activity by integrating flows over time. When a flow directly connects two stocks, it is a material flow that conserves the quantity it carries: the source stock decreases exactly as much as the target stock increases. Informational flows, such as knowledge, beliefs, and perceptions, are not conserved, as sharing information does not deplete it (Sterman 2000, 201). Because psychological systems predominantly involve informational flows, we do not discuss material flows further. This structure is the foundation of stock-and-flow models, where stocks represent the state of a system, and flows represent the processes that alter that state over time.
Stock-and-flow models are easiest to understand through a worked
example. We will create a simple model of people waiting in a queue to
introduce the syntax of sdbuildR.
We recommend constructing models incrementally and iteratively, as
this cultivates a deeper understanding of how each modification
influences system behaviour. To do so, we create a stock-and-flow model
object called sfm:
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 onlyThe object contains default simulation settings such as the total
duration, the timestep (dt) specifying the temporal
resolution of the simulation, and a solver (euler)
indicating the numerical technique used to generate output from the
model (for more details, see Karline Soetaert et
al. 2010). We update the simulation duration to ten hours:
sfm <- sfm |>
sim_settings(stop = 10, time_units = "hours", only_stocks = FALSE)We further add a model name and description to the model metadata:
sfm <- sfm |>
meta(
name = "Queue Model",
description = "A simple model of a queue of people waiting for service"
)As we build the model iteratively, the object sfm will
continuously be overwritten. sfm contains no model
variables yet. Throughout the tutorial, we use the term “variable” for
any part of the system, be that a stock, flow, constant, or auxiliary.
Though this usage may differ from other scientific fields, we here
choose to adhere to system dynamics terminology (Ford 2019; Sterman 2000). We added a stock to
the model called people. Every stock needs an initial
condition: its value at the start of the simulation. We initialize
people at 0:
sfm <- sfm |>
stock(people, eqn = 0, label = "People Waiting in Queue")Plotting the model shows its stock-and-flow diagram:
plot(sfm)To assess the model’s dynamics, we simulate it over time and visualize the resulting timeseries:
Above, we use the pipe operator |> to pass the result
of an expression to the next expression as its first argument (e.g.,
x |> f(y) is interpreted as f(x, y)). As
shown in Figure A, people remains at the same value across
the entirety of the simulation. Stocks without flows are static, as
there is no process specifying how they change.
We first added an inflow of arrivals at a constant rate
of one person per hour.
A stock with a constant inflow increases linearly:
To decrease the stock, we add an outflow of service
which removes people from the stock at a constant rate of two people per
hour.
As the outflow rate exceeds the inflow rate, the stock becomes negative.
To rectify this implausible behaviour, a naive solution may be to
include a logical statement such as
ifelse(people < 0, 0, people). However, this
computational trick would mask model misspecification. Ideally, stocks
should remain within bounds due to plausible equations and parameters.
Below, we prevent negative people by making service
proportional to the number of people in the queue:
0.5 * people. In this way, when people is
zero, the outflow is also zero.
This formulation introduces a negative feedback loop to the system (Meadows 2008): the greater the number of people, the greater the outflow. Positive feedback loops amplify change, whereas negative feedback loops bring the system back to a target state (Sterman 2000). The number of people stabilizes at a fixed level when the outflow meets the constant inflow.
Lastly, we add a second outflow from the stock representing people leaving as they are tired of waiting. This lowers the level at which the stock stabilizes.
The queue model demonstrates a principle key to understanding stock-and-flow modelling. Variables are continuous quantities that change smoothly over time. As a result, the queue model produces fractional people, and does not track individuals but rather people on aggregate. This remains a contested issue (Tafreshi et al. 2016; Franz 2022), where some treat it as a pragmatic simplification and others as a serious empirical and philosophical commitment in need of justification (Michell 1997). If variables cannot reasonably be treated as continuous quantities, other modelling frameworks may be more appropriate (e.g., Markov models for transitions between discrete states).
Constants and Auxiliaries
Stock-and-flow models can further be supplemented with two other
variable types: constants and auxiliaries. Constants are static
parameters that do not change over the time course of the simulation. In
contrast, auxiliaries are dynamic, meaning they are computed anew at
each step. They are used for intermediate computations in flow equations
or to monitor other dynamic quantities. To illustrate the difference, a
constant defined as runif(1) will be fixed to a random
number at the beginning of the simulation, whereas an equivalently
defined auxiliary will draw a new number each time step.
Constants and auxiliaries help to make equations more interpretable and easier to revise later. For instance, we may explicitly denote the service rate with a constant:
sfm <- sfm |>
constant(service_rate, eqn = .5, label = "Service Rate") |>
update(service, eqn = service_rate * people)
order_vars <- c("leave", "service", "service_rate")
plot(sfm, show_constants = TRUE, align = order_vars, order = order_vars, spacing = .5)Similarly, we can add an auxiliary to track satisfaction: the fraction of people leaving the queue who are served, rather than abandoning it out of frustration:
sfm <- sfm |>
aux(satisfaction, eqn = service / (service + leave), label = "Satisfaction")Variable types
| Characteristic | Stock | Flow | Constant | Auxiliary |
|---|---|---|---|---|
| Role in system | Defines the state of the system; accumulates the effects of flow(s) over time | Increases or decreases a stock | Specifies static quantity | Provides intermediate computations for convenience; keeps track of changing quantities |
| Varies within time horizon | ✓ | ✓ | ✗ | ✓ |
| A process taking place over time | ✗ | ✓ | ✗ | Possibly |
| Can be captured at any given moment in time | ✓ | ✗ | ✓ | ✓ |
eqn denotes
|
Initial condition | Flow rate computed at every time step | Fixed value | Value computed at every time step |
Allowed dependencies in eqn
|
Constants and (initial values of) other stocks | Any other variable | Other constants and (initial values of) stocks | Any other variable |
| Examples | Emotions, beliefs, stress, trust, resources | Coping, learning, emotion regulation | Rates, capacities, thresholds | Performance indices, ratios, sums of stocks |
The following flowchart can be used to determine a variable’s type:
Overview of package functionality
| Function | Purpose |
|---|---|
stockflow()
|
Create empty model or load template |
stock()
|
Add or modify a stock |
flow()
|
Add or modify a flow |
constant()
|
Add or modify a constant |
aux()
|
Add or modify an auxiliary |
lookup()
|
Add or modify a lookup function |
update()
|
Add or modify any variable (generic) |
simulate()
|
Simulate model |
plot()
|
Plot model diagram or simulation |
summary()
|
Run model diagnostics |
as.data.frame()
|
Get model properties in a dataframe |
sim_settings()
|
Modify simulation specifications |
meta()
|
Modify model metadata |
export_model()
|
Export models to other formats |
Global simulation variables
In some cases, it may be useful to refer to global simulation variables in the model’s equations:
| Variable | Description | Use case |
|---|---|---|
times
|
Vector with simulation times |
Use in equations, e.g., pulse(times, 5, width = dt)
|
t
|
Current time in the ODE |
Use in equations, e.g., input(t)
|
dt
|
Time step of the simulation |
Use in equations, e.g., pulse(times, 5, width = dt)
|
Simulation specifications
We may want to simulate the system over a longer time period, or with a different time step.
sfm <- sfm |>
sim_settings(
start = 0,
stop = 250,
dt = 0.001
)Simulation settings can be set directly on the model object as above,
or passed to simulate():
sim <- simulate(sfm, start = 0, stop = 250)dt refers to the time step of the simulation, which
determines how often the model’s equations are evaluated. A smaller
dt can increase the accuracy of the simulation, but also
increases computational time and the size of the resulting dataframe. To
reduce the saved output, we may save fewer timepoints, for instance,
every 0.1 days:
sfm <- sim_settings(sfm, save_by = 0.1)Or specific time points:
sfm <- sim_settings(sfm, save_times = c(1, 50, 100))Alternatively, we can specify the number of time points to save with
save_length:
sfm <- sim_settings(sfm, save_length = 100)Similarly, we may change the numerical method used to solve the
model. The default method is "euler", which is the simplest
numerical integration method. For more complex models or when higher
accuracy is needed, consider other methods like "rk4":
sfm <- sim_settings(sfm, method = "rk4")All available simulation methods can found with:
sim_methods()
#> $R
#> [1] "euler" "rk2" "rk4" "rk23bs" "ode23"
#> [6] "rk45dp6" "rk45dp7" "rk45e" "rk45f" "rk45ck"
#> [11] "rk78dp" "rk78f" "ode45" "irk3r" "irk5r"
#> [16] "irk4hh" "irk4l" "irk6kb" "irk6l" "lsoda"
#> [21] "lsodar" "lsode" "lsodes" "bdf" "bdf_d"
#> [26] "vode" "daspk" "adams" "impAdams" "impAdams_d"
#> [31] "radau"
#>
#> $Julia
#> [1] "Euler()" "ForwardEuler()" "Midpoint()" "Heun()"
#> [5] "RK4()" "BS3()" "Tsit5()" "Vern6()"
#> [9] "Vern7()" "Vern8()" "Vern9()" "Rosenbrock23()"Note that some methods may not be available in Julia and vice versa.
In case the simulation contains stochastic elements, we can set a seed to ensure that the simulation is reproducible. For example, the initial value of people could be a random number:
The seed needs to be an integer:
sfm <- sim_settings(sfm, seed = 123)The seed can also be removed to ensure variation across simulations.
sfm <- sim_settings(sfm, seed = NULL)Renaming variables
Variable names can easily be changed:
sfm <- change_name(sfm, people, new_name = customers)This will ensure that all references to people are
changed to customers.
Allowed variable names
When creating variables or changing variable names, a warning may be issued that the name was modified to be syntactically valid and unique. For example:
sfm <- change_name(sfm, customers, new_name = t)
#> Warning: A name was changed for syntactic validity or uniqueness.
#> ℹ "t" → `t_1`The name t is not usable, as this already refers to the
current time step. Similarly, names cannot contain spaces or special
characters:
sfm <- change_name(sfm, t_1, new_name = a - b)
#> Warning: A name was changed for syntactic validity or uniqueness.
#> ℹ "a - b" → `a___b`Names also cannot be duplicated:
sfm <- change_name(sfm, service, new_name = served)Removing variables
To remove a variable from the model, use discard():
sfm <- discard(sfm, served)
#> Warning: Found a lingering reference to removed variable `served`.
#> → Check equation of variable "satisfaction".Note that this cannot be undone!
Lookup functions
Lookup functions (lookup), also known as table or
graphical functions, are interpolation functions used to create custom
input-output functions, where we define the desired output (y) for a
specified input (x). They are defined by a set of x- and y-domain
points. The interpolation method defines the behaviour of the lookup
function between x-points, and the extrapolation method defines the
behaviour outside of the x-points. For example, a simple lookup function
called "graph" may look like this:
sfm <- stockflow() |>
lookup(graph,
xpts = c(0, 1, 2), ypts = c(0.5, 1, 1),
interpolation = "linear", extrapolation = "nearest"
)The function can now be used in any equation in the model like so:
sfm <- constant(sfm, x, eqn = graph(1))Custom functions
New functions can be defined such that they can be used anywhere in
the model. For example, if the logistic() function did not
exist, you could create it yourself:
sfm <- stockflow() |>
custom_func(f, eqn = function(x, slope = 1, midpoint = .5) 1 / (1 + exp(-slope * (x - midpoint))))This will create a function f() that can be used in any
equation in the model like so:
sfm <- constant(sfm, x, eqn = f(0))Documenting
To document meta-properties of the model, use meta().
For example, the model’s name, subtitle
(caption), or author. meta()
accepts any key-value pair, so custom metadata can also be added.
sfm <- meta(sfm, author = "Kyra Evers", affiliation = "University of Amsterdam")