Learning delay dynamics for multivariate stochastic processes, with application to the prediction of the growth rate of COVID-19 cases in the United States.

Delay differential equations form the underpinning of many complex dynamical systems. The forward problem of solving random differential equations with delay has received increasing attention in recent years. Motivated by the challenge to predict the COVID-19 caseload trajectories for individual states in the U.S., we target here the inverse problem. Given a sample of observed random trajectories obeying an unknown random differential equation model with delay, we use a functional data analysis framework to learn the model parameters that govern the underlying dynamics from the data. We show the existence and uniqueness of the analytical solutions of the population delay random differential equation model when one has discrete time delays in the functional concurrent regression model and also for a second scenario where one has a delay continuum or distributed delay. The latter involves a functional linear regression model with history index. The derivative of the process of interest is modeled using the process itself as predictor and also other functional predictors with predictor-specific delayed impacts. This dynamics learning approach is shown to be well suited to model the growth rate of COVID-19 for the states that are part of the U.S., by pooling information from the individual states, using the case process and concurrently observed economic and mobility data as predictors.

Introduction

The modeling of time dynamical systems is of interest in multiple scientific fields. While ordinary differential equations (ODE) have a long history, interest in delay differential equations (DDE) is more recent, albeit both have been extensively studied [1], [3], [22], [26], [61]. A DDE is a natural extension of an ODE when observed processes have an aftereffect. Unlike the situation for ODEs, the solution of a DDE depends not only on the initial state but on the entire history of the process in a time interval of length equal to the delay prior to the initial time point.

Moving beyond the basic notion of a deterministic ODE, random differential equations (RDE) [23], [46], [54], [55] and stochastic differential equations (SDE) [2] are used to accommodate probabilistic uncertainty in temporal stochastic processes. In this paper we focus on RDEs, where the random effects are manifested in the model parameters. These include coefficients, initial conditions and forcing terms that are typically smooth in time, leading to differentiable sample path solutions of the RDE, in contrast to the situation for SDEs, which have a non-differentiable forcing term, typically a Wiener process, and therefore require stochastic calculus [34].

Extensive developments in the theory of the forward problem of obtaining solutions for RDEs stand in contrast with statistical approaches, which include data-oriented methodology for the inverse problem, that is, learning the nature of the differential equation from data. This approach has been referred to as empirical dynamics [45]. The starting point is a sample of random trajectories that are viewed as independent and identically distributed (i.i.d.) realizations of an underlying smooth (continuously differentiable) stochastic process. For a smooth stochastic process , there always exists a function f such that with almost surely. This forms the basis of empirical dynamics and has led to both parametric and nonparametric modeling approaches for f [57], [66], which can be characterized as dynamics learning from data. Due to the minimal assumptions, dynamics learning covers a large number of specific ODEs that do not need to be a priori specified, however it requires that a sample of repeated realizations of the underlying stochastic process is available. Such samples form the backbone of functional data analysis [48], [58]. A pertinent example is provided by COVID-19 caseload trajectories, where one observes samples of such trajectories across geographic units such as states or countries [18].

Related but quite distinct statistical inference concerns the problem of identifying the parameters of an a priori specified dynamic system [12], [20], [37], usually just from one or very few realizations. Here the trajectories are considered non-stochastic but are typically assumed to have been measured with noise. This leads to a curve fitting problem that in simple cases can be addressed with nonlinear least squares, but for which also nonparametric and semiparametric statistical methods have been employed [38], [47]. A key difference is that in dynamics learning one requires data that can be considered as an independent and identically distributed sample of realizations of an underlying smooth stochastic process , while the parameter identification problem is often addressed given (noisy) data from one trajectory. The modeling approach that we introduce here follows the paradigm of dynamics learning, albeit in a somewhat more structured framework than empirical dynamics.

The general form of DDEs for vector functions is where is the function value or state at time t and can be either a vector of the states evaluated at discrete time delays , , so that , or alternatively an integral of the trajectory over a past period, representing a continuum of delays often called distributed delays [4], [5], [6], [8], [15], [27], [31], [35], [41]. Results on the existence and uniqueness of random differential equations with delay (RDED) are relatively recent, where [14] studied conditions for existence and uniqueness of the solutions of a RDED and [24] the specific case of a linear RDED with a forcing term.

