Do linear models with shrinkage maintain the "controlling for" property of the predictors? Suppose I run a linear model with shrinkage, say elastic net, and only 2 coefficients remain. The model then looks like
$$y = \beta_1 x_1 + \beta_2 x_2 + \epsilon$$
Can I say that $\beta_1$ is the impact of $x_1$ controlling for $x_2$, like I would if this was an ordinary least squares fit? I feel like because of all the shrinkage it's no longer an orthogonal projection so that interpretation is not accurate. 
 A: $\newcommand{\l}{\lambda}$Consider a linear model $Y = X\beta + \varepsilon$ where $X$ is full rank. We all know the standard interpretation of $\hat \beta_j$: it is the estimated effect on the response for a unit increase in $X_j$ when the other predictors are held fixed.
First I will show where in the math this interpretation comes from, and then I will apply that analysis to ridge regression (I'm not considering the lasso right now).
We know $\hat \beta = (X^T X)^{-1}X^T Y$. Consider a partition of $X$ so that we have $X = \left(Z \ \vert\  W \right)$ where $Z \in \mathbb R^n$ is the first column of $X$ and $W \in \mathbb R^{n \times (p-1)}$ is the remaining $p-1$ columns. We could shuffle $X$ before this partitioning so in effect $Z$ represents an arbitrary column. We will work out the estimated coefficient $\hat \beta_Z$ of $Z$.
From its block form we know that
$$
X^TX =
\left[
\begin{array}{c|c}
Z^TZ & Z^TW \\
\hline
W^TZ & W^TW
\end{array}
\right]
$$
By assumption $X$ is full column rank so $X^TX$ is invertible, and this means that $W$ is full column rank too so $W^TW$ is also invertible.
Let $H = W(W^TW)^{-1}W^T$ be the projection matrix onto the column space of $W$. Then by the inverse of a block diagonal matrix we have
$$
(X^TX)^{-1} = 
\left[
\begin{array}{c|c}
(Z^T(I-H)Z)^{-1} & -(Z^T(I-H)Z)^{-1}Z^TW(W^TW)^{-1} \\
\hline
C & D
\end{array}
\right].
$$
We only care about the first row of $(X^TX)^{-1}$ since that's all we'll need for $\hat \beta_Z$ so we omit the details of $C$ and $D$.
Note that $Z^T(I-H)Z = 0$ iff $Z$ lies completely in the column space of $W$, which is impossible by the assumption that $X$ is full rank.
This means that
$$
\hat \beta_Z = \frac{Z^TY}{Z^T(I-H)Z} - \frac{Z^THY}{Z^T(I-H)Z} = \frac{Z^T(I-H)Y}{Z^T(I-H)Z} = \frac{e^TY}{e^Te}
$$
where $e := (I-H)Z$ is the residual vector from regressing $Z$ on $W$. I'm using the idempotence of $H$ to conclude that $Z^T(I-H)Z = \|(I- H)Z\|^2 = \|e\|^2$. All together, we see that $\hat \beta_Z$ is equal to the result of a simple linear regression of $Y$ on the part of $Z$ that is not explained by $W$. This is where we get the interpretation that the coefficient $\hat \beta_Z$ is the effect of $Z$ after controlling for the other predictors.

Now that we have a mathematical understanding of where this comes from, we can apply this to ridge regression to see how things change. The ridge regression estimate is 
$$
\hat \beta^{(\lambda)} = (X^T X + \lambda I)^{-1}X^T Y.
$$
We can recompute $(X^TX + \lambda I)^{-1}$ in the same manner as before. Let $H_\lambda = W(W^TW + \lambda I)^{-1}W^T$ be the hat matrix for a ridge regression onto the columns of $W$. We find that
$$
(X^TX + \lambda I)^{-1} = \left[
\begin{array}{c|c}
\left(Z^T(I - H_\lambda)Z + \lambda\right)^{-1} & -\left(Z^T(I - H_\lambda)Z + \lambda\right)^{-1}Z^TW(W^TW + \lambda I)^{-1} \\
\hline
C' & D'
\end{array}
\right].
$$
where again $C'$ and $D'$ don't matter. From this we find
$$
\hat \beta_Z^{(\lambda)} = \frac{Z^T(I - H_\lambda) Y}{Z^T(I-H_\lambda)Z + \lambda}.
$$
The interpretation here isn't quite so clean because now $H_\lambda$ and in turn $I - H_\lambda$ are no longer idempotent, so $Z^T(I-H_\lambda)Z$ is not exactly the norm of the residuals of a ridge regression of $Z$ on $W$. It's still similar enough though that we could just not worry about the loss of idempotence, but an alternative interpretation is to note that this is the result of a generalized least squares simple linear ridge regression weighted with $(I-H_\lambda)^{-1}$. I can use the Woodbury matrix identity to rearrange this as
$$
(I - H_\lambda)^{-1} = \lambda^{-1}\left(WW^T + \lambda I\right)
$$
so letting $\Omega_\l = WW^T + \l I$ I have
$$
\hat\beta_Z^{(\l)} = \frac{Z^T\Omega_\l^{-1} Y}{Z^T\Omega_\l^{-1}Z + \lambda}.
$$
What I like about this form is that this is the solution to the penalized Mahalanobis distance problem
$$
\underset{b \in \mathbb R}{\text{argmin }} (Y - b Z)^T\Omega_\l^{-1}(Y - bZ) + \l b^2
$$
and the GLS part of this is what I'd get from having an error term with covariance matrix $WW^T + \l I$ which is like a similarity matrix (this is also the sort of covariance matrix that arises from a 1-way random intercepts model). 
Let $W = U\tilde DV^T$ be the "wide" SVD of $W$ so that $U$ is $n\times n$ and actually is orthonormal, and $\tilde D = {D \choose \mathbf 0}$ is $n\times p$. Then
$$
WW^T + \l I = U\left(\begin{array}{c|c}D^2 + \l I_p& \mathbf 0 \\ \hline \mathbf 0 & \l I_{n-p}\end{array}\right)U^T
$$
so for a particular $b \in \mathbb R$ under consideration I'll get a working residual vector $r = Y - b Z$ and I measure its goodness of fit via
$$
r^T\Omega_\l^{-1} r = (U^Tr)^T \left(\begin{array}{c|c}D^2 + \l I_p& \mathbf 0 \\ \hline \mathbf 0 & \l I_{n-p}\end{array}\right)^{-1} (U^Tr) \\
= \sum_{i=1}^p \frac{c_i^2}{d_i^2 + \l} + \l^{-1}\sum_{j=p+1}^n c_i^2
$$
where $c = U^T r$ is $r$ expressed w.r.t. the basis $U$. This will be small when $c$ has most of its "mass" in the first few coordinates since those have the largest $d_i^2$, so a good $b$ is not one that makes $\|Y - b Z\|^2$ small but rather one that gets $Y - b Z$ mostly into the span of the first few columns of $U$. From the low-rank approximation view of the SVD I can think of these vectors as being the ones that best approximate the column space of $W$ which is what linear regression is all about. 
Thus to summarize all this, this Mahalanobis interpretation lets us see $W$ as being "controlled for" since $Z$ explaining $Y$ actually just requires explaining the part of $Y$ that's not well within the column space of $W$. And then that coefficient is further shrunken.
You ask whether or not orthogonal projections play a role here and my many references to idempotence indicate that they do indeed matter. When $\lambda=0$ so that $I - H_\lambda$ really is an orthogonal projection, the interpretation of $(I-H)Z$ is crystal clear: we know that we can decompose $Z = HZ + (I-H)Z$, so we are just taking the part of $Z$ that is completely orthogonal to the column space of $W$. When $\lambda > 0$ we still have $Z = H_\lambda Z + (I-H_\lambda)Z$ but this doesn't correspond to an orthogonal decomposition, so $(I-H_\lambda)Z$ can lie a little bit in the column space of $W$. We still can view this as controlling for $W$, because we're controlling for the effect of $W$ as estimated with this model, but it's not as tidy as when we really are taking the part of $Z$ that is completely orthogonal to $W$ and the residuals increasingly resemble $Z$ so the less flexible our model is the less we're able to control, but it is still happening.
Hope this helps!
