Sequential Learning of Regression Models by Penalized Estimation

Abstract

When data arrive in a sequence of two or more datasets, modeling on the most recent dataset should take previous datasets into account. We specifically investigate a strategy for regression modeling when parameter estimates from previous data can be used as anchoring points, yet may not be available for all parameters, thus, covariance information cannot be reused. A procedure that updates through targeted penalized estimation, which shrinks the estimator toward a nonzero value, is presented. The parameter estimate from the previous data serves as this nonzero value when an update is sought from novel data. This naturally extends to a sequence of datasets with the same response, but potentially only partial overlap in covariates. The iteratively updated regression parameter estimator is shown to be asymptotically unbiased and consistent. The penalty parameter is chosen through constrained cross-validated log-likelihood optimization. The constraint bounds the amount of shrinkage of the updated estimator toward the current one from below. The bound aims to preserve the (updated) estimator’s goodness of fit on all-but-the-novel data. The proposed approach is compared to other regression modeling procedures. Finally, it is illustrated on an epidemiological study where the data arrive in batches with different covariate-availability and the model is refitted with the availability of a novel batch. Supplementary materials for this article are available online.

Keywords:

• Aymptotic unbiasedness
• Consistency
• Constrained cross-validation
• Generalized linear model
• Targeted ridge (${\ell }_2$) penalty
• Updating

1 Introduction

Settings where data continuously arrive and need to be incorporated into modeling have been addressed by incremental or online learning techniques (Losing, Hammer, and Wersing 2018). Transfer learning approaches have been developed for situations when there are just two datasets, and the first might only be loosely similar in structure to the second, but still informative (Pan and Yang 2010; Shilo, Rossman and Segal 2020). The idea of transfer learning has also been shown to be more generally useful, for example, for regression modeling, where it may also be more feasible to get analytical insight into algorithmic approaches (Minami et al. 2020). Some authors have also considered the combination of both challenges, that is, online transfer learning (Zhao et al. 2014). We are specifically interested in such online transfer learning for regression modeling, where covariance structure and the set of available covariates may change between datasets. In particular we will investigate a transfer learning strategy based on ridge penalized estimation.

While transfer learning for regression models allows for differences in the response as well as in the parameters when moving from dataset to dataset, the framework of transfer learning typically assumes that there is an initial model and that the corresponding structure stays the same when new data become available and need to be incorporated. Transfer learning then refers to the updating of knowledge of the model parameters when new information becomes available (Shalev-Shwartz 2012). As such it produces a sequence of parameter estimates. The latest estimate takes into account the most recent information but also exploits previously acquired data, while minimizing a certain loss function. Examples of transfer regression techniques abound. Statistical readers will be most familiar with recursive least squares (Plackett 1950) and mixed models (Verbeke and Molenbergs 2009). In this context and related, we recently showed how to perform transfer learning of the normal distribution’s precision matrix through targeted ridge penalized estimation (van Wieringen et al. 2020).

In the following, we consider settings where the response stays the same, when moving along a sequence of datasets, but the set of covariates may change as may the values they assume. For example, data of large epidemiological studies often arrive in batches. Each batch sheds light on the same phenomenon, for example, environment-disease relation, observed with possibly different covariate information. As such with the arrival of a new batch, the currently learned model of the phenomenon under study needs updating. In this article we show how updating can be achieved by means of penalized estimation.

Here we propose targeted ridge penalized estimation as an enrichment of the incremental/online transfer learning regression toolbox. The presented procedure updates the regression parameter estimate by minimization of the loss function using the novel data and with the target formed by the current parameter estimate derived from the hitherto available data. The approach is designed to not depend on covariance estimates of parameters from the previous data, as a change of available covariates will typically invalidate the usefulness of such covariance estimates. We discuss incorporation of multiple competing current estimates, if available. We study properties, for example, moments and consistency, of this estimator for transfer regression learning. We then translate all of this to the learning of the logistic regression model. The penalty parameter is chosen via constrained cross-validation to warrant learning and avoid one-off estimation. We compare transfer regression learning via targeted ridge penalized estimation to standard alternatives. Finally, we illustrate the potential of the method on an epidemiological study.

2 Methods

2.1 Model and Penalty Structure

Throughout we entertain either a linear or logistic regression model, being the most widely used variants of generalized linear models. Both are learned by means of targeted ridge penalized estimation, an extension of ridge regression (inspired by Singh, Chaubey, and Dwivedi 1986). It considers estimating the p-dimensional regression coefficient vector $\mathit{\beta }$ of the linear regression model, $\mathbf{Y}=\mathbf{X}\mathit{\beta }+\mathit{\epsilon }$ with $\mathit{\epsilon }\sim \mathcal{N}({\mathbf{0}}_n,{\sigma }_{\epsilon }^2{\mathbf{I}}_{nn})$, by minimizing the nonzero centered ridge regression loss function:(1) $$\left|\right|\mathbf{Y}-\mathbf{X}\mathit{\beta }|{|}_2^2+\lambda ||\mathit{\beta }-{\mathit{\beta }}_0|{|}_2^2,$$(1) with penalty parameter $\lambda >0$ and nonrandom shrinkage target ${\mathit{\beta }}_0$. An analytic expression for the minimizer of the above loss function exists and gives rise to the nonzero centered ridge regression estimator: $\widehat{\mathit{\beta }}(\lambda ,{\mathit{\beta }}_0)={({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}+\lambda {\mathbf{I}}_{pp})}^{-1}({\mathbf{X}}^{\mathrm{\top }}\mathbf{Y}+\lambda {\mathit{\beta }}_0)$. High-dimensionally, this expression can be evaluated computational efficiently when rewritten in terms of the singular value decomposition of X (Hastie and Tibshirani 2004). If ${\mathit{\beta }}_0={\mathbf{0}}_p$ the traditional ridge regression loss function is recovered. The motivation for this estimator becomes evident after reformulation of the penalized estimation problem above, that is, the minimization of loss (1), as a constrained estimation problem. When ${\mathit{\beta }}_0={\mathbf{0}}_p$, the regression parameter constraint is centered at the origin. But when $\lambda \cancel{=}0$ and ${\mathit{\beta }}_0\cancel{=}{\mathbf{0}}_p$ the constraint on $\mathit{\beta }$ is nonzero centered. Of particular interest is the case when an informative choice for ${\mathit{\beta }}_0$ is available. The targeted ridge estimator then becomes less biased (Singh, Chaubey, and Dwivedi 1986).

Two alternative interpretations of this estimator are relevant for the presented work. First, the nonzero centered ridge regression can be interpreted as an attempt to correct for the nonzero target. This is easily seen when using the change-of-variable $\mathit{\gamma }=\mathit{\beta }-{\mathit{\beta }}_0$. Substitution of $\mathit{\gamma }$ in the loss function (1) gives: $\left|\right|(\mathbf{Y}-\mathbf{X}{\mathit{\beta }}_0)-\mathbf{X}\mathit{\gamma }|{|}_2^2+\lambda |\left|\mathit{\gamma }\right|{|}_2^2$. In this $\mathbf{Y}-\mathbf{X}{\mathit{\beta }}_0$ is the residual vector after invoking the prior knowledge. This is regressed on X to correct for the variance in Y wrongly attributed to X by the prior knowledge as condensed in the nonzero target ${\mathit{\beta }}_0$. Secondly, from a Bayesian perspective the nonzero centered ridge penalty still corresponds (in the sense that the posterior mean coincides with the optimum of the loss function above) to a normal prior on the elements of the regression parameter $\mathit{\beta }$. But this prior now has a mean equal to the corresponding element of ${\mathit{\beta }}_0$.

2.2 Updating the Linear Regression Model

The targeted ridge regression estimator introduced in the previous section is put to use in the case where multiple datasets become available sequentially at time points $t=1,2,3,\dots$. At each time point we fit the linear regression model (now equipped with the index t to separate the datasets): ${\mathbf{Y}}_t={\mathbf{X}}_t\mathit{\beta }+{\mathit{\epsilon }}_t$ with ${\mathit{\epsilon }}_t\sim \mathcal{N}({\mathbf{0}}_{n_t},{\sigma }_{\epsilon }^2{\mathbf{I}}_{n_t,n_t})$ using knowledge from the previous case. To this end we use the nonzero centered ridge regression framework. In the above, the design matrices ${\left\{{\mathbf{X}}_t\right\}}_{t=1}^{\infty }$ may differ among datasets. For the theoretical exposé, however, the number of covariates, that is, the column dimension p of the ${\mathbf{X}}_t$’s, is fixed, while the sample sizes ${\left\{n_t\right\}}_{t=1}^{\infty }$ may differ. Moreover, the settings/values of the covariates, for example, the elements of ${({\mathbf{X}}_t)}_{i,j}$ with $i=1,\dots ,n_t$ and $j=1,\dots ,p$, are also free to vary over the t. The latter implies that the degree of collinearity among the covariates may change from dataset to dataset. Or, put mathematically, ${\mathbf{X}}_t^{\mathrm{\top }}{\mathbf{X}}_t$ need not equal ${\mathbf{X}}_{t\prime }^{\mathrm{\top }}{\mathbf{X}}_{t\prime }$ if $t\cancel{=}t\prime$. While the design matrices may thus change from dataset to dataset, the ${\mathbf{X}}_t$ are considered to be nonrandom in the remainder. Finally, the effect of these covariates on the response, however, is assumed to be fixed through time. The error variance is also assumed to be fixed through time. The latter assumption is not necessary as the theoretical results can be generalized to accommodate an error variance that differs among datasets.

