REGULARIZED REGRESSION

Ordinary Least Squares (OLS) and regularized regression methods are both approaches for estimating the model y = Xβ + ϵ. Frequentist regularized methods shrink coefficient estimates by adding a penalty term to the OLS cost function, but differ in how they do so. We examine three common frequentist regularization methods: ridge regression, LASSO, and elastic net. We also contextualize these methods against the Bayesian approach of specifying explicit priors.

Ridge regression minimizes the cost function:

⁽¹⁾ $$ \begin{array}{rl}{\displaystyle \sum_{i=1}^M}\left(\right.{y}_i\hskip0.3em -\hskip0.3em {\beta}_0\hskip0.3em -\hskip0.3em {\displaystyle \sum_{j=1}^p}{\beta}_j{x}_{ij}{\left)\right.}^2+\lambda {\displaystyle \sum_{j=1}^p}{\beta}_j^2.& \end{array} $$

The ridge cost function is the OLS cost function, which minimizes the sum of squared residuals, plus an additional penalty term $ \lambda \sum_{j=1}^p{\beta}_j^2. $ Regularization penalties shrink coefficients because any reduction in the sum of squared residuals produced by larger coefficients is partially offset by an increase in the size of the penalty term. For example, while a model with two coefficients 2 and 1 might result in a lower sum of squares compared to a model with coefficients 1.9 and 0.9, the latter has a smaller penalty term, and with a sufficiently large λ will have a smaller overall sum. Importantly, the size of the ridge penalty is proportional to the square of the coefficients, so if coefficient β ₁ is twice the magnitude of coefficient β ₂, β ₁ contributes four times as large a penalty as β ₂. Consequently, for a marginal increase in the sum of squared residuals, ridge regression prefers two coefficients of similar sizes to one large and one small coefficient.

LASSO minimizes the cost function:

(2) $$ \begin{array}{rl}{\displaystyle \sum_{i=1}^M}\left(\right.{y}_i\hskip0.3em -\hskip0.3em {\beta}_0\hskip0.3em -\hskip0.3em {\displaystyle \sum_{j=1}^p}{\beta}_j{x}_{ij}{\left)\right.}^2+\lambda {\displaystyle \sum_{j=1}^p}|{\beta}_j|.& \end{array} $$

The difference between Equation 2 and Equation 1 being that the LASSO loss function regularizes coefficients proportional to their absolute size rather than their square, and so does not produce the effect of ridge regression of preferring many small coefficients to a smaller set of larger coefficients.

The elastic net incorporates both of these penalty terms, adding an additional term α, which controls the relative weight between the two. When α = 1, the elastic net is equivalent to LASSO, and when α = 0, it is equivalent to ridge regression:

(3) $$ \begin{array}{rl}{\displaystyle \sum_{i=1}^M}\left(\right.{y}_i\hskip0.3em -\hskip0.3em {\beta}_0\hskip0.3em -\hskip0.3em {\displaystyle \sum_{j=1}^p}{\beta}_j{x}_{ij}{\left)\right.}^2+\lambda {\displaystyle \sum_{j=1}^p}(\alpha |{\beta}_j|+(1-\alpha ){\beta}_j^2).& \end{array} $$

These penalty terms encode assumptions about the likely size of the β _j effects, which can be interpreted as Bayesian priors (Hastie, Tibshirani, and Friedman Reference Hastie, Tibshirani and Friedman2009). In Bayesian terms, the ridge regression penalty expects that β _j are normally distributed around zero, which implies that coefficients are increasingly unlikely the further away from zero they are. As a Bayesian prior the LASSO penalty follows a Laplace distribution, which is sharply peaked at zero, entailing the expectation that large values of β _j are unlikely, but furthermore that many effects are exactly zero. The elastic net prior is a mixture of these two distributions. Larger λs imply stronger versions of these assumptions, or in Bayesian terms, prior distributions with smaller standard deviations and/or scale parameters.

Regularized Regressions Can Reopen Backdoor Causal Paths

In order to demonstrate the counter-intuitive effects that regularization can have, we run simulations in which X causes Y, and X and Z are correlated, but Z has no causal effect on Y, illustrated in the directed acyclic graph (DAG) (4) below.

(4)

Where X and Z are multivariate standard normal variables with a correlation ρ, and Y is a linear function of β _X X and a standard normal error term (in CF2’s data, if the response rate effect on polarization was actually zero, then X would be equivalent to time, Z to response rate, and Y to polarization):

