Data enriched linear regression

We present a linear regression method for predictions on a small data set making use of a second possibly biased data set that may be much larger. Our method fits linear regressions to the two data sets while penalizing the difference between predictions made by those two models. The resulting algorithm is a shrinkage method similar to those used in small area estimation. We find a Stein-type finding for Gaussian responses: when the model has 5 or more coefficients and 10 or more error degrees of freedom, it becomes inadmissible to use only the small data set, no matter how large the bias is. We also present both plug-in and AICc-based methods to tune our penalty parameter. Most of our results use an $L_2$ penalty, but we obtain formulas for $L_1$ penalized estimates when the model is specialized to the location setting. Ordinary Stein shrinkage provides an inadmissibility result for only 3 or more coefficients, but we find that our shrinkage method typically produces much lower squared errors in as few as 5 or 10 dimensions when the bias is small and essentially equivalent squared errors when the bias is large.


Introduction
The problem we consider here is how to combine linear regressions based on data from two sources. There is a small data set of expensive high quality observations and a possibly much larger data set with less costly observations. The big data set is thought to have similar but not identical statistical characteristics to the small one. The conditional expectation might be different there or the predictor variables might have been measured in somewhat different ways. The motivating application comes from within Google. The small data set is a panel of consumers, selected by a probability sample, who are paid to share their internet viewing data along with other data on television viewing. There is a second and potentially much larger panel, not selected by a probability sample who have opted in to the data collection process.
The goal is to make predictions for the population from which the smaller sample was drawn. If the data are identically distributed in both samples, we should simply pool them. If the big data set is completely different from the small one, then using it may not be worth the trouble.
Many settings are intermediate between these extremes: the big data set is similar but not necessarily identical to the small one. We stand to benefit from using the big data set at the risk of introducing some bias. Our goal is to glean some information from the larger data set to increase accuracy for the smaller one. The difficulty is that our best information about how the two populations are similar is in our samples from them.
The motivating problem at Google has some differences from the problem we consider here. There were two binary responses, one sample was missing one of those responses, and tree-based predictions were used. See Chen et al. (2014). This paper studies linear regression because it is more amenable to theoretical analysis, is more fundamental, and sharper statements are possible.
The linear regression method we use is a hybrid between simply pooling the two data sets and fitting separate models to them. As explained in more detail below, we apply shrinkage methods penalizing the difference between the regression coefficients for the two data sets. Both the specific penalties we use, and our tuning strategies, reflect our greater interest in the small data set. Our goal is to enrich the analysis of the smaller data set using possibly biased data from the larger one. Section 2 presents our notation and introduces L 1 and L 2 penalties on the parameter difference. Most of our results are for the L 2 penalty. For the L 2 penalty, the resulting estimate is a linear combination of the two within sample estimates. Theorem 2.1 gives a formula for the degrees of freedom of that estimate. Theorem 2.2 presents the mean squared error of the estimator and forms the basis for plug-in estimation of an oracle's value when an L 2 penalty is used. We also show how to use Stein shrinkage, shrinking the regression coefficient in the small sample towards the estimate from the large sample. Such shrinkage makes it inadmissible to ignore the large sample when there are 3 or more coefficients including the intercept.
Section 3 considers in detail the case where the regression simplifies to estimation of a population mean. In that setting, we can determine how plug-in, bootstrap and cross-validation estimates of tuning parameters behave. We get an expression for how much information the large sample can add. Theorem 3.1 gives a soft-thresholding expression for the estimate produced by L 1 penalization and equation (3.7) can be used to find the penalty parameter that an L 1 oracle would choose when the data are Gaussian.
Section 4 presents some simulated examples. We simulate the location problem for several L 2 penalty methods varying in how aggressively they use the larger sample. The L 1 oracle is outperformed by the L 2 oracle in this setting. When the bias is small, the data enrichment methods improve upon the small sample, but when the bias is large then it is best to use the small sample only. Things change when we simulate the regression model. For dimension d ≥ 5, data enrichment outperforms the small sample method in our simulations at all bias levels. We did not see such an inadmissibility outcome when we simulated cases with d ≤ 4. In our simulated examples, the data enrichment estimator performs better than plain Stein shrinkage of the small sample towards the large sample.
Section 5 presents theoretical support for our estimator. Theorem 5.1. shows that when there are 5 or more predictors and 10 or more degrees of freedom for error, then some of our data enrichment estimators make simply using the small sample inadmissible. The reduction in mean squared error is greatest when the bias is smallest, but no matter how large the bias is, we gain an improvement. The estimator we study employs a data-driven weighting of the two withinsample least squares estimators. In simulations, our plug-in estimator performed even better than the estimator from Theorem 5.1. Section 6 explains how our estimators are closer to a matrix oracle than the James-Stein estimators are, and this may explain why they outperform simple shrinkage in our simulations.
There are many statistical settings where data from one setting is used to study a different setting. They range from older methods in survey sampling, to recently developed methods for bioinformatics. Section 7 surveys some of those literatures. Section 8 has brief conclusions. The longer proofs are in Section 9 of the Appendix.
Our contributions include the following: • a new penalization method for combining data sets, • an inadmissibility result based on that method, • a comparison of L 1 and L 2 penalty oracles for the location setting, and • evidence that more aggressive shrinkage pays in high dimensions.