At the arrival of a new tth dataset, the estimates of parameters $\mathit{\beta }$ and ${\sigma }_{\epsilon }^2$ obtained from the previous $(t-1)$th dataset(s) are available. These may be used as prior knowledge in the updating of the parameter estimates from the tth dataset through:$$\begin{array}{ccc}{\widehat{\mathit{\beta }}}_t({\lambda }_t)\hfill & =\hfill & {({\mathbf{X}}_t^{\mathrm{\top }}{\mathbf{X}}_t+{\lambda }_t{\mathbf{I}}_{pp})}^{-1}[{\mathbf{X}}_t^{\mathrm{\top }}{\mathbf{Y}}_t+{\lambda }_t{\widehat{\mathit{\beta }}}_{t-1}({\lambda }_{t-1},{\widehat{\mathit{\beta }}}_{t-2})],\hfill \end{array}$$where the previous estimate of $\mathit{\beta }$ thus, plays the role of ${\mathit{\beta }}_0$. In addition, the ${\widehat{\mathit{\beta }}}_{t-1}$ and ${\widehat{\sigma }}_{\epsilon ,t-1}^2$ may be used to choose the penalty parameters by minimization of the MSE plugging in these previous estimates for $\mathit{\beta }$ and ${\sigma }_{\epsilon }^2$. Later we will choose the penalty parameter by means of cross-validation (see Section 2.4).

The updated ridge linear regression estimator is a linear combination of the observed responses from the subsequent studies. Consequently, it is normally distributed with expectation and variance (see supplementary material I for details):$$\begin{array}{c}\mathbb{E}[{\widehat{\mathit{\beta }}}_t({\lambda }_t)]=\sum _{t_h=1}^t\left\{\prod _{\tau =t_h+1}^t{\left[{\lambda }_{\tau }{({\mathbf{X}}_{\tau }^{\mathrm{\top }}{\mathbf{X}}_{\tau }+{\lambda }_{\tau }{\mathbf{I}}_{pp})}^{-1}\right]}^{I_{\{t\ge \tau \}}}\right\}\mathit{\beta }\hfill \\ -\sum _{t_h=1}^t\left\{\prod _{\tau =t_h}^t\right[{\lambda }_{\tau }{({\mathbf{X}}_{\tau }^{\mathrm{\top }}{\mathbf{X}}_{\tau }+{\lambda }_{\tau }{\mathbf{I}}_{pp})}^{-1}\left]\right\}\mathit{\beta }\hfill \\ +\left\{\prod _{\tau =1}^t\right[{\lambda }_{\tau }{({\mathbf{X}}_{\tau }^{\mathrm{\top }}{\mathbf{X}}_{\tau }+{\lambda }_{\tau }{\mathbf{I}}_{pp})}^{-1}\left]\right\}{\mathit{\beta }}_0,\hfill \\ \text{var}[{\widehat{\mathit{\beta }}}_t({\lambda }_t)]={\sigma }_{\epsilon }^2\sum _{t_h=1}^t\left[\prod _{\tau =t_h+1}^t{({\lambda }_{\tau }^2)}^{I_{\{t\ge \tau \}}}\right]\hfill \\ \left[\prod _{\tau =t_h}^t{({\mathbf{X}}_{\tau }^{\mathrm{\top }}{\mathbf{X}}_{\tau }+{\lambda }_{\tau }{\mathbf{I}}_{pp})}^{-1}\right]\hfill \\ {\mathbf{X}}_{t_h}^{\mathrm{\top }}{\mathbf{X}}_{t_h}\left[\prod _{\tau =t_h}^t{({\mathbf{X}}_{\tau }^{\mathrm{\top }}{\mathbf{X}}_{\tau }+{\lambda }_{\tau }{\mathbf{I}}_{pp})}^{-1}\right].\hfill \end{array}$$

In the above the λ_t’s have been assumed to be nonrandom (as is done by Hoerl and Kennard 1970 in the derivation of the moments). To provide some intuition, assume (i) every dataset is equipped with the same orthonormal design matrix (i.e., ${\mathbf{X}}_t^{\mathrm{\top }}{\mathbf{X}}_t={\mathbf{I}}_{pp}$) and (ii) ${\lambda }_t=\lambda$ for all t. The expectation and variance then reduce to $\mathbb{E}[{\widehat{\mathit{\beta }}}_t(\lambda )]=\mathit{\beta }+{\lambda }^t{(1+\lambda )}^{-t}({\mathit{\beta }}_0-\mathit{\beta })$ and $\text{var}[{\widehat{\mathit{\beta }}}_t(\lambda )]={\sigma }_{\epsilon }^2{(1+2\lambda )}^{-1}[1-{\lambda }^{2t}{(1+\lambda )}^{-2t}]{\mathbf{I}}_{pp}$. In particular, $\underset{t\to \infty }{\lim }\mathbb{E}[{\widehat{\mathit{\beta }}}_t(\lambda )]=\mathit{\beta }$ and $\underset{t\to \infty }{\lim }\text{var}[{\widehat{\mathit{\beta }}}_t(\lambda )]={\sigma }_{\epsilon }^2{(1+2\lambda )}^{-1}{\mathbf{I}}_{pp}$. In words, when the number of novel datasets grows large enough, the updated ridge estimator becomes unbiased. Moreover, its variance is smaller than that of the ML (maximum likelihood) regression estimator for any $\lambda >0$.

Next, we obtain more general results on the convergence of the sequence of ridge updated regression estimator and the associated linear predictor. To this end, notice that the above generated sequence of regression estimators ${\{{\widehat{\mathit{\beta }}}_t({\lambda }_t,{\widehat{\mathit{\beta }}}_{t-1})\}}_{t=1}^{\infty }$ (in which the arguments of ${\widehat{\mathit{\beta }}}_{t-1}$ have been dropped to reduce notational clutter) forms a first-order Markov chain with continuous state space ${ℝ}^p$ as:$$\begin{array}{c}{\widehat{\mathit{\beta }}}_{t+1}({\lambda }_{t+1},{\widehat{\mathit{\beta }}}_t)|{\{{\widehat{\mathit{\beta }}}_{t\prime }({\lambda }_t^{\prime },{\widehat{\mathit{\beta }}}_{t-1})\}}_{t\prime =0}^t\\ \hspace{1em}\sim {\widehat{\mathit{\beta }}}_{t+1}({\lambda }_{t+1},{\widehat{\mathit{\beta }}}_t)|{\widehat{\mathit{\beta }}}_t({\lambda }_t,{\widehat{\mathit{\beta }}}_{t-1}).\end{array}$$

Moreover, it is time-homogeneous:$$\begin{array}{ccc}\hfill & \hfill & {\widehat{\mathit{\beta }}}_{t+\tau +1}({\lambda }_{t+\tau +1},{\widehat{\mathit{\beta }}}_{t+\tau })\left|\right\{{\widehat{\mathit{\beta }}}_{t+\tau }({\lambda }_{t+\tau },{\widehat{\mathit{\beta }}}_{t+\tau -1})\hfill \\ \hfill & \hfill & \hspace{1em}=\mathit{\beta },{\lambda }_{t+\tau }=\lambda ,{\mathbf{X}}_{t+\tau +1}=\mathbf{X}\}\hfill \\ \hfill & \hfill & \hspace{1em}\sim {\widehat{\mathit{\beta }}}_{t+1}({\lambda }_{t+1},{\widehat{\mathit{\beta }}}_t)\left|\right\{{\widehat{\mathit{\beta }}}_t({\lambda }_t,{\widehat{\mathit{\beta }}}_{t-1})=\mathit{\beta },{\lambda }_t=\lambda ,{\mathbf{X}}_{t+1}=\mathbf{X}\},\hfill \end{array}$$for any $\tau \in \mathbb{N},\mathbf{X}\in {\mathcal{M}}^{n,p},\mathit{\beta }\in {ℝ}^p$, and $\lambda >0$.

The well-known theory of Markov process can now be exploited to show (under some conditions) the asymptotic unbiasedness of the ridge updated linear regression estimator and predictor (Theorem 1). The theorem’s proof, deferred – as all proofs – to the supplementary material II, requires to show that the updating process has a stationary distribution with expectation equal to $\mathit{\beta }$.

Theorem 1.

Assume the existence of an infinite sequence of studies into the linear relationship between a continuous response and a set of covariates. The data from these studies, ${\{{\mathbf{X}}_t,{\mathbf{Y}}_t\}}_{t=1}^{\infty }$, are used to fit the linear regression model by means of the updated ridge linear regression estimator, which yields the sequence of estimators ${\{{\widehat{\mathit{\beta }}}_t({\lambda }_t,{\widehat{\mathit{\beta }}}_{t-1})\}}_{t=1}^{\infty }$ which is initiated by an arbitrary, nonrandom ${\mathit{\beta }}_0$. Furthermore, let $T\in \mathbb{N}$ be sufficiently large and ${\mathbf{X}}_{\text{new}}$ be the design matrix with covariate information on novel samples for which a prediction is needed. Then,

• $\underset{t\to \infty }{\lim }{\mathbb{E}}_{{\left\{{\mathbf{Y}}_{\tau }\right\}}_{\tau =1}^{t+1}}[{\widehat{\mathit{\beta }}}_{t+1}({\lambda }_{t+1},{\widehat{\mathit{\beta }}}_t)|{\left\{{\lambda }_{\tau }\right\}}_{\tau =1}^{t+1}]=\mathit{\beta }+\mathbf{u}$ for some $\mathbf{u}\in {\cap }_{t=T}^{\infty }\text{null}({\mathbf{X}}_t)$, where $\text{null}({\mathbf{X}}_t)$ denotes the null space of the linear map induced by ${\mathbf{X}}_t$. If ${\cap }_{t=T}^{\infty }\text{null}({\mathbf{X}}_t)={\mathbf{0}}_p$, then $\mathbf{u}={\mathbf{0}}_p$.

