Stock-and-flow models present an intuitive way to formalize psychological systems as dynamic processes that unfold over time. In this vignette, we will formalize Job Demands-Resources (JD-R) theory as a stock-and-flow model. JD-R theory is a prominent framework for understanding burnout and work engagement. Note that this vignette serves as online supplemental material B 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.
Overview of System Dynamics Modelling
To develop an understanding of the system, system dynamics modelling follows a structured process (see below). The first two steps are covered in detail in the paper, and we will only cover their application to JD-R theory here.
| Step | Description | |
|---|---|---|
| 1 | Problem articulation | |
| 1(a) | Target phenomenon | Express the phenomenon to be explained as a pattern over time, which forms the reference mode throughout the modelling process. |
| 1(b) | Key variables | Select the most important variables needed to define and explain the target phenomenon. |
| 1(c) | Time horizon and time unit | Specify the time frame across which the simulation takes place and the time resolution with which variables change. |
| 2 | Dynamic hypothesis | Formulate a provisional account explaining how the target phenomenon arises endogenously from the system structure. |
| 3 | Formalization | |
| 3(a) | Formalizing variables | Formulate variables such that they can be represented as continuous quantities. |
| 3(b) | Stock-and-flow diagram | Categorize variables as constants, stocks, flows, or auxiliaries, and draw connections between variables. |
| 3(c) | Stock-and-flow model | Iteratively build and refine the stock-and-flow model to reproduce the target phenomenon. |
| 4 | Testing | Perform verification and validity tests to expose errors, misspecifications, and implausibilities in the model. |
| 5 | Application | Design and evaluate interventions to identify effective leverage points. |
Step 1. Problem Articulation
We select as our target phenomenon the development of burnout, which we characterize by the co-occurrence of gradually decreasing work engagement, increasing exhaustion, and decreasing job performance.
The key variables in JD-R theory are job demands, job resources, work engagement, exhaustion, proactive behaviour, self-undermining behaviour, and job performance. As working definitions, we follow those provided in Bakker et al. (2023), as included below.
JD-R theory primarily explains how burnout develops rather than its maintenance or recovery, suggesting a time horizon of months rather than years. We adopt six months as the time horizon for reproducing our core phenomenon. We adopt six months as the time horizon for reproducing our core phenomenon. Finally, we select a time unit of one day, as JD-R theory seems to pertain to a daily timescale (Bakker and Demerouti 2024a; Bakker et al. 2023).
| Variable | Definition |
|---|---|
| Job Demands | The physical, psychological, social, or organizational aspects of the job that require sustained physical, cognitive, and/or emotional effort and are therefore associated with certain physiological and/or psychological costs |
| Job Resources | The physical, psychological, social, or organizational aspects of the job that have motivating potential, that are functional in achieving work goals, that regulate the impact of job demands, and that stimulate learning and personal growth |
| Work Engagement | a positive, fulfilling, work-related state of mind that is characterized by vigor, dedication, and absorption. Vigor refers to high levels of energy and mental resilience while working, the willingness to invest effort in one’s work, and persistence even in the face of difficulties. Dedication implies being strongly involved in one’s work and experiencing a sense of significance, enthusiasm, and challenge. Absorption refers to being fully concentrated and happily engrossed in one’s work, whereby time passes quickly. Thus, work engagement is characterized by a high level of energy and strong identification with one’s work, whereas burnout is characterized by the opposite: a low level of energy and poor identification with one’s work |
| Exhaustion | Depletion of energy resources; also used interchangeable or as part of job strain |
| Proactive Behaviour | Also called job crafting; employees’ personal initiative to change their job demands and job resources in order to better align the design of the job with their own abilities and preferences |
| Self-Undermining Behaviour | Employees’ dysfunctional behaviors (e.g., poor communication, conflict behaviors) that create obstacles and may undermine performance. |
| Job Performance | Undefined in the literature; refers to the extent to which an individual performs well at their job (e.g., fulfilling responsibilities) |
Step 2. Dynamic Hypothesis
JD-R theory embodies the dynamic hypothesis that burnout occurs as a result of two competing feedback loops: a health impairment loop, in which excessive demands produce exhaustion, which in turn triggers self-undermining behaviour that further increases demands and depletes resources; and a motivational loop, in which resources foster work engagement, which promotes proactive behaviour that generates additional resources and lowers demands. Burnout emerges when the demands and exhaustion amplified by the health impairment loop overwhelm the resources and engagement sustained by the motivational loop. JD-R theory thus offers an endogenous explanation: it posits that self-undermining and exhaustion-induced resource depletion amplify and sustain burnout, rather than attributing it to external causes alone. In contrast, an exogenous explanation would hold rising external demands responsible without any reciprocal effect of the employee’s state on those demands. As the target phenomenon is articulated on a within-person level, the dynamic hypothesis correspondingly describes within-person dynamics.
Step 3. Formalization
Please see the paper for Step 3a and 3b; here, we only note that exhaustion was reformulated to energy.
Step 3c. Building Stock-and-Flow Models in R
Formulating a stock-and-flow diagram forms the initial step of formalization. A stock-and-flow model requires the specification of equations, parameters and initial conditions. Throughout the modelling process, our aim is to evaluate what behaviour is implied by the formalized theory, and whether it is sufficient to reproduce the target phenomenon.
For brevity, the paper only includes the formalization of the self-undermining process in JD-R theory. The complete step-by-step formalization can be found here. As a starting point, we load a base model from the model library:
sfm <- stockflow("jdr", version = 1)
print(sfm)
#>
#> ── Stock-and-Flow Model: Job Demands and Resources (JD-R) Theory ───────────────
#> 2 stocks • 2 flows
#>
#> ── Stock-Flow Structure ──
#> demands: (no flows)
#> energy: + recovery - effort
#>
#> ── Simulation Settings ──
#> Time: 0.0 to 182.0 days (dt = 0.01, save_length = 50) • euler • R • seed = 123
#> Simulation output: all variablesThe model now consists of two stocks: energy and demands. Energy has
a constant inflow rate and an outflow rate set by demands, but demands
is not connected to any flows yet. Low energy (i.e., exhaustion) should
increase job demands through self-undermining behaviour. This
qualitative relationship can be formalized with a functional
form: a mathematical representation of the effect of one variable
on another. For instance, in undermining = f(energy), the
functional form f() maps the input energy to the output
self-undermining behaviour. Self-undermining should be high when energy
is low, and low when energy is high, as energetic employees presumably
do not self-undermine.
To implement this in the model, we first try an exponential decay function:
sfm_exp <- sfm |>
constant(s_slope, eqn = 5, label = "Steep Slope") |>
constant(undermining_rate, eqn = 0.1, label = "Self-undermining Rate") |>
flow(undermining,
eqn = undermining_rate * exp(-s_slope * energy),
to = demands, label = "Self-undermining"
)
sfm_exp |>
simulate() |>
plot(main = "Self-Undermining: Exponential Decay")
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.This function reaches its maximum at zero energy, with a slope determining its steepness. As shown above, this functional form indeed produces the burnout phenomenon, but job demands grow without bounds. When a model displays implausible behaviour, it helps to revisit the theoretical propositions embedded in the equations. The exponential decay function keeps increasing demands even when energy has collapsed, implying that no energy is needed to engage in self-undermining behaviour. To formalize the notion that some energy is required to self-undermine, we may revise the model with an alternative functional form, such as a Ricker function. The Ricker is zero when energy is zero, peaks at low values of energy, and steeply drops off as energy increases (Bolker 2008). As such, the Ricker function does not further inflate demands when energy is zero, producing the burnout phenomenon with demands stabilizing at a plateau.
sfm_ricker <- sfm_exp |>
constant(location, eqn = 0.2, label = "Location") |>
constant(shape, eqn = 5, label = "Shape") |>
update(undermining,
eqn = undermining_rate * ricker(energy, location = location, shape = shape)
)
sfm_ricker |>
simulate() |>
plot(main = "Self-Undermining: Ricker Functional Form")
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.Theories rarely hinge on a single, carefully chosen initial condition. Individuals may have very different starting levels of energy, yet still experience burnout. A model is typically considered more robust if it produces the target phenomenon when initial conditions are varied within reasonable ranges. Moreover, an ill-chosen initial condition such as zero can prevent any dynamics from occurring, and different initial conditions can produce qualitatively distinct behaviours. Rather than simulating the model from a single initial condition, we can run it repeatedly in an ensemble simulation. To do so, we first initialize energy to draw from a uniform distribution bounded between 0 and 1 in both the exponential decay and Ricker model:
sfm_exp <- sfm_exp |>
update(energy, eqn = runif(1))
sfm_ricker <- sfm_ricker |>
update(energy, eqn = runif(1))Each simulation will draw a new value because eqn uses
non-standard evaluation. This means that runif(1) is stored
as an expression and not a computed value, as shown below:
as.data.frame(sfm_exp, vars = "energy", properties = "eqn")
#> type name eqn
#> 1 stock energy runif(1)To instead use a computed value, we can use the !! (bang-bang) operator to evaluate the expression and store the resulting value in the model:
sfm_exp |>
update(energy, eqn = !!runif(1)) |>
as.data.frame(vars = "energy", properties = "eqn")
#> type name eqn
#> 1 stock energy 0.904203379992396Note that we haven’t saved the above statement to the model, so the model still contains the unevaluated expression.
Both models are now ready for ensemble simulations:
sims_exp <- ensemble(sfm_exp, n = 100, save_length = 50, save_sims = TRUE)
#> Starting ensemble simulation in "Julia" with 100 simulations.
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#> ✔ Ensemble simulation completed in 3.9233 seconds.
sims_ricker <- ensemble(sfm_ricker, n = 100, save_length = 50, save_sims = TRUE)
#> Starting ensemble simulation in "Julia" with 100 simulations.
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#> ✔ Ensemble simulation completed in 2.0324 seconds.
# Plot individual simulations
which <- "sims"
# Don't show the central tendency (i.e., mean or median)
central <- "none"
plot(sims_exp,
main = "Self-undermining (exponential)",
which = which, central = central
)
plot(sims_ricker,
main = "Self-undermining (Ricker)",
which = which, central = central
)Above, we ran and visualized a hundred iterations of each model. The ensemble simulations confirm that the exponential decay model consistently produces boundless growth in demands, whereas the Ricker model robustly produces stabilized demands. In both functional forms, self-undermining can create demands even when there are no existing demands. Though producing qualitatively different dynamical behaviour, both formalizations are consistent with the verbal theory. The theory merely states that low energy produces self-undermining behaviour, but leaves underspecified how to model the effects of this relationship over time.
As ensemble simulations are often memory-intensive, only summary
statistics and not individual simulations are stored by
ensemble(). This default can be overridden with
save_sims = TRUE. Similarly, storing each time point
typically produces an unnecessarily large object, which we reduce above
by saving only fifty time points evenly spaced across the time interval
with save_length = 50. Both save_sims and
save_length can be set either on the model itself with
sim_settings() or passed to ensemble()
directly. Computational efficiency can be further improved by running
simulations in parallel, or more substantially by using a Julia backend,
as described below.
Parallel Ensemble Simulations in R
To support computationally intensive simulations,
sdbuildR enables parallelization supported by the
future package (Bengtsson
2021). The parallelization backend needs to be configured prior
to running ensemble simulations:
future::plan(future::multisession,
workers = parallelly::availableCores() - 1
)If enabled, ensemble() will distribute the ensemble
simulations across the parallel workers. Parallelization can be ended
with:
future::plan(future::sequential)Simulations with Julia
For even greater computational efficiency, simulations can be conducted in Julia. Julia is a modern, open-source programming language that reaches performance comparable to lower-level languages like C while maintaining higher-level syntax similar to R and Python (Bezanson et al. 2017). Julia is increasingly finding applications in psychology (Ernst et al. 2025; Bates et al. 2025), and offers excellent support for stock-and-flow modelling (Rackauckas 2024; Rackauckas and Nie 2017).
To enable Julia simulations in sdbuildR, first install
Julia at https://julialang.org/install/.
Next, set up the Julia environment for sdbuildR – note that
it may take up to 15-20 minutes to download all required Julia
packages:
If anything goes amiss, a detailed set-up guide is available at https://kcevers.github.io/sdbuildR/articles/julia-setup.html. Once completed, the simulation engine should be set to Julia:
sfm <- sim_settings(sfm, language = "Julia")Each (ensemble) simulation will now run in Julia. However, users do
not need to leave the R environment or interact with Julia directly, as
sdbuildR translates R to Julia code and calls Julia from R
using JuliaConnectoR package (Lenz
et al. 2022). The Julia session can be ended with
use_julia(stop = TRUE).
We now walk through the remaining formalization of JD-R theory, which is omitted from the paper for brevity. The full model is available in the model library:
An ensemble simulation shows that the model robustly produces the target phenomenon (note that we run simulations in Julia for speed):
sims <- ensemble(sfm, n = 100, save_sims = TRUE, save_length = 50)
#> Starting ensemble simulation in "Julia" with 100 simulations.
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#> ✔ Ensemble simulation completed in 2.5264 seconds.
plot(sims, vars = vars, which = "sims", central = "none")Job demands rise to an excessive degree, rapidly depleting energy and work engagement. Energy and engagement collapse to zero and do not recover, representing a permanent burnout.
Note that (ensemble) simulations can be made reproducible by setting a seed:
sfm <- sim_settings(sfm, seed = 123)As we will continue experimenting with the model, we save a copy of the current model, and we save only the main variables for computational efficiency:
vars <- c("engagement", "demands", "resources", "energy", "performance")
sfm0 <- sfm <- sim_settings(sfm, vars = vars)Complete Formalization of JD-R Theory
To follow the rationale of how this model is formalized, either a) unfold the code below or b) use the widget to walk through the steps.
Scroll through steps
Currently, we have specified the recovery inflow as a constant, which seems unlikely. Plausibly, recovery depends on the current level of energy. We revise our implementation of recovery to a Ricker function, where energy recovery is zero when energy is zero (some energy is needed to recover), recovery is highest at mid-levels of energy, and recovery drops off at higher levels of energy (where recovery is not needed).
Currently, we have specified the recovery inflow as a constant, which seems unlikely. Plausibly, recovery depends on the current level of energy. We revise our implementation of recovery to a Ricker function, where energy recovery is zero when energy is zero (some energy is needed to recover), recovery is highest at mid-levels of energy, and recovery drops off at higher levels of energy (where recovery is not needed).
sfm <- sfm_ricker |>
constant(recovery_rate, eqn = 0.3, label = "Recovery Rate") |>
flow(recovery,
eqn = recovery_rate * ricker(energy, location = 0.5, shape = shape)
)
sfm |>
simulate() |>
plot()
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.Moreover, in its current formulation, the model does not produce any energy outflow when demands are zero. However, even in the absence of demands, we presume a natural decay in energy. We thus include an explicit constant for the rate of energy decay, and revise energy’s outflow:
Moreover, in its current formulation, the model does not produce any energy outflow when demands are zero. However, even in the absence of demands, we presume a natural decay in energy. We thus include an explicit constant for the rate of energy decay, and revise energy’s outflow:
sfm <- sfm |>
constant(energy_decay_rate, eqn = 0.75, label = "Energy Decay Rate") |>
flow(effort, eqn = energy * (energy_decay_rate + demands))
summary(sfm)
#>
#> ── Stock-and-Flow Model Diagnostics ────────────────────────────────────────────
#> ✔ No problems detected!
plot(sfm)
sfm |>
simulate() |>
plot()
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.Currently, job demands can only increase, but job demands are also reduced in proportion to the current level of energy. The effect is multiplied by demands to prevent reducing demands past zero.
Currently, job demands can only increase, but job demands are also reduced in proportion to the current level of energy. The effect is multiplied by demands to prevent reducing demands past zero.
sfm <- sfm |>
update(demands, eqn = runif(1)) |>
flow(work, eqn = energy * demands, from = demands, label = "Work")
sfm |>
simulate() |>
plot()
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.In the positive feedback loop between work engagement and job resources, resources have a motivating impact on engagement and engagement drives proactive behaviour which increases resources. We first add a stock for job resources initialized at a random value. Its corresponding outflow represents the natural decay of resources over time. The outflow is buffered by energy, reflecting the idea that exhaustion drains resources. When energy is high, resources deplete at a lower rate.
In the positive feedback loop between work engagement and job resources, resources have a motivating impact on engagement and engagement drives proactive behaviour which increases resources. We first add a stock for job resources initialized at a random value. Its corresponding outflow represents the natural decay of resources over time. The outflow is buffered by energy, reflecting the idea that exhaustion drains resources. When energy is high, resources deplete at a lower rate.
sfm <- sfm |>
stock(resources, eqn = runif(1), label = "Job Resources") |>
flow(resource_decay, eqn = resource_decay_rate * resources / (1 + energy), from = resources, label = "Decay") |>
constant(resource_decay_rate, eqn = 0.1, label = "Resource Decay Rate")
sfm |>
simulate() |>
plot()
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.Similarly, we add a stock for work engagement, which is depleted by an outflow representing the natural decay of engagement over time. This outflow is also buffered by energy, as exhaustion drains engagement.
Similarly, we add a stock for work engagement, which is depleted by an outflow representing the natural decay of engagement over time. This outflow is also buffered by energy, as exhaustion drains engagement.
sfm <- sfm |>
stock(engagement, eqn = runif(1), label = "Work Engagement") |>
flow(engagement_decay, eqn = engagement_decay_rate * engagement / (1 + energy), from = engagement, label = "Dissipation") |>
constant(engagement_decay_rate, eqn = 0.2, label = "Engagement Decay Rate")
sfm |>
simulate() |>
plot()
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.With both stocks in place, we first implement the motivating effect of resources on engagement. The effect is implemented as a Hill function, which is a common way to model saturation effects in biological systems. The Hill function captures the idea that the effect of resources on engagement increases with resources but eventually saturates at high levels of resources. Additionally, according to the boost hypothesis, job demands amplify the motivating impact of resources on engagement. This effect is multiplicative according to JD-R theory. In our implementation, we further multiply the effect by energy, as we suppose energy is required to sustain motivation.
With both stocks in place, we first implement the motivating effect of resources on engagement. The effect is implemented as a Hill function, which is a common way to model saturation effects in biological systems. The Hill function captures the idea that the effect of resources on engagement increases with resources but eventually saturates at high levels of resources. Additionally, according to the boost hypothesis, job demands amplify the motivating impact of resources on engagement. This effect is multiplicative according to JD-R theory. In our implementation, we further multiply the effect by energy, as we suppose energy is required to sustain motivation.
sfm <- sfm |>
constant(m_slope, eqn = 3, label = "Medium Slope") |>
flow(motivation, eqn = motivation_rate * energy * hill(resources, m_slope) * demands, to = engagement, label = "Motivation") |>
constant(motivation_rate, eqn = 0.3, label = "Motivation Rate")
sfm |>
simulate() |>
plot()
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.Next, we implement the beneficial effect of engagement on resources. Engaged employees are more likely to engage in proactive behavior, which increases job resources. We implement this effect as a Hill function, which is a common way to model saturation effects in biological systems. The Hill function captures the idea that the effect of engagement on proactive behavior increases with engagement but eventually saturates at high levels of engagement.
Next, we implement the beneficial effect of engagement on resources. Engaged employees are more likely to engage in proactive behavior, which increases job resources. We implement this effect as a Hill function, which is a common way to model saturation effects in biological systems. The Hill function captures the idea that the effect of engagement on proactive behavior increases with engagement but eventually saturates at high levels of engagement.
sfm <- sfm |>
flow(proactive, eqn = proactive_rate * hill(engagement, m_slope), to = resources, label = "Proactive behaviour") |>
constant(proactive_rate, eqn = 0.2, label = "Proactive Behaviour Rate")
sfm |>
simulate() |>
plot()
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.As of now, job resources are only cultivated through proactive behaviour, and job demands are only created by self-undermining behaviour. However, JD-R theory implies that both job resources and demands are partially exogenously driven. Both are defined as ‘physical, psychological, social, or organizational aspects of the job’. We capture these exogenous drivers with a new inflow to job resources and demands. Here, exogenously provided resources drop off exponentially at higher levels of resources, reflecting the idea that employers are less likely to provide resources when employees are already well-resourced:
As of now, job resources are only cultivated through proactive behaviour, and job demands are only created by self-undermining behaviour. However, JD-R theory implies that both job resources and demands are partially exogenously driven. Both are defined as “physical, psychological, social, or organizational aspects of the job”. We capture these exogenous drivers with a new inflow to job resources and demands. Here, exogenously provided resources drop off exponentially at higher levels of resources, reflecting the idea that employers are less likely to provide resources when employees are already well-resourced:
sfm <- sfm |>
flow(exo_resources, eqn = exo_resource_rate * exp(-s_slope * resources), to = resources, label = "Exogenous support") |>
constant(exo_resource_rate, eqn = 0.1, label = "New resource rate")
sfm |>
simulate() |>
plot()
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.In the absence of more information about the functional form of exogenous demands, we use the same functional form as for exogenous resources. Exogenously provided demands drop off exponentially at higher levels of demands, reflecting the idea that employers are less likely to assign new tasks when employees are already at high workloads:
In the absence of more information about the functional form of exogenous demands, we use the same functional form as for exogenous resources. Exogenously provided demands drop off exponentially at higher levels of demands, reflecting the idea that employers are less likely to assign new tasks when employees are already at high workloads:
sfm <- sfm |>
flow(exo_demands, eqn = exo_demand_rate * exp(-s_slope * demands), to = demands, label = "Exogenous tasks") |>
constant(exo_demand_rate, eqn = 0.3, label = "New task rate")
sfm |>
simulate() |>
plot()
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.We next implement the buffer hypothesis, which states that though job demands decrease energy, this effect is buffered by resources. In other words, when resources are high, demands are less energy-depleting.
We next implement the buffer hypothesis, which states that though job demands decrease energy, this effect is buffered by resources. In other words, when resources are high, demands are less energy-depleting.
sfm <- sfm |>
flow(effort,
eqn = energy * (energy_decay_rate + demands) / (1 + resources)
)
sfm |>
simulate() |>
plot()
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.Furthermore, the outflow from job demands is influenced by
engagement: engaged employees ‘optimize’ job demands, which we interpret
to mean that they are more effective at reducing demands. We specify
this by multiplying energy by (1 + engagement), which
captures the idea that engagement is not necessary for reducing
demands (whereas * engagement would turn the outflow from
demands to zero when engagement is zero).
Furthermore, the outflow from job demands is influenced by
engagement: engaged employees “optimize” job demands, which we interpret
to mean that they are more effective at reducing demands. We specify
this by multiplying energy by (1 + engagement), which
captures the idea that engagement is not necessary for reducing
demands (whereas * engagement would turn the outflow from
demands to zero when engagement is zero).
sfm <- sfm |>
flow(work, eqn = energy * demands * (1 + engagement))
sfm |>
simulate() |>
plot()
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.Lastly, we formalize job performance as the sum of engagement and energy, which is consistent with the idea that both engagement and exhaustion are key determinants of performance.
Lastly, we formalize job performance as the sum of engagement and energy, which is consistent with the idea that both engagement and exhaustion are key determinants of performance.
sfm <- sfm |>
aux(performance, eqn = engagement + energy, label = "Job Performance")
sfm |>
simulate() |>
plot()
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
sfm <- sim_settings(sfm, vars = c("engagement", "demands", "resources", "energy", "performance"))Currently, we have specified the recovery inflow as a constant, which seems unlikely. Plausibly, recovery depends on the current level of energy. We revise our implementation of recovery to a Ricker function, where energy recovery is zero when energy is zero (some energy is needed to recover), recovery is highest at mid-levels of energy, and recovery drops off at higher levels of energy (where recovery is not needed).
Moreover, in its current formulation, the model does not produce any energy outflow when demands are zero. However, even in the absence of demands, we presume a natural decay in energy. We thus include an explicit constant for the rate of energy decay, and revise energy's outflow:
Currently, job demands can only increase, but job demands are also reduced in proportion to the current level of energy. The effect is multiplied by demands to prevent reducing demands past zero.
In the positive feedback loop between work engagement and job resources, resources have a motivating impact on engagement and engagement drives proactive behaviour which increases resources. We first add a stock for job resources initialized at a random value. Its corresponding outflow represents the natural decay of resources over time. The outflow is buffered by energy, reflecting the idea that exhaustion drains resources. When energy is high, resources deplete at a lower rate.
Similarly, we add a stock for work engagement, which is depleted by an outflow representing the natural decay of engagement over time. This outflow is also buffered by energy, as exhaustion drains engagement.
With both stocks in place, we first implement the motivating effect of resources on engagement. The effect is implemented as a Hill function, which is a common way to model saturation effects in biological systems. The Hill function captures the idea that the effect of resources on engagement increases with resources but eventually saturates at high levels of resources. Additionally, according to the boost hypothesis, job demands amplify the motivating impact of resources on engagement. This effect is multiplicative according to JD-R theory. In our implementation, we further multiply the effect by energy, as we suppose energy is required to sustain motivation.
Next, we implement the beneficial effect of engagement on resources. Engaged employees are more likely to engage in proactive behavior, which increases job resources. We implement this effect as a Hill function, which is a common way to model saturation effects in biological systems. The Hill function captures the idea that the effect of engagement on proactive behavior increases with engagement but eventually saturates at high levels of engagement.
As of now, job resources are only cultivated through proactive behaviour, and job demands are only created by self-undermining behaviour. However, JD-R theory implies that both job resources and demands are partially exogenously driven. Both are defined as 'physical, psychological, social, or organizational aspects of the job'. We capture these exogenous drivers with a new inflow to job resources and demands. Here, exogenously provided resources drop off exponentially at higher levels of resources, reflecting the idea that employers are less likely to provide resources when employees are already well-resourced:
In the absence of more information about the functional form of exogenous demands, we use the same functional form as for exogenous resources. Exogenously provided demands drop off exponentially at higher levels of demands, reflecting the idea that employers are less likely to assign new tasks when employees are already at high workloads:
We next implement the buffer hypothesis, which states that though job demands decrease energy, this effect is buffered by resources. In other words, when resources are high, demands are less energy-depleting.
sfm <- sfm |>
flow(effort,
eqn = energy * (energy_decay_rate + demands) / (1 + resources)
)Furthermore, the outflow from job demands is influenced by engagement: engaged employees 'optimize' job demands, which we interpret to mean that they are more effective at reducing demands. We specify this by multiplying energy by `(1 + engagement)`, which captures the idea that engagement is not *necessary* for reducing demands (whereas `* engagement` would turn the outflow from demands to zero when engagement is zero).
sfm <- sfm |>
flow(work, eqn = energy * demands * (1 + engagement))Lastly, we formalize job performance as the sum of engagement and energy, which is consistent with the idea that both engagement and exhaustion are key determinants of performance.
sfm <- sfm |>
aux(performance, eqn = engagement + energy, label = "Job Performance")Step 4. Testing
Verification Tests
Verification involves assessing whether our model behaves as we intended it to. Although it has been designed to do so, an increasingly complex model can generate unexpected behaviours. A model should conform to known real-world behaviours, physical limitations and logical constraints (Sterman 2000, 871). Such tests are called ‘reality checks’ in the system dynamics literature (Ventana Systems 2022; Peterson and Eberlein 1994), which is a more appropriate term for business and engineering than psychological systems. For example, in a population model, setting birth rates to zero should result in no new people being born. Similarly, daily work hours should never exceed 24 hours, the severity of a headache cannot become negative, and income cannot grow to infinity. Any violation indicates the model needs to be reformulated. To implement such verification tests in a stock-and-flow model, we make use of unit tests, a concept from software engineering (Duggan 2016, Ch. 6; Fowler and Beck 2019; Martin 2012). A unit test compares the behaviour of a small aspect of the model (i.e., a unit) to an explicitly formulated expectation. In JD-R theory, we may for example check that self-undermining never becomes negative, as behaviours should be strictly zero or positive:
As undermining refers to the entire timeseries of
self-undermining behaviour, we use all() to check that all
time points are equal to or above zero. To assess whether our
expectation holds, we run verify(), which simulates the
model and checks all unit tests on its output:
verify(sfm)
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#>
#>
#> ── Stock-and-Flow Unit Test Results ────────────────────────────────────────────
#>
#> 1/1 test passed.
#>
#> ✔ 1. undermining is at least 0 (for all values)The model passes the test, increasing our confidence in its
plausibility. A test label
("undermining is at least 0 (for all values)") has been
automatically generated based on the test’s expectation, but may also be
customized by passing a label. To expose more unrealistic
behaviours, the model can be subjected to extreme conditions.
Extreme values, such as zero, negative, or infinite values, tend to
reveal equation errors more readily than variations within plausible
ranges (Peterson and Eberlein 1994; Barlas
1996). For example, when job demands start and remain at zero, no
tasks ever enter the system, such that there is no work to perform well
on. In this scenario, job performance should be low. In
unit_test(), we can set the conditions under which an
expectation should hold with condition, which should be
specified as a named list with only constants or initial values of
stocks. Here, we expect the last value of job performance to be low when
demands start at zero and its inflow rates at zero:
conditions <- list(demands = 0, exo_demand_rate = 0, undermining_rate = 0)
sfm <- unit_test(sfm,
expr = tail(performance, 1) < 0.1,
conditions = conditions
)We run our tests:
verify(sfm)
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#>
#>
#> ── Stock-and-Flow Unit Test Results ────────────────────────────────────────────
#>
#> 1/2 tests passed.
#>
#> ✔ 1. undermining is at least 0 (for all values)
#> ✖ 2. the last 1 value of performance is less than 0.1 (demands = 0,
#> exo_demand_rate = 0, undermining_rate = 0)
#> Expected: TRUE Actual: FALSEOur test failed. To understand why, we may selectively plot the failed test:
sfm |>
verify() |>
plot(status = "fail")
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.Despite demands being zero throughout the simulation, performance
reaches high levels. Our formal model enables us to exactly pinpoint the
reason for this implausible behaviour, namely in performance’s
eqn.
as.data.frame(sfm, vars = performance, properties = "eqn")
#> type name eqn
#> 1 aux performance engagement + energyAs JD-R theory merely states that performance is increased by engagement and decreased by exhaustion, we have defined performance simply as the sum of engagement and energy. Demands are thus not directly necessary to perform well. This may be rectified by simply revising job performance to also depend on demands, such that all tests pass:
sfm <- update(sfm, performance, eqn = demands * (engagement + energy))
verify(sfm)
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#>
#>
#> ── Stock-and-Flow Unit Test Results ────────────────────────────────────────────
#>
#> 2/2 tests passed.
#>
#> ✔ 1. undermining is at least 0 (for all values)
#> ✔ 2. the last 1 value of performance is less than 0.1 (demands = 0,
#> exo_demand_rate = 0, undermining_rate = 0)We add this modification to our saved model as well:
sfm0 <- update(sfm0, performance, eqn = demands * (engagement + energy))Note that this is identical to the version stored in the model
library, which can be loaded using stockflow():
sfm0 <- sfm <- stockflow("jdr", version = 3)We include some additional verification tests checks that are omitted from the paper for brevity. For more options and details on unit testing, see the unit tests vignette.
sfm <- sfm |>
# Job demands and energy should negatively correlate
unit_test(expr = cor(demands, energy) < -.2) |>
# Job resources and energy should positive correlate
unit_test(expr = cor(resources, energy) > .2) |>
# Job performance and work engagement should positively correlate
unit_test(expr = cor(engagement, performance) > .2) |>
# Job performance and energy should positively correlate
unit_test(expr = cor(performance, energy) > .2) |>
# Behaviours (e.g., the flows proactive and work) should always be positive
unit_test(expr = all(proactive >= 0)) |>
unit_test(expr = all(work >= 0))
verify(sfm)
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#>
#>
#> ── Stock-and-Flow Unit Test Results ────────────────────────────────────────────
#>
#> 6/6 tests passed.
#>
#> ✔ 1. the correlation between demands and energy is less than -0.2
#> ✔ 2. the correlation between resources and energy is greater than 0.2
#> ✔ 3. the correlation between engagement and performance is greater than 0.2
#> ✔ 4. the correlation between performance and energy is greater than 0.2
#> ✔ 5. proactive is at least 0 (for all values)
#> ✔ 6. work is at least 0 (for all values)When motivation is zero, motivation should be zero at all time points:
sfm <- unit_test(sfm,
expr = all(motivation == 0),
conditions = list(motivation_rate = 0)
)
verify(sfm)
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#>
#>
#> ── Stock-and-Flow Unit Test Results ────────────────────────────────────────────
#>
#> 7/7 tests passed.
#>
#> ✔ 1. the correlation between demands and energy is less than -0.2
#> ✔ 2. the correlation between resources and energy is greater than 0.2
#> ✔ 3. the correlation between engagement and performance is greater than 0.2
#> ✔ 4. the correlation between performance and energy is greater than 0.2
#> ✔ 5. proactive is at least 0 (for all values)
#> ✔ 6. work is at least 0 (for all values)
#> ✔ 7. motivation is equal to 0 (for all values) (motivation_rate = 0)As an extreme condition test, we can initialize all stocks at zero. In this case, only job demands rises.
sfm <- unit_test(sfm,
expr = all(is.finite(c(engagement, resources, energy, demands))),
conditions = list(engagement = 0, resources = 0, energy = 0, demands = 0)
)
result <- verify(sfm)
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
print(result)
#>
#> ── Stock-and-Flow Unit Test Results ────────────────────────────────────────────
#> 8/8 tests passed.
#> ✔ 1. the correlation between demands and energy is less than -0.2✔ 2. the correlation between resources and energy is greater than 0.2✔ 3. the correlation between engagement and performance is greater than 0.2✔ 4. the correlation between performance and energy is greater than 0.2✔ 5. proactive is at least 0 (for all values)✔ 6. work is at least 0 (for all values)✔ 7. motivation is equal to 0 (for all values) (motivation_rate = 0)✔ 8. is.finite([engagement, resources, energy, demands]) (for all values)
#> (engagement = 0, resources = 0, energy = 0, demands = 0)As another robustness check, we initialize all stocks at high values, which the model is able to handle.
sfm <- unit_test(sfm,
expr = all(is.finite(c(engagement, resources, energy, demands))),
conditions = list(engagement = 5, resources = 5, energy = 5, demands = 5)
)
result <- verify(sfm)
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
print(result)
#>
#> ── Stock-and-Flow Unit Test Results ────────────────────────────────────────────
#> 9/9 tests passed.
#> ✔ 1. the correlation between demands and energy is less than -0.2✔ 2. the correlation between resources and energy is greater than 0.2✔ 3. the correlation between engagement and performance is greater than 0.2✔ 4. the correlation between performance and energy is greater than 0.2✔ 5. proactive is at least 0 (for all values)✔ 6. work is at least 0 (for all values)✔ 7. motivation is equal to 0 (for all values) (motivation_rate = 0)✔ 8. is.finite([engagement, resources, energy, demands]) (for all values)
#> (engagement = 0, resources = 0, energy = 0, demands = 0)✔ 9. is.finite([engagement, resources, energy, demands]) (for all values)
#> (engagement = 5, resources = 5, energy = 5, demands = 5)Show all unit tests:
unit_tests(sfm)
#>
#> ── Stock-and-Flow Unit Tests ───────────────────────────────────────────────────
#> 9 tests • 9/9 active • 3/9 include conditions
#> • 1. the correlation between demands and energy is less than -0.2
#> `cor(demands, energy) < -0.2`
#> • 2. the correlation between resources and energy is greater than 0.2
#> `cor(resources, energy) > 0.2`
#> • 3. the correlation between engagement and performance is greater than 0.2
#> `cor(engagement, performance) > 0.2`
#> • 4. the correlation between performance and energy is greater than 0.2
#> `cor(performance, energy) > 0.2`
#> • 5. proactive is at least 0 (for all values)
#> `all(proactive >= 0)`
#> • 6. work is at least 0 (for all values)
#> `all(work >= 0)`
#> • 7. motivation is equal to 0 (for all values) (motivation_rate = 0)
#> `all(motivation == 0)`
#> Conditions: motivation_rate = 0
#> • 8. is.finite([engagement, resources, energy, demands]) (for all values)
#> (engagement = 0, resources = 0, energy = 0, demands = 0)
#> `all(is.finite(c(engagement, resources, energy, demands)))`
#> Conditions: engagement = 0, resources = 0, energy = 0, demands = 0
#> • 9. is.finite([engagement, resources, energy, demands]) (for all values)
#> (engagement = 5, resources = 5, energy = 5, demands = 5)
#> `all(is.finite(c(engagement, resources, energy, demands)))`
#> Conditions: engagement = 5, resources = 5, energy = 5, demands = 5Plot all tests:
sfm |>
verify() |>
plot()
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.Uncertainty
To illustrate the impact of aleatory uncertainty in our JD-R model, we substitute the deterministic formulation of demand influx for a stochastic process. This more closely aligns with JD-R theory, which posits job demands may fluctuate rapidly on a daily basis due to environmental volatility (Bakker and Demerouti 2024b; Downes et al. 2021).
sfm <- sfm0 |>
constant(demand_mean, eqn = .5) |>
constant(demand_sigma, eqn = 10) |>
constant(exo_demand_rate, eqn = .25) |>
aux(D_deterministic, eqn = exo_demand_rate * (demand_mean - demands)) |>
aux(D_stochastic, eqn = demand_sigma * sqrt(demands) * rnorm(1) * sqrt(dt)) |>
update(exo_demands, eqn = D_deterministic + D_stochastic)
sim <- simulate(sfm, seed = 1)
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
plot(sim)As shown above, job demands now exhibit autocorrelated fluctuations, which in turn create variability in job performance and energy. Intrinsic variability in one variable may thus propagate to other variables, even when the latter are strictly deterministically formulated.
Removing Variables
Theories tend to inflate over time, acquiring more assumptions, constructs, and interactions in an attempt to gain explanatory breadth and depth (Haslam 2016; Meehl 1990a; Smid 2023). Though this may indeed improve a theory’s explanatory power, it can stand in direct opposition to the principle of parsimony (Keas 2018). Formal models enable a direct comparison between the predictions of a more extensive versus a simpler version of the theory [also known as a perturbation analysis; Weisberg (2013), p. 89]. In JD-R theory, proactive and self-undermining behaviours are newer additions to the original theory proposed in 2001 (Bakker et al. 2023), and we may wonder what their contribution is to the model’s behaviour. We could remove self-undermining behaviour, and compare these two models:
sfm2 <- discard(sfm0, undermining)Note that discard() does not modify equations, such that
references to the removed variable may still linger in the model. A
quick check confirms that there are no undefined variables or other
common mistakes in the model specification:
summary(sfm2)
#>
#> ── Stock-and-Flow Model Diagnostics ────────────────────────────────────────────
#> ✔ No problems detected!
vars <- c("engagement", "demands", "resources", "energy")
n <- 100
sims1 <- ensemble(sfm0, n = n)
#> Starting ensemble simulation in "Julia" with 100 simulations.
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#> ✔ Ensemble simulation completed in 1.4494 seconds.
plot(sims1, vars = vars)
sims2 <- ensemble(sfm2, n = n)
#> Starting ensemble simulation in "Julia" with 100 simulations.
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#> ✔ Ensemble simulation completed in 1.7894 seconds.
plot(sims2, vars = vars, main = "Without Self-Undermining Behaviour")As presented above, in our implementation of JD-R theory, removing self-undermining behaviour has little effect on the model’s behaviour. Self-undermining behaviour is thus less empirically relevant: it is not required for burnout to occur. More generally, systematically removing components reveals which aspects of the theory are necessary for the target phenomenon and which are theoretically redundant, directly assessing empirical relevance (Dongen et al. 2025).
Challenging the Model Boundary
Endogenous and exogenous variables can be distinguished by looking at their dependencies:
dependencies(sfm0, name = performance)
#> $performance
#> [1] "demands" "engagement" "energy"By reversing the dependencies, we obtain which variables depend on job performance:
dependencies(sfm0, name = performance, reverse = TRUE)
#> $performance
#> character(0)Job performance has no effect on any part of the system. In other words, it is merely an outcome variable, illustrating a case of open-loop thinking. However, it seems plausible that engagement does not only increase job performance, but that performance itself contributes to engagement. To represent this idea, we add a new inflow to engagement that grows with performance:
sfm2 <- sfm0 |>
constant(performance_effect,
eqn = .225,
label = "Effect Job Performance on Work Engagement"
) |>
flow(pride,
eqn = performance_effect * hill(performance, slope = 2, midpoint = 0.5),
to = engagement
)
sim <- simulate(sfm0)
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
sim2 <- simulate(sfm2)
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
plot(sim)
plot(sim2)As shown above, this feedback loop saves the system from burning out. Expanding the model boundary thus has the potential to change model behaviour.
Step 5. Application
Exploring Model-Implied Phenomena
Parameters (i.e., constants) are often of greater theoretical
interest than initial conditions, as their variation can for instance
represent individual differences, contextual factors, or uncertainty in
the theory’s assumptions. By systematically varying parameters, we can
map how variation in the model’s inputs relates to variation in its
outputs. In other words, we learn exactly what behaviour is implied by
the theory. Parameters can be varied in ensemble simulations by
redefining their eqn to draw from a distribution, or by
passing ensemble() a set of values in
conditions. For example, we can simulate two values of
motivation_rate, the rate at which resources increase
engagement.
Before running the ensemble simulation, ensure that the simulation output is reduced, for example by only saving 50 evenly spaced time points:
sfm <- sfm0 |> sim_settings(save_length = 50)
# Define values to vary
conditions <- list(motivation_rate = c(0.2, 0.7))
# Retain individual simulations
sfm <- sim_settings(sfm, save_sims = TRUE)
# Generate ensemble
n <- 100
sims <- ensemble(sfm, n = n, conditions = conditions)
#> Starting ensemble simulation in "Julia" with 200 simulations in total.
#> ℹ 2 conditions x 100 simulations per condition.
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#> ✔ Ensemble simulation completed in 2.1648 seconds.“What If” Scenarios: Developing Interventions
Our ensemble simulations revealed that boosting the rate at which
resources increase engagement could be an effective intervention to
prevent burnout. For example, we may imagine that this parameter could
be targeted by a training that teaches employees to make better use of
their existing resources (Bakker and Van
Wingerden 2021). To develop an intuition of the system’s response
to increasing motivation_rate, we can simulate an idealized
intervention which is active for a particular time period [more
realistic implementations can be explored at later modelling stages;
Sterman (2000), p. 914].
To implement this, motivation_rate first needs to be
converted from a constant to a stock, as it should increase over the
time course of the simulation:
sfm <- change_type(sfm0, motivation_rate, new_type = "stock")Next, we create a pulse function that is 1 for a period of two weeks
and 0 otherwise. The use of input and interpolation functions should be
preferred over using if-statements. Floating-point precision errors
introduce small numerical inaccuracies in the solver. As a result, hard
logical conditions like if (t == 0.5) can yield
unpredictable results, where the condition may fail to occur at all. In
contrast, interpolation functions make the model more robust to
numerical errors. In pulse(), we set the starting time of
the intervention to 21 days and its duration to 14 days. Additionally,
we pass the global variable times as its first argument,
which specifies the simulation time vector. Other types of external
inputs can be created with the step(), ramp(),
and pulse() functions.
sfm <- sfm |>
constant(start, eqn = 21) |>
constant(duration, eqn = 14) |>
constant(intervention, eqn = pulse(times, start, width = duration)) |>
flow(intervention_effect,
eqn = 0.05 * intervention(t),
to = motivation_rate, label = "Intervention Effect"
)Run a single simulation:
sfm |>
simulate() |>
plot()
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.Run an ensemble simulation:
sfm <- sim_settings(sfm, save_sims = TRUE)
n <- 100
sims <- ensemble(sfm, n = n)
#> Starting ensemble simulation in "Julia" with 100 simulations.
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#> ✔ Ensemble simulation completed in 2.8231 seconds.
plot(sims,
which = "sims", sim = 1:n, alpha = .25, central = "none",
vars = c("engagement", "demands", "motivation_rate")
)To compute the effectiveness of the intervention, we increase the simulation length and only save the last timepoint.
n <- 1000
sfm <- sim_settings(sfm,
# Only save engagement at the last time point
vars = "engagement", stop = 1000, save_times = 1000
)
sims <- ensemble(sfm, n = n)
#> Starting ensemble simulation in "Julia" with 1000 simulations.
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#> ✔ Ensemble simulation completed in 32.8425 seconds.
get_prop <- function(sims) {
as.data.frame(sims,
direction = "wide",
which = "sims"
)$engagement |>
round() |>
table() |>
proportions() |>
as.data.frame() |>
setNames(c("Engagement", "Proportion"))
}
get_prop(sims)
#> Engagement Proportion
#> 1 0 0.784
#> 2 1 0.216A shorter intervention is not as effective.
sfm2 <- update(sfm, duration, eqn = 7)
sims2 <- ensemble(sfm2, n = n)
#> Starting ensemble simulation in "Julia" with 1000 simulations.
#> ℹ Julia session is already initialized
#> ℹ In case of issues, run `use_julia(restart = TRUE)` to restart Julia.
#> ✔ Ensemble simulation completed in 33.0858 seconds.
get_prop(sims2)
#> Engagement Proportion
#> 1 0 0.878
#> 2 1 0.122Informing Experimental and Statistical Design
Finally, formal models are powerful tools for supporting
theoretically informed experimental and statistical design. The
derivation chain from theory to empirical test involves a multitude of
decisions that the theory itself does not constrain. Though
meta-analyses may quantify the impact of such decisions, they do not
resolve whether discrepant findings reflect mere design artefacts or
genuine challenges to the theory (Meehl
1990b). By contrast, a formally specified theory predicts what
discrepancies are implied by the theory itself. As an
illustration, we use our JD-R model to predict the results of a
cross-lagged panel model, a widely used analysis in the JD-R literature
(Upadyaya et al. 2016; Hakanen et al. 2008;
Sorjonen et al. 2024). We assess the cross-lagged relationship
between job demands and work engagement with lavaan (Rosseel 2012) on an ensemble dataset (n =
10,000), sampled at day 20 (wave 1) and three months later at day
110 (wave 2).
As shown below, our model implies that engagement and demands have strong autoregressive effects. Counter-intuitively, higher demands lead to higher future work engagement. Engagement has little effect on future demands. Furthermore, our JD-R model allows us to assess how these effects depend on the time between waves. For instance, as shown below, though the autoregressive effects stay consistent, the effect of engagement on future demands flips in sign as the time between waves increases. If provided with only a verbal theory, there would be no principled basis for anticipating this lag dependence.
As shown below, our model implies that engagement has a strong positive effect on itself, whereas demands have little effect on future demands. Higher demands lead to lower future work engagement, but higher engagement increases future demands. Furthermore, our JD-R model allows us to assess how these effects depend on the time between waves. For instance, as shown below, though the autoregressive effect of job demands is initially negative, it flips in sign as the time between waves increases. If provided with only a verbal theory, there would be no principled basis for anticipating this lag dependence.
Stop Julia session:
use_julia(stop = TRUE)
#> ✔ Closed Julia session.Session Information
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
#> [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
#> [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
#> [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] lavaan_0.7-2 bslib_0.12.0 kableExtra_1.4.1 sdbuildR_2.2.3
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.6 xfun_0.60 ggplot2_4.0.3
#> [4] htmlwidgets_1.6.4 visNetwork_2.1.4 quadprog_1.5-8
#> [7] vctrs_0.7.3 tools_4.6.1 crosstalk_1.2.2
#> [10] generics_0.1.4 parallel_4.6.1 stats4_4.6.1
#> [13] fansi_1.0.7 tibble_3.3.1 pkgconfig_2.0.3
#> [16] data.table_1.18.6.1 RColorBrewer_1.1-3 S7_0.2.2
#> [19] desc_1.4.3 lifecycle_1.0.5 compiler_4.6.1
#> [22] farver_2.1.2 stringr_1.6.0 textshaping_1.0.5
#> [25] mnormt_2.1.2 codetools_0.2-20 httpuv_1.6.17
#> [28] htmltools_0.5.9 sass_0.4.10 yaml_2.3.12
#> [31] plotly_4.12.1 pillar_1.11.1 pkgdown_2.2.1
#> [34] later_1.4.8 jquerylib_0.1.4 tidyr_1.3.2
#> [37] MASS_7.3-65 cachem_1.1.0 mime_0.13
#> [40] tidyselect_1.2.1 digest_0.6.39 stringi_1.8.9
#> [43] dplyr_1.2.1 purrr_1.2.2 fastmap_1.2.0
#> [46] grid_4.6.1 cli_3.6.6 magrittr_2.0.5
#> [49] DiagrammeR_1.0.12 pbivnorm_0.6.0 withr_3.0.3
#> [52] scales_1.4.0 promises_1.5.0 rmarkdown_2.31
#> [55] httr_1.4.8 igraph_2.3.3 otel_0.2.0
#> [58] deSolve_1.42 JuliaConnectoR_1.1.5 ragg_1.5.2
#> [61] memoise_2.0.1 shiny_1.14.0 evaluate_1.0.5
#> [64] knitr_1.51 viridisLite_0.4.3 rlang_1.3.0
#> [67] downlit_0.4.5 Rcpp_1.1.2 xtable_1.8-8
#> [70] glue_1.8.1 xml2_1.6.0 svglite_2.2.2
#> [73] rstudioapi_0.19.0 jsonlite_2.0.0 R6_2.6.1
#> [76] systemfonts_1.3.2 fs_2.1.0