Data enriched regression
Consider linear regression with a response Y ∈ R and predictors X ∈ R d . The model for the small data set is for a parameter β ∈ R d and independent errors ε i with mean 0 and variance σ 2 S . Now suppose that the data in the big data set follow where γ ∈ R d is a bias parameter and ε i are independent with mean 0 and variance σ 2 B . The sample sizes are n in the small sample and N in the big sample.
There are several kinds of departures of interest. It could be, for instance, that the overall level of Y is different in S than in B but that the trends are similar. That is, perhaps only the intercept component of γ is nonzero. More generally, the effects of some but not all of the components in X may differ in the two samples. One could apply hypothesis testing to each component of γ but that is unattractive as the number of scenarios to test for grows as 2 d .
Let X S ∈ R n×d and X B ∈ R N ×d have rows made of vectors X i for i ∈ S and i ∈ B respectively. Similarly, let Y S ∈ R n and Y B ∈ R N be corresponding vectors of response values. We use V S = X T S X S and V B = X T B X B .

Partial pooling via shrinkage and weighting
Our primary approach is to pool the data but put a shrinkage penalty on γ. We estimate β and γ by minimizing where λ ∈ [0, ∞] and P (γ) ≥ 0 is a penalty function. There are several reasonable choices for the penalty function, including γ 2 2 , X S γ 2 2 , γ 1 , and X S γ 1 .
For each of these penalties, setting λ = 0 leads to separate fitsβ andβ +γ in the two data sets. Similarly, taking λ = ∞ constrainsγ = 0 and amounts to pooling the samples. We will see that varying λ shifts the relative weight applied to the two samples. In many applications one will want to regularize β as well, but in this paper we only penalize γ. The criterion (2.1) does not account for different variances in the two samples. If we knew the variance ratio we might multiply the sum over i ∈ B by τ = σ 2 S /σ 2 B . A large value of τ has the consequence of increasing the weight on the B sample. Choosing τ is largely confounded with choosing λ because λ also adjusts the relative weight of the two samples. We use τ = 1, but in a context with some prior knowledge about the magnitude of σ 2 S /σ 2 B another value could be used. Our inadmissibility result does not depend on knowing the correct τ .
The L 1 penalties have an advantage in interpretation because they identify which parameters or which specific observations might be differentially affected. The quadratic penalties are more analytically tractable, so we focus most of this paper on them.
Both quadratic penalties can be expressed as X T γ 2 2 for a matrix X T . The rows of X T represent a hypothetical target population of N T items for prediction. The matrix V T = X T T X T is then proportional to the matrix of mean squares and mean cross-products for predictors in the target population.
If we want to remove the pooling effect from one of the coefficients, such as the intercept term, then the corresponding column of X T should contain all zeros. We can also constrain γ j = 0 (by dropping its corresponding predictor) in order to enforce exact pooling on the j'th coefficient.
A second, closely related approach is to fitβ S by minimizing i∈S (Y i −X i β) 2 , fitβ B by minimizing i∈B (Y i − X i β) 2 , and then estimate β bŷ for some 0 ≤ ω ≤ 1. In some special cases the estimates indexed by the weighting parameter ω ∈ [n/(n + N ), 1] are a relabeling of the penalty-based estimates indexed by the parameter λ ∈ [0, ∞]. In other cases, the two families of estimates differ. The weighting approach allows simpler tuning methods. Although we find in simulations that the penalization method is superior, we can prove stronger results about the weighting approach.
Given two values of λ we consider the larger one to be more 'aggressive' in that it makes more use of the big sample bringing with it the risk of more bias in return for a variance reduction. Similarly, aggressive estimators correspond to small weights ω on the small target sample.