(5) $$ \begin{array}{rl}\begin{array}{rll}\left[\begin{array}{l}\qquad X\qquad \\ {}\qquad Z\qquad \end{array}\right]& \sim \mathcal{N}\left(\right.\left[\begin{array}{l}\qquad 0\qquad \\ {}\qquad 0\qquad \end{array}\right],\left[\begin{array}{ll}\qquad 1\qquad & \qquad \rho \qquad \\ {}\qquad \rho \qquad & \qquad 1\qquad \end{array}\right]\left)\right.,& \\ {}& \\ {}Y& \sim {\beta}_XX+\mathcal{N}(0,1).\\ {}\end{array}& \end{array} $$

From this DGP, we simulate ten thousand datasets with β _X ∼ U(0, 1) and ρ = −0.9 (approximately the correlation between time and response rate in CF2’s data).

We then fit the following model to the simulated data:

(6) $$ \begin{array}{rl}{y}_i={\beta}_0+{\beta}_X{X}_i+{\beta}_Z{Z}_i+{\epsilon}_i& \end{array} $$

using four estimators: OLS, ridge regression, LASSO, and elastic net (with α = 0.5, an equal mix of the two penalties). For each regularized estimator, we set four λ penalties: 0.00615, 0.154, 0.615, and 2.46, which correspond approximately to a weakly informative prior for the average simulated case, and priors that are 20%, 10%, and 5% the width of the standard deviation of that weakly informative prior for ridge regression.

Figure 1 shows the average coefficient estimates for β _X and β _Z for each estimator, which illustrates the interaction between effect size and λ. We first consider bias (the average signed difference between the estimate and simulated coefficient). The OLS estimate is unbiased for both β _X and β _Z across the range of the simulated β _X terms. With a penalty equivalent to a weakly informative prior (λ = 0.00615, row A), the regularized estimates are only marginally biased. With a larger penalty (λ = 0.154, row B), the bias in the estimated β _x is substantially larger. For ridge regression, this bias increases linearly as the simulated β _X increases. For LASSO (and to a lesser extent elastic net), the level of bias flattens off with sufficiently large β _X values. As λ gets larger, the bias in the β _X estimates also increases (rows C and D). These biases are strongest for LASSO—which with the highest penalty almost always shrinks the estimates to zero—and smallest for ridge regression, with elastic net in between.

Figure 1. Estimated OLS, Ridge, LASSO, and Elastic Net Coefficients for Simulated Values of $ {\beta}_X $ and $ \lambda $ (Local Linear Smooths)

Note: $ {\beta}_Z $ always simulated as zero.

That regularized regression induces bias and shrinks β _X toward zero is well understood—this is the price these estimators pay for reducing variance. What is less appreciated is that regularized estimators also bias β _Z away from zero. This occurs because some of the X → Y effect is left unaccounted for by the biased β _X estimate, and so is transferred to Z via the Z ↔ X → Y backdoor path (for a demonstration, see Section S3 of the Supplementary Material). The β _Z bias is most pronounced for ridge regression: as the λ penalty increases, the bias in β _Z initially increases (compare the ridge regression lines in rows B/C). As λ increases further, the ridge β _Z estimates are themselves regularized, and move back toward zero. Although this means the absolute bias in β _Z is reduced with the largest λ terms, the effect of β _Z relative to β _X continues to increase, with estimators with large penalties assigning almost equal values to β _X and β _Z . For LASSO, the bias in β _Z is only noticeable at the smaller levels of λ (i.e., row B); with higher penalties β _Z is almost always pulled to zero. Elastic net is again between the ridge and LASSO estimates.

Making Statistical Inferences with Regularized Regressions

To make inferences about underlying population parameters (e.g., the true effect of response rate on estimates of mass polarization), we typically use estimates of sampling variance. With unbiased estimators this is straightforward; in repeated samples, we would expect 95% of our confidence intervals (CIs) to overlap the population value. With regularized estimators, this is more complicated because (1) the sampling variance of some regularized estimators is poorly defined, and (2) even where it is well defined (such as for ridge regression) the resulting CIs may have no overlap with the underlying population parameter.

A particular challenge for the LASSO and elastic net estimators is that “we still do not have a general, statistically valid method of obtaining standard errors of LASSO estimates” (Kyung et al. Reference Kyung, Gill, Ghosh and Casella2010, 377) even using bootstrap methods. LASSO and elastic net are therefore poorly suited to inference.

This problem does not apply to ridge regression, where the variance of ridge estimates is

(7) $$ \begin{array}{rl}Var({\widehat{\beta}}_{ridge}|X)={\sigma}^2{({X}^TX+\lambda I)}^{-1}{X}^TX{({X}^TX+\lambda I)}^{-1}.& \end{array} $$