• $\underset{t\to \infty }{\lim }{\mathbb{E}}_{{\left\{{\mathbf{Y}}_{\tau }\right\}}_{\tau =1}^{t+1}}[{\mathbf{X}}_{\text{new}}{\widehat{\mathit{\beta }}}_{t+1}({\lambda }_{t+1},{\widehat{\mathit{\beta }}}_t)|{\left\{{\lambda }_{\tau }\right\}}_{\tau =1}^{t+1}]={\mathbf{X}}_{\text{new}}\mathit{\beta }$.

A few notes on this theorem are to be pointed out. First, the theorem describes asymptotic behavior. But it does so not in the traditional sense, in which the sample size goes by one at a time. Here that may be if we include a new study with a degenerated sample size. But as most studies are expected to have a sample size exceeding one, the limit is reached in larger step sizes. Moreover, the data of the new study is not added to the data of the previous study, rather the former is weighed against a summary obtained from the latter: the current updated ridge linear regression estimator.

Another, more important note concerns the role of the penalty parameters ${\left\{{\lambda }_t\right\}}_{t=1}^{\infty }$ in Theorem 1. While in the evaluation of the asymptotic expectation they are conditioned on, their specific values are irrelevant, as the proof only requires the fact that ${\lambda }_t>0$ for all t. In practice, we choose the penalty parameters by means of cross-validation or some other means involving data. The penalty parameters are thus, random variables themselves. This may cause different limiting behavior, although Theorem 1 indicates that their actual realizations are irrelevant. Indeed, simulations show that, even with penalty parameters chosen by means of cross-validation, the estimates converge to the true parameters.

In addition to asymptotic unbiasedness the updated ridge linear regression estimator, we show it to be consistent (Theorem 2), although under assumptions on the penalty parameter sequence.

Theorem 2.

Assume the existence of an infinite sequence of studies on the linear relationship between a continuous response and a set of covariates. The data from these studies, ${\{{\mathbf{X}}_t,{\mathbf{Y}}_t\}}_{t=1}^{\infty }$, are used to fit the linear regression model by means of the updated ridge linear regression estimator, which yields the sequence of estimators ${\{{\widehat{\mathit{\beta }}}_t({\lambda }_t,{\widehat{\mathit{\beta }}}_{t-1})\}}_{t=1}^{\infty }$ which is initiated by an arbitrary, nonrandom ${\mathit{\beta }}_0$. Furthermore, let $T\in \mathbb{N}$ be sufficiently large and ${\cap }_{t=T}^{\infty }\text{null}({\mathbf{X}}_t)={\mathbf{0}}_p$. Assume that the penalty parameter sequence ${\left\{{\lambda }_t\right\}}_{t=1}^{\infty }$ is chosen such that $\underset{t\to \infty }{\lim }{\sigma }_{\epsilon }^{2}p{d}_1^2({\mathbf{X}}_t){\lambda }_t^{-2}=0$ with $d_1({\mathbf{X}}_t)$ the largest singular value of ${\mathbf{X}}_t$. Then, for every c > 0:$$\begin{array}{c}\underset{t\to \infty }{\lim }P\left[\right||{\widehat{\mathit{\beta }}}_t({\lambda }_t,{\widehat{\mathit{\beta }}}_{t-1})-\mathit{\beta }||\ge c|{\left\{{\lambda }_{\tau }\right\}}_{\tau =1}^t]\to 0,\\ \underset{t\to \infty }{\lim }P\left[\right||{\mathbf{X}}_{\text{new}}{\widehat{\mathit{\beta }}}_t({\lambda }_t,{\widehat{\mathit{\beta }}}_{t-1})-{\mathbf{X}}_{\text{new}}\mathit{\beta }||\ge c|{\left\{{\lambda }_{\tau }\right\}}_{\tau =1}^t]\to 0.\end{array}$$

The proposed updating scheme can be seen as a frequentist analogue of Bayesian updating (Berger 2013). In the Bayesian narrative, the ridge penalty corresponds to a normal prior. The resulting posterior of $\mathit{\beta }$ is also a normal distribution. When updating, this normal posterior then serves as (normal) prior for the next update of $\mathit{\beta }$. The prior for the t + 1th update then becomes ${\mathit{\beta }}_{t+1}|{\sigma }^2\sim \mathcal{N}[{\mathit{\beta }}_t({\lambda }_t,{\mathit{\beta }}_{t-1}),{\sigma }^2{\lambda }_{t+1}^{-1}{\mathbf{I}}_{pp}]$, which yields a posterior mean $\mathbb{E}[{\mathit{\beta }}_{t+1}|{\mathbf{Y}}_{t+1},{\mathbf{X}}_{t+1},{\sigma }^2,{\mathit{\beta }}_t({\lambda }_t,{\beta }_{t-1})]$ coinciding with the frequentist estimator ${\widehat{\mathit{\beta }}}_{t+1}({\lambda }_{t+1},{\mathit{\beta }}_t)$.

Intuition (behind the proofs) of Theorems 1 and 2 is best formed in the low-dimensional case with all design matrices ${\mathbf{X}}_t$ of full rank. Furthermore, denote the ML regression estimator of tth dataset by ${\widehat{\mathit{\beta }}}_t={({\mathbf{X}}_t^{\mathrm{\top }}{\mathbf{X}}_t)}^{-1}{\mathbf{X}}_t^{\mathrm{\top }}{\mathbf{Y}}_t$. Then, note that ${({\mathbf{X}}_t^{\mathrm{\top }}{\mathbf{X}}_t+{\lambda }_t{\mathbf{I}}_{pp})}^{-1}{\mathbf{X}}_t^{\mathrm{\top }}{\mathbf{Y}}_t={\widehat{\mathit{\beta }}}_t-{\lambda }_t{({\mathbf{X}}_t^{\mathrm{\top }}{\mathbf{X}}_t+{\lambda }_t{\mathbf{I}}_{pp})}^{-1}{\widehat{\mathit{\beta }}}_t$. The tth update ridge regression estimator can then be expressed by either of the two expression below:$$\begin{array}{c}{\widehat{\mathit{\beta }}}_t({\lambda }_t,{\widehat{\mathit{\beta }}}_{t-1})={\widehat{\mathit{\beta }}}_t+\sum _{\tau =0}^{t-1}{(-1)}^{\tau +1}\\ \hspace{1em}\left\{\prod _{{\tau }^{\prime }=0}^{\tau }\right[{\lambda }_{t-{\tau }^{\prime }}{({\mathbf{X}}_{t-{\tau }^{\prime }}^{\mathrm{\top }}{\mathbf{X}}_{t-{\tau }^{\prime }}+{\lambda }_{t-{\tau }^{\prime }}{\mathbf{I}}_{pp})}^{-1}\left]\right\}\\ \hspace{1em}({\widehat{\mathit{\beta }}}_{t-\tau }-{\widehat{\mathit{\beta }}}_{t-\tau -1}),\\ {\widehat{\mathit{\beta }}}_t({\lambda }_t,{\widehat{\mathit{\beta }}}_{t-1})=\sum _{\tau =0}^t\left\{\prod _{{\tau }^{\prime }=0}^{\tau }\right[{({\mathbf{X}}_{t-{\tau }^{\prime }}^{\mathrm{\top }}{\mathbf{X}}_{t-{\tau }^{\prime }}+{\lambda }_{t-{\tau }^{\prime }}{\mathbf{I}}_{pp})}^{-1}\\ \hspace{1em}{\mathbf{X}}_{t-\tau \prime }^{\mathrm{\top }}{\mathbf{X}}_{t-\tau \prime }\left]\right\}{\widehat{\mathit{\beta }}}_{t-\tau }\\ =\sum _{\tau =0}^t{\mathbf{W}}_{t-\tau \prime }{\widehat{\mathit{\beta }}}_{t-\tau },\end{array}$$with ${\widehat{\mathit{\beta }}}_0={\mathit{\beta }}_0$ and ${\mathbf{W}}_t\le {\mathbf{W}}_{t-1}$ if $t\prime <t$. The result of Theorem 1 becomes insightful after taking the expectation of the former expression of the preceeding display, in which the weighted sum of the ML regression estimator differences vanishes for large t. The latter expression of the same display can loosely be interpreted as a weighted average, which—under an appropriate regularization scheme (see Theorem 2)—has a vanishing variance as $t\to \infty$.

Finally, the target need not be formed by the most recently updated ridge regression estimate. It may be replaced by alternative estimates obtained from the preceding dataset(s). For instance, the ridge regression estimate obtained from all preceding datasets may play the role of the target. These targets may even be simultaneously available prior to the estimation of the regression parameter from the current data. One would then preferably have the current data choose among these targets. This can be done straightforwardly within the current framework. Hereto replace the targeted ridge penalty by what could be considered a mixture of such penalties. For $G\in \mathbb{N}$ such targets ${\mathit{\beta }}_{0,t,1},\dots ,{\mathit{\beta }}_{0,t,G}$, this yields the loss function: $\left|\right|{\mathbf{Y}}_t-{\mathbf{X}}_t\mathit{\beta }|{|}_2^2+{\lambda }_t{\sum }_{g=1}^G{\alpha }_{t,g}||\mathit{\beta }-{\mathit{\beta }}_{0,t,g}|{|}_2^2$, where all ${\alpha }_{t,g}\ge 0$ and ${\sum }_{g=1}^G{\alpha }_{t,g}=1$. This loss is minimized with respect to $\mathit{\beta }$ by$$\begin{array}{ccc}\widehat{\mathit{\beta }}({\lambda }_t,{\alpha }_{t,1},\dots ,{\alpha }_{t,G})\hfill & =\hfill & {({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}+{\lambda }_t{\mathbf{I}}_{pp})}^{-1}\hfill \\ \hfill & \hfill & ({\mathbf{X}}^{\mathrm{\top }}\mathbf{Y}+{\lambda }_t{\sum }_{g=1}^G{\alpha }_{t,g}{\mathit{\beta }}_{0,t,g}).\hfill \end{array}$$

