Intraindividual time-varying dynamic network of affects: linear autoregressive mixed-effects models for ecological momentary assessment

An interesting recent development in emotion research and clinical psychology is the discovery that affective states can be modeled as a network of temporally interacting moods or emotions. Additionally, external factors like stressors or treatments can influence the mood network by amplifying or dampening the activation of specific moods. Researchers have turned to multilevel autoregressive models to fit these affective networks using intensive longitudinal data gathered through ecological momentary assessment. Nonetheless, a more comprehensive examination of the performance of such models is warranted. In our study, we focus on simple directed intraindividual networks consisting of two interconnected mood nodes that mutually enhance or dampen each other. We also introduce a node representing external factors that affect both mood nodes unidirectionally. Importantly, we disregard the potential effects of a current mood/emotion on the perception of external factors. We then formalize the mathematical representation of such networks by exogenous linear autoregressive mixed-effects models. In this representation, the autoregressive coefficients signify the interactions between moods, while external factors are incorporated as exogenous covariates. We let the autoregressive and exogenous coefficients in the model have fixed and random components. Depending on the analysis, this leads to networks with variable structures over reasonable time units, such as days or weeks, which are captured by the variability of random effects. Furthermore, the fixed-effects parameters encapsulate a subject-specific network structure. Leveraging the well-established theoretical and computational foundation of linear mixed-effects models, we transform the autoregressive formulation to a classical one and utilize the existing methods and tools. To validate our approach, we perform simulations assuming our model as the true data-generating process. By manipulating a predefined set of parameters, we investigate the reliability and feasibility of our approach across varying numbers of observations, levels of noise intensity, compliance rates, and scalability to higher dimensions. Our findings underscore the challenges associated with estimating individualized parameters in the context of common longitudinal designs, where the required number of observations may often be unattainable. Moreover, our study highlights the sensitivity of autoregressive mixed-effect models to noise levels and the difficulty of scaling due to the substantial number of parameters.


INTRODUCTION
We consider simple directed networks comprising two moods and one external factors node.The two moods interact mutually by enhancing (blue) or dampening (red) each other under the unidirectional influence of external factors.
Assuming intensive longitudinal data through Ecological Momentary Assessment (EMA), we formalize the mathematical representation of these networks, Figure S1, by an Exogenous Linear Autoregressive Mixed-effects model (LARMEx).The autoregressive coefficients represent the interactions and the external factors are treated as exogenous covariates.We let every parameter in the model have fixed and random components aiming at networks that are allowed to have variable structures over reasonable units of time like days or weeks depending on the study design.Harnessing the rich theoretical and computational literature that furnishes classic linear mixed-effects models, we transform this formulation to a classical one and use the existing methods and tools.We assume our model is the true data generating process and simulate data using a predefined set of parameters.Then we investigate the performance and feasibility of this approach in delivering reliable estimates for different choices of the number of observations and the intensity of noise.

MATHEMATICAL REPRESENTATION
Linear autoregressive (AR) processes are used to model stochastic dynamical systems of discrete data through difference equations.That is, the current state of the system, s t , is assumed to depend linearly on the preceding states, s t−1 , . . ., s t−p and a constant term, β 0 , analogous to the intercept in linear regression.All the unknown factors which might affect the system are encapsulated in the error term ϵ t and ϵ t is usually assumed to be white noise with zero mean (Hamilton, 1994), In this formulation, s is usually supposed to be a time series comprising sufficient observations over time, which makes it possible to detect its temporal dynamics.The autoregressive coefficients β i determine how the current value of s is affected by its preceding values by p steps back in time.In the simplest form of an AR process, the dependence on previous states is only taken one step back and the so-called Lag(1) processes are formulated by s t = β 1 s t−1 + β 0 + ϵ t .In the absence of noise, this process has a fixed point of If such a process reaches its fixed point it will stay there forever.Otherwise, such a point determines the asymptotic behavior of the process in the long run.In this work we only consider stable systems that will eventually converge to fixed points and require that |β 1 | < 1 for a single variable Lag(1) process.
The immediate generalization of AR processes to accommodate multiple states is a common practice whenever one needs to study the dynamics of several variables simultaneously.Exogenous covariates might also be included in a linear fashion to account for the effects of known external factors.
LARMEx models the mean response as a combination of trait-like characteristics of an individual as fixed effects (FE) and state-like features varying over days (weeks) as random effects (RE).Needless to say, existing techniques allow for the estimation of fixed effects and the variances of random effects, and the values of random effects are predicted using these estimations.This way, the model accounts for the dependency in different levels of data, e.g. during a day, and at the same time compensates for the lack of statistical power due to few observations on the first level (Funatogawa and Funatogawa, 2018;Fitzmaurice et al., 2011).
We build a dynamic affect network comprising two interacting nodes (e.g., positive and negative mood in simplification) and an exogenous factor (E) representing an emotionally meaningful event or situation that may be considered as the input for the dynamical system.