We must be clear what this variance captures: how much a ridge estimate $ {\widehat{\beta}}_{ridge} $ is expected to differ across repeated samples (when λ = 0, this is the same as OLS). By design, the λ penalty decreases variance. However, λ also increases bias (indeed, it reduces variance by increasing bias). This has the pernicious effect of shrinking the ridge CIs at the same time as increasing the distance between the ridge estimates and the underlying parameter of interest. As bias grows, CI coverage falls because more of the discrepancy between $ {\widehat{\beta}}_{ridge} $ and β is due to bias (which the CIs do not account for) rather than sampling variation (which the CIs do account for). If the bias is large enough, CIs routinely exclude the true value of β. Narrow ridge regression CIs do not necessarily indicate a precise estimate of the parameter of interest, but simply that ridge regression would produce similar estimates with new samples. The resulting CIs may lead to misleading inferences about population parameters.

We can again demonstrate this problem using simulations. Returning to the simulated data we used in Figure 1, we calculate CIs for OLS and ridge regression. For each simulation, we record whether the CI overlaps the simulated parameters for β _X and β _Z , illustrated in Figure 2. The OLS CIs have the appropriate coverage (95%) across all simulated values of β _X . For the ridge CIs, however, we again see an interaction between λ and the size of β _X . When λ is small and β _X is small, the CIs still have coverage close to the 95% benchmark. As β _X gets larger the coverage rate begins to fall for both β _X and β _Z . This problem is exacerbated by larger values of λ—the coverage rate falls precipitously and approaches zero. At the highest levels of λ and β _X , this means that the CIs for β _X and β _Z never overlap the population values for these parameters. In all but the most favorable of circumstances, ridge regression is simply not an appropriate tool for statistical inference.

Figure 2. Estimated OLS and Ridge CI Coverage for Simulated Values of $ {\beta}_X $ and $ \lambda $ (Local Logistic Smooths)

Note: $ {\beta}_Z $ always simulated as zero.

AN EXAMPLE OF RIDGE REGRESSION ESTIMATES IMPLYING UNREALISTIC PRIORS ABOUT THE EFFECT OF TIME

In light of the issues we discuss above, we turn to CF2’s analysis. CF2 analyze the relationship between survey response/cooperation rate and partisan polarization for survey questions asked by Pew between 2004 and 2018. They operationalize polarization as the absolute Cohen’s D of mean differences between Republican and Democrat identifiers’ responses. They estimate models for six issue areas: economy, energy, immigration, civil rights, welfare, and foreign policy, specified as

(8) $$ \begin{array}{rl}{y}_{it}={\beta}_0+{\beta}_1r{r}_{it}+{\beta}_2congres{s}_t+{\beta}_3yea{r}_t+{\epsilon}_{it,}& \end{array} $$

where rr _it is survey i’s response rate at time t, congress _t is congressional polarization at time t, and year _t is the year in which the survey was conducted. For the contact and cooperation rate versions of these models, the β ₁rr _it term is replaced with β ₄contact _it and β ₅cooperation _it . Equation 8, which adds a year term to an earlier model estimated by Cavari and Freedman (Reference Cavari and Freedman2018, henceforth CF1) was suggested by Mellon and Prosser (Reference Mellon and Prosser2021, hereafter MP) in a critique of CF1 as one of several ways they could deal with the problem of time trends.Footnote ¹ CF2 estimate their models using ridge regression (unlike MP, who use OLS) with λ chosen by the KM4 method. CF2 reports statistically significant negative effects of response and cooperation rates (i.e., lower response rates increase the estimated level of polarization) for three areas: economy, energy, and immigration.

CF2’s KM4 ridge regression estimates imply strong priors that are between 1.5% and 5.9% the width of a weakly informative prior. Are these priors plausible when considering time trends? A priori, there are good reasons to think not. Many social processes can be modeled as random walks or similar processes. Independent random walks will tend to correlate with time (Granger and Newbold Reference Granger and Newbold1974) because their variance increases as a function of time (see Section S2 of the Supplementary Material).

Moreover, we can show that many social science variables are correlated with time in practice. We calculate the annual mean of binary, ordinal and continuous variables asked in 15 or more years of the GSS (Smith et al. Reference Smith, Davern, Freese and Morgan2019) and estimate a linear model of each variable on time using OLS. Figure 3 shows a histogram of the standardized coefficients for the effect of time on each variable. The distribution has a far higher dispersion than CF2’s priors imply, with a standard deviation of 0.657 and heavy tails close to 1 and −1.

Figure 3. Two Hundred Forty Nine Standardized Coefficients of Time on GSS Variable Annual Averages