The targets are effectively averaged weightedly to form a novel target. The weights of this average are unknown. Like the penalty parameter λ_t, the weights are tuning parameters and are to be chosen prior to the estimation, for example, through a procedure like cross-validation (cf. Section 2.4). The multi-targeted ridge regression estimator is not accompanied by theoretical results, although similar asymptotic results may be obtained. Instead, we show in simulation that the estimator shrinks toward the preferred target (see supplementary material III).

2.3 Updating the Logistic Regression Model

The results of the previous section carry over to generalized linear models. This is exemplified here for the logistic regression model. Hereto the same experimental set-up as for the linear regression model is assumed to apply, now with a binary response: $Y_i\in \{0,1\}$ for $i=1,\dots ,n$. The probability of a ‘1’ – and thereby that of a ‘0’ – is modeled by the logistic link function: $P(Y_i=1)=\hspace{0.17em}\text{exp}\hspace{0.17em}({\mathbf{X}}_{i,\ast }\mathit{\beta })/[1+\hspace{0.17em}\text{exp}\hspace{0.17em}({\mathbf{X}}_{i,\ast }\mathit{\beta })]$ for every i. Again the regression parameter $\mathit{\beta }$ is estimated by maximization of the log-likelihood augmented with a nonzero centered ridge penalty:$$\begin{array}{ccc}\widehat{\mathit{\beta }}(\lambda )\hfill & =\hfill & \arg \underset{\mathit{\beta }\in {ℝ}^p}{\max }{\sum }_{i=1}^n\{Y_i{\mathbf{X}}_{i,\ast }\mathit{\beta }-\hspace{0.17em}\text{log}\hspace{0.17em}[1+\hspace{0.17em}\text{exp}\hspace{0.17em}({\mathbf{X}}_{i,\ast }\mathit{\beta })\left]\right\}\hfill \\ \hfill & \hfill & -\lambda \left|\right|\mathit{\beta }-{\mathit{\beta }}_0|{|}_2^2.\hfill \end{array}$$

The corresponding estimating equation is(2) $${\mathbf{X}}^{\mathrm{\top }}[\mathbf{Y}-{\overrightarrow{\mathbf{g}}}^{-1}(\mathbf{X};\mathit{\beta })]-\lambda (\mathit{\beta }-{\mathit{\beta }}_0)={\mathbf{0}}_p,$$(2) where ${\overrightarrow{\mathbf{g}}}^{-1}(\mathbf{X};\mathit{\beta })={[g^{-1}({\mathbf{X}}_{1,\ast };\mathit{\beta }),\dots ,g^{-1}({\mathbf{X}}_{n,\ast };\mathit{\beta })]}^{\mathrm{\top }}$ with link function $g(·;·)$ defined by $g^{-1}(·;·)=\hspace{0.17em}\text{exp}\hspace{0.17em}(·;·)/[1+\hspace{0.17em}\text{exp}\hspace{0.17em}(·;·)]$ and ${\mathbf{X}}_{i,\ast }$ denoting the ith row of the design matrix. The root of Equation (2)(2) $${\mathbf{X}}^{\mathrm{\top }}[\mathbf{Y}-{\overrightarrow{\mathbf{g}}}^{-1}(\mathbf{X};\mathit{\beta })]-\lambda (\mathit{\beta }-{\mathit{\beta }}_0)={\mathbf{0}}_p,$$(2) is found by means of a Newton–Raphson algorithm, reformulated as an iteratively reweighted least squares algorithm. This requires a minor modification from that presented in Schaefer, Roi, and Wolfe (1984) and Le Cessie and van Houwelingen (1992) for the zero-centered ridge logistic regression estimator. Starting from an initial guess ${\widehat{\mathit{\beta }}}^{(0)}$, the estimate is updated by$$\begin{array}{ccc}{\widehat{\mathit{\beta }}}^{(k+1)}(\lambda )\hfill & =\hfill & {({\mathbf{X}}^{\mathrm{\top }}\mathbf{W}\mathbf{X}+\lambda {\mathbf{I}}_{pp})}^{-1}({\mathbf{X}}^{\mathrm{\top }}\mathbf{W}\mathbf{Z}+\lambda {\mathit{\beta }}_0)\hfill \\ \hfill & =\hfill & [{\lambda }^{-1}{\mathbf{I}}_{pp}-{\lambda }^{-2}{\mathbf{X}}^{\mathrm{\top }}{({\mathbf{W}}^{-1}+{\lambda }^{-1}\mathbf{X}{\mathbf{X}}^{\mathrm{\top }})}^{-1}\mathbf{X}]\hfill \\ \hfill & \hfill & ({\mathbf{X}}^{\mathrm{\top }}\mathbf{W}\mathbf{Z}+\lambda {\mathit{\beta }}_0),\hfill \end{array}$$ where W diagonal with ${(\mathbf{W})}_{ii}=\hspace{0.17em}\text{exp}\hspace{0.17em}({\mathbf{X}}_i{\widehat{\mathit{\beta }}}^{(k)}){[1+\hspace{0.17em}\text{exp}\hspace{0.17em}({\mathbf{X}}_i{\widehat{\mathit{\beta }}}^{(k)})]}^{-2}$ and the adjusted response variable $\mathbf{Z}=\{\mathbf{X}{\widehat{\mathit{\beta }}}^{(k)}+{\mathbf{W}}^{-1}[\mathbf{Y}-{\overrightarrow{\mathbf{g}}}^{-1}(\mathbf{X};{\mathit{\beta }}^{(k)})]\}$. The last expression in the preceding display is obtained by means of the Woodbury matrix identity and to be used in the high-dimensional case for computational efficiency. In particular, the computationally most demanding task, the evaluation of $\mathbf{X}{\mathbf{X}}^{\mathrm{\top }}$, can be performed prior to the iterations. This efficiency gain could also have been obtained by means of singular value decomposition (Eilers et al. 2001). The recursive formula above is applied until convergence, which yields the desired nonzero centered ridge logistic regression estimate. The inclusion of a nonzero target ${\mathit{\beta }}_0$ generally increases the computing time of the iteratively reweighted least squares algorithm (see supplementary material IV). For informative targets, this increase is relatively small. Off targets, however, may yield larger computing times.

We now assume data ${\{{\mathbf{X}}_t,{\mathbf{Y}}_t\}}_{t=1}^{\infty }$ from a sequence of studies on the same logistic regression model, that is, with a common regression parameter $\mathit{\beta }$ but possibly different design matrices, are available. From each study, we estimate the parameter by means of the nonzero centered ridge logistic regression estimator with the estimate of the previous study as shrinkage target. The resulting estimate ${\widehat{\mathit{\beta }}}_{t+1}[{\lambda }_{t+1},{\widehat{\mathit{\beta }}}_t({\lambda }_t,{\widehat{\mathit{\beta }}}_{t-1})]$ is that $\mathit{\beta }$ which solves$$\begin{array}{ccc}{\mathbf{X}}_{t+1}^{\mathrm{\top }}[{\mathbf{Y}}_{t+1}-{\overrightarrow{\mathbf{g}}}^{-1}({\mathbf{X}}_{t+1};\mathit{\beta })]-{\lambda }_{t+1}[\mathit{\beta }-{\widehat{\mathit{\beta }}}_t({\lambda }_t,{\widehat{\mathit{\beta }}}_{t-1})]\hfill & =\hfill & {\mathbf{0}}_p.\hfill \end{array}$$

This yields a sequence of estimators ${\{{\widehat{\mathit{\beta }}}_t({\lambda }_t,{\widehat{\mathit{\beta }}}_{t-1})\}}_{t=2}^{\infty }$, in which (again) the arguments of the ${\widehat{\mathit{\beta }}}_{t-1}$ have been dropped to reduce notational clutter. The t + 1th estimator ${\widehat{\mathit{\beta }}}_{t+1}({\lambda }_{t+1},{\widehat{\mathit{\beta }}}_t)$ of this sequence relates, after convergence of the iteratively reweighted least squares algorithm, to its predecessor by:$$\begin{array}{ccc}{\widehat{\mathit{\beta }}}_{t+1}({\lambda }_{t+1},{\widehat{\mathit{\beta }}}_t)\hfill & =\hfill & {({\mathbf{X}}_{t+1}^{\mathrm{\top }}{\mathbf{W}}_{t+1}{\mathbf{X}}_{t+1}+{\lambda }_{t+1}{\mathbf{I}}_{pp})}^{-1}\hfill \\ \hfill & \hfill & [{\mathbf{X}}_{t+1}^{\mathrm{\top }}{\mathbf{W}}_{t+1}{\mathbf{Z}}_{t+1}+{\lambda }_{t+1}{\widehat{\mathit{\beta }}}_t({\lambda }_t,{\widehat{\mathit{\beta }}}_{t-1})],\hfill \end{array}$$with ${\mathbf{W}}_{t+1}$ and ${\mathbf{Z}}_{t+1}$ now involving ${\mathbf{X}}_{t+1}$ and ${\mathbf{Y}}_{t+1}$.

The thus defined sequence of updated logistic regression estimators ${\{{\widehat{\mathit{\beta }}}_t({\lambda }_t,{\widehat{\mathit{\beta }}}_{t-1})\}}_{t=2}^{\infty }$ can again be seen as being generated by a first-order Markov chain with ${ℝ}^p$ as its state space. As before, we exploit this fact to show asymptotic properties of the updated ridge logistic regression estimator. Theorem 3 states the asymptotic unbiasedness and consistency of the updated estimator and predictor.