In this paper we propose dynamics learning from a sample of multivariate functional data, where the model generating the observed derivative trajectories is assumed to be a RDE with a delay component. The randomness in the RDE is included in the forcing term, part of which is explained by additional covariates. These are stochastic processes with individual delay components, which in the motivating COVID-19 application correspond to trajectories of mobility and economic activity. The unexplained remainder is a drift process, which appears as in equation (2).

We utilize tools from the theory of RDE to provide a foundation for the proposed population models and establish that all the models described in Section 2 correspond to RDED with unique solutions, which depend on the model parameters. The goal is then to learn these unknown parameters, which govern the dynamics encapsulated in the population RDED, by pooling information across the sample of observed trajectories. The presence of a large number of covariates and delay components to choose from gives rise to the challenge of model selection. To address this we propose an initial pruning step as described in Section 4.2 followed by a backfitting step to optimize over a set of delay components as outlined in Section 4.3. The proposed approach differs from existing optimization strategies [41], [59], [60], [65] and statistical methodology [36] for parameter estimation that has been previously deployed to identify a DDE based on noisy data observed for a single non-stochastic trajectory. To our knowledge, dynamics learning where one has samples of stochastic processes has so far not been explored for the case of an underlying RDED, even for the case of a one-dimensional stochastic process. We furthermore demonstrate in this paper that such an approach is well suited for modeling the growth rates of COVID-19.

We consider both discrete and distributed delay models and establish existence and uniqueness of the solutions in the sense [14] in Section 3. For the case of discrete delays we harness functional concurrent regression models [13], [33], [33], [51], [53] and for distributed delays the functional history-index model [40], [52], which incorporates a range of recent past values of the process. For dynamics learning of RDEDs, we adopt a two stage procedure. In a first step we utilize functional linear regression [17], [42], [63] with history index to learn the distributed delay, where the regression parameter function then corresponds to a history index function for the process of interest. In a second step the resulting linear predictor, which is the inner product of the history index function and the predictor process of interest, is used as a predictor along with additional covariate processes with covariate-specific delays in a concurrent model to fit the derivative process.

We apply this model to predict the time-dynamic growth rate of COVID-19 for individual states in the United States, where the sample of case trajectories of all states provides the sample of trajectories that is the starting point for the proposed dynamics learning. Predictor processes that we consider include the cumulative case (caseload) process, daily economic indicators and changes in mobility patterns. Modeling the time evolution dynamics of the COVID-19 pandemic is of great importance to understand and interpret the underlying associations as well as for deploying resources and formulating policies in the face of great uncertainty [7], [11], [32]. The proposed methodology is also of interest to assess the dynamics of many other empirically observed complex multivariate stochastic processes for which samples of observed trajectories are available. We show by means of leave-one-out predictions that by employing dynamics learning for the proposed RDED model with delay components one obtains considerably more accurate time-dynamic growth rate predictions for COVID-19 caseload curves compared to models without delays, demonstrating the importance of considering the inclusion of lags when modeling the growth rate of COVID-19.

Proposed models

Let denote a multivariate stochastic process where is a continuously differentiable process of interest, is a vector function of additional covariates, and is a time window of interest. Consider a RDED with discrete delays, where g corresponds to the initial condition stochastic process, the , are discrete delays, and are smooth functions whose regularity will be specified below in Section 3. In the above, is a random drift process that is independent of .

While the inclusion of discrete delays is an extension of the classical functional concurrent regression model [48], the functional linear regression model with history index [40], [52] might be a better choice when the derivative trajectories depend not only on the predictor process value at a single past instant but on the entire continuum in the recent past. We model these distributed delays as where g again is an initial condition process. For the purpose of illustration and technical derivations, we assume that is a univariate process in (4); the corresponding multivariate generalization is straightforward.

For COVID-19 modeling a hybrid model that combines distributed delay in with discrete delays for all other predictors turned out to be particularly suitable. In this practically relevant variant of a RDED one postulates that for all , with an initial condition .

We remark here that the RDED models proposed above are linear in the model parameters, although general RDED models determined by nonlinear functions are possible and might provide greater modeling flexibility, at the cost of additional modeling complexity. For example, it is not difficult to develop population models for quadratic and polynomial versions, in analogy to functional polynomial models [62]. Modeling with linear RDEDs has the advantages that such models are easy to apply and implement, come with excellent interpretability of the model parameters, and exhibit good empirical performance in our application. For the rest of the manuscript we therefore only consider linear RDEDs and leave the development of nonlinear RDEDs for future research.

