# How can I get information about coefficient std. error, t value and p value for Regularized Linear Model

In R we can get a lot of information from summary model.

> summary(lm(mpg~wt+hp,data=mtcars))

Call:
lm(formula = mpg ~ wt + hp, data = mtcars)

Residuals:
Min     1Q Median     3Q    Max
-3.941 -1.600 -0.182  1.050  5.854

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 37.22727    1.59879  23.285  < 2e-16 ***
wt          -3.87783    0.63273  -6.129 1.12e-06 ***
hp          -0.03177    0.00903  -3.519  0.00145 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.593 on 29 degrees of freedom
Multiple R-squared:  0.8268,    Adjusted R-squared:  0.8148
F-statistic: 69.21 on 2 and 29 DF,  p-value: 9.109e-12


Are Std. Error, t value P value still meaningful in regularized case? and How can we get the information in addition to "coefficient estimate" in ridge or lasso regression?

• I don't know if you get notified about my edits, but I've added a substantial bit more on ridge regression in GLMs. – jld Nov 3 '16 at 20:21

${\large \text{Section 1: Ridge regression with } \lambda \text{ known}}$

Linear Regression:

Under the assumption that $\lambda$ is known (which is a huge assumption) we can work out the distribution of the ridge regression estimator $\hat \beta_\lambda$ in an interpretable way.

We know that $\hat \beta_\lambda = (X^T X + \lambda I)^{-1} X^T Y$. Taking the SVD of $X$ we have $X = UDV^T$ so $$\hat \beta_\lambda = (V D^2 V^T + \lambda I)^{-1} V D U^T Y = V (D^2 + \lambda I)^{-1}D U^T Y$$

$$= V (D^2 + \lambda I)^{-1}D \left[DV^TVD^{-1} \right] U^T Y$$

$$= V (D^2 + \lambda I)^{-1}D^2V^T \left[VD^{-1} U^T Y\right]$$

$$= V (D^2 + \lambda I)^{-1}D^2V^T \hat \beta.$$

Letting $A_\lambda = V (D^2 + \lambda I)^{-1}D^2V^T$ we can see that we have $$\hat \beta_\lambda = A_\lambda \hat \beta \sim \mathcal N(A_\lambda \beta, \sigma^2 A_\lambda (X^T X)^{-1} A_\lambda).$$ This means that $A_\lambda$ captures both the increase in bias that we've incurred and the corresponding reduction in variance, and we can see how it does so by shrinking $\beta$ along the smaller principal components.

Working out the variance we have $$A_\lambda (X^T X)^{-1} A_\lambda = V (D^2 + \lambda I)^{-2} D^2 V^T$$

so

$$\hat \beta_\lambda \sim \mathcal N\left(V (D^2 + \lambda I)^{-1}D^2V^T \beta, \sigma^2V (D^2 + \lambda I)^{-2} D^2 V^T\right).$$

From this we can work out the usual tests since everything except $\beta$ and $\sigma^2$ is known in the above equation. It is important to note the role that the bias of $\hat \beta_\lambda$ plays in these tests, though, as described here for instance.

GLMs:

Here I'll discuss how to adapt the standard GLM procedure to handle an $L_2$ penalty on $\beta$. For a given $\lambda$ this will allow us to get $\hat \beta_\lambda$, and if we ignore uncertainty in $\lambda$ (i.e. assume it is known) we can then make inference about $\beta$. This is basically straight from McCullagh and Nelder's GLM book, so I'm going to be pretty light on details and just highlight how the addition of the penalty changes things.

We will maximize the penalized log-likelihood $l$ by Fisher scoring (which is equivalent to Newton-Raphson if the canonical link is used). This means that we iterate $$b^{(m+1)} = b^{(m)} - E\left( \frac{\partial l^2}{\partial \beta \partial \beta^T} \right)^{-1}_{(m)} \left(\frac{\partial l}{\partial \beta} \right)_{(m)}.$$