Theorem 3.

Assume an infinite sequence of studies into the generalized linear relationship between a binary response and a set of covariates. The data from these studies, ${\{{\mathbf{X}}_t,{\mathbf{Y}}_t\}}_{t=1}^{\infty }$, are used to fit the logistic regression model by means of the updated ridge logistic regression estimator, which yields the sequence of estimators ${\{{\widehat{\mathit{\beta }}}_t({\lambda }_t,{\widehat{\mathit{\beta }}}_{t-1})\}}_{t=1}^{\infty }$ which is initiated by an arbitrary, nonrandom ${\mathit{\beta }}_0$. Let $T\in \mathbb{N}$ be sufficiently large. Then:

• $\underset{t\to \infty }{\lim }\mathbb{E}[{\widehat{\mathit{\beta }}}_{t+1}({\lambda }_{t+1},{\widehat{\mathit{\beta }}}_t)]=\mathit{\beta }+\mathbf{u}$ for some $\mathbf{u}\in {\cap }_{t=T}^{\infty }\text{null}({\mathbf{X}}_t)$. If ${\cap }_{t=T}^{\infty }\text{null}({\mathbf{X}}_t)={\mathbf{0}}_p$, then $\mathbf{u}={\mathbf{0}}_p$.

• if ${\cap }_{t=T}^{\infty }\text{null}({\mathbf{X}}_t)={\mathbf{0}}_p$ and ${\left\{{\lambda }_t\right\}}_{t=1}^{\infty }$ such that $\underset{t\to \infty }{\lim }2p^{1/2}|d_1({\mathbf{X}}_t)|{\lambda }_t=0$ with $d_1({\mathbf{X}}_t)$ the largest singular value of ${\mathbf{X}}_t$, for every c > 0:$$\begin{array}{c}\underset{t\to \infty }{\lim }P\left[\right||{\widehat{\mathit{\beta }}}_t({\lambda }_t,{\widehat{\mathit{\beta }}}_{t-1})-\mathit{\beta }||\ge c|{\left\{{\lambda }_{\tau }\right\}}_{\tau =1}^t]\to 0,\\ \underset{t\to \infty }{\lim }P\left[\right||{\mathbf{X}}_{\text{new}}{\widehat{\mathit{\beta }}}_t({\lambda }_t,{\widehat{\mathit{\beta }}}_{t-1})-{\mathbf{X}}_{\text{new}}\mathit{\beta }||\ge c|{\left\{{\lambda }_{\tau }\right\}}_{\tau =1}^t]\to 0.\end{array}$$

The updated ridge logistic regression estimator too is a frequentist analogue of a Bayesian updating scheme. Now the analogy is not exact but approximate. This is due to the fact that the normal prior is not conjugate for the logistic likelihood and, consequently, the posterior is not a well-known and characterized distribution. Nonetheless, the latter may be approximated in a Laplacian manner by a normal distribution (cf. Bishop 2006), for which the Bernstein–von Mises theorem (van der Vaart 1998) provides conditions that ensure its quality. The mean of this approximating normal corresponds to the mode of posterior (i.e., the updated ridge logistic regression estimator), while its variance relates to the curvature of the posterior at its mode. This normal-like posterior then serves as the prior when updating the knowledge on the logistic regression parameter at the arrival of a novel dataset.

The targeted ridge logistic regression estimator can also be modified to accommodate multiple targets acquired from different estimators. Again, this amounts to averaging the targets weightedly.

2.4 Choice of the Tuning Parameter

We present informative choices for the tuning parameters, ${\left\{{\lambda }_t\right\}}_{t=1}^{\infty }$ and ${\mathit{\beta }}_0$, of updated ridge estimators. We choose the ridge penalty parameter of a novel dataset through a form of cross-validation (see Arlot and Celisse 2010 for a survey). Cross-validation chooses the penalty parameter that yields the best performance of the model on novel data. In absence of novel data part of data, referred to as the test or left-out data, is set aside to serve as such. The remaining data, called the training or left-in data, are used to learn the model that is evaluated on the test data. The thus acquired performance depends on the selection of the test data, and may accidentally yield an overly optimistic performance. To remove the effect of the particular the data are split into training and test data repetitively, say, K times, thus, giving rise to K-fold cross-validation. Splitting may be done randomly but is done here in a stratified manner. The latter procedure ensures approximate equally sized test datasets and the single occurrence of each sample in a test dataset. The performance averaged over the test datasets is considered to be representative for novel data. The chosen penalty parameter optimizes this performance. The model’s performance is usually measured by its fit on the test data. For the linear regression model, this amounts to the cross-validated sum-of-squared $K^{-1}{\sum }_{k=1}^K\left|\right|{\mathbf{Y}}_t^{(k)}-{\mathbf{X}}_t^{(k)}{\widehat{\mathit{\beta }}}_t^{(-k)}({\lambda }_t)|{|}_2^2$. In this loss, ${\mathbf{Y}}_t^{(k)}$ and ${\mathbf{X}}_t^{(k)}$ are formed by subsetting ${\mathbf{Y}}_t$ and ${\mathbf{X}}_t$ to the kth test data, while ${\widehat{\mathit{\beta }}}_t^{(-k)}$ is the ridge regression estimator obtained from all but the kth test data. In case of logistic regression, the cross-validated likelihood $K^{-1}{\sum }_{k=1}^K\mathcal{L}[{\mathbf{Y}}_t^{(k)},{\mathbf{X}}_t^{(k)},{\widehat{\mathit{\beta }}}_t^{(-k)}({\lambda }_t)]$ is used as a performance measure. The optimal λ_t is found by minimization/maximization of the appropriate criterion using standard machinery. In the remainder, we use K = 10-fold cross-validation, which has been found to be a good compromise when trying to balance the bias and variance of the cross-validation criterion (Friedman, Hastie, and Tibshirani 2001).

The proposed cross-validation procedure does not take into account the datasets prior to the tth one. Unforeseen circumstances may have led to a realization of the tth dataset that is not representative. For instance, the response-covariate relation may be absent, that is, $\mathit{\beta }={\mathbf{0}}_p$, in this particular dataset. Cross-validation, which optimizes the predictive performance in the tth dataset, is likely to select a rather small λ_t that will not shrink toward the updated ridge regression estimator obtained from the previous t – 1 datasets. As such the information from these prior datasets accumulated in the updated ridge regression estimator is ignored. With the tth dataset being the odd-one-out, estimation of the regression parameter from the t + 1th dataset starts anew with an uninformative shrinkage target. To safeguard against such phenomena, we still optimize the cross-validated performance measure with respect to the penalty parameter but with a constraint on the vanishing of the fit on the preceding datasets. The minimization takes place over the set $\mathcal{D}$ defined as$$\begin{array}{c}\{{\lambda }_t>0:(1-f_t)K^{-1}{\sum }_{k=1}^K{\sum }_{\tau =1}^{t-1}||{\mathbf{Y}}_{\tau }-{\mathbf{X}}_{\tau }{\widehat{\mathit{\beta }}}_t^{(-k)}({\lambda }_t)|{|}_2^2\\ \hspace{1em}\le {\sum }_{\tau =1}^{t-1}\left|\right|{\mathbf{Y}}_{\tau }-{\mathbf{X}}_{\tau }{\widehat{\mathit{\beta }}}_{t-1}({\lambda }_{t-1})\left|{|}_2^2\right\}\end{array}$$with $f_t=n_t/({\sum }_{\tau =1}^t{n}_{\tau })$ denoting the sample size fraction of the tth dataset in relation to the sample size accumulated over all datasets up to t. In the optimization, we enforce the constraint using a barrier function. Constraining the choice of λ_t to $\mathcal{D}$ ensures that the fit of previously acquired dataset does not vanish more than the fraction f_t. In the presence of multiple targets the constraint simply demarcates a viable subspace of the penalty parameters λ_t and the ${\alpha }_{t,g}$’s. Finally, a similar constraint can be conceived for the choice of the ridge logistic regression estimator’s penalty parameter through cross-validation. The sum-of-squares are then to be replaced by the (negative) log-likelihoods.

To provide some intuition behind the proposed constrained cross-validation, we study the t = 2-case where the constraint, based on the estimate obtained from the first dataset, bounds the penalty parameter of the regression parameter estimate of the second dataset. Hereto we evaluate, using the moments of quadratic forms, the expectation of the quantities that determine the constraint$$\begin{array}{ccc}\hfill & \hfill & \mathbb{E}\left[\right||{\mathbf{Y}}_1-{\mathbf{X}}_1{\widehat{\mathit{\beta }}}_1({\lambda }_1)|{|}_2^2\left|{\lambda }_1\right]\hfill \\ \hfill & \hfill & \hspace{1em}=A_1+B_1:={\lambda }_1^2{({\mathit{\beta }}_0-\mathit{\beta })}^{\mathrm{\top }}{({\mathbf{X}}_1^{\mathrm{\top }}{\mathbf{X}}_1+{\lambda }_1{\mathbf{I}}_{pp})}^{-1}\hfill \\ \hfill & \hfill & \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}{\mathbf{X}}_1^{\mathrm{\top }}{\mathbf{X}}_1{({\mathbf{X}}_1^{\mathrm{\top }}{\mathbf{X}}_1+{\lambda }_1{\mathbf{I}}_{pp})}^{-1}({\mathit{\beta }}_0-\mathit{\beta })\hfill \\ \hfill & \hfill & \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}+{\sigma }^2\text{tr}\left\{{[{\mathbf{I}}_{n_1{n}_1}-{\mathbf{X}}_1{({\mathbf{X}}_1^{\mathrm{\top }}{\mathbf{X}}_1+{\lambda }_1{\mathbf{I}}_{pp})}^{-1}{\mathbf{X}}_1^{\mathrm{\top }}]}^2\right\},\hfill \end{array}$$