Existence and uniqueness of solutions

Functional concurrent model with discrete delays

Here we discuss the existence and uniqueness of the solutions of the RDEDs corresponding to the concurrent model (3) and the history-index functional linear model in (4), as well as the hybrid model in (5) used in the data application.

Consider a complete probability space , where is the σ-algebra on Ω and P is a probability measure. Denote the space of random variables with finite moment by , that is, for if . A sequence of real-valued random variables is said to converge to a random variable Y in the moment if as . We define a process to be moment continuous if for any . The notions of -moment differentiability and -moment Riemann integrability are defined similarly (see [54] Chapter 4, pages 92 and 100).

We say that a stochastic process is a solution of an RDED in the sense on the interval if X is -moment differentiable on I, -moment continuous on , and X satisfies the RDED including the corresponding initial condition [14]. We require the following assumptions on the coefficient functions and predictor processes in model (3).

) and are continuous on .

is -moment continuous on .

The predictor processes satisfy and are -moment continuous on , for .

Calatayud et al. [14] showed the existence and uniqueness of the solution of the following general form of a random delay differential equation with discrete delay , given by and established the following result.

If r satisfies the Lipschitz condition , for , , and a k with , then the random delay differential equation (6) has a unique solution in the sense.

The concurrent model in (3) is a special case of the random delay differential equation (6) with a specific choice of r. A key step is to show that the function f defined in (3) is -moment continuous. The existence and uniqueness of the solutions of (3) can then be obtained following similar arguments to those in the proof of Theorem 2.2 in [14]. We summarize this result in Proposition 1. All proofs are in the Appendix.

When, under the regularity conditions(A0)–(A2)and assuming that the coefficientsare continuous on, the functional concurrent DDE model given by(3)has a unique solution in thesense provided that the initial conditionis-moment continuous.

Proposition 1 requires the -moment continuity of on for the uniqueness of the solution in the sense. If almost surely, then there exists a version of the -solution which solves the RDED in the sample path sense, see Theorem 2.3 in [14] and [16]. Note that in (3), we do not include the process itself as one of the predictors; included is only the delayed version , where . Hence solutions in the sample path sense do not require g to be differentiable almost surely. For , see Remark 3.

The explicit sample path solution in the special case of a linear RDED has been derived in [14], [16], [24]. The linear RDED often includes random coefficients and a random time-varying forcing term. In (3), the terms involving the predictor processes , the drift process , and the intercept coefficient function form the constituents of the random forcing term; the varying coefficient functions are non-random. Using the method of steps as outlined in the proof of Theorem 2.3 in [14], we can provide a solution to (3) in the sample path sense as follows. For , (3) is equivalent to an ordinary RDE given by where f is defined as in (3). Then, the solution of (7) is given by Next, for , Solving (8) leads to Repeating this argument until the entire domain is covered one obtains the sample path solution. Due to the existence and uniqueness of the solution of (3) in the sense and by Remark 1, the sample path solution obtained as above is unique.

In applications, the delay parameter can take the value zero and when in (3) one arrives at an ordinary RDE without delay. For this case, existence and uniqueness of solutions has been well studied (see, e.g., Theorem 5.1.1 in [54] and [55]). The sample path solution in this case can be obtained via the integration factor method and is given by where with f defined as in (3). Here consists of the intercept term , the predictor processes along with the corresponding coefficient functions , and the drift process .

Functional linear model with history index for distributed delay

We start with the simplest history index model with only one predictor process , characterized as a first order RDED described in (4). We prove the existence of a unique solution of (4) by adopting similar arguments as those behind Proposition 2.1 and Theorem 2.2 of [14]. The extension to multivariate predictor processes is straightforward. Observe that one can rewrite model (4) as where with the initial condition .

Assume that is -moment continuous on and define such that for a -moment continuous process and , , Here is a random functional whose first argument is the trajectory and whose second argument is . The following proposition with proof in Appendix B provides the continuity of in t in the sense. For the random process X on and the time variable , the history index model (9) can be expressed as For the existence of a unique solution of (11) we need the following assumptions, in addition to (A0)–(A2):The following result characterizes the solution of (11), where the integral in condition (b) of the following Proposition is defined in the -moment Riemann integral sense [for definition of -moment Riemann integrability see [54], Chapter 4 page 100].

Supposeis-moment continuous for. Thenas defined in(10)is continuous in t in thesense on the domain.

