# Analysis Documentation

A Naive Theory-driven Analysis of Intra-individual Dynamics

Author

Yongfu Liao

Published

August 12, 2023

## 1 Logistic Growth

The ordinary differential equation in (1.1) expresses the rate of logistic growth of P. To compute the trajectory of the growth of P over time, the initial value of P needs to be first determined. This is often termed the initial condition of a dynamical system and we denote the initial condition as P_0 here.

\begin{equation} \frac{d P}{dt} = r P (1 - \frac{P}{K}) \end{equation} \tag{1.1}

With a given P_0, we can integrate equation (1.1) over time to arrive at any value P_t where t>0, as shown in (1.2):

\begin{equation} P_t = P_0 + \int_{0}^{t} r P (1 - \frac{P}{K}) dt \end{equation} \tag{1.2}

Computationally, equation (1.2) can be re-expressed in forward Euler form as (1.3), where \Delta t \rightarrow 0. In practice, it would be computationally prohibitive for \Delta t to approach zero. Thus a small value of \Delta t is chosen to give a good enough approximation to the P_t values. In addition, more complicated but also more accurate procedures such as the Runge-Kutta 4 (RK4) method are often used instead of the forward Euler method to achieve better numerical stability.

\begin{cases} P_\text{t=0} &= P_0 \\ P_{t + \Delta t} &= P_t + r P_t (1 - \frac{P_t}{K}) \Delta t \\ \end{cases} \tag{1.3}

In R, expression (1.3) can be implemented with the numerical integrators provided in the `deSolve` package. The code below uses `deSolve` to compute the trajectory of P_t over 49 days and reproduces Figure 4 in the main article.

Show the code

    library(deSolve)   # ODE numerical solvers
    source("utils.R")  # helper functions

    # Define ordinary differential equation for logistic population growth
    lpg = function(Time, State, Pars) {
        with( as.list(c(State, Pars)), {
            dP = r * P * (1 - P/K)
            return( list(c(dP)) )
        })
    }

    # ODE solver
    ODE = function( r=.3, K=10, P0=.01 ) {
        pars = c( r = r, K = K)     # ODE parameters
        init = c( P = P0 )          # Initial condition
        times = seq(0, 49, by=.01)  # Time points to compute values
        out = ode( init, times, lpg, pars )   # Integrate
        out
    }

    # Plot
    par( mfrow= c(1, 2) )

    plot(1, type="n", xlab = "Day", ylab = "P", ylim=c(0,10), xlim=c(0,50))
    lines( ODE(r=.2, K=10, P0=.01), col = col.alpha(2,1), lwd=2)
    lines( ODE(r=.3, K=10, P0=.01), col = col.alpha(2,1), lwd=3)
    lines( ODE(r=.4, K=10, P0=.01), col = col.alpha(2,1), lwd=4)
    legend(30, 2.5, c("r = .2  K = 10",
                      "r = .3  K = 10",
                      "r = .4  K = 10"), col=2, lwd=2:4)

    plot(1, type="n", xlab = "Day", ylab = "P", ylim=c(0,10), xlim=c(0,50))
    abline(h=c(10,8,6), col="grey", lty="dashed", lwd=2)
    lines( ODE(r=.3, K=10, P0=.01), col = 4, lwd=4)
    lines( ODE(r=.3, K= 8, P0=.01), col = 4, lwd=3)
    lines( ODE(r=.3, K= 6, P0=.01), col = 4, lwd=2)
    legend(30, 2.5, c("r = .3  K = 10",
                      "r = .3  K =  8",
                      "r = .3  K =  6"), col=4, lwd=4:2 )

## 2 Gaussian Processes

Gaussian processes are used to model the covariations between a set of categories as a function of the distances between them. At the core of Gaussian processes is a kernel function that defines how the covariance is computed from the distance. There are multiple popular choices for formulating the kernel function. Here, we select the form specified in (2.1). Figure 1 plots three kernel functions, showing how the covariance decays differentially as the distance increases according to different parameter values of \rho.