and$$\begin{array}{ccc}\hfill & \hfill & \mathbb{E}\left[\right||{\mathbf{Y}}_1-{\mathbf{X}}_1{\widehat{\mathit{\beta }}}_2^{(-k)}({\lambda }_2)|{|}_2^2|{\lambda }_1,{\lambda }_2]\hfill \\ \hfill & \hfill & \hspace{1em}=A_2^{(-k)}+B_2^{(-k)}+C_2^{(-k)}\hfill \\ \hfill & \hfill & \hspace{1em}:={\lambda }_1^2{\lambda }_2^2{({\mathit{\beta }}_0-\mathit{\beta })}^{\mathrm{\top }}{({\mathbf{X}}_1^{\mathrm{\top }}{\mathbf{X}}_1+{\lambda }_1{\mathbf{I}}_{pp})}^{-1}\hfill \\ \hfill & \hfill & \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}{({\left[{\mathbf{X}}_2^{(-k)}\right]}^{\mathrm{\top }}{\mathbf{X}}_2^{(-k)}+{\lambda }_2{\mathbf{I}}_{pp})}^{-1}{\mathbf{X}}_1^{\mathrm{\top }}{\mathbf{X}}_1\hfill \\ \hfill & \hfill & \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}{\{{\left[{\mathbf{X}}_2^{(-k)}\right]}^{\mathrm{\top }}{\mathbf{X}}_2^{(-k)}+{\lambda }_2{\mathbf{I}}_{pp}\}}^{-1}{({\mathbf{X}}_1^{\mathrm{\top }}{\mathbf{X}}_1+{\lambda }_1{\mathbf{I}}_{pp})}^{-1}({\mathit{\beta }}_0-\mathit{\beta })\hfill \\ \hfill & \hfill & \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}+{\sigma }^2\text{tr}\left\{\right[{\mathbf{I}}_{n_1,n_1}-{\lambda }_2{\mathbf{X}}_1{({\mathbf{X}}_1^{\mathrm{\top }}{\mathbf{X}}_1+{\lambda }_1{\mathbf{I}}_{pp})}^{-1}\hfill \\ \hfill & \hfill & \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}{({\left[{\mathbf{X}}_2^{(-k)}\right]}^{\mathrm{\top }}{\mathbf{X}}_2^{(-k)}+{\lambda }_2{\mathbf{I}}_{pp})}^{-1}{\mathbf{X}}_1^{\mathrm{\top }}]\hfill \\ \hfill & \hfill & \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}[{\mathbf{I}}_{n_1,n_1}-{\lambda }_2{\mathbf{X}}_1{({\mathbf{X}}_1^{\mathrm{\top }}{\mathbf{X}}_1+{\lambda }_1{\mathbf{I}}_{pp})}^{-1}\hfill \\ \hfill & \hfill & \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}{({\left[{\mathbf{X}}_2^{(-k)}\right]}^{\mathrm{\top }}{\mathbf{X}}_2^{(-k)}+{\lambda }_2{\mathbf{I}}_{pp})}^{-1}{\mathbf{X}}_1^{\mathrm{\top }}{]}^{\mathrm{\top }}\}\hfill \\ \hfill & \hfill & \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}+{\sigma }^2\text{tr}\{{\mathbf{X}}_2^{(-k)}{({\left[{\mathbf{X}}_2^{(-k)}\right]}^{\mathrm{\top }}{\mathbf{X}}_2^{(-k)}+{\lambda }_2{\mathbf{I}}_{pp})}^{-1}{\mathbf{X}}_1^{\mathrm{\top }}{\mathbf{X}}_1\hfill \\ \hfill & \hfill & \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}{({\left[{\mathbf{X}}_2^{(-k)}\right]}^{\mathrm{\top }}{\mathbf{X}}_2^{(-k)}+{\lambda }_2{\mathbf{I}}_{pp})}^{-1}{\left[{\mathbf{X}}_2^{(-k)}\right]}^{\mathrm{\top }}\}.\hfill \end{array}$$

We can now conclude that the set $\mathcal{D}$ is nonempty. Hereto set ${\lambda }_2=\infty$. Then $A_2^{(-k)}=A_1,B_2^{(-k)}=B_1$, and $C_2=0$. Consequently, the two expectations above are equal. This is no surprise as ${\widehat{\mathit{\beta }}}_2^{(-k)}({\lambda }_2)$ is shrunken toward $\widehat{\mathit{\beta }}({\lambda }_1)$ if ${\lambda }_2=\infty$. Furthermore, the expressions of the two preceding displays reveal the relationship between λ₁ and λ₂. While this relationship is not trivial, some conclusions may be deduced. For instance,

1. If n > p, the target ${\mathit{\beta }}_0$ is off, and λ₁ is (then naturally) close to zero. This implies $A_1\approx 0\approx A_2^{(-k)}$, while simultaneously, $B_2^{(-k)}\ge B_1$ and $C_2^{(-k)}\ge 0$. In particular, $B_2^{(-k)}$ and C₂ are strictly decreasing in λ₂. Hence, λ₂ has to be large and shrinks to the unbiased estimate of $\mathit{\beta }$ obtained from the first dataset.

2. If the target ${\mathit{\beta }}_0$ is close to $\mathit{\beta }$ and λ₁ large. Then, $A_1\approx A_2^{(-k)}$ and $B_1\approx B_2^{(-k)}$. To ensure $C_2^{(-k)}$ is close to zero the penalty parameter λ₂ has to be large too. Hence, ${\widehat{\mathit{\beta }}}_2^{(-k)}({\lambda }_2)$ is shrunken considerably (and correctly) to $\mathit{\beta }$ as well.

In both cases, the ${\widehat{\mathit{\beta }}}_1({\lambda }_1)$ is thus, a good “warm start.” In the first case it shields off the misinformed initial target ${\mathit{\beta }}_0$ and provides the unbiased estimate instead, while in the second one it propagates the well-informed initial target ${\mathit{\beta }}_0$. The scenario’s conclusions hinge upon the assumption that f₂ is small. If not, the multiplication factor $1-f_2$ in the constraint provides leverage and may allow λ₂ to deviate from these two extreme scenario’s. In particular, if n₁ is small, this leverage is considerable, hardly enforcing a constraint on the choice of λ₂.

Additionally, we have studied the behavior of constrained cross-validation in simulation (see supplementary material V). These simulations reveal that constrained cross-validation is beneficial (in the sense of yielding a reduction of the estimator’s loss) if (i) the novel data, on which the estimator is fitted, are from a different distribution then the “historic data,” while (ii) “historic data” are drawn from the correct model, and (iii) the target is informative. If either the target is uninformative or the novel data are from the correct model, the constrained cross-validation shows neither a gain over nor does it fall behind regular cross-validation. In all, constrained cross-validation is thus, only a safeguard against “outlying” novel data, outlying in the sense that it deviates from the historic data.

An informative choice on ${\mathit{\beta }}_0$ that initiates the updating may be obtained from profound knowledge of the relationship described by the linear regression model. Usually, however, such knowledge is at best present in tacit form. One way out could be literature providing univariate estimates of (some of) the β_j’s that constitute $\mathit{\beta }$. An alternative would “sacrifice” the first dataset to obtain an initial estimate with which the updating is initiated. Traditional ridge regression (with a zero target) or, following literature, univariate estimates from covariate-wise simple linear regression could be used to arrive at this initial estimate. Due to the shrinkage to zero, the former tends to underestimate the elements of $\mathit{\beta }$, whereas the latter has the opposite result (as confounding covariates are omitted). Alternatively, an estimate from other transfer regression learning procedures—preferably an unbiased one—can be used to initiate the updating.

Initialization is required not only at the initial time point. With the onset of the big data age, more and more information is registered over time. Consequently, “early” datasets comprise a limited number of covariates, while more recent ones comprise many additional ones. Or, it may be too costly to measure certain covariates at each instance but the budget allows for their measurement at (say) every other study. Hence, in practice not every covariate will be present in each dataset. A pragmatic choice of the target is then to use the latest estimates element-wise.

3 Comparisons

We conduct two comparisons. The first illustrates that updating is beneficial compared to de novo estimation. Secondly, we contrast updating, both in papyro and in silico, to estimation from the pooled data.

3.1 Regular Ridge Regression

In a small simulation study, we contrast the behavior of the updated ridge linear regression estimator to its traditional ridge counterpart. First, we generate an initial dataset from $\mathit{\beta }$ and estimate with the traditional ridge regression estimator in combination with a cross-validated penalty parameter. This estimate serves as the ${\mathit{\beta }}_0$ to initiate the updated ridge regression. Now sequentially, we draw another 25 datasets. From each dataset, we estimate $\mathit{\beta }$ using both the traditional and updated ridge estimators with cross-validation (without any constraint but positivity) for penalty parameter determination. For the latter, the previous updated ridge regression estimate is used as target. We have repeated this process a hundred times.

We have drawn all datasets as follows. Throughout we have set the jth element of the vector of regression coefficients $\mathit{\beta }$ as ${\beta }_j=(j-50)/20$ for $j=0,1,2,\dots ,100$. Each element of the $(n\times p)$-dimensional tth design matrix ${\mathbf{X}}_t$ with n = 25 is drawn from the standard normal distribution. We then form the response through ${\mathbf{Y}}_t={\mathbf{X}}_t\mathit{\beta }+{\mathit{\epsilon }}_t$ with each element of ${\mathit{\epsilon }}_t$ sampled from $\mathcal{N}(0,0.04)$. Hence, we generate ${\mathbf{X}}_t$, ${\mathit{\epsilon }}_t$, and consequently ${\mathbf{Y}}_t$, anew at each t.