γ is continuous on and is continuous on .

on and is -moment continuous on .

The processis asolution of(11)if and only if for,

X is-moment continuous on

for each

for each.

The following result establishes the existence and uniqueness of a solution for (11). The proof extends the classical Picard theorem for deterministic ODE to RDED in the sense, via the Banach fixed-point theorem, following a similar line of arguments as in [14] (see Appendix B). For the existence and uniqueness of the solution of the hybrid model in (5), which is the most promising of the models considered for modeling the growth rate of COVID-19, the following corollary applies.

Assume the regularity conditions(A0)–(A2)and(B0)–(B1). Then forwe have thatand a solution of(4)exists and is unique in thesense.

Suppose the regularity conditions(A1),(A2)and(B1)hold. Ifand the coefficientsare continuous on, and γ is continuous on, then the model in(5)has a uniquesolution.

Learning dynamics from samples of multivariate stochastic processes and practical considerations for COVID-19 caseload modeling

Obtaining derivatives

In data-driven differential equation modeling, given a sample of observed processes, a necessary first step is the estimation of derivatives. In applications such as the COVID-19 case trajectories the available data are usually noise-contaminated, and therefore some care is needed to obtain viable derivative estimates. The situation for the observed data for one of the trajectories is reflected by the nonparametric regression model where are the grid points where the process X is measured and are the measurement errors, which typically are considered to be independent mean zero random variables with finite variance.

In our data application, the available data correspond to the cumulative case counts of COVID-19, which we consider to be noisy realizations of a smooth underlying function that is observed on a daily grid. For the practical estimation of the derivatives one has various choices that all require a tuning parameter. Since in preliminary studies local polynomial fitting yielded the smallest leave-one-out prediction error relative to the difference quotients, which serve as nearly unbiased targets, we opted to obtain derivatives by applying local polynomial regression [28], [44] as per Appendix C, where further details are provided.

Target model and initial lag and covariate process selection

For modeling the COVID-19 caseload process X for the states of the U.S. we adopt the RDED model (5), where the time delay of the effect of the process X at time is distributed on the interval , and are discrete time delays or time lags for the covariate or predictor processes . Our starting point is that the data represent n independent realizations of the stochastic process , where each realization (corresponding to a state) is indexed by i and n is the number of states included in the analysis.

Given , to estimate the history index weight function in the historical model without covariates, we preform a function-to-scalar linear regression for each grid point , i.e., This can be implemented using function FLM in the R package fdapace [19]. With estimates of for each k in hand, the estimate of the history index weight function, for all , is obtained using a two dimensional kernel smoother.

For each of the covariates , we conduct an initial lag selection based on functional varying-coefficient regression including the corresponding variable as a single predictor, considering the models where . Given a lag τ, for each subject i and time point , we obtain the leave-one-out prediction where and are estimated by fitting model (14) by ordinary least squares, excluding the data from subject i. The lags are chosen to minimize the leave-one-out prediction error across subjects,

With initial choices of the lags and the estimate obtained from model (13), for each , variable selection is then performed using LASSO [56], where . This can be implemented using the R package ncvreg [9], [10], where the tuning parameter λ is chosen by leave-one-out cross-validation. Considering all K time points simultaneously, we obtain proportions of the number of time points where each variable is selected, i.e.,

Subsequently, selection among the time-varying predictors other than the process X is performed by applying a threshold on these proportions . This threshold is chosen by leave-one-out cross-validation,

Lag selection by backfitting

Let denote the set of selected predictor processes obtained using the method described in Section 4.2. Updated lag parameters are then chosen for the variables in through an iterative backfitting algorithm such that the total leave-one-out prediction error is minimized. Given the iteration of the lag vector, one obtains the update for the lag by minimizing the leave-one-out integrated mean squared prediction error, i.e., where for

Here the functions represent the fitted coefficients corresponding to the iteration of the lag, with the observation held out for prediction. This process is iterated over until the lag vector converges. Lags are initialized at the vector , where are as described in Section 4.2 and the stopping rule is such that the algorithm terminates if no lags have changed after a complete cycle of updates for all predictors. After the final lags are estimated, we lightly smooth the estimated coefficients with local linear smoothing, which is a standard procedure in concurrent functional regression modeling [25], [29]. Specifically, the initial estimates are smoothed to obtain the effect functions presented in Fig. 3, Fig. 4, defined as where with representing the Epanechnikov kernel with a bandwidth of days.