To show how different these priors are from CF2’s, we can consider how likely we are to observe a 0.5 or larger magnitude unit-scale coefficient. The weakest of CF2’s priors implies the probability of observing a 0.5 or larger unit-scale coefficient (considerably smaller than the mean OLS-year estimate on their data) is 2.8 × 10⁻¹⁷. By contrast, a N(0, 0.657) prior gives a probability of 0.447 of observing an effect size magnitude of 0.5 or larger.

CF2’s narrow priors mean they strongly shrink the year coefficient in their models, which risks reopening the backdoor causal paths that the inclusion of a time control is designed to block.

AN EXAMPLE OF HOW A LARGE SHRINKAGE PARAMETER AFFECTS EMPIRICAL RESULTS

To demonstrate the impact of CF2’s regularization approach, we reestimate their models using a variety of approaches: their original KM4 ridge regressions, OLS, ridge regression with λ chosen by CV, a weakly informative N(0, 1) unit-scaled prior, a prior informed by the GSS distribution of time effect sizes N(0, 0.657), LASSO with λ chosen by CV, and elastic net (α = 0.5) with λ again chosen by CV. These estimators show the impact of different regularization approaches, but they all share a common approach to controlling for time (a linear trend), which may not fully account for time-series problems. Given this, we benchmark them against a more sophisticated repeated cross-section model suggested by Lebo and Weber (Reference Lebo and Weber2015) (ARIMA-MLM—we elaborate on this approach in Section S1 of the Supplementary Material). These estimates for the key response and cooperation rate variables are shown in Figure 4.Footnote ²

Figure 4. Estimated Effects of Response and Cooperation Rates from Replication of CF2’s Analysis Using Six Estimators: OLS, Ridge Regression ( $ \lambda $ Chosen Using KM4), Ridge Regression ( $ \lambda $ Chosen by CV), LASSO ( $ \lambda $ Chosen by CV), Elastic Net ( $ \lambda $ Chosen by CV), Bayesian Regression (Weakly Informative Priors), Bayesian Regression (Informed Priors Based on the GSS Distribution of Time Trends), and ARIMA-MLM

Note: CIs omitted from LASSO and elastic net estimates for the reasons discussed in Section “Making Statistical Inferences with Regularized Regressions.” CIs for the ridge (KM4) estimates are drawn, but are smaller than the size of point estimate marker. See Appendix R6 for full tables.

The results in Figure 4 present a striking pattern—six of the seven estimators of the linear-time-trend model produce similar results, and one does not. With the exception of the economy issue area, CF2’s KM4 results are considerably closer to zero than either the OLS estimate or the other regularized estimates. All of the unregularized estimates, and nearly all of the other regularized estimates apart from KM4, have estimates in the opposite direction to those reported by CF2. Those estimates that CF2 report as negative have swapped directions compared to their unregularized equivalents—the same effect we demonstrated earlier in our simulations—and are driven not by the data, but by high λ penalties.Footnote ³ Only in one other case—ridge regression with λ chosen by CV for the economy—are the estimates of response and cooperation rate effects negative. Even these estimates do not support an effect of response rates on mass polarization: the CIs for both overlap zero.

Appendix R3 presents further simulations that examine the likelihood of false positives using CF2’s approach if the true effect of response/cooperation rates was zero. The results are stark. We observe false positives more than 98% of the time for the effect of response/cooperation rates for energy, immigration, and foreign policy and more than 83% for the economy. In other words, CF2’s results are what we would expect to see if there were no true effect of response/cooperation rate on mass polarization.

Should we conclude that the effects of response/cooperation rates on mass polarization are null or even positive (i.e., the opposite of CF2’s claims)? While there are plausible theoretical models that lower response rates could deflate measures of mass polarization (see Appendix R4), the evidence is insufficient to make that claim. While some OLS estimates of the linear-time-trend model suggest positive effects, the ARIMA-MLM approach (and several other alternate model specifications) never find significant effects in either direction, suggesting that the linear-time-trend specification does not fully account for structural properties such as auto-correlation. We believe that ARIMA-MLM is a better model specification for removing spurious time effects, so prioritize its result over the linear specification used by CF2.

We also refrain from claiming a precise null effect. As Section S4 of the Supplementary Material shows, the data are underpowered even if the linear time trend specification is correct. We suspect that matched individual-level and administrative data (e.g., Clinton, Lapinski, and Trussler Reference Clinton, Lapinski and Trussler2022) may be a more fruitful avenue for exploring this question. Our claim is merely that the apparent evidence of an effect of response rates on mass polarization is driven by an inappropriate estimator and that the data are insufficient to test for effects of the claimed size.