Penalized Estimation and Forecasting of Multiple Subject Intensive Longitudinal Data

Abstract

Intensive longitudinal data (ILD) is an increasingly common data type in the social and behavioral sciences. Despite the many benefits these data provide, little work has been dedicated to realize the potential such data hold for forecasting dynamic processes at the individual level. To address this gap in the literature, we present the multi-VAR framework, a novel methodological approach allowing for penalized estimation of ILD collected from multiple individuals. Importantly, our approach estimates models for all individuals simultaneously and is capable of adaptively adjusting to the amount of heterogeneity present across individual dynamic processes. To accomplish this, we propose a novel proximal gradient descent algorithm for solving the multi-VAR problem and prove the consistency of the recovered transition matrices. We evaluate the forecasting performance of our method in comparison with a number of benchmark methods and provide an illustrative example involving the day-to-day emotional experiences of 16 individuals over an 11-week period.

Appendix

In this technical appendix, we discuss some theoretical aspect of LASSO estimation in the multi-VAR setting, namely concerning its consistency and sparsistency.

1.1 Consistency

Consistency of LASSO estimation for single (stable) VAR models was established in the seminal paper by Basu and Michailidis (2015a, b), building upon such results in the regression setting by Loh and Wainwright (2012a, b). In the multi-VAR setting, the model is inherently unidentifiable. It could be that the LASSO solution is consistent for some particular $\varvec{\mu }^{*}$, $\varvec{\Delta }^{*}_{k}$ in the model (7), or over a subset of such identifications, but this problem still appears largely unresolved. Some related result though can be found in the discussion on sparsistency below following Ollier and Viallon (2017). Here, we shall discuss a weaker form of consistency of $\hat{\mathbf {B}}_{k}=\hat{\varvec{\mu }}+\hat{\varvec{\Delta }}_{k}$ for $\mathbf {B}_k^{*}$. The arguments are quite straightforward and shed some light on the problem, and also seemingly were not made in the related literature yet.

We first describe the basic result for a single VAR model expressed in the regression form (5), and then turn to a multi-VAR model. We index the model quantities with subscript k or superscript (k), $k=1,\ldots ,K$, representing the individual models in the multi-VAR setting. After expanding the quadratic term of the objective function (6), the estimation equation can be rewritten as in Basu and Michailidis (2015a) in terms of the quantities

