price nbath nbed sqft city pets
1 1199 2 3 1200 Columbia yes
2 1199 2 3 1450 Columbia yes
3 1299 2 3 1950 Columbia yes
4 1175 2 3 1312 Columbia yes
5 1349 2 4 2015 Columbia yes
6 1225 2 3 1336 Columbia yesSTAT 516 Lec 04
Multiple linear regression (part 2/2)
Karl Gregory
2026-10-26
SC apartments data
As in part 1/2, consider these data from “Apartment for Rent Classified” (2019)
There are \(n = 214\) observations.
Setup
Consider data \((\mathbf{x}_1,Y_1),\dots,(\mathbf{x}_n,Y_n)\), with each \(\mathbf{x}_i = (x_{i1},\dots,x_{ip})^T\).
The multiple linear regression model is
\[ Y_i = \beta_0 + x_{i1}\beta_1 + \dots + x_{ip}\beta_p + \varepsilon_i, \quad i = 1,\dots,n, \] where
- \(\mathbf{x}_1,\dots,\mathbf{x}_n\) are vectors in \(\mathbb{R}^p\) of covariate or predictor values.
- \(Y_1,\dots,Y_n\) are the response values
- \(\beta_0, \beta_1, \dots, \beta_p\) are the regression coefficients.
- \(\varepsilon_1,\dots,\varepsilon_n\) are iid \(\text{Normal}(0,\sigma^2)\) error terms.
- \(\sigma^2\) is the error term variance.
Goals in multiple linear regression
In part 1/2, we addressed these goals:
- Estimate the regression coefficients \(\beta_0\) and \(\beta_1,\dots,\beta_p\).
- Estimate the error term variance \(\sigma^2\).
- Perform inference on \(\beta_1,\dots,\beta_p\).
- Build a CI for \(\beta_0 + \beta_1 x_{\operatorname{new},1} + \dots + \beta_p x_{\operatorname{new},p}\) at any \(\mathbf{x}_{\operatorname{new}}\).
- Build a prediction interval for \(Y\) at any \(\mathbf{x}_{\operatorname{new}}\).
- Decompose the variation in \(Y\) into (sums of) sums of squares.
- Check whether the model assumptions are satisfied.
- Identify outliers and understand their effects.
- Include interaction effects.
Goals in multiple linear regression (cont)
In part 2/2 we focus on these:
- Test for significance of a subset of covariates.
- Understand how correlations among the covariates affect inferences.
- Do variable selection.
Review of F distributions
For \(W_1\sim \chi_{\nu_1}^2(\phi)\), \(W_2 \sim \chi_{\nu_2}^2\) independent, \(\displaystyle R = \frac{W_1 / \nu_1}{W_2 / \nu_2} \sim F_{\nu_1,\nu_2}(\phi)\).
\(F_{\nu_1,\nu_2}(\phi)\) denotes the \(F\) distribution with
- numerator degrees of freedom \(\nu_1\)
- denominator degrees of freedom \(\nu_2\)
- noncentrality parameter \(\phi \geq 0\)
If \(\phi > 0\) the distribution is a non-central F distribution.
When \(\phi = 0\) we just write \(F_{\nu_1,\nu_2}\) to denote the “central” F distribution.
We will encounter ratios of sums of squares which have \(F\) distributions.
Plot of some F distribution pdfs
nu1 <- c(1,2,3,5,5,5,50,50)
nu2 <- c(3,3,3,10,10,10,50,50)
phi <- c(0,0,0,0,4,8,0,4)
f <- seq(.01,4,length=200)
dfmat <- matrix(0,length(f),200)
for(j in 1:length(nu1)){
dfmat[j,] <- df(f,df1 = nu1[j],df2=nu2[j],ncp=phi[j])
}
lab <- paste("(df1,df2,phi) = (",
apply(cbind(nu1,nu2,phi),1,paste,collapse = ","),
")",sep="")Plot of some F distribution pdfs
The overall F-test
We may wish to test whether any covariates are important, that is \[ \text{$H_0$: $\beta_j = 0$ for all $j = 1,\dots,p$}. \]
The overall F-test of significance is carried out as follows:
Fit the model with all the covariates and obtain the value \[ F_{\operatorname{test}} = \frac{\operatorname{MS}_{\operatorname{Reg}}}{\operatorname{MS}_{\operatorname{Error}}} \left( = \frac{\operatorname{SS}_{\operatorname{Reg}}/ p}{\operatorname{SS}_{\operatorname{Error}}/ (n - (p+1))} \right) \]
Reject \(H_0\) at \(\alpha\) if \(F_{\operatorname{test}} > F_{p,n - (p + 1),\alpha}\).
Obtain p-value is \(P(F > F_{\operatorname{test}})\), where \(F \sim F_{p,n-(p+1)}\).
This test statistic and p-value are reported by summary() on lm().
Exercise
Show that the test statistic of the overall F test can be written \[ F_{\operatorname{test}} = \frac{\operatorname{MS}_{\operatorname{Reg}}}{\operatorname{MS}_{\operatorname{Error}}} = \frac{(n - (p+1))}{p}\frac{R^2}{1 - R^2}, \] where \(R^2\) is the coefficient of determination.
SC apartments
Call:
lm(formula = log(price) ~ city + nbath + nbed + log(sqft) + pets +
log(sqft):city, data = scapts)
Residuals:
Min 1Q Median 3Q Max
-0.582 -0.133 -0.004 0.093 0.684
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.88590 0.68527 4.21 3.8e-05 ***
cityColumbia 1.52500 0.75905 2.01 0.0459 *
cityGreenville 0.87340 0.88640 0.99 0.3256
cityRock Hill 0.33480 1.36863 0.24 0.8070
nbath -0.13511 0.04508 -3.00 0.0031 **
nbed 0.07974 0.03470 2.30 0.0226 *
log(sqft) 0.62871 0.10617 5.92 1.3e-08 ***
petsyes -0.00804 0.03089 -0.26 0.7949
cityColumbia:log(sqft) -0.25861 0.11112 -2.33 0.0209 *
cityGreenville:log(sqft) -0.14616 0.13059 -1.12 0.2644
cityRock Hill:log(sqft) -0.07234 0.19835 -0.36 0.7157
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.207 on 203 degrees of freedom
Multiple R-squared: 0.43, Adjusted R-squared: 0.402
F-statistic: 15.3 on 10 and 203 DF, p-value: <2e-16Exercise
Suppose you fit a regression model with \(p = 10\) predictors on a data set with \(214\) observations, and you obtain \(\hat \sigma = 0.207\) and \(R^2 = 0.43\). It turns out that these pieces of information are sufficient to determine all the values in the ANOVA table. See if you can fill in all the values.
| Source | Df | SS | MS | F value | p-value |
|---|---|---|---|---|---|
| Reg | \(p\) | \(\operatorname{SS}_{\operatorname{Reg}}\) | \(\operatorname{MS}_{\operatorname{Reg}}\) | \(F_{\operatorname{test}}\) | \(P(F > F_{\operatorname{test}})\) |
| Error | \(n-(p+1)\) | \(\operatorname{SS}_{\operatorname{Error}}\) | \(\operatorname{MS}_{\operatorname{Error}}\) | ||
| Total | \(n-1\) | \(\operatorname{SS}_{\operatorname{Tot}}\) |
Testing for significance of a subset of covariates
Consider testing the significance of a subset \(D \subset \{1,\dots,p\}\) of covariates: \[ \text{$H_0$: $\beta_j = 0$ for all $j \in D$.} \] Use the full-reduced model F-test:
- Let \(s\) be the number of covariates in \(D\) and compute \[ F_{\operatorname{test}} =\frac{(\operatorname{SS}_{\operatorname{Error}}(\text{Reduced}) - \operatorname{SS}_{\operatorname{Error}}(\text{Full}))/s}{\operatorname{SS}_{\operatorname{Error}}(\text{Full})/(n -(p+1))}, \]
- “Full” is the model with all \(p\) covariates.
- “Reduced” is the model after dropping the covariates in \(D\).
Reject \(H_0\) at \(\alpha\) if \(F_{\operatorname{test}} > F_{s,n-(p+1),\alpha}\).
Obtain p-value as \(P(F > F_{\operatorname{test}})\), where \(F \sim F_{s,n-(p+1)}\).
SC apartments
Check whether the number of bedrooms or bathrooms contributes significantly to the rent.
That is, test \(H_0\): \(\beta_{\operatorname{nbath}} = \beta_{\operatorname{nbed}} = 0\).
Fit full and reduced model and obtain \(\operatorname{SS}_{\operatorname{Error}}\):
SC apartments
Full model: \[\begin{align*} \widehat{\log(\texttt{price})} = 2.886 &+ (1.525) \cdot \texttt{Columbia} + (0.873) \cdot \texttt{Greenville} \\ &+ (0.335) \cdot \texttt{Rock Hill} + (-0.135) \cdot \texttt{nbath} \\ &+ (0.08) \cdot \texttt{nbed} + (0.629) \log(\texttt{sqft}) \\ &+ (-0.008) \cdot \texttt{pets} \\ &+ (-0.259) \cdot \texttt{Columbia} \times \log(\texttt{sqft}) \\ &+ (-0.146) \cdot \texttt{Greenville} \times \log(\texttt{sqft}) \\ &+ (-0.072) \cdot \texttt{Rock Hill} \times \log(\texttt{sqft}). \end{align*}\]
Reduced model (without \(\texttt{nbath}\) and \(\texttt{nbed}\)): \[\begin{align*} \widehat{\log(\texttt{price})} = 3.15 &+ (0.892) \cdot \texttt{Columbia} + (0.872) \cdot \texttt{Greenville} \\ &+ (0.117) \cdot \texttt{Rock Hill} + (0.581) \log(\texttt{sqft}) \\ &+ (-0.007) \cdot \texttt{pets} \\ &+ (-0.165) \cdot \texttt{Columbia} \times \log(\texttt{sqft}) \\ &+ (-0.147) \cdot \texttt{Greenville} \times \log(\texttt{sqft}) \\ &+ (-0.044) \cdot \texttt{Rock Hill} \times \log(\texttt{sqft}). \end{align*}\]
We obtain \(\operatorname{SS}_{\operatorname{Error}}(\text{Full}) = 8.69\) and \(\operatorname{SS}_{\operatorname{Error}}(\text{Reduced}) = 9.125\).
SC apartments
Compute the test statistic for the full-reduced model F test.
We obtain \(F_{\operatorname{test}} = 5.078\) and \(F_{s,n-(p+1),0.05} = 3.04\), and the p-value is \(0.007\).
SC apartments
Check whether the city and its interactions are useful in the model.
SC apartments
Full model: \[\begin{align*} \widehat{\log(\texttt{price})} = 2.886 &+ (1.525) \cdot \texttt{Columbia} + (0.873) \cdot \texttt{Greenville} \\ &+ (0.335) \cdot \texttt{Rock Hill} + (-0.135) \cdot \texttt{nbath} \\ &+ (0.08) \cdot \texttt{nbed} + (0.629) \log(\texttt{sqft}) + (-0.008) \cdot \texttt{pets} \\ &+ (-0.259) \cdot \texttt{Columbia} \times \log(\texttt{sqft}) \\ &+ (-0.146) \cdot \texttt{Greenville} \times \log(\texttt{sqft}) \\ &+ (-0.072) \cdot \texttt{Rock Hill} \times \log(\texttt{sqft}). \end{align*}\]
Reduced model (without the city variables and their interactions): \[\begin{align*} \widehat{\log(\texttt{price})} = 3.832 &+ (-0.112) \cdot \texttt{nbath} + (0.042) \cdot \texttt{nbed} \\ &+ (0.483) \log(\texttt{sqft}) + (-0.06) \cdot \texttt{pets} \end{align*}\]
We obtain \(\operatorname{SS}_{\operatorname{Error}}(\text{Full}) = 8.69\) and \(\operatorname{SS}_{\operatorname{Error}}(\text{Reduced}) = 10.586\).
SC apartments
Compute the test statistic for the full-reduced model F test.
We obtain \(F_{\operatorname{test}} = 7.381\) and \(F_{s,n-(p+1),0.05} = 2.143\), and the p-value is \(3.723\times 10^{-7}\).
Full-reduced model F test for a single covariate
If we test \(H_0\): \(\beta_j = 0\) for a single covariate using the full-reduced model F test, then:
The test statistic \(F_{\operatorname{test}}\) will be equal to the square of the test statistic \(T_{\operatorname{test}}\) for testing \(H_0\): \(\beta_j = 0\) in the full model.
The p-value will be the same as that of the test based on \(T_{\operatorname{test}}\). This is due to the fact that \(T \sim t_\nu\) implies \(T^2 \sim F_{1,\nu}\).
SC apartments
Use the full-reduced model F test to test whether the pets variable has an influence on the price.
lm_red <- lm(log(price) ~ city + nbath + nbed + log(sqft) +
log(sqft):city, data = scapts)
lm_full <- lm(log(price) ~ city + nbath + nbed + log(sqft) + pets
+ log(sqft):city, data = scapts)
SSE_red <- sum(lm_red$residuals^2)
SSE_full <- sum(lm_full$residuals^2)
s <- 1 # significance of a single covariate being tested
Ftest <- (SSE_red - SSE_full)/s / ( SSE_full / (n - (p + 1)))
pval <- 1-pf(Ftest,s,n-(p+1))We obtain \(\sqrt{F} = 0.26 = |T_{\operatorname{test}}|\). The p value is \(0.795\). Compare this to that given in the summary() output.
Effect of correlations among the covariates
From before \(\operatorname{Var}\hat \beta_j = \sigma^2 \Omega_{jj}/n\). An alternate expression gives \[ \operatorname{Var}\hat \beta_j = \frac{1}{1 - R^2_j} \frac{\sigma^2}{\sum_{i=1}^n(x_{ji} - \bar x_j)^2}, \] where \(R^2_j\) is the \(R^2\) from regressing \(x_j\) on the other covariates.
So multicollinearity of \(x_j\) with the other covariates “inflates” \(\operatorname{Var}\hat \beta_j\):
Makes confidence intervals for \(\beta_j\) wider.
Makes tests of \(H_0\): \(\beta_j = 0\) vs \(H_1\): \(\beta_j \neq 0\) less powerful.
Call \(\displaystyle \frac{1}{1- R^2_j}\) the variance inflation factor (VIF) for \(x_j\), \(j = 1,\dots,p\).
SC apartments
Add to the data set a spurious predictor highly correlated with the square footage.
Check the effect of this on our inferences for \(\beta_{\log(\operatorname{sqft})}\).
SC apartments
SC apartments
2.5 % 97.5 %
(Intercept) 1.5347 4.2371
cityColumbia 0.0284 3.0216
cityGreenville -0.8743 2.6211
cityRock Hill -2.3638 3.0334
nbath -0.2240 -0.0462
nbed 0.0113 0.1482
log(sqft) 0.4194 0.8380
petsyes -0.0689 0.0529
cityColumbia:log(sqft) -0.4777 -0.0395
cityGreenville:log(sqft) -0.4036 0.1113
cityRock Hill:log(sqft) -0.4634 0.3187SC apartments
2.5 % 97.5 %
(Intercept) 1.5277 4.2407
cityColumbia 0.0250 3.0271
cityGreenville -0.8921 2.6512
cityRock Hill -2.3715 3.0395
nbath -0.2242 -0.0457
nbed 0.0111 0.1483
log(sqft) 0.1214 1.1149
petsyes -0.0692 0.0535
log(x) -0.4489 0.4705
cityColumbia:log(sqft) -0.4785 -0.0390
cityGreenville:log(sqft) -0.4081 0.1139
cityRock Hill:log(sqft) -0.4643 0.3199The width of the CI for \(\beta_{\operatorname{log(sqft)}}\) was \(0.419\).
With the new covariate the width of the CI for \(\beta_{\operatorname{log(sqft)}}\) becomes \(0.994\).
So including log(x) in the model makes our estimation of \(\beta_{\operatorname{log(sqft)}}\) less accurate.
SC apartments
Call:
lm(formula = log(price) ~ city + nbath + nbed + log(sqft) + pets +
log(sqft):city, data = scapts)
Residuals:
Min 1Q Median 3Q Max
-0.582 -0.133 -0.004 0.093 0.684
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.88590 0.68527 4.21 3.8e-05 ***
cityColumbia 1.52500 0.75905 2.01 0.0459 *
cityGreenville 0.87340 0.88640 0.99 0.3256
cityRock Hill 0.33480 1.36863 0.24 0.8070
nbath -0.13511 0.04508 -3.00 0.0031 **
nbed 0.07974 0.03470 2.30 0.0226 *
log(sqft) 0.62871 0.10617 5.92 1.3e-08 ***
petsyes -0.00804 0.03089 -0.26 0.7949
cityColumbia:log(sqft) -0.25861 0.11112 -2.33 0.0209 *
cityGreenville:log(sqft) -0.14616 0.13059 -1.12 0.2644
cityRock Hill:log(sqft) -0.07234 0.19835 -0.36 0.7157
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.207 on 203 degrees of freedom
Multiple R-squared: 0.43, Adjusted R-squared: 0.402
F-statistic: 15.3 on 10 and 203 DF, p-value: <2e-16The p-value for log(sqft) is very small.
SC apartments
Call:
lm(formula = log(price) ~ city + nbath + nbed + log(sqft) + pets +
log(sqft):city + log(x), data = scapts)
Residuals:
Min 1Q Median 3Q Max
-0.5816 -0.1335 -0.0037 0.0931 0.6840
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.8842 0.6879 4.19 4.1e-05 ***
cityColumbia 1.5261 0.7613 2.00 0.0463 *
cityGreenville 0.8796 0.8985 0.98 0.3288
cityRock Hill 0.3340 1.3721 0.24 0.8079
nbath -0.1350 0.0453 -2.98 0.0032 **
nbed 0.0797 0.0348 2.29 0.0230 *
log(sqft) 0.6181 0.2519 2.45 0.0150 *
petsyes -0.0079 0.0311 -0.25 0.7999
log(x) 0.0108 0.2332 0.05 0.9631
cityColumbia:log(sqft) -0.2588 0.1115 -2.32 0.0212 *
cityGreenville:log(sqft) -0.1471 0.1324 -1.11 0.2679
cityRock Hill:log(sqft) -0.0722 0.1988 -0.36 0.7168
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.207 on 202 degrees of freedom
Multiple R-squared: 0.43, Adjusted R-squared: 0.399
F-statistic: 13.9 on 11 and 202 DF, p-value: <2e-16The p-value for log(sqft) is not nearly as small when log(x) is included!
Getting VIFs with vif() from the car package
We can use the R package car from Fox and Weisberg (2019).
First time must install the package with install.package("car").
GVIF Df GVIF^(1/(2*Df))
city 1.89e+08 3 23.97
nbath 2.97e+00 1 1.72
nbed 3.97e+00 1 1.99
log(sqft) 6.18e+00 1 2.49
pets 1.19e+00 1 1.09
city:log(sqft) 1.97e+08 3 24.12 GVIF Df GVIF^(1/(2*Df))
city 1.95e+08 3 24.07
nbath 2.98e+00 1 1.73
nbed 3.97e+00 1 1.99
log(sqft) 3.46e+01 1 5.89
pets 1.20e+00 1 1.10
log(x) 3.01e+01 1 5.49
city:log(sqft) 2.02e+08 3 24.23Note the change in VIF for log(sqft) due to including x in the model!
Variable selection
Sometimes the number of potentially important predictors is quite large.
Large \(p\) tends to increase the VIFs, leading to low power.
So we may wish to discard some predictors. We briefly discuss:
- Best subset selection with Mallow’s \(C(p)\)
- Forward and backward stepwise selection with AIC
- LASSO selection
And most importantly:
- The dangers of naïve post-selection inference!!
Best subset selection with Mallow’s \(C_p\)
Given \(q\) available covariates, there are \(2^q\) possible subset models (why?).
Mallow’s \(C_p\) can be used to compare subset models: Let \[ C_p = (n - (p + 1))\left[\frac{\operatorname{MS}_{\operatorname{Error}}(\operatorname{subset})}{\operatorname{MS}_{\operatorname{Error}}(\operatorname{all})} - 1\right] + (p + 1), \] where
- \(p\) is the number of predictors in the subset model.
- \(\operatorname{MS}_{\operatorname{Error}}(\operatorname{subset})\) is the \(\operatorname{MS}_{\operatorname{Error}}\) of the subset model.
- \(\operatorname{MS}_{\operatorname{Error}}(\operatorname{all})\) is the \(\operatorname{MS}_{\operatorname{Error}}\) of the model with all the covariates.
If the subset model is adequate, \(\operatorname{MS}_{\operatorname{Error}}(\operatorname{subset})\) estimates the same target as \(\operatorname{MS}_{\operatorname{Error}}(\operatorname{all})\), so the first term should be small and \(C_p \approx p+1\).
Can look at \(C_p\) values for all subset models of each size \(p=0,1,2,\dots,q\)
Want smallest model such that \(C_p \approx p+1\).
SC apartments
Can run the following to identify best model of each size and obtain the Mallow’s \(C_p\) statistic for each of these models.
赤池 弘次 (あかいけひろつぐ)
Introduced Akaike’s Information Criterion (AIC).
Akaike’s Information Criterion (AIC) for comparing models
For a given model, i.e. set of covariates, AIC is defined as
\[ \operatorname{AIC}= 2(p+1) - 2~\underbrace{\ell(\hat \sigma^2, \hat \beta_0,\hat \beta_1,\dots,\hat \beta_p)}_{\text{log-likelihood}}. \] The log-likelihood is the log of the joint pdf of the data (STAT 512).
AIC can be used to compare several models for the same data.
The “best” model is the one which minimizes AIC.
Stepwise selection based on AIC
Stepwise selection:
- Backward: Begin with all the predictors and remove one at a time.
- Forward: Begin with no predictors and add one at a time.
In each step remove/add predictor to get largest decrease in AIC.
If a decrease in AIC is not possible, stop.
SC apartments
Below we run backward selection with the step() function:
Start: AIC=-664
log(price) ~ city + nbath + nbed + log(sqft) + pets + log(sqft):city
Df Sum of Sq RSS AIC
- pets 1 0.003 8.69 -666
- city:log(sqft) 3 0.244 8.93 -664
<none> 8.69 -664
- nbed 1 0.226 8.92 -660
- nbath 1 0.385 9.07 -656
Step: AIC=-666
log(price) ~ city + nbath + nbed + log(sqft) + city:log(sqft)
Df Sum of Sq RSS AIC
<none> 8.69 -666
- city:log(sqft) 3 0.257 8.95 -665
- nbed 1 0.224 8.92 -662
- nbath 1 0.386 9.08 -658SC apartments
The chosen model has all the covariates except \(\texttt{pets}\).
The fitted model is:
\[\begin{align*} \widehat{\log(\texttt{price})} = 2.886 &+ (1.525) \cdot \texttt{Columbia} + (0.873) \cdot \texttt{Greenville} \\ &+ (0.335) \cdot \texttt{Rock Hill} + (-0.135) \cdot \texttt{nbath} \\ &+ (0.08) \cdot \texttt{nbed} + (0.629) \log(\texttt{sqft}) \\ &+ (-0.008) \cdot \texttt{Columbia} \times \log(\texttt{sqft}) \\ &+ (-0.259) \cdot \texttt{Greenville} \times \log(\texttt{sqft}) \\ &+ (-0.146) \cdot \texttt{Rock Hill} \times \log(\texttt{sqft}). \end{align*}\]
SC apartments
Here we run forward selection with the step() function:
Start: AIC=-563
log(price) ~ 1
Df Sum of Sq RSS AIC
+ log(sqft) 1 4.16 11.1 -630
+ nbed 1 2.73 12.5 -604
+ nbath 1 1.31 13.9 -581
+ city 3 1.03 14.2 -572
<none> 15.2 -563
+ pets 1 0.02 15.2 -562
Step: AIC=-630
log(price) ~ log(sqft)
Df Sum of Sq RSS AIC
+ city 3 1.834 9.25 -662
+ nbath 1 0.267 10.82 -633
+ pets 1 0.194 10.89 -631
<none> 11.08 -630
+ nbed 1 0.000 11.08 -628
Step: AIC=-662
log(price) ~ log(sqft) + city
Df Sum of Sq RSS AIC
+ nbath 1 0.1529 9.10 -664
<none> 9.25 -662
+ nbed 1 0.0358 9.21 -661
+ pets 1 0.0101 9.24 -660
+ city:log(sqft) 3 0.1228 9.13 -659
Step: AIC=-664
log(price) ~ log(sqft) + city + nbath
Df Sum of Sq RSS AIC
+ nbed 1 0.1471 8.95 -665
<none> 9.10 -664
+ city:log(sqft) 3 0.1804 8.92 -662
+ pets 1 0.0077 9.09 -662
Step: AIC=-665
log(price) ~ log(sqft) + city + nbath + nbed
Df Sum of Sq RSS AIC
+ city:log(sqft) 3 0.2571 8.69 -666
<none> 8.95 -665
+ pets 1 0.0163 8.93 -664
Step: AIC=-666
log(price) ~ log(sqft) + city + nbath + nbed + log(sqft):city
Df Sum of Sq RSS AIC
<none> 8.69 -666
+ pets 1 0.0029 8.69 -664SC apartments
For these data forward and backward stepwise selection resulted in the same model.
But this need not be the case!
LASSO selection
The LASSO estimators \(\hat \beta_0^L,\hat \beta_1^L,\dots,\hat \beta_p^L\) are obtained by minimizing \[ Q_\lambda(b_0,b_1,\dots,b_p) = \sum_{i=1}^n(Y_i - (b_0 + b_1 x_{i1} + \dots + b_p x_{ip}))^2 + \lambda \sum_{j=1}^p|b_j|, \] where \(\lambda > 0\) is a tuning parameter.
The penalty term \(\lambda \sum_{j=1}^p|b_j|\) can cause \(\hat \beta^L_j = 0\) for some \(j\).
For \(\lambda\) large enough, all the \(\hat \beta^L_j\) will be equal to zero.
So LASSO performs variable selection and estimation simultaneously.
Drawback: Hard to build CIs based on \(\hat \beta_0^L,\hat \beta_1^L,\dots,\hat \beta_p^L\).
The dangers of post-selection inference
It is dangerous to:
- Ask the data which variables to include in the model.
- Use the same afterwards to perform inference on the coefficients.
Illustration:
Add 50 spurious predictors to the SC apartment prices data.
See how many we find to be significant.
SC apartments with spurious covariates added
Call:
lm(formula = log(price) ~ ., data = scaptsX)
Residuals:
Min 1Q Median 3Q Max
-0.4069 -0.1059 -0.0098 0.0918 0.5983
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.72e+00 5.73e-02 117.41 < 2e-16 ***
nbath -7.57e-02 4.37e-02 -1.73 0.08541 .
nbed 3.52e-02 4.37e-02 0.81 0.42173
sqft 4.61e-04 9.51e-05 4.85 2.9e-06 ***
cityColumbia -2.70e-01 4.41e-02 -6.12 7.4e-09 ***
cityGreenville -1.45e-01 3.99e-02 -3.65 0.00036 ***
cityRock Hill -1.16e-01 5.21e-02 -2.22 0.02801 *
petsyes -1.27e-04 3.30e-02 0.00 0.99692
x1 -4.26e-02 1.63e-02 -2.61 0.00988 **
x2 -6.73e-03 1.74e-02 -0.39 0.69872
x3 -1.34e-02 1.67e-02 -0.80 0.42602
x4 1.80e-02 1.63e-02 1.10 0.27188
x5 3.44e-02 1.57e-02 2.19 0.02968 *
x6 -4.50e-02 1.55e-02 -2.91 0.00412 **
x7 -5.93e-03 1.66e-02 -0.36 0.72086
x8 2.53e-03 1.68e-02 0.15 0.88078
x9 -2.73e-02 1.56e-02 -1.75 0.08167 .
x10 8.99e-03 1.71e-02 0.53 0.59902
x11 2.28e-02 1.51e-02 1.51 0.13291
x12 1.50e-03 1.67e-02 0.09 0.92848
x13 2.63e-03 1.53e-02 0.17 0.86330
x14 1.52e-02 1.80e-02 0.84 0.40091
x15 2.07e-02 1.59e-02 1.30 0.19416
x16 1.62e-02 1.62e-02 1.00 0.32077
x17 -2.21e-02 1.66e-02 -1.34 0.18368
x18 7.99e-04 1.57e-02 0.05 0.95951
x19 2.09e-02 1.75e-02 1.20 0.23334
x20 -8.96e-03 1.56e-02 -0.57 0.56701
x21 -2.82e-02 1.60e-02 -1.76 0.07985 .
x22 6.12e-03 1.68e-02 0.36 0.71668
x23 9.89e-03 1.53e-02 0.65 0.51923
x24 1.78e-03 1.57e-02 0.11 0.90977
x25 5.01e-03 1.47e-02 0.34 0.73381
x26 6.96e-04 1.59e-02 0.04 0.96513
x27 -2.46e-02 1.66e-02 -1.48 0.14005
x28 -2.10e-02 1.50e-02 -1.40 0.16293
x29 2.93e-02 1.60e-02 1.83 0.06888 .
x30 -1.50e-02 1.72e-02 -0.87 0.38512
x31 -2.98e-02 1.57e-02 -1.89 0.06018 .
x32 2.45e-03 1.50e-02 0.16 0.87082
x33 -2.56e-02 1.51e-02 -1.70 0.09152 .
x34 2.53e-02 1.60e-02 1.58 0.11566
x35 -1.67e-02 1.54e-02 -1.08 0.27967
x36 8.31e-03 1.55e-02 0.53 0.59346
x37 1.06e-02 1.53e-02 0.69 0.48948
x38 -2.86e-02 1.61e-02 -1.78 0.07704 .
x39 9.53e-03 1.71e-02 0.56 0.57831
x40 1.55e-02 1.52e-02 1.02 0.31080
x41 1.03e-02 1.61e-02 0.64 0.52228
x42 5.70e-03 1.54e-02 0.37 0.71096
x43 1.62e-03 1.61e-02 0.10 0.92012
x44 5.53e-03 1.64e-02 0.34 0.73577
x45 2.04e-02 1.60e-02 1.28 0.20373
x46 -1.29e-02 1.55e-02 -0.83 0.40503
x47 1.23e-02 1.58e-02 0.78 0.43867
x48 -4.58e-02 1.62e-02 -2.83 0.00526 **
x49 6.31e-05 1.39e-02 0.00 0.99638
x50 -7.86e-05 1.64e-02 0.00 0.99619
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.201 on 156 degrees of freedom
Multiple R-squared: 0.589, Adjusted R-squared: 0.438
F-statistic: 3.92 on 57 and 156 DF, p-value: 8.63e-12SC apartments with spurious covariates added
We reject \(H_0\): \(\beta_j = 0\) at \(\alpha = 0.05\) for \(4\) of the spurious predictors.
So the Type I error rate was \(4/50 = 0.08\).
Now do backwards stepwise selection to throw some variables away.
Then see how many of the spurious predictors we find “significant”.
SC apartments with spurious covariates added
Call:
lm(formula = log(price) ~ nbath + sqft + city + x1 + x5 + x6 +
x9 + x11 + x17 + x21 + x27 + x29 + x31 + x33 + x34 + x38 +
x45 + x48, data = scaptsX)
Residuals:
Min 1Q Median 3Q Max
-0.4359 -0.1020 -0.0111 0.0935 0.6932
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.68e+00 4.66e-02 143.47 < 2e-16 ***
nbath -6.78e-02 3.49e-02 -1.95 0.0532 .
sqft 5.52e-04 5.57e-05 9.90 < 2e-16 ***
cityColumbia -2.78e-01 3.88e-02 -7.16 1.6e-11 ***
cityGreenville -1.55e-01 3.43e-02 -4.53 1.0e-05 ***
cityRock Hill -1.24e-01 4.44e-02 -2.78 0.0059 **
x1 -3.43e-02 1.41e-02 -2.44 0.0156 *
x5 3.75e-02 1.34e-02 2.79 0.0057 **
x6 -3.67e-02 1.35e-02 -2.71 0.0073 **
x9 -2.68e-02 1.32e-02 -2.02 0.0443 *
x11 2.91e-02 1.31e-02 2.23 0.0272 *
x17 -2.68e-02 1.41e-02 -1.90 0.0591 .
x21 -2.80e-02 1.38e-02 -2.02 0.0447 *
x27 -2.40e-02 1.39e-02 -1.73 0.0859 .
x29 2.66e-02 1.38e-02 1.93 0.0552 .
x31 -2.03e-02 1.36e-02 -1.49 0.1384
x33 -2.04e-02 1.31e-02 -1.55 0.1220
x34 2.32e-02 1.40e-02 1.66 0.0985 .
x38 -2.22e-02 1.38e-02 -1.60 0.1103
x45 2.03e-02 1.37e-02 1.49 0.1391
x48 -3.28e-02 1.42e-02 -2.31 0.0218 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.19 on 193 degrees of freedom
Multiple R-squared: 0.543, Adjusted R-squared: 0.495
F-statistic: 11.5 on 20 and 193 DF, p-value: <2e-16SC apartments with spurious covariates added
Backwards stepwise selection keeps \(15\) of the \(50\) spurios predictors.
Among these \(15\), we reject \(H_0\): \(\beta_j = 0\) at \(\alpha = 0.05\) for \(7\) of them.
So the post-selection Type I error rate was \(7 / 15= 0.467\) .
WARNING: Selecting variables and then getting p-values in the selected model often leads to astonishingly inflated Type I error rates.
References
“Apartment for Rent Classified.” 2019. UCI Machine Learning Repository.
Fox, John, and Sanford Weisberg. 2019. An R Companion to Applied Regression. Third. Thousand Oaks CA: Sage. https://socialsciences.mcmaster.ca/jfox/Books/Companion/.