Fig. 3

Predicted growth rates by the proposed RDED (5) (red, dashed) and by a RDE without lags (blue, dotted) for US States, Alabama–Montana. Observed rates are the black solid curves. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

Fig. 4

Predicted growth rates by the proposed RDED (5) (red, dashed) and by a RDE without lags (blue, dotted) for US States, Nebraska–Wisconsin. Observed rates are the black solid curves.

Random differential equation with delay fitting and prediction for COVID 19 case trajectories in the United States

Data description

We obtained daily confirmed cases across states in the United States from the COVID-19 Data Repository by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University, which is publicly available at https://github.com/CSSEGISandData/COVID-19 and was accessed on September 19, 2020. The positive and negative COVID-19 test results per day and state were obtained from the COVID Tracking Project and are publicly available at https://COVIDtracking.com/ (accessed on September 19, 2020). The former variable refers to the number of people with confirmed or probable case of COVID-19 (see https://COVIDtracking.com/about-data/data-definitions for more details) and the latter refers to the total number of people with a completed negative PCR test. We used effective positivity rates (EPR), defined as the ratio of effective positive results (positive test or probable case) to the effective total tests (positive or probable case plus negative test results).

We also obtained features from the Opportunity Insights Economic Tracker Data [21], which is publicly available at https://github.com/OpportunityInsights/EconomicTracker (accessed on September 19, 2020). This database contains daily and state level information for several economic activity indicators such as credit and debit card expenditure across different types of activities, some of them stratified by zip codes with low/medium/high income level. We refer to [21] for more details. Additionally, we downloaded Google mobility data, publicly available at https://www.google.com/COVID19/mobility/ (accessed on September 19, 2020). These data include indicators of mobility pattern changes in different areas such as parks, residential locations, retail locations, among others; and information such as the percent change in revenue as well as total number of small businesses for high/low income consumers and across different economic activities; see Table 1 for a complete listing.

Table 1

List of all variables considered (see Data Description Section for further details).

Item	Feature description
1	cases per million
2	ratio of effective positive results to effective total tests
3	ratio of effective positive results to state population in 2019
4	individuals who are currently hospitalized with COVID-19
5	spending in all categories
6	spending in arts, entertainment, and recreation categories
7	spending in accommodation and food service categories
8	spending in general merchandise stores and apparel and accessories categories
9	spending in grocery and food store categories
10	spending in health care and social assistance categories
11	spending in transportation and warehousing categories
12	spending among high (top quartile) income ZIP codes in all categories
13	spending among low (bottom quartiles) income ZIP codes in all categories
14	spending among middle (middle two quartiles) income ZIP codes in all categories
15	time spent at retail and recreation locations
16	time spent at grocery and pharmacy locations
17	time spent at parks
18	time spent at inside transit stations
19	time spent at work places
20	time spent at residential locations
21	time spent outside of residential locations
22	number of small businesses open
23	number of small businesses open among high(top quartile) income ZIP codes
24	number of small businesses open among low(bottom quartile) income ZIP codes
25	number of small businesses open among middle(middle two quartiles) income ZIP codes
26	number of small businesses open in transportation
27	number of small businesses open in education and health services
28	number of small businesses open in leisure and hospitality
29	net revenue for small businesses
30	net revenue for small businesses among high(top quartile) income ZIP codes
31	net revenue for small businesses among low(bottom quartile) income ZIP codes
32	net revenue for small businesses among middle(middle two quartiles) income ZIP codes
33	net revenue for small businesses in transportation
34	net revenue for small businesses in education and health services
35	net revenue for small businesses in leisure and hospitality

We briefly describe a few features below and refer to the Economic Tracker source [21] for further details:

‘spend all’: Seasonally adjusted credit and debit card spending relative to the period Jan 4 to Jan 31, 2020, in all merchant categories, reported as a seven day moving average formed by using the values for the current and the previous six days.

‘spend tws’: Seasonally adjusted credit and debit card spending relative to the period Jan 4 to Jan 31, 2020, in the category transportation and warehousing, and also reported as a seven day moving average.

‘spend all inchigh’: Seasonally adjusted credit and debit card spending by consumers living in high median income (above 75% quantile) zip codes, and relative to the period Jan 4 to Jan 31, 2020, in all merchant categories, again reported as a seven day moving average.