k = \eta ~ \text{exp}( -\rho \text{d}^\text{2} ) \tag{2.1}

Show the code

    library(latex2exp)
    curve( 2 * exp( - 0.2 * x^2 ), from=0, to=9.5, lwd=2,
           ylab = "K (covariance)", xlab = "distance" )
    curve( 2 * exp( - 0.1 * x^2 ),  add=T, lwd=2, col=2 )
    curve( 2 * exp( - 0.05 * x^2 ), add=T, lwd=2, col=4 )
    legend( 6.4, 2, TeX( paste0("\\eta = 2,  \\rho = ",c(.2,.1,.05)) ),
            lwd=3, col=c(1,2,4) )

Figure 1: Gaussian process kernel

With a chosen kernel function, the covariance matrix between any set of categories can be computed, given their pairwise distances. In the main article, we apply the kernel function in (2.1) to time series observations, defining the difference in time between two observations as the distance between them. We can therefore arrive at a covariance matrix specifying the covariance between any pair of observations, and sample the full set of time series observations from the multivariate normal with the aforementioned covariance matrix. Formally, this is expressed in (2.2):

\begin{aligned} \begin{bmatrix} N_1 \\ N_2 \\ \vdots \\ N_T \end{bmatrix} &\sim \text{MVNormal}( \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \end{bmatrix}, \textbf{K} ) \\ k_{i,j} &= \eta ~ \text{exp}( -\rho \text{d}^\text{2}_{i,j} ) \end{aligned} \tag{2.2}

The code below reproduces Figure 3 in the main article, in which four independent sets of time series are drawn from a common multivariate normal, whose covariance matrix is generated by a kernel with the parameter \eta=2 and \rho=.15.

Show the code

    # Simulate time series with Gaussian process
    GP = function(max_cov, rate, distMat, delta=0) {  # Gaussian kernel function
        m = max_cov * exp( - rate * distMat^2 )
        diag(m) = diag(m) + delta
        m
    }
    n_days = 49
    obs_times = seq(0, n_days, length=150)
    obs_dist = dist(obs_times, diag=T, upper=T) |> as.matrix()
    K = GP( max_cov=2, rate=.15, distMat=obs_dist )

    set.seed(10)
    plot(1, type="n", xlim=c(0,50), ylim=c(-5,5.2), xlab="Day", ylab = "N")
    for ( i in 2:5 ) {
        N = MASS::mvrnorm( 1, rep(0,length(obs_times)), K )
        lines( obs_times, N,  col=col.alpha(i,.3), lwd=2 )
        points( obs_times, N, col=col.alpha(i,.8), pch=c(15,17,23,19)[i-1] )
    }
    legend(-1.1, 5.3, paste("draw",1:4), col = 2:5, pch=c(15,17,23,19),
           lty=1, lwd=2, cex=1.1 )

### Linking P to N

As described in the main article, the mean of the multivariate normal in (2.2) can be used to model the long-term impact of P on N, as formulated in Eq. (4.1) in the main article and listed below:

\begin{align*} &\begin{cases} P_\text{t=0} ~~ = P_0 \\ P_{t + \Delta t} = P_t + r P_t (1 - \frac{P_t}{K}) \Delta t \\ \end{cases} \\ & \begin{bmatrix} N_1 \\ N_2 \\ \vdots \\ N_T \end{bmatrix} \sim \text{MVNormal}( - \gamma \begin{bmatrix} P^1 \\ P^2 \\ \vdots \\ P^T \end{bmatrix}, \textbf{K} ) \tag{4.1} \\ & ~~~~ k_{i,j} = \eta ~ \text{exp}( -\rho \text{d}^\text{2}_{i,j} ) \\ \end{align*}

The code below implements (4.1) and reproduces Figure 7 in the main article.