Mood Networks
Let m i,t = [m 1 , m 2 ] T i,t be the 2 × 1 vector corresponding to the two mood values from the ith measurement day at any observation occasion t = 1, 2, . . ., n i where n i is the number of observations per day and T denotes the transpose of a matrix.Note that n i might differ for every i despite the fact that it is planned to collect the same amount of data every day.In the absence of external driving forces like pleasant or unpleasant events and noise, the mood network (Figure S2a) is represented by a first-order autoregressive system in which B ar i = β ar + b ar i are 2 × 2 matrices of auto-and cross-lagged regression coefficients, and B c i = β c + b c i are 2 × 1 vectors of constant terms corresponding to the two moods.Every coefficient matrix and constant term is a sum of fixed-and random-effects terms representing the subject specific, β, and day specific, b i , components.The random effects are assumed to be normally distributed with the mean zero and a certain variancecovariance, b ∼ N(0, G).Different structures of G are used to account for possible scenarios representing the dependencies amongst random effects.We restrict ourselves to a compound symmetric (CS) type where all random effects share the same variance, the diagonal elements of G, and the same covariance, off-diagonal terms.Furthermore, autoregressive and exogenous components are assumed to be independent.For more structures we refer the reader to (Funatogawa and Funatogawa, 2018) and references therein.
To avoid the complexity of restricting moods to take only positive values, we assume that they are centered and can take any value.Nevertheless, we set the initial values, constant terms and the range of exogenous factors in a way that for reasonably low noise levels, most of the trajectories remain in the interval [−1, 1].The dynamics represented by Equation S2 is asymptotically stable given that all the eigenvalues of B ar i have modulus less than one, Figure S2b (Lütkepohl, 2005).The asymptotic behavior is determined by the autoregressive coefficients and constant terms.Note that the coefficients take positive and negative values so the moods interact through activating and suppressing each other, although in a more complex networks there might be only positive or negative interactions.We also assume that the constant term has a zero fixed-effect component to reduce the number of parameters.In dealing with real data we put no constraint on the coefficients to account for any possibility.In Figure S2b, depending on the magnitudes of coefficients, different patterns of decay are seen.For larger values, variables take a longer time to relax toward the asymptotic level.

Network with Exogenous Factors
Assuming that a collection of external factors, E, acts on mood nodes and its effect could vary over days, the network of moods as discussed above could be extended in two ways (Figure S3a).On the one hand, we could assume that on every day a subject is exposed to a baseline level of input, e.g.stressors, which might be different for every day.This has already been treated by adding a constant term, B c i = β c + b c i as fixed and random intercepts for day i, Equation S2.On the other hand, one could assume that E has a timedependent component as well, which is measured along with the moods.For instance, a day could be marked by various unpleasant events or alternating positive and negative events.This term enters as a time-varying covariate to the formulation in Equation S2 and forms an exogenous autoregressive process of In this formulation, B e i could be considered as a measure of contemporaneous reactivity to external events in the long run and short periods of time, like days.The term B e i e i,t + B c i determines how the asymptotic values of moods are affected by external amplifiers.Since E varies over time, the asymptotic levels are changing from time to time and, therefore, a fluctuating pattern emerges over the mood trajectories.In this way, the differences in exposure to stressors, for instance, are accounted for.The stressor is measured longitudinally alongside moods.In Figure S3b below, we included E as an impulse, i.e., it has a nonzero value at t = 5 and is zero elsewhere.This is done only for presentation purposes, and in practice, it might take any value at different points.Opposite perturbations in the mood levels are visible compared to Figure S2b.In a more realistic scenario, one would introduce a positive random value at each time point with a non-negative value as the fixed effect and generate random effects from a truncated normal distribution to keep the elements of e i,t positive.If this was the case an overall upward (downward) profile would be visible in data indicating the exciting (inhibiting) effect of the stressors on negative (positive) mood.