‘revenue all’: Seasonally adjusted percent change in net revenue for small businesses and indexed to the period Jan 4 to Jan 31, 2020, reported as a seven day moving average.

‘revenue ss70’: Seasonally adjusted percent change in net revenue for small businesses in leisure and hospitality and indexed to the period Jan 4 to Jan 31, 2020, reported as a seven day moving average.

‘merchant inchigh’: Seasonally adjusted percent change in the number of small businesses open in high median income (above 75% quantile) zip codes and indexed to the period Jan 4 to Jan 31, 2020, reported as a seven day moving average.

‘gps parks’: Corresponds to the time spent at parks relative to a baseline which is defined as the median value, for the corresponding day of the week, during the period Jan 3–Feb 6, 2020.

We employ local linear regression to slightly smooth the feature trajectories; this also allows to impute missing observations for some states and days. For this step, a Gaussian kernel with bandwidth 1.5 days was employed and we utilized the R package fdapace version 0.5.5 for the computational implementation [19]. The effective test positivity rate was smoothed analogously. We use local quadratic regression [28] to obtain the derivatives , implemented by using the R package locfit version 1.5-9.4 [39]. For data analysis we consider states with population at least 1 million as the states with smaller populations had less COVID-19 cases, and their inclusion would therefore increase the overall noise level in the data and negatively impact the analysis. This is a pragmatic choice and is not a limitation of the proposed method. In the end, we included the data from 44 states (excluding Alaska, Delaware, North and South Dakota, Vermont and Wyoming).

As a result of the LASSO variable selection described in Section 4.3, which was implemented with the R package ncvreg version 3.12.0 [9] with threshold in (16), the features utilized in the final model are as follows, where the selected lags are listed in Table 2 .

Table 2

Predictors and their corresponding lags (see Data Description Section for further details).

Predictor	variable	lag
Park mobility activity	U1(t)	0
Number of people currently hospitalized	U2(t)	14
Effective test positivity rate	U3(t)	7
Credit/debit card spending on arts, entertainment, and recreation	U4(t)	3

, the total cases per million at time t for state i, .

, the relative time spent at park areas at time t in state i.

, the number of hospitalized people at time t in state i.

, the “positive testing rate” defined as the ratio of effective positive results (positive test or probable case) to the effective total tests (positive or probable case plus negative test results) at time t for state i.

, the credit/debit card relative spending in the arts, entertainment, and recreation category (see ‘spend aer’ predictor in Table 1) at time t in state i.

Model fitting and results

In order to study the delay time dynamics of COVID-19 in the United States, we focus on the state-specific time-varying cumulative confirmed cases per million people, , for , which is referred to as the case process in the following. Here, the time domain that we consider starts on April 5 and ends on August 12 in 2020, where data are recorded once per day, on an equidistant grid for 130 days. We estimate the history coefficient function γ in the RDED model (5) as described in Section 4.2 with , and adopt the functional varying-coefficient regression in model (5) with other variables as listed in Table 1. The estimated coefficient surface and some of its cross-sections are displayed in Fig. 1 .

Fig. 1

(a) Estimated coefficient surface in the RDED model (5), where t corresponds to number of days since April 5, 2020 and s takes values from 0 to τ0. (b) t-sliced cross sections for t = 1,26,51,76,101,126 of as functions of s.

For each fixed t, , the estimated history index in model (5) is found to change from positive near current time to negative with increasing lags. This means that the recent past of the process has a positive association with the current derivative while days further into the past are negatively associated. The right panel shows that the shape of changes from more or less quadratic to linear from the early days of COVID-19 to recent times, indicating that in April 2020 the contrast between the effects of near-current and past caseload on the current derivative was more pronounced than it is more recently.

Predictor process effects on growth rate

The results are visualized in Fig. 2 . States which are suffering from faster spread of the virus are characterized by higher growth rates. The historically weighted integrated case count per million has a positive effect on the growth rate. This reflects dynamic explosive behavior [45], so that states will tend to have caseloads that move even further away from the mean caseload per million across states as time moves forward, which could be either moving further above the average across states or further below it. Thus the positive effect of the history index reflects increasing variability in the caseloads across states. Potential reasons for this are the distinct policies that states implement in terms of mobility, business and social restrictions, as well as cultural acceptance of such restrictions. Spatial infectivity dynamics with waves of infections moving spatially across the U.S. also might play a major role.

Fig. 2