Quadratic penalties and degrees of freedom
The quadratic penalty takes the form P (γ) = X T γ 2 2 = γ T V T γ for a matrix X T ∈ R r×d and V T = X T T X T ∈ R d×d . The value r is d or n in the examples above and could take other values in different contexts. Our criterion becomes (2.2) Here and below x means the Euclidean norm x 2 . Given the penalty matrix X T and a value for λ, the penalized sum of squares (2.2) is minimized byβ λ andγ λ satisfying To avoid uninteresting complications we suppose that the matrix X T X is invertible. The representation (2.3) also underlies a convenient computational approach to fittingβ λ andγ λ using r rows of pseudo-data just as one does in ridge regression.
The estimateβ λ can be written in terms ofβ as the next lemma shows. Lemma 2.1. Let X S , X B , and X T in (2.2) all have rank d. Then for any λ ≥ 0, the minimizersβ andγ of (2.2) satisfŷ Proof. The normal equations of (2.2) are Solving the second equation forγ, plugging the result into the first and solving forβ, yields the result with This expression for W λ simplifies as given and simplifies further when The remaining challenge in model fitting is to choose a value of λ. Because we are only interested in making predictions for the S data, not the B data, the ideal value of λ is one that optimizes the prediction error on sample S. One reasonable approach is to use cross-validation by holding out a portion of sample S and predicting the held-out values from a model fit to the held-in ones as well as the entire B sample. One may apply either leave-one-out cross-validation or more general K-fold cross-validation. In the latter case, sample S is split into K nearly equally sized parts and predictions based on sample B and K − 1 parts of sample S are used for the K'th held-out fold of sample S.
We prefer to use criteria such as AIC, AICc, or BIC in order to avoid the cost and complexity of cross-validation. To compute AIC and alternatives, we need to measure the degrees of freedom used in fitting the model. We define the degrees of freedom to be whereŶ S = X Sβλ . This is the formula of Ye (1998) and Efron (2004) adapted to our setting where the focus is only on predictions for the S data. We will see later that the resulting AIC type estimates based on the degrees of freedom perform similarly to our focused cross-validation described above.
Theorem 2.1. For data enriched regression the degrees of freedom given at (2.5) satisfies df(λ) = tr(W λ ) where W λ is given in Lemma 2.1.
where ν 1 , . . . , ν d are the eigenvalues of Proof. Section 9.1 in the Appendix.
With a notion of degrees of freedom customized to the data enrichment context we can now define the corresponding criteria such as AIC(λ) = n log(σ 2 S (λ)) + n 1 + 2df(λ) n and AICc(λ) = n log(σ 2 S (λ)) + n 1 + The AIC is more appropriate than BIC here since our goal is prediction accuracy, not model selection. We prefer the AICc criterion of Hurvich and Tsai (1989) because it is more conservative as the degrees of freedom become large compared to the sample size.
Next we illustrate some special cases of the degrees of freedom formula in Theorem 2.1. First, suppose that λ = 0, so that there is no penalization on γ. Then df(0) = tr(I) = d as is appropriate for regression on sample S only.
We can easily see that the degrees of freedom are monotone decreasing in λ. As λ → ∞ the degrees of freedom drop to df(∞) = d j=1 ν j /(1 + ν j ). This can be much smaller than d. For instance if V S = nΣ and V B = N Σ for some positive definite Σ ∈ R d×d , then all ν j = n/N and so df(∞) = d/(1 + N/n) ≤ dn/N .
Monotonicity of the degrees of freedom makes it easy to search for the value λ which delivers a desired degrees of freedom. We have found it useful to investigate λ over a numerical grid corresponding to degrees of freedom decreasing from d by an amount ∆ (such as 0.25) to the smallest such value above df(∞). It is easy to adjoin λ = ∞ (sample pooling) to this list as well.

Predictive mean square errors
Here we develop an oracle's choice for λ and a corresponding plug-in estimate. We work in the case where V S = V T and we assume that V S has full rank. Given λ, the predictive mean square error is E( X S (β − β) 2 ).
We will use the matrices V 1/2 S and M from Theorem 2.1 and the eigendecomposition M = U DU T where the j'th column of U is u j and D = diag(ν j ).
Theorem 2.2. The predictive mean square error of the data enrichment estimator is The first term in (2.9) is a variance term. It equals dσ 2 S when λ = 0 but for λ > 0 it is reduced due to the use of the big sample. The second term represents the error, both bias squared and variance, introduced by the big sample.

A plug-in method
A natural choice of λ is to minimize the predictive mean square error, which must be estimated. We propose a plug-in method that replaces the unknown parameters σ 2 S and κ 2 j from Theorem 2.2 by sample estimates. For estimatesσ 2 S andκ 2 j we chooseλ (2.10) From the sample data we takeσ 2 where the positive part operation on a symmetric matrix preserves its eigenvectors but replaces any negative eigenvalues by 0. Similar results can be obtained with This estimator is somewhat simpler but (2.11) has the advantage of being at least as large asσ 2 B V −1 B while (2.12) can degenerate to 0.

James-Stein shrinkage estimators
For background on these estimators, see Efron and Morris (1973b). We shrink towards a target vector, to get better estimators of θ S = V 1/2 S β. To make use of the big data set we shrinkθ S towardŝ We consider two shrinkerŝ Each of these makesθ S inadmissible in squared error loss as an estimate of θ S , when d ≥ 3. The squared error loss on the θ scale is (2.14) When d ≥ 3 and our quadratic loss is based on V S , we can makeβ S inadmissible by shrinkage, so long as d ≥ 3. Copas (1983) found that ordinary least squares regression is inadmissible when d ≥ 4. Stein (1960) also obtained an inadmissibility result for regresssion, but under stronger conditions than Copas needs. Copas (1983) applies no shrinkage to the intercept but shrinks the rest of the coefficient vector towards zero. In this problem it is reasonable to shrink the entire coefficient vector as the big data set supplies a nonzero default intercept.