In Figure 1 the 5%, 50%, and 95% quantiles of the traditional and updated ridge estimates of β_j with $j\in \{0,30,50,70,100\}$ are plotted against t. These quantiles of the traditional ridge estimates of these elements of $\mathit{\beta }$ are constant over t (left panel of Figure 1). Those of the updated ridge estimates (right panel of Figure 1) clearly improve as t increases. The improvement is 2-fold: (i) they become less biased, and (ii) the distance between the 5% and 95% quantiles vanishes.

Fig. 1 The top panels show the $(5\%,95\%$-quantile intervals of the traditional (left) and updated (right) ridge estimates of β_j with $j\in \{0,30,50,70,100\}$ plotted against t. The solid, colored line inside these intervals is the corresponding 50% quantile. The dotted, gray lines are the true values of the β_j’s. The bottom panels show the mean squared errors over hundred runs of the mixed model’s ML, traditional (based on all data) and updated ridge regression estimators (initiated with ${\mathit{\beta }}_0={\mathbf{0}}_p$ and ${\mathit{\beta }}_0=\mathit{\beta }$) plotted against t. All datasets are sampled from the same regression model, but in the right panel every tenth dataset is sampled from an empty model.

We have conducted additional simulations with the elements of design matrix sampled from more intricate distribution (see supplementary material VI). They reveal that a varying location of the covariates does not affect the performance of the updated ridge regression estimator. The introduction of collinearity among the covariates, however, does, as it slows down the convergence of the updated ridge regression estimator to the true parameter values.

3.2 Mixed Model

A natural and widespread alternative to the proposed approach would be to employ a mixed model that includes a dataset-specific random effect of the covariates. The mixed model is then refitted at the arrival of a new dataset, each time to an enlarged dataset encompassing the newly arrived data and all those that preceded it. This may be computationally demanding but is considered a mere practical nuisance for the moment. To make matters more precise, we assume, for the dataset formed from all those up to the tth one, the mixed model ${\mathbb{Y}}_t={\mathbb{X}}_t\mathit{\beta }+{\mathbb{Z}}_t{\mathbb{G}}_t+{\mathrm{e}}_t$ with:

• the vector ${\mathbb{Y}}_t={({\mathbf{Y}}_1^{\mathrm{\top }},\dots ,{\mathbf{Y}}_t^{\mathrm{\top }})}^{\mathrm{\top }}$ of length ${\tilde{n}}_t$-dimension where ${\tilde{n}}_t=n_1+\cdots +n_t$,

• the $({\tilde{n}}_t\times p)$-dimensional matrix ${\mathbb{X}}_t={({\mathbf{X}}_1^{\mathrm{\top }},\dots ,{\mathbf{X}}_t^{\mathrm{\top }})}^{\mathrm{\top }}$,

• the $({\tilde{n}}_t\times tp)$-dimensional matrix ${\mathbb{Z}}_t$ constructed as a t × t block matrix with $\text{diag}({\mathbb{Z}}_t)=({\mathbf{X}}_1,\dots ,{\mathbf{X}}_t)$ and off-diagonal blocks all equaling the zero matrix of appropriate dimensions,

• the vector ${\mathbb{G}}_t={({\mathit{\gamma }}_1^{\mathrm{\top }},\dots ,{\mathit{\gamma }}_t^{\mathrm{\top }})}^{\mathrm{\top }}$ of length $tp$ with the dataset-specific random effects, and

• the error vector ${\mathrm{e}}_t={({\mathit{\epsilon }}_1^{\mathrm{\top }},\dots ,{\mathit{\epsilon }}_t^{\mathrm{\top }})}^{\mathrm{\top }}$ of length ${\tilde{n}}_t$.

We assume the random effects and errors to be independent and obey the multivariate normal laws, ${\mathit{\gamma }}_t\sim \mathcal{N}({\mathbf{0}}_p,{\sigma }_{\gamma }^2{\mathbf{I}}_{pp})$ and ${\mathit{\epsilon }}_t\sim \mathcal{N}({\mathbf{0}}_n,{\sigma }_{\epsilon }^2{\mathbf{I}}_{nn})$ for all t, and with zero covariance across datasets, for example, $\text{cov}({\mathit{\epsilon }}_t,{\mathit{\epsilon }}_{t\prime })={\mathbf{0}}_{nn}$ for $t\cancel{=}t\prime$.

The maximum likelihood estimator of the above formulated mixed model’s fixed regression parameter $\mathit{\beta }$ given the variance parameters ${\sigma }_{\epsilon }^2$ and ${\sigma }_{\delta }^2$ is (see Bates and DebRoy 2004):$$\begin{array}{ccc}{\widehat{\mathit{\beta }}}_t^{\mathrm{(}\text{me}\mathrm{)}}\hfill & =\hfill & {\left[{\mathbb{X}}^{\mathrm{\top }}{(\xi \mathbb{Z}{\mathbb{Z}}^{\mathrm{\top }}+{\mathbf{I}}_{{\tilde{n}}_t,{\tilde{n}}_t})}^{-1}\mathbb{X}\right]}^{-1}{\mathbb{X}}^{\mathrm{\top }}{(\xi \mathbb{Z}{\mathbb{Z}}^{\mathrm{\top }}+{\mathbf{I}}_{{\tilde{n}}_t,{\tilde{n}}_t})}^{-1}\mathbb{Y},\hfill \end{array}$$with its first two moments:$$\begin{array}{ccc}\mathbb{E}\left[{\widehat{\mathit{\beta }}}_t^{\mathrm{(}\text{me}\mathrm{)}}\right]\hfill & =\hfill & {\left[{\mathbb{X}}^{\mathrm{\top }}{(\xi \mathbb{Z}{\mathbb{Z}}^{\mathrm{\top }}+{\mathbf{I}}_{{\tilde{n}}_t,{\tilde{n}}_t})}^{-1}\mathbb{X}\right]}^{-1}{\mathbb{X}}^{\mathrm{\top }}{(\xi \mathbb{Z}{\mathbb{Z}}^{\mathrm{\top }}+{\mathbf{I}}_{{\tilde{n}}_t,{\tilde{n}}_t})}^{-1}{\mathbb{X}}_t\mathit{\beta },\hfill \\ \text{var}\left[{\widehat{\mathit{\beta }}}_t^{\mathrm{(}\text{me}\mathrm{)}}\right]\hfill & =\hfill & {\left[{\mathbb{X}}^{\mathrm{\top }}{(\xi \mathbb{Z}{\mathbb{Z}}^{\mathrm{\top }}+{\mathbf{I}}_{{\tilde{n}}_t,{\tilde{n}}_t})}^{-1}\mathbb{X}\right]}^{-1}{\mathbb{X}}^{\mathrm{\top }}{(\xi \mathbb{Z}{\mathbb{Z}}^{\mathrm{\top }}+{\mathbf{I}}_{{\tilde{n}}_t,{\tilde{n}}_t})}^{-1}\hfill \\ \hfill & \hfill & ({\sigma }_{\epsilon }^2{\mathbf{I}}_{{\tilde{n}}_t,{\tilde{n}}_t}+{\sigma }_{\gamma }^2\mathbb{Z}{\mathbb{Z}}^{\mathrm{\top }}){(\xi \mathbb{Z}{\mathbb{Z}}^{\mathrm{\top }}+{\mathbf{I}}_{{\tilde{n}}_t,{\tilde{n}}_t})}^{-1}\hfill \\ \hfill & \hfill & \mathbb{X}{\left[{\mathbb{X}}^{\mathrm{\top }}{(\xi \mathbb{Z}{\mathbb{Z}}^{\mathrm{\top }}+{\mathbf{I}}_{{\tilde{n}}_t,{\tilde{n}}_t})}^{-1}\mathbb{X}\right]}^{-1},\hfill \end{array}$$ in which $\xi ={\sigma }_{\epsilon }^{-2}{\sigma }_{\gamma }^2$.

We now compare the updated ridge regression estimator and the maximum likelihood estimator of the mixed model’s fixed effect parameter with respect to their mean squared error. Theorem 4 does so for a sequence of studies with orthonormal design matrices.

Theorem 4

(Mean squared error of mixed versus updated estimator). Assume the existence of an infinite sequence of studies into the linear relationship between a continuous response and a set of covariates. The data from these studies, ${\{{\mathbf{X}}_t,{\mathbf{Y}}_t\}}_{t=1}^{\infty }$, are used to fit the linear regression model by means of the updated ridge linear regression estimator, which yields the sequence of estimators ${\{{\widehat{\mathit{\beta }}}_t({\lambda }_t,{\widehat{\mathit{\beta }}}_{t-1})\}}_{t=1}^{\infty }$ which is initiated by an arbitrary, nonrandom ${\mathit{\beta }}_0$. Let ${\mathbf{X}}_t$ be orthonormal and ${\lambda }_t>{\sigma }_{\epsilon }{({\sigma }_{\epsilon }^2+{\sigma }_{\gamma }^2)}^{-1/2}{2}^{t/2}{T}^{1/2}$ for $1\le t\le T$. Then, when initiated by ${\widehat{\mathit{\beta }}}_1({\lambda }_1)={\widehat{\mathit{\beta }}}_1^{\mathrm{(}\text{me}\mathrm{)}}$, the updated ridge regression estimator outperforms (in the mean squared error sense), the maximum likelihood estimator of the mixed model’s fixed effects parameter as given above. Hence, MSE$[{\widehat{\mathit{\beta }}}_T({\lambda }_T)|{\left\{{\lambda }_t\right\}}_{t=1}^T]<\text{MSE}[{\widehat{\mathit{\beta }}}_T^{\mathrm{(}\text{me}\mathrm{)}}]$.

This theorem states that the updated linear regression estimator may, when initiated appropriately and equipped with the right penalty parameter scheme, outperform in the MSE sense the maximum likelihood estimator of the mixed model’s fixed effects parameter. Theorem 4) also holds if ${\sigma }_{\gamma }^2=0$. The reasons for this are 2-fold. The first prerequisite is an unbiased initiation of the sequence of updated ridge regression estimators. The updated estimators are then unbiased too. As Theorem 4 is concerned with the MSE, and the mixed model’s ML estimator is unbiased too, this leaves to compare the variances. The condition on the sequence of penalty parameters then warrants that the updated ridge regression estimator’s variance shrinks faster than that of the mixed model’s ML estimator as t increases.