Unknown Effects as Noise
So far we considered only deterministic models where the dynamics of the system is fully known once the parameters are fixed.However, our knowledge of real systems is usually imperfect either under the influence of measurement errors or lack of information about unknown influences on the dynamics.To account for all the uncertainties due to unknown factors, but not measurement errors, we add a white noise element, ϵ i,t , to Equation S3.Treating measurement errors adds another layer of complexity to the model and hence is avoided here for the sake of simplicity.Then the full model of mood networks in the presence of external factors as an exogenous linear autoregressive mixed-effects model becomes in which ϵ i,t ∼ N(0, σ 2 I 2×2 ).The choice of a diagonal form for ϵ i,t guarantees the model parameters be identifiable.However, the total noise in the system depends also on the covariance matrix of the random effects, G in Equation ??.Given that in practice there is no restriction on G, any potential contemporaneous correlation between different mood variables, as empirical EMA data usually exhibit (Epskamp, 2020), could occur as well.
Taken E into account, there might be an increasing variance and divergent profiles up-and downwards.Adding this noise to the system results in Figure S4.

DATA GENERATING PROCESS
Let m i,t = [m 1 , m 2 ] T i,t be the 2 × 1 vector corresponding to the mood values from the ith day at any observation occasion t = 1, 2, . . ., n i where n i is the number of observations per day and T denotes the transpose of a matrix.
Assuming that a collection of external factors (E) act on mood nodes and their effect could vary between days, the evolution of m i,t is represented by a LARMEx model as in which β and b represent the fixed and random effects (FE and RE henceforth).The temporal (Granger-causal) connections between the moods are governed by the lag(1) autoregressive term, (β ar + b ar i )m i,t−1 .The remaining terms represent the contemporaneous effects of E and the constants (intercepts).For simplicity, we assume no connections from moods to E.

GENERATING A SET OF KNOWN PARAMETERS
For simulating data with this model, we need to specify the parameters.To achieve replicability, one can predefine a random number generator (RNG).In Julia, we implement this using the following struct.
These matrices are generated using keyword arguments to the constructor of parLARMEx.
It is also possible to initialize these with custom matrices or random values by passing [].
The variance-covariance of the random effects is built using b_var to be This is not the final variance-covariance of the RE, because one needs to control for the stability of the autoregressive process.To this end, we sample a numger of nSamp = 20000 sets from a multivariate normal distribution of MVN(0, G) and retain nL2Max = 100 sets for which the absoloute eigenvalues of β ar + b ar are less than one.These are stored in b_ar, b_e, b_c, and the variance-covariance of this latter sample as b_cov.
This process is accomplished by calling the function genRE_CS(rng,G,B_AR,maxSbj,nSamp) inside the constructor which updates the fields corresponding to the random effects and their variance-covariance structure.In what follows we import these and other function definitions from the file helperSim.jlon GitHub, and then construct the parameters using a replicable RNG: [1]: include("helperSim.jl"); [2]: rng = MersenneTwister(1984); par = parLARMEx(rng=rng);

GENERATING DATA
We assume that the data generating process has a two-level structure.The multiple observations are made on level I and these are nested in level II units.For EMA data, level I could be the daily observations, as many as nObsL1, which are nested in single days with a total number of nL2.These characteristics together with the varianc-covariance of noise, sigma, initial values for days, M0, and the vector of known exogenous factors, E, are stored in another struct.The following code snippet setts up the configuration of the desired data set.The RNG is carried over to perovide replicablity.Finally, genData() gives the simulated data along with the noiseless data as dataframes and the signal-to-noise-ratio (SNR).

