Is there a generalization of the "within" transformation to partial pooling?

Question

The random effects estimator for panel data can be cast as a penalized verion of the "fixed-effects" estimator from econometrics. In both cases, the model is $$ y_{it} = \mathbf{D}_{i}\alpha_i + \mathbf{X}_{it}\beta + \epsilon $$

Where $\mathbf{D}_{i}$ is a matrix of dummies indicating the cross-sectional unit, and $\alpha_i$ are the respective intercepts.

The random effects estimator used in mixed models is equivalent to ridge regression: $$ \hat\alpha = \left(\mathbf{D}^T\mathbf{D}+ \lambda\mathbf{I}\right)^{-1}\mathbf{D}^Ty $$

(ignoring the X's). Obviously the smoothing parameter needs to be selected. Ridge tends to use cross-validation, but likelihood-based methods also exist. These are implemented for example in mgcv, using the syntax

m <- gam(y ~ s(id, bs = "re"), method = "REML")

If I were to specify the smoothing parameter to zero via

m <- gam(y ~ s(id, bs = "re", sp = 0), method = "REML")

then the result would be identical to an OLS regression on the dummy variables, which is to say that $\lambda = 0$.

Now, fixed effects regressions in econometrics are commonly fitted using the "within" transformation (or more generally, the method of alternating projections). Basically, the data are de-meaned: $$ y_{it} - \bar{y}_i = \left(\mathbf{D}_{i} - \bar{\mathbf{D}}_{i}\right)\alpha_i + \left(\mathbf{X}_{it} - \bar{\mathbf{X}}_{i}\right)\beta + \epsilon $$

The individual intercepts disappear, and the model can be fit by OLS. This is much simpler to compute than approaches that actually evaluate $\left(\mathbf{D}^T\mathbf{D}+ \lambda\mathbf{I}\right)^{-1}$, because the dimension of that matrix tends to be huge.

So here are my questions.

First: For a given smoothing parameter $\lambda > 0$, is there a faster way to compute $[\hat\alpha, \hat\beta]$ than penalized OLS?? I'd like to avoid the matrix inversion especially.

Second: For a given $\lambda$, is there some $f(., \lambda)$ such that $$ y_{it} - f(y_{i}, \lambda) = \left(\mathbf{D}_{i} - f(\mathbf{D}_{i}, \lambda)\right)\alpha_i + \left(\mathbf{X}_{it} - f(\mathbf{X}_{i}, \lambda)\right)\beta + \epsilon $$ yields estimates identical to the ridge regression when fit with OLS?

The plm package uses quasi-demeaning as $f(, \theta) = f(\bar{x}_i\theta)$, where $\theta = 1−[\sigma_u^2 /(\sigma_u^2 +T \sigma_e^2 )]^.5$, where the variances are estimated in the first step of FGLS and $T$ is the number of time periods. But what if I have an exogenously-specified $\lambda$? Is there a one-to-one relationship between $\lambda$ and $\theta$?

Answer 1

First answer: Fastest way might be either: method of alternating projections to find $\hat{\beta}$ and then use the group means to find $\hat{\alpha}$ or you can use an iterative solver like gradient descent / stochastic gradient descent, conjugate gradient solvers, or most suited for your problem a sparse conjugate gradient solver such as LSMR or LSQR for each of the systems of equations.

Second answer: When I implemented fixed effects and ridge regression in my software I assumed users would like to get unbiased estimates of the fixed effects for absorbing to the model matrix (i.e., believe would be the right approach to solve endogeneity at least in panel data context). I first absorb the effects and then estimate the ridge regression (using standardized model matrix after absorbing fixed effects). I am unsure of a function for the second question, but invite you to reflect on why you want to penalize the fixed effects. Remember that, at least in panel data, true fixed effects are not even consistently estimated.