Theorem 4 assumes an orthonormal design matrix for all studies. This is rather restrictive from a practical perspective. In principle, however, a result as in Theorem 4 can be obtained for ${\mathbf{X}}_t=\mathbf{X}$ with X of sufficient rank. The proof follows that of Theorem 4, but the mathematics is more cumbersome and the result less insightful.

Theorem 4 assumes the variance parameters ${\sigma }_{\epsilon }^2$ and ${\sigma }_{\delta }^2$ known. In practice, these parameters too need to be estimated from data. This is done in the simulation study that follows.

We compare the performance of the presented updated ridge regression estimator to the maximum likelihood estimator of the mixed model’s fixed effect parameter in two simulations. For reference, the comparison also includes the traditional ridge regression estimator using all available data. The first simulation investigates the result of Theorem 4 under more realistic assumptions. The second simulation reveals a better performance of the updated ridge regression estimator than the maximum likelihood estimator of the mixed model’s fixed effect parameter, when for a particular study in the sequence accidentally data from a different model are sampled. This accidentally erroneous sampling from a different model will introduce bias in the latter estimator, as it cannot reduce the influence of the data from this study. The updated ridge regression estimator, however, decides by the size of the penalty parameter how it weighs the current data against the regression estimate obtained from the preceding studies. In principle, it may favor the latter, effectively ignoring the data from the erroneously sampled study. This may cause a small reduction in efficiency compared to its competitor, but this is expected to be outweighed by the avoided bias. Both simulations are modified from that presented in Section 3.1. The settings are identical except for the fact that ${\sigma }_{\epsilon }^2=1$. Moreover, in a second scenario the tth dataset with $t\in \{t:t\hspace{1em}\mathrm{mod}10=0\}$ is sampled from an empty model, that is, ${\mathbf{Y}}_t={\mathit{\epsilon }}_t$. With the arrival of each new dataset the parameter $\mathit{\beta }$ is estimated using the ML estimator of the mixed model’s fixed parameter (from that t when the accumulated sample size exceeds the dimension). But also using the updated ridge regression estimator initiated uninformatively with ${\mathit{\beta }}_0={\mathbf{0}}_p$ and informatively with ${\mathit{\beta }}_0=\mathit{\beta }$. We choose its penalty parameter through constrained cross-validation to prevent the deterioration of the fit on the preceding datasets. The quadratic loss, for example, $\left|\right|{\widehat{\mathit{\beta }}}_t[{\lambda }_t,{\widehat{\mathit{\beta }}}_{t-1}({\lambda }_{t-1},{\widehat{\mathit{\beta }}}_{t-2})]-\mathit{\beta }|{|}_2^2$, of the mixed model’s fixed parameter ML, traditional (using all data) and updated ridge regression estimators are plotted against t (bottom panels of Figure 1). In the first scenario with all datasets drawn from the same linear regression model the maximum likelihood estimator of the mixed model’s fixed parameter performs better than the uninformatively initiated updated ridge estimator, but worse than the informatively initiated one. The same holds for the traditional ridge regression estimator based on all data, which in turn performs slightly better than the mixed model as it involves fewer—only fixed—parameters. For small t, this picture is also seen in the second scenario with every tenth dataset sampled from an empty model. But the maximum likelihood estimator of the mixed model’s fixed parameter suffers most (relatively) from the “outlying” datasets, and is for larger t even overtaken by the uninformatively initiated updated ridge regression estimator. Again, the traditional ridge regression estimator based on all data shows similar behavior as the mixed model, but roles between them are now reversed due to the flexibility of the latter by virtue of its random effects.

4 Application

The fingertips study registers annual information on a range of public health indicators for the counties of England (https://fingertips.phe.org. uk/). Following van Schaik et al. (2019), who advocate the use of the fingertips data for the illustration of novel statistical techniques, we apply the proposed transfer learning procedure to relate the counties’ suicide rate to the other health indicators. Building on the R-script of van Schaik et al. (2019), we have downloaded the data using the fingertipsR-package (Fox and Flowers 2009). Preprocessing of the data amounts to removal of cases without full information on the health indicators registered. This is done per year as the set of registered health indicators varies over the years. In particular, the size of this set tends to increase over time. Moreover, we have zero-centered the data year-wise. The resulting data comprise the years 2008–2016. Over these years the number of health indicators ranges from 1 to 23, from a total of 23 different ones, while the number of counties with full case information ranges from 57 to 159 over the years (refer to the supplementary material VII for more details).

Analysis of the data commences with the first year available, that is, 2008. We have fitted the linear model with the targeted ridge regression estimator with a zero target, that is, ${\mathit{\beta }}_0={\mathbf{0}}_p$, and chosen the penalty parameter by cross-validation. For subsequent years, we have used the same estimator, but have formed the target element-wise from previous estimates: the jth element of the target is taken from the most recent estimate. If the jth covariate was not present in the data of the preceding year, we have taken the corresponding estimate from the year before that, and so on, until it is taken from the initial target for the year 2008. Moreover, we have chosen the penalty parameter by unrestricted and constrained cross-validation (as discussed in Section 2.4). Finally, for the purpose of reference, we have also estimated the regression parameter year-wise with the maximum likelihood regression estimator.

The trajectories of the resulting estimates are shown in Figure 2 and the supplementary material VII. The main takeaway of the plots is the volatile behavior of the ML regression estimates. In part this is due to a different sample size and, probably more important, the varying set of health indicators registered each year. The updated ridge regression estimator, irrespectively of the employed cross-validation method, yields much smoother trajectories. Hence, providing last year’s parameter estimate as a suggestion for the estimation of the current year’s one appears to harness against the involvement of a different set of covariates. Indirectly, the usefulness of this suggestion for the learning of the updating of the estimate is also evident by the increase of both the un- and constrained cross-validated penalty parameters over the years (not shown). This indicates that the suggestion is becoming more relevant as years go by. A convenient consequence of the updated ridge regression estimator’s smooth trajectories is the consistency of the sign of the regression parameter’s estimates over the years, most obvious from the estimators’ trajectories for two elements of the regression parameter (provided in the supplementary material VII), facilitating a sensible interpretation. Of course, there is no free lunch. The updated ridge regression estimators exhibit a small deterioration in the fit as can be witnessed from the residuals (see supplementary material VII), although this appears to be negligible.

Fig. 2 The panels show the trajectories of the ML (left) and the updated ridge regression, with its penalty parameter chosen via constrained LOOCV, estimates. Each trajectory represents a single covariate. The presence of a health indicator in the data of a particular year is evident from a symbol on its trajectory at the corresponding year. The symbol is omitted in years that the health indicator was not registered.

5 Conclusion

We presented methodology for continuously learning regression models from studies executed over large periods of time with data coming in bit-by-bit as they run their course. It could be considered a frequentist analogy to Bayesian updating and comprises sequential targeted ridge penalized estimation, shrinking toward the latest estimate, when novel data become available. At each update, the penalty parameter represents the extent to which the latest estimate yields an adequate model of the novel data. The iteratively updated estimator and its associated linear predictor have been shown to be asymptotically unbiased and consistent. The penalty parameter is chosen through cross-validation in which the search domain of the penalty parameter is constrained to ensure that the newly learned estimate causes little to no vanishing of the fit of the historic data in comparison to the previous estimate. In a comparison with other standard statistical transfer learning procedures, situations were identified where these procedures were outperformed by the proposed method. Finally, the proposed method was illustrated on data from an epidemiology study, showing its promise.

The current incremental/online transfer learning proposal employs a single penalty parameter per update. Hence, all elements of the novel estimate are shrunken in the same amount to the current one. However, it may be that certain elements of the current estimate provide a better fit than others on the novel data. Then, to differentiate the shrinkage among the novel estimate’s elements, the single-parameter penalty is to be replaced by a multi-parameter one: ${[\mathit{\beta }-{\widehat{\mathit{\beta }}}_{t-1}({\mathbf{\Lambda }}_{t-1},{\widehat{\mathit{\beta }}}_{t-1})]}^{\mathrm{\top }}{\mathbf{\Lambda }}_t[\mathit{\beta }-{\widehat{\mathit{\beta }}}_{t-1}({\mathbf{\Lambda }}_{t-1},{\widehat{\mathit{\beta }}}_{t-1})]$, the penalty parameter ${\mathbf{\Lambda }}_t$ is now a matrix. This matrix may be unstructured, leaving one with the problem of selecting $\frac{1}{2}p(p+1)$ penalty parameters. More practically, it may be diagonal and possibly parameterized by a few scalar penalty parameters, thereby shrinking groups of elements similarly. Alternatively, a parameterization in line with generalized ridge regression (Hoerl and Kennard 1970; Hemmerle 1975; Lawless 1981), only shrinking differentially in the p canonical directions, may be considered. Either choice provides greater flexibility in the information propagated between updates.

Supplementary Materials

Proofs: Proofs of theoretical results. (.pdf-file).

Illustration: R code for the analysis (.r-files).

Illustration: Plots (.pdf-file).