Letting $A_{(m)} = - E\left( \frac{\partial l^2}{\partial \beta \partial \beta^T} \right)_{(m)}$ we have

$$A_{(m)} b^{(m+1)} = A_{(m)} b^{(m)} + \left(\frac{\partial l}{\partial \beta} \right)_{(m)}.$$

In the aforementioned McCullagh and Nelder, they show that for the unpenalized case $A_{(m)} = X^T W^{(m)} X$, with $W^{(m)}$ being the diagonal weight matrix, so for our case we have $A_{(m)} = X^T W^{(m)} X + \lambda I$.

Again in the unpenalized case, we can show that $\left(\frac{\partial l}{\partial \beta} \right)_{(m)} = X^T W^{(m)} G^{(m)} (Y - \mu^{(m)})$ for $G^{(m)}$ diagonal with $G^{(m)}_{ii} = g'(\mu^{(m)}_i)$. Now, with the penalty, we have $\left(\frac{\partial l}{\partial \beta} \right)_{(m)} = X^T W^{(m)} G^{(m)} (Y - \mu^{(m)}) - \lambda b^{(m)}$.

Putting this together, we have $$(X^T W^{(m)} X + \lambda I) b^{(m+1)} = (X^T W^{(m)} X + \lambda I) b^{(m)} + X^T W^{(m)} G^{(m)} (Y - \mu^{(m)}) - \lambda b^{(m)}$$ $$= X^T W^{(m)} (X b^{(m)} + G^{(m)}(Y - \mu^{(m)})) = X^T W^{(m)} (g(\mu^{(m)}) + G^{(m)}(Y - \mu^{(m)}))$$

so we can see that this exactly corresponds to performing IRLS on the same working response as in the unpenalized case, except we simply use the weighted ridge regression estimator at each step rather than the unpenalized WLS estimator.

As for the variance and distribution of $\hat \beta_\lambda$, we can just get this from the final IRLS fit, so $$Var(\hat \beta_\lambda) = (X^T W X + \lambda I)^{-1} X^T W X (X^T W X + \lambda I)^{-1} = A I(\beta) A := \Sigma_\lambda$$

where $I(\beta)$ is the observed information coming from the actual likelihood of $\beta$ (which is not penalized).

Since linear weighted ridge regression is just WLS on a modified data set, and WLS is an MLE, we have that $\hat \beta_\lambda$ is asymptotically normal with variance $\Sigma_\lambda$ and we can base hypothesis tests off of this. As before, the bias of $\hat \beta_\lambda$ needs to be remembered.

${\large \text{Section 2: } \lambda \text{ unknown}}$

More realistically, if we do not want to assume that $\lambda$ is known then the main techniques available to us (to the best of my knowledge) are the bootstrap and a Bayesian analysis. Since ridge regression corresponds to a normal prior on $\beta$, and the normal distribution is its own conjugate prior, we have that the MAP and posterior mean are the same, i.e. $\hat \beta_\lambda$ is both the MAP and the posterior mean. This should make Bayesian inference very tractable. For the bootstrap, we still have the bias issue so it should be used with care. More details to come.

${\large \text{Section 3: The lasso}}$

For the lasso things get tricky but I believe that this paper by Lockhart, Taylor, and the Tibshiranis addresses your question (although I must confess I have only skimmed it -- as I have time I'll flesh this section out more).

• How do you account for uncertainty in the value of the parameter $\lambda$? – whuber Oct 28 '16 at 21:52
• How about logistic regression? – Haitao Du Oct 28 '16 at 21:54
• Good points both. I'll follow up in a bit. – jld Oct 28 '16 at 22:04
• @hxd1011 If we assume that $\lambda$ is known I wonder if we could use the linear mode result on the last iteration of IRLS – jld Oct 29 '16 at 0:00
• @whuber you're very welcome; that was a very good point. Overall this has been a great excuse to really dig into all of this stuff. – jld Nov 2 '16 at 17:54