Time dynamics of COVID-19

We apply tools from functional data analysis to model cumulative trajectories of COVID-19 cases across countries, establishing a framework for quantifying and comparing cases and deaths across countries longitudinally. It emerges that a country’s trajectory during an initial first month “priming period” largely determines how the situation unfolds subsequently. We also propose a method for forecasting case counts, which takes advantage of the common, latent information in the entire sample of curves, instead of just the history of a single country. Our framework facilitates to quantify the effects of demographic covariates and social mobility on doubling rates and case fatality rates through a time-varying regression model. Decreased workplace mobility is associated with lower doubling rates with a roughly 2 week delay, and case fatality rates exhibit a positive feedback pattern.

S.1 Details on data pre-processing Table 1 presents information such as the first up-crossing time of 20 confirmed cases for each country in the study. The cross-sectional ranks and integrated ranks are estimated as per Section S.6.
The Google community mobility trends data measure how the frequency and length of stay at different locations have changed relative to a baseline level prior to the COVID-19 pandemic. The baseline activity level is defined as the median activity value for the corresponding day of the week in the period of Jan. 3 to Feb. 6, 2020 during which time most countries (excluding China) did not implement any distancing efforts. Available categories include retail, grocery, park, transit, workplace, and residential categories.
Since these categories are highly collinear we consider just workplace mobility as an index of social activity (see Section S.7).  Table 1: First up-crossing times during 2020 of 20+ confirmed cases and dates after the 67 days time window for each of the 64 countries in the study along with their cross-sectional ranks at t = 5 and t = 65 days, as well as the integrated rank. Figure 1 shows the imputed curves for the confirmed cases through the FPCA method with K = 2 eigenfunctions for the countries Brazil, Chile, India, Italy, Switzerland, United Kingdom, and the United States. The first 2 eigenfunctions explain 97% of the total variation and the quality of the fits suggests that including just two components in the Karhunen-Loève expansion works well.

S.3 Functional Concurrent Regression
For an observed response curve Y (t), functional predictor X(t) and baseline covariate U , the functional concurrent regression model is given by where β 0 (t), β 1 (t) and β 2 (t) are smooth coefficient functions and (t) is a zero mean Gaussian process. Various estimation techniques are available for the intercept function β 0 (t) and the slope functions β 1 (t) and β 2 (t) for both densely and sparsely observed functional data [1,2,3].
For assessing the goodness of fit, one can use the dynamic coefficient of determination .
Larger values of R 2 (t) indicate that a larger fraction of the variability in the response Y (t) is explained by the linear model in X(t) and U . A positive-valued slope function β 1 (t) indicates a positive association between Y (t) and X(t) at time t, while a negative-valued slope implies a negative association, and this also applies to the scalar covariates' dynamic associations with the response Y at time t.
In some cases there may be a time lag in the association between the functional response Y (t) and the functional predictor X(t), for t > ∆ where ∆ ≥ 0 denotes the lag, which is usually unknown.
When one does not have prior knowledge about ∆, its value can be selected by optimizing a data-adaptive criterion, as follows: Let Y i (t j ) denote the response process of the i th subject observed at the j th time point t j and X i (t j ) be the corresponding predictor process observed at the same time point. Let U i be the baseline covariate associated with the i th subject. LetŶ ∆ i (t j ) be the fitted value for Y i (t j ) obtained after fitting the functional concurrent regression model based on ∆ in (1) using all but the i th subject's observations.
Then the mean-normalized leave-one-out prediction error P error (∆) is defined as for situations where the response curves Y i are on different scales, which is the case in our application. The optimal ∆ is chosen aŝ where I ∆ is the set of potential candidates for ∆.

S.4 Historical Functional Linear Model
Prediction of doubling rates is based on the historical functional linear regression model with scalar response Y and functional predictors C(s), W (s), with t − 13 ≤ s ≤ t − 1, where is a mean-zero error term. The functions β 1 (s), β 2 (s) are then estimated by representing them in the eigenfunctions of C(s), W (s) respectively, which are assumed to form a basis of the function space L 2 . Writing ξ C,k , φ C,k (s), ξ W,k , φ W,k (s) to denote functional principal component scores and eigenfunctions of C(s) and W (s), and letting we truncate the sum at K included terms so that the model becomes The β 1,k can then be estimated by simply fitting a linear regression model with predictors ξ C,k , ξ W,k , k = 1, 2, . . . , K. The evaluation of the model is the same as in Section S.3, where we used R 2 (t) = 1 − var( (t)) var(Y (t)) .

S.5 Empirical Dynamics
The foundation of FDA relies on the assumption that the observed sample of trajectories is generated from an underlying smooth and square integrable process. We also assume that derivatives exist that can then be utilized to study the underlying dynamics of the process [4,5,6]. For both functional and longitudinal data, empirical dynamics [7] provides a principled approach to learn and quantify the dynamics. The motivation for empirical dynamics is the fact that for a differentiable Gaussian process Y (t) with mean µ(t) one where β(t) = d dt log[var{Y (t)}] is a smooth dynamic coefficient function and Z(t) is a random drift process independent of Y (t). This leads to the linear model formulation that can be analyzed under the functional concurrent regression framework described in Section S.3. For non-Gaussian processes, model (4) yields effective approximations in the least squares sense, where one can use the coefficient of determination to gauge the fraction of variance in d dt Y (t) explained linearly by Y (t). On domains where R 2 (t) is large, the linear term β(t)(Y (t) − µ(t)) in model (4) provides a systematic approach to study simultaneously the empirical dynamics of Y (t) and the dependence on other baseline or functional covariates. The lag ∆ can either be set to zero for a truly concurrent approach, or may be selected using a data adaptive criterion as described in Section S.3 to model a relationship that involves a time delay.
For analyzing empirical dynamics, one needs to estimate derivatives given the observations Y i (t j ), with notation borrowed from Section S.3. For our applications, we choose to estimate derivatives by local quadratic smoothing [8,9] using the Epanechnikov kernel and a bandwidth of 2 days.

S.6 Rank dynamics
For functional data, cross-sectional ranks and their temporal dynamics may be investigated through the rank processes and summary statistics for rank dynamics as per [10].
For a generic stochastic process Y : T → R, our starting point is the cross-sectional distribution P (Y (t) ≤ z) =: F t (z), for each t ∈ T . Without loss of generality, we consider T = [0, 1]. The rank processes S i associated with trajectories Y i are then given by To summarize the overall rank of each trajectory and quantify how each rank trajectory varies with time, we consider two summary measures for rank processes, individual-specific integrated rank ρ i and rank volatility ν i , defined as respectively. For the COVID-19 data, as described before, the caseload trajectories C(t) are observed on the same time grid T = [0, 66] for each country. A straightforward estimate for the rank processes in (6) is then given by replacing the cross-sectional distribution function F t by its empirical counterpart, i.e., Hence, individual-specific integrated rank ρ i and rank volatility ν i in (7) can be estimated by plugging in the estimated rank processes in (8) and taking numerical integration. The results are shown in Figure 2.  (7) for empirical ranks during the first 67 days since exposure.