Show the code

    # Integrated form of logistic growth
    P_growth = function(t, r, K, P0)
        K / ( 1 + exp(-r*t)*(K-P0)/P0 )  # "K" is carrying capacity here
    P = P_growth(obs_times, r=.3, K=10, P0=.01)
    gamma = .3

    set.seed(10)
    plot(1, type="n", xlim=c(0,50), ylim=c(-7,10), xlab="Day", ylab = "")
    lines(obs_times, P, lwd=3, col=col.alpha(2))
    for ( i in 2:5 ) {
        N = MASS::mvrnorm( 1, -gamma * P, K )  # "K" is the covariance matrix here
        lines( obs_times, N, col=col.alpha(i,.3), lwd=2 )
        points( obs_times, N, col=col.alpha(i,.8), pch=c(15,17,23,19)[i-1] )
    }
    legend(-1.1, 10, paste("draw",1:4), col = 2:5, pch=c(15,17,23,19),
           lty=1, lwd=2, cex=1.1 )

## 3 Simulation

The simulation as described in the main article is conducted in two steps. First, the 98 observations of the N-process are simulated with R (R Core Team, 2022), as documented in the `get_params()` function in `model/utils.R`. Subsequently, using `cmdstanr` (Gabry, Češnovar, & Johnson, 2023) as the interface, the parameter values, as well as the simulated N-process, are passed in as data to Stan (Carpenter et al., 2017) to conduct the remaining simulation. The code for the second part of the simulation is available in `model/sim.stan`. We reproduce the simulation and Figure 4 of the main article below.

Show the code

    source("utils.R")
    set.seed(1024)

    # Set parameters for simulation
    #   see main article for a description of the parameters
    Params = get_params(
        n_days = 49,
        n_obs_times = 98,
        # P-process
        P0     =  .01 ,
        K      =   10 ,
        r      =  .3 ,
        # N-process
        max_cov =  2,
        rate    = .2,
        # Measurement of P-process
        a1     =  .78 ,
        b1     =  .744 ,
        c1     =   3.14 ,
        s1     =  .5 ,
        # Measurement of N-process
        a2     =  .6 ,
        c2     =   0 ,
        s2     =  .5 ,
        N      =  NULL,
        obs_times = NULL
    )

    # Run simulation
    sim = simulate(Params, "sim.stan")

    Running MCMC with 1 chain...

    Chain 1 Iteration: 1 / 1 [100%]  (Sampling)
    Chain 1 finished in 0.0 seconds.

Show the code

    # Plot Figure 4
    idx_obs = seq(1, sim$params$n_obs_times, by=1) |>
        as.integer() |> unique()
    with( c(sim$sim, sim$params) , {
        plot(1, type="n", ylim=c(-6, 11), xlim=c( 0, 49 ),
             xlab="Day", ylab="Measurement")
        # Process curve
        lines(obs_times, P, col=col.alpha(2,.8), lwd=2, lty=2)
        lines(obs_times, N+6.5, col=col.alpha(1,.65), lwd=2, lty=2 )
        # Measurement latent score
        lines( obs_times, M1, col=col.alpha(2,.5), lwd=2 )
        # Measurement observed score
        points( obs_times[idx_obs], P_meas[idx_obs], pch=19, col=col.alpha(2,.8) )
        legend(27,-1, c( TeX("Aversive Environmental Impact ($N_t$)"),
                         TeX("Cognitive Reappraisal ($P_t$)"),
                         TeX("Measurement (latent, $M_t$) "),
                         TeX("Measurement (observed, $P_t^{obs}$)")
                        ),
               col = c(col.alpha(1,.65), col.alpha(2,.8), col.alpha(2,.8), col.alpha(2,.8)),
               lty = c(2,2,1,NA),
               pch = c(NA,NA,NA,19),
               lwd = 3)
    })

## 4 Inference