Concurrent effect functions for predictor processes U in the RDED (5).

In terms of the effects of covariate processes, our results indicate that higher park mobility during late spring and summer (May onwards) is associated with a decrease in growth rates. For each percent increase in park mobility, the growth rate decreased by as much as 10 cases per million per day, with the effect reaching its maximum in mid-summer. This indicates negligible risk and potential benefits from outdoors recreation activities.

The other covariate processes that were selected by LASSO as predictors in the RDED model (5) are patient counts at hospitals, the testing positivity rate, and credit card spending on recreational activities, including dining at restaurants. All of these were found to have positive effects on the growth rate throughout nearly the entire time domain. These effects were included in model (5) with covariate process-specific lags, as per Table 2. The number currently hospitalized was found to be most predictive for the case growth rate when including a two-week lag, whereas for test positivity rate a one-week lag was found to be most predictive. This may suggest that individuals who are infected by individuals with less severe symptoms (i.e. a non-hospitalized COVID-positive person) tend to learn of their disease one week after their infector did.

The two-week difference for hospitalized individuals is likely a reflection of the time needed to develop serious symptoms which warrant admission for medical care. The effect of recreational spending on viral spread was much more immediate; the optimal lag in this case is only 3 days. This suggests that states with higher recreational economies during the pandemic see an almost immediate uptick in growth rate, with cases increasing by as much as 25 infections per million per day for each percent of additional recreational spending.

Predictive performance of the random differential equation with delay model

Figs. 3 and 4 display the fitted growth rates provided by the RDED model (5), with covariate specific lags, and also the results of fitting a differential equation model where the covariates do not have individual lags, i.e. a concurrent RDE model given by (5) in which . Fitted curves are the result of leave-one-out predictions for each state, where the data from the left-out state were not used for fitting the model and then the state's data are plugged into the fitted model to obtain the predicted case growth rate. These predictions are compared against the actual observed growth rates as obtained by local polynomial regression (see Section 4.1 and Appendix C). Visual inspection (as well as measures of goodness-of-fit, such as integrated mean squared prediction error) indicates that the incorporation of delays with predictor specific lags in the model substantially improves the fits. This provides compelling evidence for modeling COVID-19 growth rates with RDEDs that include covariate-specific lags. A vanilla RDE without predictor delays does not incorporate historical information beyond the previous day and struggles to explain deviations in phase (as in, e.g. Illinois). The problem of modeling phase variation is instead foisted upon the coefficient functions, which are ill-equipped to reflect it, and thus at times the straight RDE model without predictor delays falters in terms of amplitude modeling as well (e.g. for Colorado). In contrast, the proposed RDED incorporates past information through the included delays and thus reflects phase variability in a natural way, especially through the inclusion of distributed delays via the history index. Including a two-week history of case counts through the history index function γ in model (5) substantially enhances the flexibility of the model. The estimates of coefficient functions are no longer uncontrollably altered by the effects of unaccounted lags, and as a result the fitted model performs much better.

Model-benchmarked evaluation of performance

A comparison of observed vs. RDED model-predicted growth rates allows one to assess the performance of states in terms of controlling spread. A state whose observed growth rate falls below the predicted curve corresponds to a negative residual and suggests better-than-expected viral control, and vice versa for observed rates which fall above the predicted curve. Fig. 5, Fig. 6 display the residual curves for the fitted RDED model (5), where positive residuals correspond to the number of excess new cases per day relative to the predictions obtained from fitting model (5), as described above, including the two-week history index and other covariate processes with their respective lags. These residuals reflect the unexplained residual stochastic process Z that is an integral and principally unpredictable part of the RDED model (5) on which our predictions are based.

Fig. 5

Residuals for the predictions from the RDED (5) for US States, Alabama–Montana.

Fig. 6

Residuals for US States, Nebraska–Wisconsin.

Major outbreaks, which can be interpreted as an unforeseen amount of excess cases, are revealed through positive spikes in the residual plots. For example, New York's residuals leap up to as high as 170 excess cases per million per day during a spike in Spring. Arizona's summer outbreak is captured in the residuals as well, subsequently course-correcting with a streak of fewer-than-expected new cases around day 100. These spikes also suggest that in the two weeks prior to the uptick, there were substantially more cases in the state population than the numbers reported in the JHU dataset.