$$\begin{aligned} \widehat{\Gamma }_{k}=\frac{1}{N}\mathbf {Z}^{(k)'}\mathbf {Z}^{(k)}=\frac{1}{N}(\mathrm {I}_{d}\otimes \mathcal {X}^{(k)}\mathcal {X}^{(k)'}),\quad \hat{\gamma }_{k}=\frac{1}{N}\mathbf {Z}^{(k)'}\mathbf {Y}^{(k)}. \end{aligned}$$

(31)

Estimation consistency is proved under the following two conditions on these quantities:

• Restricted eigenvalue condition: The matrix $\widehat{\Gamma }_{k}$ is said to satisfy this condition with parameters $\alpha _{k},\tau _{k}>0$, if

  $$\begin{aligned} \beta _{k}'\widehat{\Gamma }_{k}\beta _{k}\ge \alpha _{k}\Vert \beta _{k}\Vert _{2}^{2}-\tau _{k}\Vert \beta _{k}\Vert _{1}^{2},\quad \beta _{k}\in \mathbb {R}^{q}, \end{aligned}$$

  (32)

  with $q=pd^{2}$.

• Deviation condition: This condition is satisfied if

  $$\begin{aligned} \Vert \hat{\gamma }_{k}-\widehat{\Gamma }_{k}\mathbf {B}_{k}^{*}\Vert _{\infty } \le Q_{k}(\mathbf {B}_{k}^{*},\Sigma _{k,\varepsilon }) \sqrt{\frac{\log q}{N}}, \end{aligned}$$

  (33)

  for a deterministic function $Q_{k}$.

Let $s_{k}=\Vert \mathbf {B}_{k}^{*}\Vert _{0}$ denote the sparsity of the model. Under the conditions above and assuming $s_{k}\tau _{k}\le \alpha _{k}/32$, Proposition 4.1 of Basu and Michailidis (2015a) states that any solution $\hat{\mathbf {B}}$ of (6) satisfies: for any $\lambda \ge 4Q_{k}(\mathbf {B}_{k}^{*},\Sigma _{k,\varepsilon })\sqrt{\frac{\log q}{N}}$,

$$\begin{aligned} \Vert \hat{\mathbf {B}}_{k}-\mathbf {B}_{k}^{*}\Vert _{1}\le \frac{64s_{k}\lambda }{\alpha _{k}},\quad \Vert \hat{\mathbf {B}}_{k}-\mathbf {B}_{k}^{*}\Vert _{2}\le \frac{16\sqrt{s_{k}}\lambda }{\alpha _{k}}. \end{aligned}$$

(34)

Additionally, a result on the support of thresholded estimators of $\hat{\mathbf {B}}_{k}$ is also available. The consistency results in (34) apply to generic LASSO estimators as long as the quantities $\widehat{\Gamma }_{k},\hat{\gamma }_{k}$ satisfy the restricted eigenvalue and deviation conditions.

Among the key contributions of Basu and Michailidis (2015a) are their results (Propositions 4.2 and 4.3) proving that $\widehat{\Gamma }_{k}$ and $\hat{\gamma }_{k}$ satisfy the restricted eigenvalue and deviation conditions with high enough probabilities, and expressing the various parameters involved in the conditions $(\alpha _{k},\tau _{k},Q_{k}(\mathbf {B}_{k}^{*},\Sigma _{k,\varepsilon }))$ in terms of the VAR model parameters. Furthermore, in the restricted eigenvalue condition, $\tau _{k}$ can be chosen so that $s_{k}\tau _{k}\le \alpha _{k}/32$. We also note that the right-hand side of the inequalities (34) are expected to be negligible for small $\lambda$ and hence small $\log q/N$. The case when the logarithm of the dimension compares to the sample size through this way is the typical LASSO scenario.

In the multi-VAR setting, the optimization problem (9) can be expressed through the objective function

$$\begin{aligned} -\sum _{k=1}^{K}2\mathbf {B}_k'\hat{\gamma }_k + \sum _{k=1}^{K}\mathbf {B}_k'\widehat{\Gamma }_{k}\mathbf {B}_k + \lambda _{1}\Vert \varvec{\mu }\Vert _1 + \sum _{k=1}^{K}\lambda _{2,k}\Vert \mathbf {B}_{k}-\varvec{\mu }\Vert _1. \end{aligned}$$

(35)

A consistency bound for the minimizer $\hat{\mathbf {B}}_{k}$ of (35) can still be obtained similarly as for single VAR models if one is willing to make the assumption

$$\begin{aligned} \Vert \hat{\varvec{\mu }}\Vert _0 \le s_0. \end{aligned}$$

(36)

The constraint (36) could be imposed while optimizing (35) or choosing $\lambda _1$ appropriately large, or inferred to hold (with high enough probability) from sparsistency result, if available. Indeed, under (36), a consistency bound can be derived easily as in the proof of Proposition 3.3 in Basu and Michailidis (2015a, b). That is, observe first that

$$\begin{aligned}&-\sum _{k=1}^{K}2\hat{\mathbf {B}}_k'\hat{\gamma }_k + \sum _{k=1}^{K}\hat{\mathbf {B}}_k'\widehat{\Gamma }_{k}\hat{\mathbf {B}}_k + \lambda _{1}\Vert \hat{\varvec{\mu }}\Vert _1 + \sum _{k=1}^{K}\lambda _{2,k}\Vert \hat{\mathbf {B}}_{k}-\hat{\varvec{\mu }}\Vert _1\\&\quad \le -\sum _{k=1}^{K}2\mathbf {B}_k^{*'}\hat{\gamma }_k + \sum _{k=1}^{K}\mathbf {B}_k^{*'}\widehat{\Gamma }_{k}\mathbf {B}_k^{*} + \lambda _{1}\Vert \hat{\varvec{\mu }}\Vert _1 + \sum _{k=1}^{K}\lambda _{2,k}\Vert \mathbf {B}_{k}^{*}-\hat{\varvec{\mu }}\Vert _1 \end{aligned}$$

and rearranging the terms and setting $\mathbf {v}_k=\hat{\mathbf {B}}_k-\mathbf {B}_k^*$, we deduce

$$\begin{aligned} \sum _{k=1}^{K}\mathbf {v}_k \widehat{\Gamma }_k \mathbf {v}_k \le \sum _{k=1}^{K}2\mathbf {v}_k' (\hat{\gamma }_k-\widehat{\Gamma }_k \mathbf {B}_k^{*})+\sum _{k=1}^{K}\lambda _{2,k}\left( \Vert \mathbf {B}_k^{*}-\hat{\varvec{\mu }}\Vert _1 - \Vert \mathbf {B}_{k}^{*} - \hat{\varvec{\mu }} + \mathbf {v}_k\Vert _1 \right) . \end{aligned}$$

With $\hat{J}_{k}=\text {supp}\{ \mathbf {B}_k^{*} - \hat{\varvec{\mu }}$ being the index support of $\mathbf {B}_k^{*} - \hat{\varvec{\mu }}\}$, repeating the argument in Basu and Michailidis (2015a, b), we get

$$\begin{aligned} 0 \le \sum _{k=1}^{K} \mathbf {v}_k' \widehat{\Gamma }_k \mathbf {v}_k \le \sum _{k=1}^{K}\left( \frac{3\lambda _{2,k}}{2}\Vert (\mathbf {v}_{k})_{\hat{J}_k}\Vert _1 -\frac{\lambda _{2,k}}{2}\Vert (\mathbf {v}_{k})_{\hat{J}_k^{c}}\Vert _1 \right) \end{aligned}$$

(37)

as long as $\lambda _{2,k}\ge 4Q_k(\mathbf {B}_k^{*},\Sigma _{k,\varepsilon })\sqrt{\frac{\log q}{N}}$ (with the function $Q_k$ from the deviation condition), where $(\cdot )_{\hat{J}}$ and $(\cdot )_{\hat{J}^c}$ denote restrictions to the index sets $\hat{J}$ and $\hat{J}^c$, respectively. Then,

$$\begin{aligned} \sum _{k=1}^{K}\lambda _{2,k}\Vert (\mathbf {v}_k)_{\hat{J}_k^{c}}\Vert _1 \le 3\sum _{k=1}^{K}\lambda _{2,k}\Vert (\mathbf {v}_k)_{\hat{J}_k}\Vert _1 \end{aligned}$$

and one also has

$$\begin{aligned}&\sum _{k=1}^{K}\lambda _{2,k}\Vert \mathbf {v}_k\Vert _1 \le 4\sum _{k=1}^{K}\lambda _{2,k}\Vert (\mathbf {v}_k)_{\hat{J}_k}\Vert _1 \nonumber \\&\quad \le 4\sum _{k=1}^{K}\lambda _{2,k}(s_0+s_k)^{1/2}\Vert \mathbf {v}_k\Vert _2 \le 4 \sqrt{\sum _{k=1}^{K}\lambda _{2,k}^{2}(s_0+s_k)}\Vert \mathbf {v}\Vert _2, \end{aligned}$$

(38)

by Cauchy–Schwarz inequality (twice) and the fact that $|\text {supp}\{\hat{J}_{K}\}|\le s_0+s_k$, where $\Vert v\Vert _2^2=\sum _{k=1}^{K}\Vert v_k\Vert _2^2$. Similarly, by the restricted eigenvalue condition (32) for each model and assuming $s_k\tau _k \le \alpha _k/32$, we have

$$\begin{aligned} \sum _{k=1}^{K} \mathbf {v}_k'\widehat{\Gamma }_k \mathbf {v}_k \ge \sum _{k=1}^{K}\frac{\alpha _k}{2}\Vert \mathbf {v}_k\Vert _2^2 \ge \frac{\min \{\alpha _k\}}{2}\Vert \mathbf {v}\Vert _2^2. \end{aligned}$$

(39)

A combination of (37)–(39) yields, e.g.,

$$\begin{aligned} \frac{\min \{\alpha _k\}}{2}\Vert \mathbf {v}\Vert _2^2 \le 6\sqrt{\sum _{k=1}^{K}\lambda _{2,k}^2(s_0+s_k)}\Vert \mathbf {v}\Vert _2 \end{aligned}$$

(40)

or

$$\begin{aligned} \Vert \mathbf {v}\Vert _2 \le \frac{12\sqrt{\sum _{k=1}^{K}\lambda _{2,k}^2(s_0+s_k)}}{\min \{\alpha _k\}}. \end{aligned}$$

(41)

This is the multi-VAR analogue of the second consistency bound in (34). One can similarly obtain a bound on $\Vert \mathbf {v}\Vert _1$ analogous to the first one in (34).

1.2 Sparsistency

We comment here briefly on the possibility of recovering the supports of $\varvec{\mu }^{*}$ and $\varvec{\Delta }_{k}^{*}$. The same issue of (non)identifiability is fundamental here as well. Some result nevertheless are available in the literature for special cases. Assuming effectively that $s\lambda _1/\lambda _{2,k}=cK^{1/2}$, Ollier and Viallon (2017) gave conditions for identifiability and sparsistency with the limiting common parameter of interest $\varvec{\mu }^{*}$ defined as the entrywise median of $\varvec{B}_k^{*}$. Their approach goes through verifying a particular well-known irrepresentability condition on a design matrix. It could in principle be adapted to the multi-VAR context but the value of this effort might be questionable. First, irrepresentability conditions are quite restrictive and difficult to verify, and as a result, adaptive LASSO versions are advocated for. The setting where the limiting parameter of interest is necessarily related to the median could also be viewed restrictive.