Two statistical models were constructed following the description of the main article. The first model, which infers the parameter values from a single set of time series data, was fitted three times—the first was fitted with 49 observations and the second and third with 98 observations. The simulated observations for the first and second fit have a higher measurement error (\sigma = 0.5) whereas the simulated observations for the third fit have a lower measurement error (\sigma = 0.1). The second model, extended from the first model to include two parallel time series observations, was fitted with a total of 98 observations, with 49 observations for each time series. See `model/fit.R` for the code to fit the models.

### MCMC Diagnostics

Four MCMC chains were used for each fit, with 150 collected samples and 300 warmup samples for each chain. All fits, except for the fit with a lower measurement error, proceeded smoothly with no divergent transitions and maximum tree depth reached and exhibited well mixing of the chains, as indicated by the trace plots1. We list the diagnostic summaries for the four fits below:

Show the code

    m1.49    = readRDS("data/m1-49.RDS")     # 1st model with 49 obs.
    m1.98    = readRDS("data/m1-98.RDS")     # 1st model with 98 obs.
    m1.98.sd = readRDS("data/m1-98-sd.RDS")  # 1st model with 98 obs. & low measurement error
    m2       = readRDS("data/m2-49.RDS")     # 2nd model with 49 obs.

    for ( i in 1:4 ) {
        m = c(m1.49, m1.98, m1.98.sd, m2)[[i]]
        n = c("m1.49", "m1.98", "m1.98.sd", "m2")[i]
        cat("-------------\n")
        cat("[", n, "]\n")
        print( m$diagnostic_summary() )
    }

    -------------
    [ m1.49 ]
    $num_divergent
    [1] 0 0 0 0

    $num_max_treedepth
    [1] 0 0 0 0

    $ebfmi
    [1] 0.9107932 0.7029512 0.9116294 0.8402846

    -------------
    [ m1.98 ]
    $num_divergent
    [1] 0 0 0 0

    $num_max_treedepth
    [1] 0 0 0 0

    $ebfmi
    [1] 1.0852703 1.0762497 0.9705682 0.7346965

    -------------
    [ m1.98.sd ]

    Warning: 226 of 600 (38.0%) transitions hit the maximum treedepth limit of 10.
    See https://mc-stan.org/misc/warnings for details.

    $num_divergent
    [1] 0 0 0 0

    $num_max_treedepth
    [1]   0 106   0 120

    $ebfmi
    [1] 0.9060043 0.9156679 1.0691940 0.8909156

    -------------
    [ m2 ]
    $num_divergent
    [1] 0 0 0 0

    $num_max_treedepth
    [1] 0 0 0 0

    $ebfmi
    [1] 0.5914937 0.6104637 0.8680307 0.5642768

The third model `m1.98.sd` has a high percentage of MCMC sampler reaching maximum treedepth. This signals that the MCMC sampler cannot efficiently explore the posterior distribution, which can arise for many reasons and is idiosyncratic to the constructed model. Given that all other models have no difficulty for the sampler and the contradictory fact that imposing stronger constraints on the model (through more accurate observations) leads to a difficult posterior for the sampler, we infer that the task we wish the model to accomplish—inferring and separating two latent dynamics (the P- and N-process) from a single time series—is inherently hard or even impossible since the data and the assumed data-generating process do not strongly identify such latent dynamics.

### Posterior Predictions

In Bayesian inference, samples drawn from the parameters’ posterior distribution (i.e., posterior samples) are available after the model is fitted to the data, which can then be used for generating predictions that are compared to empirical data. The posterior distribution can be thought of as the fitted model’s belief of the parameters’ true values. To generate predictions according to the posterior, we simply pass in the posterior samples as data to the simulation described in Section 3. For each of our fitted models, 600 sets of posterior samples are available (150 sets \times 4 chains). The size of a set of posterior samples corresponds to the number of parameters in the model. To avoid cluttering the plot, we draw only 20 sets of posterior samples for plotting. See `model/post_predict_cache.R` for the code to collect posterior samples. We reproduce Figure 5 in the main article below.