The consistency of a state's handling of the virus can be seen through the volatility of the residual curve. States like Colorado and Oregon exhibit only minor deviations from zero, which suggests their new cases are quite predictable given a two-week history of the state case counts and covariate processes. Some states, particularly in the Northeast, including New York, New Jersey, and Rhode Island, only exhibit high volatility during the early days of the pandemic before achieving a period of stability. On the other hand, states such as Alabama and Louisiana exhibit less-predictable case counts throughout.

Computation times and code availability

The codes for reproducing the above analysis are available at https://osf.io/w48zd/?view_only=7cd3e2c686a74fd086e4eae8534214c0. The run times of the different steps in the analysis using R version 3.6.3 (2020-02-29) running under CentOS Linux 7 (Core) on x86_64-pc-linux-gnu (64-bit) platform are summarized in Table 3 .

Table 3

Run time for the different steps in the data analysis.

Step	Description	Number of Cores	Time in Seconds
1	preprocessing	single core	9.475
2	initial lag selection	8 cores in parallel	422.505
3	history index function estimation	single core	3.414
4	lasso	8 cores in parallel	31.316
5	backfitting	single core	2565.715
6	smoothing	single core	0.012
	Total		3032.437

Discussion and concluding remarks

In this paper we propose dynamics learning exemplified by learning a RDED from a sample of observed trajectories, which is motivated by modeling COVID-19 daily new cases from the previous behavior of the case trajectory itself through a distributed delay or history index component and additional covariate process components. RDEDs form the driving mechanism of real life processes that include a feedback loop into the future. Depending on the application and the exact nature of the feedback, one might incorporate discrete delays or distributed delays. While discrete delays borrow information from some isolated times in the past, a distributed delay takes into account the entire continuum in the recent past and turns out to be closely related to a functional history index model that has been considered in the statistics literature.

In the current framework, we have not included inference for the estimated model coefficient functions and the lags. It is left for future research to investigate the convergence of these estimated parameters to their population targets and to obtain guarantees on the convergence rates. Some words of caution are in order. Any method that includes a time delay component suffers from the limitation that a larger set of possible lags means one has to sacrifice a part of the data that are available for the analysis, due to the initialization. We have assumed that the trajectories corresponding to the different states are i.i.d. realizations, which may not hold if there is significant spatial correlation based on geographic and demographic similarities across the states. It would be interesting to study adaptations of the proposed methodology when the samples are correlated instead of being independent. From an application point of view, a different and perhaps more refined perspective could potentially be gained by applying the proposed method at a finer level to county level data.

Establishing the existence and uniqueness of the solutions of RDEDs has found increasing interest in recent years [14], [24]. We contribute here to the forward problem of obtaining solutions of an RDED by targeting the inverse problem analytically and combine this with functional data analysis approaches to estimate the parameter functions from an observed sample of multivariate stochastic processes for both discrete and distributed delay setups. Our approach builds on and contributes to empirical dynamics and the literature on inferring dynamics and differential equations from functional data, and extends it to models that include time delays.

A key feature of the proposed approach is the statistical estimation of the model parameters in the population RDED, which has inbuilt model selection steps that aid in estimating the optimal lags and also to select the important predictors through the lasso pruning steps, with majority voting across all time points. This reduces the number of predictors in the model and thus its complexity prior to the optimization step of the lags that are associated with the predictors. To handle this complex lag optimization step, instead of a brute force search method, we adopted a backfitting algorithm, which significantly simplifies the optimization step and also leads to faster convergence. We provide a complete toolkit that starts with a sample of multivariate trajectories and extracts from the observed data an estimated RDED that best explains the dynamics inherent in the data, along with automatic variable selection and delay estimation.

The application of our method to the modeling of COVID 19 growth rates in the United States demonstrates excellent performance in terms of predictions of the trajectories. The predictions from the random differential equation model without predictor lags grossly overestimates the COVID-19 growth rates for significant periods of time in New York, New Jersey, Pennsylvania, Massachusetts, Connecticut, Virginia, West Virginia, New Hampshire, California, Texas, Colorado, Oregon, Washington, Maine and Hawaii but underestimates the growth rates for Alabama, Arizona, Arkansas, Florida, Louisiana, Nebraska, Nevada, South Carolina, Tennessee and Wisconsin; see Fig. 3, Fig. 4. This highlights the importance of incorporating time lags in differential equation modeling of functional data, when one has a sample of realized stochastic processes, specially when the underlying processes tend to have an aftereffect.