TRANSFORMING TO CLASSICAL MIXED-EFFECTS MODEL
By stacking mood values, parameters and covariates for every day separately, and forming the design matrices for fixed and random effects, X and Z, data from a single day takes the following matrix form which is equivalent to a linear mixed effects formulation.
In this formulation for a network of two moods and one external nodes one has and In its general form, the above equation for k moods, Y i is a k(n i − 1) × 1 vector of mood values, β is a (k 2 + 2k) × 1 vector of fixed effects, b i is a (k 2 + 2k) × 1 vector of random effects, X i and Z i are k(n i − 1) × (k 2 + 2k) design matrices.The random effects b i and residuals ϵ i are assumed to be independent with a multivariate normal (MVN) distribution of This step is done by the following function which gives another data frame suitable for performing the estimation.Here, it is possible to introduce missing values to the analysis, e.g.miss = .2for a 20% of missing values at random.By default, it is assumed that there is no missing values, miss = 0.The code is available in helperSim.jlon GitHub.

SETTING UP THE FORMULA FOR FITTING
The prepared data has more columns representing the response variable M, and the predictors including the connections in the network M11, M12, M21, M22, E1, E2, as well as the constant terms C1, C2.It also contains the level I time points tL1, and the level II identifiers idL2.We provide a function in helperSim.jl,on GitHub, setFormula(fitData) which constructs the formula suitable to feed in the Julia package MixedModels.jl.This function is restricted only to a dataframe which has the exact columns as shown above.

FITTING
We use the MixedModels.jlpackage in Julia to perform the estimation.It is similar to the lme4 in R but exhibits faster performance in our case.This can be done using the GitHub file, simBootstrap.jl.It uses the following function to generate data repeatedly in a loop.loopSim (csvDir,M0_max,sigma,n,P) This function creates a folder called bootstrap in the current working directory and for each noise intensity, sigma, saves the true and estimated parameters of each simulation in different folders, re and reh, based on the number of days as CSV files.Fixed effects and the variances of random effects are collected in two separate CSV files.
The code could take some time, depending on the computational power.Using the estimated parameters one can construct bootstrap confidence intervals for fixed effects and the variance of random effects as well as the relative errors of the estimations.For details regarding these analyses, please refer to the main text.We presented an extension of linear mixed-effects models by adding an autoregressive component to model the network representation of mental state.Specifically, we assumed a simple network of two causally interacting moods under the unidirectional influence of an external-factors node.This framework is suitable for intensive longitudinal data.

Figure S1 :
Figure S1: Nodes M1 and M2 interact with each other leading to temporal correlation, meaning that the value of the source node (variable) at time t has an impact on the value of its target at time t + 1. Looped arrows indicate autoregressive effects.Exogenous factors influence moods at the same time in a one-way fashion.Blue and red arrows indicate activation and suppression respectively.Different thicknesses stand for distinct interaction strengths.Dashed arrows indicate contemporaneous effects while solid ones are indicative of possible temporal causalities.

Figure S2 :
Figure S2: Tow interacting nodes.Network of two nodes M1 and M2, e.g.positive and negative mood, without external forces (a), realizations for 5 simulated days of 10 observation each (b).Data is generated using Equation S2 with initial values for different days sampled uniformly from the interval [−0.5, 0.5].

Figure S3 :
Figure S3: Tow interacting nodes with external factors.Network of two nodes M1 and M2, e.g.positive and negative mood, with a fixed amount of negative external input E acting on both at time point 5 (a), realizations for 5 simulated days of 10 observation each (b).An analogous model could be postulated for positive events or situations as inputs.

Figure S4 :
Figure S4: Tow interacting nodes with external factors and noise.Realizations from a network of two nodes M1 and M2, e.g.positive and negative mood, for 5 simulated days of 10 observation each.External factors act on both nodes at time point 5 in the presence of noise.