Show the code

    source("post_predict.R", encoding = "UTF-8")

    Running MCMC with 1 chain...

    Chain 1 Iteration: 1 / 1 [100%]  (Sampling)
    Chain 1 finished in 0.0 seconds.

    Running MCMC with 1 chain...

    Chain 1 Iteration: 1 / 1 [100%]  (Sampling)
    Chain 1 finished in 0.0 seconds.

    Running MCMC with 1 chain...

    Chain 1 Iteration: 1 / 1 [100%]  (Sampling)
    Chain 1 finished in 0.0 seconds.

    Running MCMC with 1 chain...

    Chain 1 Iteration: 1 / 1 [100%]  (Sampling)
    Chain 1 finished in 0.0 seconds.

## References

Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., … Riddell, A. (2017). Stan: A Probabilistic Programming Language. Journal of Statistical Software, 76, 1. https://doi.org/10.18637/jss.v076.i01
Gabry, J., Češnovar, R., & Johnson, A. (2023). Cmdstanr: R interface to ’CmdStan’ [Manual].
R Core Team. (2022). R: A language and environment for statistical computing [Manual]. Vienna, Austria: R Foundation for Statistical Computing. Retrieved from https://www.R-project.org/

## Computing Environment

All packages except `stom` are available from CRAN. To install `stom`, simply run:

    install.packages("remotes")  # if not installed
    remotes::install_github("liao961120/stom")

### Used Packages

    cmdstanr::cmdstan_version()

    [1] "2.32.1"

    sessionInfo()

    R version 4.1.3 (2022-03-10)
    Platform: x86_64-w64-mingw32/x64 (64-bit)
    Running under: Windows 10 x64 (build 19045)

    Matrix products: default

    locale:
    [1] C
    system code page: 950

    attached base packages:
    [1] stats     graphics  grDevices utils     datasets  methods   base

    other attached packages:
    [1] latex2exp_0.9.6 deSolve_1.35

    loaded via a namespace (and not attached):
     [1] compiler_4.1.3       pillar_1.9.0         cmdstanr_0.6.0
     [4] tools_4.1.3          digest_0.6.31        jsonlite_1.8.4
     [7] evaluate_0.21        lifecycle_1.0.3      tibble_3.2.1
    [10] checkmate_2.1.0      gtable_0.3.3         pkgconfig_2.0.3
    [13] rlang_1.1.0          cli_3.6.1            rstudioapi_0.14
    [16] parallel_4.1.3       yaml_2.3.7           stom_0.0.0.9000
    [19] xfun_0.39            fastmap_1.1.1        withr_2.5.0
    [22] dplyr_1.1.2          stringr_1.5.0        knitr_1.42
    [25] generics_0.1.3       vctrs_0.6.1          htmlwidgets_1.6.2
    [28] systemfonts_1.0.4    tidyselect_1.2.0     grid_4.1.3
    [31] data.table_1.14.8    svglite_2.1.1        glue_1.6.2
    [34] R6_2.5.1             processx_3.8.1       fansi_1.0.4
    [37] distributional_0.3.2 rmarkdown_2.21       tensorA_0.36.2
    [40] farver_2.1.1         ggplot2_3.4.2        posterior_1.4.1
    [43] magrittr_2.0.3       ps_1.7.5             backports_1.4.1
    [46] scales_1.2.1         htmltools_0.5.5      MASS_7.3-55
    [49] abind_1.4-5          colorspace_2.1-0     utf8_1.2.3
    [52] stringi_1.7.12       munsell_0.5.0

## Footnotes

1. To inspect the trace plots, run the code in `model/diagnose.R` for each of the fitted models by modifying the path to the fitted model (`model/data/m*.RDS`) at the beginning of the script.↩︎