The location model
The simplest instance of our problem is the location model where X S is a column of n ones and X B is a column of N ones. Then the vector β is simply a scalar intercept that we call µ and the vector γ is a scalar mean difference that we call δ. The response values in the small data set are Y i = µ + ε i while those in the big data set are Y i = (µ + δ) + ε i . In the location family we lose no generality taking the quadratic penalty to be λδ 2 .
The quadratic criterion is i∈S ( Choosing a value for ω corresponds to choosing . The degrees of freedom in this case reduce to df(λ) = ω, which ranges from df(0) = 1 down to df(∞) = n/(n + N ).

Oracle estimator of ω
The mean square error ofμ(ω) is The mean square optimal value of ω (available to an oracle) is Pooling the data corresponds to ω pool = n/(N +n) and makesμ equal the pooled meanȲ P ≡ (nȲ S + NȲ B )/(n + N ). Ignoring the large data set corresponds to ω S = 1. Here ω pool ≤ ω orcl ≤ ω S . The mean squared error reduction for the oracle is after some algebra. If δ = 0, then as min(n, N ) → ∞ we find ω orcl → 1 and the optimal ω corresponds to simply using the small sample and ignoring the large one. If instead δ = 0 and N → ∞ for finite n, then the effective sample size for data enrichment may be defined using (3.1) as The mean squared error from data enrichment with n observations in the small sample, using the oracle's choice of λ, matches that of n IID observations from the small sample. We effectively gain up to σ 2 S /δ 2 observations worth of information. This is an upper bound on the gain because we will have to estimate λ.
Equation (3.2) shows that the benefit from data enrichment is a small sample phenomenon. The effect is additive not multiplicative on the small sample size n. As a result, more valuable gains are expected in small samples. In some of the motivating examples we have found the most meaningful improvements from data enrichment on disaggregated data sets, such as specific groups of consumers.
If we bootstrap the S and B samples independently M times and choose ω to then the minimizing value tends to ω plug as M → ∞. Thus bootstrap methods give an approach analogous to plug-in methods, when no simple plug-in formula exists. We can also determine the effects of cross-validation in the location setting, and arrive at an estimate of ω that we can use without actually cross-validating.
Consider splitting the small sample into K parts that are held out one by one in turn. The K − 1 retained parts are used to estimate µ and then the squared error is judged on the held-out part. That is If n is a multiple of K and we average over all of the K-fold sample splits we might use, then an analysis in Section 9.3 shows that K-fold cross-validation chooses a weighting centered around .
( 3.3) Cross-validation allows ω < 0. This can arise when the bias is small and then sampling alone makes the held-out part of the small sample appear negatively correlated with the held-in part. The effect can appear with any K. We replace any ω cv,K < n/(n + N ) by n/(n + N ).
Leave-one-out cross-validation has K = n (and r = 1) so it chooses a weight centered around ω cv,n = [δ 2 Smaller K, such as choosing K = 10 versus n, tend to make ω cv,K smaller resulting in less weight onȲ S . In other words, 10-fold cross-validation makes more aggressive use of the large sample than does leave-one-out.
Remark 1. The cross-validation estimates do not make use ofσ 2 B because the large sample is held fixed. They are in this sense conditional on the large sample. Our oracle takes account of the randomness in set B, so it is not conditional. One can define a conditional oracle without difficulty, but we omit the details. Neither the bootstrap nor the plug-in methods are conditional, as they approximate our oracle. Taking ω bapi as a representor of unconditional methods and ω cv,n as a representor of conditional ones, we see that the latter has a larger denominator while they both have the same numerator, at least whenδ 2 0 >σ 2 S /n. This suggests that conditional methods are more aggressive and we will see this in the simulation results.
With an L 1 penalty on δ we find from Theorem 3.1 that That is, the estimator movesȲ S towardsȲ B by an amount λ/n except that it will not move past the pooled averageȲ P . The optimal choice of λ is not available in closed form.

An L 1 oracle
The L 2 oracle depends only on moments of the data. The L 1 case proves to be more complicated, depending also on quantiles of the error distribution. To investigate L 1 penalization, we suppose that the errors are Gaussian. Then we can compute E((μ(λ) − µ) 2 ) for the L 1 penalization by a lengthy expression broken into several steps below. That expression is not simple to interpret. But we can use it to numerically find the best value of λ for an oracle using the L 1 penalty. That then allows us to compare the L 1 and L 2 oracles in Section 4.1. Let D =Ȳ B −Ȳ S and then define Lemma 9.1 in the Appendix gives these identities: Section 9.5 of the Appendix shows that Substituting the quantities above into (3.7) yields a computable expression for the loss in the L 1 penalized case.

Numerical examples
We have simulated some special cases of the data enrichment problem. First we simulate the pure location problem which has d = 1. Then we consider the regression problem with varying d.

Location
We simulated Gaussian data for the location problem. The large sample had N = 1000 observations and the small sample had n = 100 observations: Our data had µ = 0 and σ 2 S = σ 2 B = 1. We define the relative bias as We investigated a range of relative bias values. It is only a small simplification to take σ 2 S = σ 2 B . Doubling σ 2 B has a very similar effect to halving N . Equal variances might have given a slight relative advantage to a hypothesis testing method as described below.
The accuracy of our estimates is judged by the relative mean squared error E((μ−µ) 2 )/(σ 2 S /n). Simply takingμ =Ȳ S attains a relative mean squared error of 1. Figure 1 plots relative mean squared error versus relative bias for a collection of estimators, with the results averaged over 10,000 simulated data sets. We used the small sample only method as a control variate.
The solid curve in Figure 1 shows the oracle's value. It lies strictly below the horizontal S-only line. None of the competing curves lie strictly below that line. None can becauseȲ S is an admissible estimator for d = 1 (Stein, 1956). The second lowest curve in Figure 1 is for the oracle using the L 1 version of the penalty. The L 1 penalized oracle is not as effective as the L 2 oracle and it is also more difficult to approximate. The highest observed predictive MSEs come from a method of simply pooling the two samples. That method is very successful when the relative bias is near zero but has an MSE that becomes unbounded as the relative bias increases. Now we discuss methods that use the data to decide whether to use the small sample only, pool the samples or choose an amount of shrinkage. We may list them in order of their worst case performance. From top (worst) to bottom (best) in Figure 1 they are: hypothesis testing, 5-fold cross-validation, 10-fold cross-validation, AICc, leave-one-out cross-validation, and then the simple plugin method which is minimax among this set of choices. AICc and leave-one-out are very close. Our cross-validation estimators used ω = max(ω cv,K , n/(n + N )) where ω cv,K is given by (3.3).
The hypothesis testing method is based on a two-sample t-test of whether δ = 0. If the test is rejected at α = 0.05, then only the small sample data is used. If the test is not rejected, then the two samples are pooled. That test was based on σ 2 B = σ 2 S which may give hypothesis testing a slight advantage in this setting (but it still performed poorly).
The AICc method performs virtually identically to leave-one-out cross-validation over the whole range of relative biases.
None of these methods makes any other one inadmissible: each pair of curves crosses. The methods that do best at large relative biases tend to do worst at relative bias near 0 and vice versa. The exception is hypothesis testing. Compared to the others it does not benefit fully from low relative bias but it recovers the quickest as the bias increases. Of these methods hypothesis testing is best at the highest relative bias, K-fold cross-validation with small K is best at the lowest relative bias, and the plug-in method is best in between.
Aggressive methods will do better at low bias but worse at high bias. What we see in this simulation is that K-fold cross-validation is the most aggressive followed by leave-one-out and AICc and that plug-in is least aggressive. These findings confirm what we saw in the formulas from Section 3. Hypothesis testing does not quite fit into this spectrum: its worst case performance is much worse than the most aggressive methods yet it fails to fully benefit from pooling when the bias is smallest. Unlike aggressive methods it does very well at high bias.

Regression
We simulated our data enrichment method for the following scenario. The small sample had n = 1000 observations and the large sample had N = 10,000. The true β was taken to be 0. This is no loss of generality because we are not shrinking β towards 0. The value of γ was taken uniformly on the unit sphere in d dimensions and then multiplied by a scale factor that we varied.
We considered d = 2, 4, 5 and 10. All of our examples included an intercept column of 1s in both X S and X B . The other d − 1 predictors were sampled from a Gaussian distribution with covariance C S or C B , respectively.
In one simulation we took C S and C B to be independent Wishart(I, d − 1, d − 1) random matrices. In the other simulation, they were sampled as a spiked covariance model (Johnstone, 2001). There C S = I d−1 + ρuu T and C B = I d−1 + ρvv T where u and v are independently and uniformly sampled from the unit sphere in R d−1 and ρ ≥ 0 is a parameter that measures the lack of proportionality between covariances. We chose ρ = d so that the sample specific portion of the variance has comparable magnitude to the common part.
The variance in the small sample was σ 2 S = 1. To model the lower quality of the large sample we used σ 2 B = 2. We scaled the results so that regression using sample S only yields a mean squared error of 1. We computed the risk of an L 2 oracle, as well as sampling errors when λ is estimated by the plug-in formula, by our bias-adjusted plug-in formula and via AICc. In addition we considered the simple weighted combination ωβ S + (1 − ω)β B with ω chosen by the plug-in formula. To optimize (2.10) over λ we used the optimize function in R which is based on golden section search (Brent, 1973). We also included a shrinkage estimator. Because our simulated runs all had β S = 0 it is not reasonable to include shrinkage ofβ S towards zero in the comparison; we cannot in practice shrink towards the truth. Instead, we used the positive part Stein shrinkage estimate (2.13) shrinkingβ S towardsβ B but not past it. That shrinkage requires an estimateσ 2 S of σ 2 S . We used the true value, σ 2 S = 1, giving the shrinkage estimator a slight advantage. We did not include hypothesis testing in this example, because there are 2 d possible ways to decide which parameters to pool and which to estimate separately.
We simulated each of the two covariance models with each of the four dimensions 10,000 times. For each method we averaged the squared prediction errors (β − β) T V S (β − β) and then divided those mean squared errors by the one for using the small sample only. Figures 2 and 3 show the results.
At small bias levels pooling the samples is almost as good as the oracle. But the loss for pooling samples grows without bound when the bias increases. For d = 2, shrinkage amounts to using the small sample only but for d > 2 it performs universally better than the small sample. When comparing methods we see that the curves usually cross. Methods that are best at low bias tend not to be best at high bias. Note however that there is a lot to gain at low bias while the methods differ only a little at high bias. As a result the more aggressive methods making greater use of the large data are more likely to yield a big improvement.
The weighting estimator generally performs better than the shrinkage estimator in that it offers a meaningful improvement at low bias costing a minor relative loss at high bias. We analyze that estimator in Section 5. The plug-in estimator also is generally better than shrinkage. The AICc estimator is generally better than both of those. We do not show the bias adjusted plug-in estimators. Their performance is almost identical to AICc (always within about 0.015). Of those, the one using Θ bapi was consistently at least as good as Θ bapi and sometimes a little better.
The greatest gains are at or near zero bias. Table 1 shows the quadratic losses at δ = 0 normalized by the loss attained by the oracle. Pooling is almost Table 1 This table shows the quadratic loss (2.14) normalized by that of the oracle when δ = 0 (the no bias condition). The methods are described in the text. as good as the oracle in this case but we rule it out because of its extreme bad performance when the bias is large. Some of our new estimators yield very much reduced squared error compared to the shrinkage estimator. For example three of the new methods' squared errors are just less than half that of the shrinkage estimator for d = 10 and the Wishart covariances.

Inadmissibility
Section 4 gives empirical support for our proposal. Several of the estimators perform better than ordinary shrinkage. In this section we provide some theoretical support. We provide a data enriched estimator that makes least squares on the small sample inadmissible. The estimator is derived for the proportional design case but inadmissibility holds even when V B = X T B X B is not proportional to V S = X T S X S . The inadmissibility is with respect to a loss function E( X T (β − β) 2 ) where V T = X T T X T is proportional to V S . To motivate the estimator, suppose for the moment that V B = N Σ, V S = nΣ and V T = Σ for a positive definite matrix Σ. Then the weighting matrix W λ in Lemma 2.1 simplifies to W λ = ωI where ω = (N + nλ)/(N + nλ + N λ). As a resultβ = ωβ S + (1 − ω)β B and we can find and estimate an oracle's value for ω.
We show below that the resulting estimator ofβ with estimated ω dominateŝ β S (making it inadmissible) under mild conditions that do not require V B ∝ V S . We do need the model degrees of freedom to be at least 5, and it will suffice to have the error degrees of freedom in the small sample regression be at least 10. The result also requires a Gaussian assumption in order to use a lemma of Stein's. Write . This error is minimized by the oracle's parameter value .
When V S = nΣ and V B = N Σ, we find The plug-in estimator iŝ To allow a later bias adjustment, we generalize this plug-in estimator. Let h(σ 2 B ) be any nonnegative measurable function ofσ 2 Here are the conditions under whichβ S is made inadmissible by the data enrichment estimator.
Theorem 5.1. Let X S ∈ R n×d and X B ∈ R N ×d be fixed matrices with X T S X S = nΣ and X T holds for any matrix X T with X T T X T = Σ and anyω =ω plug,h given by (5.2). Proof. Section 9.7 in the Appendix.
The condition on m can be relaxed at the expense of a more complicated statement. From the details in the proof, it suffices to have d ≥ 5 and m(1 − 4/d) ≥ 2.
Because E(γ T Σγ) > γ T Σγ we find thatω plug is biased upwards, making it conservative. In the proportional design case we find that the bias is dσ 2 S /n + dσ 2 B /N . That motivates a bias adjustment, replacingγ T Σγ byγ T Σγ − dσ 2 where values below n/(n + N ) get rounded up. This bias-adjusted estimate of ω is not covered by Theorem 5.1. Subtracting onlyσ 2 which corresponds to taking h(σ 2 B ) ≡ 0 in equation (5.2). Data enrichment witĥ ω given by (5.4) makesβ S inadmissible whether or not the motivating covariance proportionality holds.

A matrix oracle
In this section we look for an explanation of how data enrichment might be more accurate than Stein shrinkage. We generalize our estimator tô and then find the optimal matrix W .
It is interesting that when we are free to choose the entire d × d matrix W , then the optimal choice is the same for all weighting matrices V T .
The penalized least squares criterion (2.1) leads to a matrix weighting of the two within-sample estimators. The weighting matrix W λ is in a one dimensional family indexed by 0 ≤ λ ≤ ∞. The optimal W from (6.1) is not generally in that family.
Both W λ from criterion (2.1) and W from (6.1) trade off bias and variance, through the appearance γγ T , σ 2 S , and σ 2 B , which for (2.1) appear in the formula for the optimal λ. The advantage of working with W λ instead of W is that W λ yields a one parameter family of candidate weighting matrices to search over.
When V S and V B are both proportional to the same positive definite matrix V T , then the data enrichment oracle chooses W = ωI d where which mimicks the form of the optimal W in equation (6.1), replacing numerator and denominator by traces after multiplying both by V T .
The James-Stein shrinker chooses W = ω JS I d where If we approximateγ T V Sγ by its expectation tr((γγ T + σ 2 .
The presence of 2σ 2 S in the numerator leads the James-Stein approach to make less aggressive use of the big data set than data enrichment does. We believe that this is why the James-Stein method did not perform well in our simulations.

Related literatures
There are many disjoint literatures that study problems like the one we have presented. They do not seem to have been compared before, the literatures seem to be mostly unaware of each other, and there is a surprisingly large variety of problem contexts. Some quite similar sounding problems turn out to differ on critically important details. We give a brief summary of those topics here.
The key ingredient in our problem is that we care more about the small sample than the large one. Were that not the case, we could simply pool all the data and fit a model with indicator variables picking out one or indeed many different special subsets of interest. Without some kind of regularization, that approach ends up being similar to taking λ = 0 and hence does not borrow strength.
The closest match to our problem setting comes from small area estimation in survey sampling. The monograph by Rao (2003) is a comprehensive treatment of that work and Ghosh and Rao (1994) provide a compact summary. In that context the large sample may be census data from the entire country and the small sample (called the small area) may be a single county or a demographically defined subset. Every county or demographic group may be taken to be the small sample in its turn. The composite estimator (Rao, 2003, Chapter 4.3) is a weighted sum of estimators from small and large samples. The estimates being combined may be more complicated than regressions, involving for example ratio estimates. The emphasis is usually on scalar quantities such as small area means or totals, instead of the regression coefficients we consider. One particularly useful model (Ghosh and Rao, 1994, equation (4.2)) allows the small areas to share regression coefficients apart from an area specific intercept. Then BLUP estimation methods lead to shrinkage estimators similar to ours.
Our methods and results are similar to empirical Bayes methods, drawing heavily on ideas of Charles Stein. A Stein-like result also holds for multiple regression in the context of just one sample. We mentioned already the regression shrinkers of Copas (1983) and Stein (1960). Efron and Morris (1973a) find that the Stein effect for shrinking to a common mean takes place at dimension 4 and George (1986) finds that the effect takes place at dimension 3+q when shrinking means towards a q-dimensional linear manifold.
A similar problem to ours is addressed by Chen and Chen (2000). Like us, they have (X, Y ) pairs of both high and low quality. In their setting both high and low quality pairs are defined for the same set of individuals. Their given sample has all of the low quality data and the high quality data are available only on a simple random sample of the subjects. Boonstra et al. (2013a) consider a genomics problem where there are both low and high quality versions of X, from two different technical platforms, but all data share the same Y . All observations have the low quality X's while a subset have both high and low quality X measurements. They take a Bayesian approach. Boonstra et al. (2013b) handle the same problem via shrinkage estimates. A crucial difference in our setting, is that the subjects are completely different in our two samples; no (X, Y ) pair in one data set comes from the same person as an (X, Y ) pair in the other data set. Mukherjee and Chatterjee (2008) use shrinkage methods to blend two estimators. One is a case-control estimate of a log odds ratio. The other is a case-only estimator, derived under an assumption of gene-environment independence. They also derive and employ a plug-in estimator. Their target parameter is scalar so no Stein effect could be expected. Chen et al. (2009) address the same issue via L 1 and L 2 shrinkage based methods, and give some asymptotic covariances.
In chemometrics, a calibration transfer problem (Feudale et al., 2002) comes up when one wants to adjust a model to new spectral hardware. There may be a regression model linking near-infrared spectroscopy data to a property of some sample material. The transfer problem comes up for data from a new machine. Sometimes one can simply run a selection of samples through both machines but in other cases that is not possible, perhaps because one machine is remote (Woody et al., 2004). Their primary and secondary instruments correspond to our small and big samples respectively. Their emphasis is on transfering either principal components regression or partial least squares models, not the plain regressions we consider here.
A common problem in marketing is data fusion, also known as statistical matching. Variables (X, Y ) are measured in one sample while variables (X, Z) are measured in another. There may or may not be a third sample with some measured triples (X, Y, Z). The goal in data fusion is to use all of the data to form a large synthetic data set of (X, Y, Z) values, perhaps by imputing missing Z for the (X, Y ) sample and/or missing Y for the (X, Z) sample. When there is no (X, Y, Z) sample, some untestable assumptions must be made about the joint distribution, because it cannot be recovered from its bivariate margins. The text by D' Orazio et al. (2006) gives a comprehensive summary of what can and cannot be done. Many of the approaches are based on methods for handling missing data (Little and Rubin, 2009).
Medicine and epidemiology among other fields use meta-analysis (Borenstein et al., 2009). In that setting there are (X, Y ) data sets from numerous environ-ments, no one of which is necessarily of primary importance.
Our problem is an instance of what machine learning researchers call domain adaptation. They may have fit a model to a large data set (the 'source') and then wish to adapt that model to a smaller specialized data set (the 'target'). This is especially common in natural language processing. NIPS 2011 included a special session on domain adaptation. In their motivating problems there are typically a very large number of features (e.g., one per unique word appearing in a set of documents). They also pay special attention to problems where many of the data points do not have a measured response. Quite often a computer can gather high dimensional X while a human rater is necessary to produce Y . Daumé (2009) surveys various wrapper strategies, such as fitting a model to weighted combinations of the data sets, deriving features from the reference data set to use in the target one and so on. Cortes and Mohri (2011) consider domain adaptation for kernel-based regularization algorithms, including kernel ridge regression, support vector machines (SVMs), or support vector regression (SVR). They prove pointwise loss guarantees depending on the discrepancy distance between the empirical source and target distributions, and demonstrate the power of the approach on a number of experiments using kernel ridge regression. We have given conditions under which adaptation is always beneficial.
A related term in machine learning is concept drift (Widmer and Kubat, 1996). There a prediction method may become out of date as time goes on. The term drift suggests that slow continual changes are anticipated, but they also consider that there may be hidden contexts (latent variables in statistical teminology) affecting some of the data.

Conclusions
We have studied a middle ground between pooling a large data set into a smaller target one and ignoring it completely. Looking at the left side of Figures 1, 2 and 3 we see that in the low bias cases the more aggressive methods have a clear advantage. Fortune favors the bold. Pooling is the boldest and wins the most when bias is small. But pooling has unbounded risk as bias increases. That is, misfortune also favors the bold. Our shrinkage methods provide a compromise. In higher dimensional settings of Figures 2 and 3, we see that AICc and bias adjusted plug-in gain a lot of efficiency when the bias is low. When the bias is high, they are squeezed into a narrow band between the oracle performance and that ofβ S which ignores the big data set. As a result, the new methods show large improvements compared to shrinkage when the bias small but only lose a little when the bias is large.

Appendix I: proofs
This appendix presents proofs of the results in this article. They are grouped into sections by topic, with some technical supporting lemmas separated into their own sections.

Proof of Theorem 2.1
Proof. First Next with X T = X S , and M = V We place V 1/2 S V −1/2 S between these factors and absorb them left and right. Then we reverse the order of the factors and repeat the process, yielding tr(W λ ) = tr(I + λM + λI) −1 (I + λM ).
Writing M = U diag(ν 1 , . . . , ν d )U T for an orthogonal matrix U and simplifying yields the result.

Derivation of equation (3.3)
We suppose for simplicity that n = rK for an integer r, so the K folds have equal size. In that caseȲ S,−k = (nȲ S − rȲ S,k )/(n − r). Now After some algebra, the numerator of (9.1) is and the denominator is The only quantity in ω cv which depends on the specific K-way partition used isσ 2 S,K . If the groupings are chosen by sampling without replacement, then under this sampling, using the finite population correction for simple random sampling, where s 2 S = σ 2 S n/(n − 1). This simplifies to Replacingσ 2 S,K in ω cv by its expectation yields (3.3).

Derivation of equation (3.7)
Let f = n/(n + N ) be the fraction of the pooled data coming from the small sample and F = 1−f be the fraction from the large sample. Define D =Ȳ B −Ȳ S and c = λ/(nF ) = λ(1/n + 1/N ). Then We replaceȲ B andȲ S by linear combinations of D and another variable H chosen to be statistically independent of D.
In terms of these independent variables we havê where c 0 = α/(α + 1) − f , and c 1 = 1/(α + 1). Without loss of generality, µ = 0. Then D ∼ N (δ, σ 2 S /n + σ 2 B /N ) independently of H ∼ N (δ, α 2 σ 2 S /n + σ 2 B /N ). After some algebra, In addition to first and second moments of D and H, we need some expectations of functions of D involving the sign function and some indicators. They are given by Lemma 9.1 below. Let ϕ and Φ be the probabilty density and cumulative distribution functions respectively of the N (0, 1) distribution.

Supporting lemmas for inadmissibility
In this section we first recall Stein's Lemma. Then we prove two technical lemmas used in the proof of Theorem 5.1.
Lemma 9.3. Let η ∼ N (0, I d ), b ∈ R d , and let A > 0 and B > 0 be constants. Let Then Proof. First, By Stein's lemma (Lemma 9.2), we have and thus after collecting terms.
Note thatβ S = β + (X T S X S ) −1 X T S ε S andβ B = β + γ + (X T B X B ) −1 X T B ε B . It is convenient to define η S = Σ 1/2 (X T S X S ) −1 X T S ε S and η B = Σ 1/2 (X T B X B ) −1 X T B ε B .
The quantities A and B are, after conditioning, the constants that appear in technical Lemma 9.3. Similarly C and D introduced below match the constants used in Lemma 9.4.
With these substitutions and some algebra, We now apply the two technical lemmas from Section 9.6. Since η * S is independent of (b, A, B) and Q ∼ χ 2 (m) , by Lemma 9.3, we have where C = 2 − 4/d and D = (m/d)( b − η * S 2 + dnN −1σ2 B /σ 2 S ). Now suppose that d ≥ 5. Then C ≥ 2 − 4/5 > 1 and so conditionally on η S , η B , andσ 2 B , the requirements of Lemma 9.4 are satisfied by C, D and Q. Therefore For the general plug-inω plug,h we replace dσ 2 B /N above by h(σ 2 B ). This quantity depends onσ 2 B and is independent ofσ 2 S , η B and η S . It appears within B where we need it to be non-negative in order to apply Lemma 9.3. It also appears within D which becomes (m/d)( b − η * S 2 + nh(σ 2 B )/σ 2 S ). Even when we take var(h(σ 2 B )) = 0 we still get var(D) > 0 and so the first inequality in (9.11) is still strict. Now suppose that V B = N Σ. The distributions of η S ,σ 2 S andσ 2 B remain unchanged but now η B ∼ N 0, Σ 1/2 V −1 B Σ 1/2 σ 2 B independently of the others. The changed distribution of η B does not affect the application of Lemma 9.3 because that lemma is invoked conditionally on η B . Similarly, Lemma 9.4 is applied conditionally on η B . The changed distribution of η B changes the distribution of D but we can still apply (9.11).
The expectation at (9.9) is negative when d = 4, as can be verified by a one dimensional quadrature. For this reason, the inadmissibilty result requires d > 4. 9.8. Proof of Theorem 6.1 Recall thatβ(W ) = Wβ S + (I − W )β B . The first two moments ofβ(W ) are E(β(W )) = β + (I − W )γ, and var(β(W )) = σ 2 The loss is L(W ) = E((β(W ) − β) T V T (β(W ) − β)) and We will use two rules from Brookes (2011) for matrix differentials. If A, B and C don't depend on the matrix X then the differential of tr(XA) and tr(AX) are both AdX, and, the differential of tr(AXBX T C) is BX T CA + B T X T A T C T times dX.
The differential of L(W ) when W changes is times dW . Let W * be the hypothesized optimal matrix given at (6.1). It is symmetric, as are V S , V B and V T . We may therefore write the differential at that matrix as where G = γγ T . This differential vanishes, showing that W * satisfyies first order conditions. The differential of this differential is 2(G+σ 2 S V S +σ 2 B V B W )V T which is positive definite, and so W * must be a minimum.