Lecture 6 Live Session: Model Building and Variable Selection
R code, explanations and interpretation for both case studies (Credit and Finance)
1 How to use this document
This document follows the Lecture 6 live-session slides in order. Each section has three parts:
- The R code from the slides.
- What the code does: an explanation of each important line.
- Interpretation: what the output means, both statistically and for a business decision.
Before you knit. Save Credit.csv and
Finance.csv in the same folder as this .Rmd
file. In RStudio, use Session → Set Working Directory → To Source
File Location.
Packages used. Install them once with the code
below. The chunk is set to eval = FALSE, so it does not run
each time you knit.
install.packages(c("dplyr", "broom", "purrr", "MASS", "car",
"leaps", "olsrr", "lmtest"))library(dplyr) # select(), the %>% pipe
library(broom) # glance(): one-row model summaries
library(purrr) # map_dfr(): apply a function to several models
library(MASS) # stepAIC(): AIC-based stepwise selection
library(car) # vif(): variance inflation factors
library(leaps) # regsubsets(): best-subsets selection
library(olsrr) # p-value-based selection, standardized coefficients, partial correlations
library(lmtest) # bptest(): Breusch-Pagan test for constant varianceMASS and dplyr both contain a function
called select(). Because MASS is loaded after
dplyr, this document always writes
dplyr::select() to make clear which one is meant.
The workflow used in both case studies
| Step | What we do |
|---|---|
| Frame | State the goal, the outcome variable and the candidate predictor variables |
| Explore | Correlation matrix; look for predictor variables that are near-duplicates |
| Fit | Fit the full model with every candidate predictor variable |
| Select | Backward, forward and stepwise selection |
| Best subsets | Exhaustive search for the best model of each size |
| Compare | AIC, BIC, Mallows’ Cp, adjusted R² |
| Diagnose | VIF and residual checks |
| Justify | Defend the final model in business terms |
The golden rule: selection methods propose a model; the analyst decides. Finish with diagnostics and a justification, not with one algorithm’s printout.
Three ways to choose predictor variables
- Stepwise family (backward, forward, both directions). These are greedy searches: each makes the best single move available at each step (add or drop one predictor variable). Backward and forward selection never reconsider a move; stepwise selection can add back a predictor variable it removed earlier, but it still makes one greedy move at a time. The final model is the best one reachable along that path, not necessarily the best model overall.
- Best subsets. Searches every combination of predictor variables and reports the best model of each size. It is exhaustive, not greedy.
- The criteria. Models are compared with adjusted R² (higher is better) and with AIC, BIC and Mallows’ Cp (lower is better). Each criterion rewards fit and penalizes model size. BIC penalizes size the most, so it tends to choose the smallest models.
2 Case Study 1: Credit
Goal. Using the Credit data, decide which variables belong in a model for Balance (average credit-card balance, in dollars), and justify the choice to a manager.
- n = 400 customers.
- Outcome variable:
Balance(continuous). - Numeric predictor variables:
Income,Limit,Rating,Cards,Age,Education. - Categorical predictor variables:
Gender,Student,Married,Ethnicity. - That gives 10 candidate predictor variables.
LimitandRatingoverlap heavily, so watch for that.
2.1 Explore: read the data and scan the correlations
credit <- read.csv("Credit.csv", row.names = 1, strip.white = TRUE)
str(credit)'data.frame': 400 obs. of 11 variables:
$ Income : num 14.9 106 104.6 148.9 55.9 ...
$ Limit : int 3606 6645 7075 9504 4897 8047 3388 7114 3300 6819 ...
$ Rating : int 283 483 514 681 357 569 259 512 266 491 ...
$ Cards : int 2 3 4 3 2 4 2 2 5 3 ...
$ Age : int 34 82 71 36 68 77 37 87 66 41 ...
$ Education: int 11 15 11 11 16 10 12 9 13 19 ...
$ Gender : chr "Male" "Female" "Male" "Female" ...
$ Student : chr "No" "Yes" "No" "No" ...
$ Married : chr "Yes" "Yes" "No" "No" ...
$ Ethnicity: chr "Caucasian" "Asian" "Asian" "Asian" ...
$ Balance : int 333 903 580 964 331 1151 203 872 279 1350 ...What the code does
read.csv()reads the file into a data frame calledcredit.row.names = 1tells R that the first column (an unnamed column of row numbers 1–400) holds row labels, not data. R uses it for row names, so it does not become a variable.strip.white = TRUEremoves extra spaces around unquoted text values.str()lists every variable with its type. We see 400 rows and 11 variables: the outcome variable plus the 10 candidates.
Correction to the slide code. The slide follows
read.csv(..., row.names = 1)withcredit <- credit[ , -1] # drop row numbers. Becauserow.names = 1has already removed the row-number column, that line deletes the first real variable, Income. It has been removed here. (Usecredit[ , -1]only when you read the file withoutrow.names = 1.)
# Scan correlations among the numeric variables
credit %>%
dplyr::select(where(is.numeric)) %>%
cor() %>%
round(2) Income Limit Rating Cards Age Education Balance
Income 1.00 0.79 0.79 -0.02 0.18 -0.03 0.46
Limit 0.79 1.00 1.00 0.01 0.10 -0.02 0.86
Rating 0.79 1.00 1.00 0.05 0.10 -0.03 0.86
Cards -0.02 0.01 0.05 1.00 0.04 -0.05 0.09
Age 0.18 0.10 0.10 0.04 1.00 0.00 0.00
Education -0.03 -0.02 -0.03 -0.05 0.00 1.00 -0.01
Balance 0.46 0.86 0.86 0.09 0.00 -0.01 1.00What the code does
%>%(the pipe) passes the result of one step into the next.dplyr::select(where(is.numeric))keeps only the numeric columns, becausecor()cannot use text or factor columns.cor()computes the Pearson correlation, r, for every pair of numeric variables, andround(2)shows two decimals.
Interpretation
| Pair | r | Meaning |
|---|---|---|
| Limit and Rating | 0.997 (shows as 1.00) | Near-duplicates. They carry almost the same information. |
| Limit and Balance | 0.86 | Strong positive linear association |
| Rating and Balance | 0.86 | About the same as Limit, as expected |
| Income and Balance | 0.46 | Moderate positive linear association |
A model that keeps both Limit and Rating will have multicollinearity: two predictor variables are so strongly related that the model cannot separate their effects. Their coefficients become unstable, their standard errors are inflated, and the variance inflation factor (VIF) will be very large. We flag this now and resolve it at the diagnosis step. With two decimals, r = .997 prints as 1.00; the next chunk uses three decimals to show the true value.
2.2 Correlation matrix including the binary predictor variables
cor() needs numbers. A binary factor
(one with two categories) can be entered as a 0/1 indicator
variable. Pearson r between a numeric variable and a
0/1 indicator is called the point-biserial
correlation.
# Categorical variables as factors
credit$Gender <- factor(trimws(credit$Gender))
credit$Student <- factor(credit$Student)
credit$Married <- factor(credit$Married)
credit$Ethnicity <- factor(credit$Ethnicity)
# Binary factors as 0/1 indicators (point-biserial r)
credit$Female <- ifelse(credit$Gender == "Female", 1, 0)
credit$Student_Yes <- ifelse(credit$Student == "Yes", 1, 0)
credit$Married_Yes <- ifelse(credit$Married == "Yes", 1, 0)
cor_vars <- c("Balance", "Income", "Limit", "Rating",
"Cards", "Age", "Education", "Female",
"Student_Yes", "Married_Yes")
cor_matrix <- cor(credit[ , cor_vars])
round(cor_matrix, 3) Balance Income Limit Rating Cards Age Education Female
Balance 1.000 0.464 0.862 0.864 0.086 0.002 -0.008 0.021
Income 0.464 1.000 0.792 0.791 -0.018 0.175 -0.028 -0.011
Limit 0.862 0.792 1.000 0.997 0.010 0.101 -0.024 0.009
Rating 0.864 0.791 0.997 1.000 0.053 0.103 -0.030 0.009
Cards 0.086 -0.018 0.010 0.053 1.000 0.043 -0.051 -0.023
Age 0.002 0.175 0.101 0.103 0.043 1.000 0.004 0.004
Education -0.008 -0.028 -0.024 -0.030 -0.051 0.004 1.000 -0.005
Female 0.021 -0.011 0.009 0.009 -0.023 0.004 -0.005 1.000
Student_Yes 0.259 0.020 -0.006 -0.002 -0.026 -0.030 0.072 0.055
Married_Yes -0.006 0.036 0.031 0.037 -0.010 -0.073 0.049 0.012
Student_Yes Married_Yes
Balance 0.259 -0.006
Income 0.020 0.036
Limit -0.006 0.031
Rating -0.002 0.037
Cards -0.026 -0.010
Age -0.030 -0.073
Education 0.072 0.049
Female 0.055 0.012
Student_Yes 1.000 -0.077
Married_Yes -0.077 1.000What the code does
factor()tells R to treat a variable as categorical. In a regression, R then creates indicator variables for it automatically.trimws()removes any leading or trailing spaces. Some copies of Credit.csv store the value as" Male"with a leading space. In this copy,strip.white = TRUEhas already removed it, but in a copy where the value is in quotation marks it does not.trimws()guarantees the two levels are"Female"and"Male"either way.ifelse(condition, 1, 0)creates each 0/1 indicator. Each indicator is named after the category coded 1:Female = 1means female;Student_Yes = 1andMarried_Yes = 1mean yes.cor_varslists the columns to correlate, andcredit[ , cor_vars]keeps only those columns.Ethnicityis left out. It has three unordered categories (African American, Asian, Caucasian), so no single numeric coding is meaningful.
Interpretation
- Correlations with Balance (first row):
- Rating (r = .864) and Limit (r = .862) have strong positive linear associations with Balance.
- Income (r = .464) has a moderate positive association.
- Cards (.086), Age (.002) and Education (−.008) show essentially no linear association.
- Binary predictor variables:
Student_Yes(r = .259): students (coded 1) tend to carry higher balances than non-students.Female(.021) andMarried_Yes(−.006) show essentially no linear association with Balance.- The coding sets the sign. Coding 1 = Male instead would flip the sign but not the size.
- Among the predictor variables:
- Limit and Rating correlate at r = .997, and both correlate about .79 with Income.
- Entering Limit and Rating together creates severe multicollinearity.
- Correlation describes linear association. It does not establish causation.
2.3 Fit the full model
full <- lm(Balance ~ Income + Limit + Rating + Cards + Age + Education +
Gender + Student + Married + Ethnicity, data = credit)
summary(full)
Call:
lm(formula = Balance ~ Income + Limit + Rating + Cards + Age +
Education + Gender + Student + Married + Ethnicity, data = credit)
Residuals:
Min 1Q Median 3Q Max
-161.64 -77.70 -13.49 53.98 318.20
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -489.86112 35.80118 -13.683 < 2e-16 ***
Income -7.80310 0.23423 -33.314 < 2e-16 ***
Limit 0.19091 0.03278 5.824 1.21e-08 ***
Rating 1.13653 0.49089 2.315 0.0211 *
Cards 17.72448 4.34103 4.083 5.40e-05 ***
Age -0.61391 0.29399 -2.088 0.0374 *
Education -1.09886 1.59795 -0.688 0.4921
GenderMale 10.65325 9.91400 1.075 0.2832
StudentYes 425.74736 16.72258 25.459 < 2e-16 ***
MarriedYes -8.53390 10.36287 -0.824 0.4107
EthnicityAsian 16.80418 14.11906 1.190 0.2347
EthnicityCaucasian 10.10703 12.20992 0.828 0.4083
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 98.79 on 388 degrees of freedom
Multiple R-squared: 0.9551, Adjusted R-squared: 0.9538
F-statistic: 750.3 on 11 and 388 DF, p-value: < 2.2e-16glance(full)# A tibble: 1 × 12
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 0.955 0.954 98.8 750. 1.11e-253 11 -2399. 4823. 4875.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>What the code does
lm()fits a multiple linear regression. The formula reads “Balance is predicted by Income, Limit, …, Ethnicity”.- R turns each factor into indicator variables. Gender, Student and
Married each get one (
GenderMale,StudentYes,MarriedYes). Ethnicity gets two (EthnicityAsian,EthnicityCaucasian), with African American as the reference category. summary()prints the coefficients, their t tests, R² and the overall F test.glance()(from broom) returns the same fit statistics as one row: R², adjusted R², AIC, BIC and others.- AIC scales differ between functions.
glance()andAIC()report AIC = 4823.4 for this model, butstepAIC()(next section) reports 3686.2 for the same model, because the functions use different additive constants. Only differences in AIC within one function’s output are meaningful, so never compare AIC values across functions.
Correction to the slide code. The slide writes
lm(Balance ~ ., data = credit). The dot means “every other column”. By this pointcreditalso contains the three 0/1 indicators created for the correlation matrix, so the dot would enter each binary variable twice, once as the factor and once as the indicator. R cannot estimate both and reports three coefficients asNA. Writing the 10 predictor variables out avoids this.
Interpretation
- R² = .955: the 10 candidates together explain about 95.5% of the variance in Balance.
- Six predictor variables are individually significant at α = .05: Income, Limit, Rating, Cards, Age (p = .037) and Student (p < .001).
- Four are not significant: Education (p = .492), Gender (p = .283), Married (p = .411) and Ethnicity (neither of its two indicators is significant; both p > .23).
- Reading two coefficients:
StudentYes= 425.75: holding the other predictor variables constant, students carry balances about $426 higher than non-students.Income= −7.80: holding the others constant, each additional unit of income (incomes are in thousands of dollars) is associated with a balance about $7.80 lower.- Note the sign of Income. On its own, Income correlates positively with Balance (r = .46), but its coefficient here is negative. Higher-income customers tend to have higher credit limits (r = .79), which goes with higher balances. Once Limit is held constant, higher income is associated with lower balances. A predictor variable’s role can change once other predictor variables are controlled.
- Collinearity warning: the full model keeps both Limit and Rating. Their individual coefficients are unreliable because of the collinearity, which the selection methods and diagnostics must deal with.
2.4 Selection 1–3: backward, forward and stepwise (AIC)
stepAIC() compares models using AIC
(Akaike Information Criterion): model fit, penalized for the number of
estimated coefficients. Lower AIC is better.
null <- lm(Balance ~ 1, data = credit)
summary(null)
Call:
lm(formula = Balance ~ 1, data = credit)
Residuals:
Min 1Q Median 3Q Max
-520.02 -451.26 -60.51 342.98 1478.98
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 520.01 22.99 22.62 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 459.8 on 399 degrees of freedomWhat the code does. Balance ~ 1 is the
intercept-only (null) model: it has no predictor
variables and predicts the mean balance for everyone. Forward selection
starts here. Its single coefficient, 520.01 (printed rounded), is the
sample mean of Balance: about $520.
2.4.1 Backward selection
backward <- stepAIC(full, direction = "backward") # drop termsStart: AIC=3686.22
Balance ~ Income + Limit + Rating + Cards + Age + Education +
Gender + Student + Married + Ethnicity
Df Sum of Sq RSS AIC
- Ethnicity 2 14084 3800814 3683.7
- Education 1 4615 3791345 3684.7
- Married 1 6619 3793349 3684.9
- Gender 1 11269 3798000 3685.4
<none> 3786730 3686.2
- Age 1 42558 3829288 3688.7
- Rating 1 52314 3839044 3689.7
- Cards 1 162702 3949432 3701.0
- Limit 1 331050 4117780 3717.7
- Student 1 6326012 10112742 4077.1
- Income 1 10831162 14617892 4224.5
Step: AIC=3683.7
Balance ~ Income + Limit + Rating + Cards + Age + Education +
Gender + Student + Married
Df Sum of Sq RSS AIC
- Married 1 4545 3805359 3682.2
- Education 1 4757 3805572 3682.2
- Gender 1 10760 3811574 3682.8
<none> 3800814 3683.7
- Age 1 45650 3846464 3686.5
- Rating 1 49473 3850287 3686.9
- Cards 1 166806 3967621 3698.9
- Limit 1 340250 4141064 3716.0
- Student 1 6372573 10173387 4075.5
- Income 1 10838891 14639705 4221.1
Step: AIC=3682.18
Balance ~ Income + Limit + Rating + Cards + Age + Education +
Gender + Student
Df Sum of Sq RSS AIC
- Education 1 5399 3810759 3680.7
- Gender 1 11019 3816378 3681.3
<none> 3805359 3682.2
- Age 1 43545 3848904 3684.7
- Rating 1 46929 3852289 3685.1
- Cards 1 170729 3976088 3697.7
- Limit 1 352112 4157472 3715.6
- Student 1 6461978 10267338 4077.2
- Income 1 10860901 14666261 4219.8
Step: AIC=3680.75
Balance ~ Income + Limit + Rating + Cards + Age + Gender + Student
Df Sum of Sq RSS AIC
- Gender 1 10861 3821620 3679.9
<none> 3810759 3680.7
- Age 1 43983 3854741 3683.3
- Rating 1 49520 3860279 3683.9
- Cards 1 170898 3981656 3696.3
- Limit 1 347654 4158413 3713.7
- Student 1 6472072 10282831 4075.8
- Income 1 10855780 14666539 4217.8
Step: AIC=3679.89
Balance ~ Income + Limit + Rating + Cards + Age + Student
Df Sum of Sq RSS AIC
<none> 3821620 3679.9
- Age 1 44472 3866091 3682.5
- Rating 1 49263 3870883 3683.0
- Cards 1 172930 3994549 3695.6
- Limit 1 347898 4169518 3712.7
- Student 1 6462599 10284218 4073.9
- Income 1 10845018 14666638 4215.8summary(backward)
Call:
lm(formula = Balance ~ Income + Limit + Rating + Cards + Age +
Student, data = credit)
Residuals:
Min 1Q Median 3Q Max
-170.00 -77.85 -11.84 56.87 313.52
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -493.73419 24.82476 -19.889 < 2e-16 ***
Income -7.79508 0.23342 -33.395 < 2e-16 ***
Limit 0.19369 0.03238 5.981 4.98e-09 ***
Rating 1.09119 0.48480 2.251 0.0250 *
Cards 18.21190 4.31865 4.217 3.08e-05 ***
Age -0.62406 0.29182 -2.139 0.0331 *
StudentYes 425.60994 16.50956 25.780 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 98.61 on 393 degrees of freedom
Multiple R-squared: 0.9547, Adjusted R-squared: 0.954
F-statistic: 1380 on 6 and 393 DF, p-value: < 2.2e-16What the code does
- Starting from the full model,
stepAIC()tries removing each predictor variable one at a time, and removes the one whose removal lowers AIC the most. - It repeats this until no removal lowers AIC.
- Each block of printed output is one step.
How to read each step
- The heading line gives the current model and its AIC.
- Each row is one possible move:
-means drop that predictor variable. - The table is sorted from lowest to highest AIC.
<none>means “make no change”. It is the benchmark.- Df is the number of coefficients the move removes. Ethnicity has 2, because its three levels use two coefficients.
- Sum of Sq is the increase in the residual sum of squares from dropping that predictor variable.
- RSS is the residual sum of squares after the move.
- AIC is the model’s AIC after the move.
- Rule: take the top move if its AIC is below
<none>; otherwise stop.
Interpretation, step by step
| Step | Model AIC | Move taken | Why |
|---|---|---|---|
| 1 | 3686.2 | Drop Ethnicity | Lowest AIC (3683.7). Dropping Education, Married or Gender would also lower AIC, but only one move is made per step. |
| 2 | 3683.7 | Drop Married | Lowest AIC (3682.2). Education nearly ties; Married’s RSS is slightly smaller. |
| 3 | 3682.2 | Drop Education | Lowest AIC (3680.7). Dropping Gender also helps, but less. |
| 4 | 3680.7 | Drop Gender | Lowest AIC (3679.9), and the only move below
<none> |
| 5 | 3679.9 | Stop | <none> is the top row. Dropping any remaining
predictor variable raises AIC. |
Importance in the AIC table compares two numbers
from the same step: the AIC after dropping a predictor variable minus
the AIC of <none>.
- At step 3, dropping Income would raise AIC by 4219.8 − 3682.2 = +537.6, and dropping Student by 4077.2 − 3682.2 = +395.0. These two carry information no other predictor variable supplies. (These differences use the rounded AIC values shown in the output.)
- Dropping Rating raises AIC by only 3685.1 − 3682.2 = +2.9, because Limit (r = .997) carries almost the same information.
Final backward model: Income, Limit, Rating, Cards, Age and Student (AIC = 3679.9). Backward selection never adds a predictor variable back: once removed, a variable stays out.
2.4.2 Forward selection
forward <- stepAIC(null, direction = "forward",
scope = list(lower = formula(null), upper = formula(full))) # add termsStart: AIC=4905.56
Balance ~ 1
Df Sum of Sq RSS AIC
+ Rating 1 62904790 21435122 4359.6
+ Limit 1 62624255 21715657 4364.8
+ Income 1 18131167 66208745 4810.7
+ Student 1 5658372 78681540 4879.8
+ Cards 1 630416 83709496 4904.6
<none> 84339912 4905.6
+ Gender 1 38892 84301020 4907.4
+ Education 1 5481 84334431 4907.5
+ Married 1 2715 84337197 4907.5
+ Age 1 284 84339628 4907.6
+ Ethnicity 2 18454 84321458 4909.5
Step: AIC=4359.63
Balance ~ Rating
Df Sum of Sq RSS AIC
+ Income 1 10902581 10532541 4077.4
+ Student 1 5735163 15699959 4237.1
+ Age 1 649110 20786012 4349.3
+ Cards 1 138580 21296542 4359.0
+ Married 1 118209 21316913 4359.4
<none> 21435122 4359.6
+ Education 1 27243 21407879 4361.1
+ Gender 1 16065 21419057 4361.3
+ Limit 1 7960 21427162 4361.5
+ Ethnicity 2 51100 21384022 4362.7
Step: AIC=4077.41
Balance ~ Rating + Income
Df Sum of Sq RSS AIC
+ Student 1 6305322 4227219 3714.2
+ Married 1 95068 10437473 4075.8
+ Limit 1 94545 10437996 4075.8
+ Age 1 90286 10442255 4076.0
<none> 10532541 4077.4
+ Education 1 20819 10511722 4078.6
+ Ethnicity 2 67040 10465501 4078.9
+ Cards 1 2094 10530447 4079.3
+ Gender 1 948 10531593 4079.4
Step: AIC=3714.24
Balance ~ Rating + Income + Student
Df Sum of Sq RSS AIC
+ Limit 1 194718 4032502 3697.4
+ Age 1 44620 4182600 3712.0
<none> 4227219 3714.2
+ Married 1 13083 4214137 3715.0
+ Gender 1 12168 4215051 3715.1
+ Cards 1 10608 4216611 3715.2
+ Education 1 1400 4225820 3716.1
+ Ethnicity 2 22322 4204897 3716.1
Step: AIC=3697.37
Balance ~ Rating + Income + Student + Limit
Df Sum of Sq RSS AIC
+ Cards 1 166410 3866091 3682.5
+ Age 1 37952 3994549 3695.6
<none> 4032502 3697.4
+ Gender 1 13345 4019157 3698.0
+ Married 1 6660 4025842 3698.7
+ Education 1 5795 4026707 3698.8
+ Ethnicity 2 17704 4014797 3699.6
Step: AIC=3682.52
Balance ~ Rating + Income + Student + Limit + Cards
Df Sum of Sq RSS AIC
+ Age 1 44472 3821620 3679.9
<none> 3866091 3682.5
+ Gender 1 11350 3854741 3683.3
+ Education 1 5672 3860419 3683.9
+ Married 1 3121 3862970 3684.2
+ Ethnicity 2 14756 3851335 3685.0
Step: AIC=3679.89
Balance ~ Rating + Income + Student + Limit + Cards + Age
Df Sum of Sq RSS AIC
<none> 3821620 3679.9
+ Gender 1 10860.9 3810759 3680.7
+ Married 1 5450.6 3816169 3681.3
+ Education 1 5241.7 3816378 3681.3
+ Ethnicity 2 11517.3 3810102 3682.7summary(forward)
Call:
lm(formula = Balance ~ Rating + Income + Student + Limit + Cards +
Age, data = credit)
Residuals:
Min 1Q Median 3Q Max
-170.00 -77.85 -11.84 56.87 313.52
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -493.73419 24.82476 -19.889 < 2e-16 ***
Rating 1.09119 0.48480 2.251 0.0250 *
Income -7.79508 0.23342 -33.395 < 2e-16 ***
StudentYes 425.60994 16.50956 25.780 < 2e-16 ***
Limit 0.19369 0.03238 5.981 4.98e-09 ***
Cards 18.21190 4.31865 4.217 3.08e-05 ***
Age -0.62406 0.29182 -2.139 0.0331 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 98.61 on 393 degrees of freedom
Multiple R-squared: 0.9547, Adjusted R-squared: 0.954
F-statistic: 1380 on 6 and 393 DF, p-value: < 2.2e-16What the code does
- Forward selection starts from the null model and, at each step, adds the single predictor variable that lowers AIC the most.
scopesets the search range. The smallest model allowed is the null model; the largest is the full model, so the 10 candidates are the variables that can be added.- Only
+rows appear, because forward selection never removes a variable. For a+row, Sum of Sq is the reduction in RSS from adding the predictor variable.
Interpretation, step by step
| Step | Current AIC | Variable added | New AIC | Change |
|---|---|---|---|---|
| 1 | 4905.6 | Rating | 4359.6 | −546.0 |
| 2 | 4359.6 | Income | 4077.4 | −282.2 |
| 3 | 4077.4 | Student | 3714.2 | −363.2 |
| 4 | 3714.2 | Limit | 3697.4 | −16.8 |
| 5 | 3697.4 | Cards | 3682.5 | −14.9 |
| 6 | 3682.5 | Age | 3679.9 | −2.6 |
| 7 | 3679.9 | Stop | Every addition raises AIC; the smallest increase, adding Gender, is +0.8. |
The changes are computed from the rounded AIC values printed in the output, so a few differ by 0.1 from the unrounded values (for example, −545.9 rather than −546.0 at step 1).
Key teaching point: Limit. A predictor variable’s value depends on what is already in the model.
- Step 1: Limit almost ties Rating (AIC 4364.8 against 4359.6).
- Step 2: with Rating in the model, adding Limit would raise AIC (4361.5), because it repeats Rating’s information.
- Step 3: once Income is in the model, adding Limit would lower AIC slightly (4075.8 against 4077.4), but adding Student lowers it far more, so Student enters first.
- Step 4: Limit enters, because it still adds a little unique information beyond Rating, Income and Student.
Final forward model: Rating, Income, Student, Limit, Cards and Age. These are the same six predictor variables as backward selection, added in a different order.
2.4.3 Stepwise selection (both directions)
stepwise <- stepAIC(full, direction = "both") # add or dropStart: AIC=3686.22
Balance ~ Income + Limit + Rating + Cards + Age + Education +
Gender + Student + Married + Ethnicity
Df Sum of Sq RSS AIC
- Ethnicity 2 14084 3800814 3683.7
- Education 1 4615 3791345 3684.7
- Married 1 6619 3793349 3684.9
- Gender 1 11269 3798000 3685.4
<none> 3786730 3686.2
- Age 1 42558 3829288 3688.7
- Rating 1 52314 3839044 3689.7
- Cards 1 162702 3949432 3701.0
- Limit 1 331050 4117780 3717.7
- Student 1 6326012 10112742 4077.1
- Income 1 10831162 14617892 4224.5
Step: AIC=3683.7
Balance ~ Income + Limit + Rating + Cards + Age + Education +
Gender + Student + Married
Df Sum of Sq RSS AIC
- Married 1 4545 3805359 3682.2
- Education 1 4757 3805572 3682.2
- Gender 1 10760 3811574 3682.8
<none> 3800814 3683.7
+ Ethnicity 2 14084 3786730 3686.2
- Age 1 45650 3846464 3686.5
- Rating 1 49473 3850287 3686.9
- Cards 1 166806 3967621 3698.9
- Limit 1 340250 4141064 3716.0
- Student 1 6372573 10173387 4075.5
- Income 1 10838891 14639705 4221.1
Step: AIC=3682.18
Balance ~ Income + Limit + Rating + Cards + Age + Education +
Gender + Student
Df Sum of Sq RSS AIC
- Education 1 5399 3810759 3680.7
- Gender 1 11019 3816378 3681.3
<none> 3805359 3682.2
+ Married 1 4545 3800814 3683.7
- Age 1 43545 3848904 3684.7
+ Ethnicity 2 12011 3793349 3684.9
- Rating 1 46929 3852289 3685.1
- Cards 1 170729 3976088 3697.7
- Limit 1 352112 4157472 3715.6
- Student 1 6461978 10267338 4077.2
- Income 1 10860901 14666261 4219.8
Step: AIC=3680.75
Balance ~ Income + Limit + Rating + Cards + Age + Gender + Student
Df Sum of Sq RSS AIC
- Gender 1 10861 3821620 3679.9
<none> 3810759 3680.7
+ Education 1 5399 3805359 3682.2
+ Married 1 5187 3805572 3682.2
- Age 1 43983 3854741 3683.3
+ Ethnicity 2 12029 3798730 3683.5
- Rating 1 49520 3860279 3683.9
- Cards 1 170898 3981656 3696.3
- Limit 1 347654 4158413 3713.7
- Student 1 6472072 10282831 4075.8
- Income 1 10855780 14666539 4217.8
Step: AIC=3679.89
Balance ~ Income + Limit + Rating + Cards + Age + Student
Df Sum of Sq RSS AIC
<none> 3821620 3679.9
+ Gender 1 10861 3810759 3680.7
+ Married 1 5451 3816169 3681.3
+ Education 1 5242 3816378 3681.3
- Age 1 44472 3866091 3682.5
+ Ethnicity 2 11517 3810102 3682.7
- Rating 1 49263 3870883 3683.0
- Cards 1 172930 3994549 3695.6
- Limit 1 347898 4169518 3712.7
- Student 1 6462599 10284218 4073.9
- Income 1 10845018 14666638 4215.8summary(stepwise)
Call:
lm(formula = Balance ~ Income + Limit + Rating + Cards + Age +
Student, data = credit)
Residuals:
Min 1Q Median 3Q Max
-170.00 -77.85 -11.84 56.87 313.52
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -493.73419 24.82476 -19.889 < 2e-16 ***
Income -7.79508 0.23342 -33.395 < 2e-16 ***
Limit 0.19369 0.03238 5.981 4.98e-09 ***
Rating 1.09119 0.48480 2.251 0.0250 *
Cards 18.21190 4.31865 4.217 3.08e-05 ***
Age -0.62406 0.29182 -2.139 0.0331 *
StudentYes 425.60994 16.50956 25.780 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 98.61 on 393 degrees of freedom
Multiple R-squared: 0.9547, Adjusted R-squared: 0.954
F-statistic: 1380 on 6 and 393 DF, p-value: < 2.2e-16What the code does. With
direction = "both", every step considers both dropping a
predictor variable that is in the model (- rows) and adding
back one that was removed earlier (+ rows).
Interpretation
- On these data the path matches backward selection: drop Ethnicity, Married, Education, then Gender, then stop.
- Why the paths match: at every step the search also
checked whether adding back an earlier removal would help, and it never
did. For example, at step 3 adding back Married (AIC 3683.7) or
Ethnicity (3684.9) would have raised AIC above
<none>(3682.2). - Final model: Income, Limit, Rating, Cards, Age and Student (AIC = 3679.9).
Summary of the three AIC searches
| Method | Selected predictor variables | Number |
|---|---|---|
| Backward (AIC) | Income, Limit, Rating, Cards, Age, Student | 6 |
| Forward (AIC) | the same six | 6 |
| Stepwise (AIC) | the same six | 6 |
All three agree, but all three keep both Limit and Rating. The AIC penalty for one extra coefficient is small, and Rating still lowers AIC slightly (by about 3 points), so a greedy AIC search does not flag the redundancy. Agreement among methods is evidence that the result is stable, not proof that it is the best model.
2.5 Diagnose and justify the final model
Final_VS <- lm(Balance ~ Income + Limit + Rating + Cards + Age + Student, data = credit)
summary(Final_VS)
Call:
lm(formula = Balance ~ Income + Limit + Rating + Cards + Age +
Student, data = credit)
Residuals:
Min 1Q Median 3Q Max
-170.00 -77.85 -11.84 56.87 313.52
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -493.73419 24.82476 -19.889 < 2e-16 ***
Income -7.79508 0.23342 -33.395 < 2e-16 ***
Limit 0.19369 0.03238 5.981 4.98e-09 ***
Rating 1.09119 0.48480 2.251 0.0250 *
Cards 18.21190 4.31865 4.217 3.08e-05 ***
Age -0.62406 0.29182 -2.139 0.0331 *
StudentYes 425.60994 16.50956 25.780 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 98.61 on 393 degrees of freedom
Multiple R-squared: 0.9547, Adjusted R-squared: 0.954
F-statistic: 1380 on 6 and 393 DF, p-value: < 2.2e-16vif(Final_VS) Income Limit Rating Cards Age Student
2.776906 229.238479 230.869514 1.439007 1.039696 1.009064 What the code does
Final_VSrefits the six-predictor model chosen by the stepwise searches.vif()computes each predictor variable’s variance inflation factor: how much the variance of its coefficient is inflated by its correlation with the other predictor variables.- A VIF of 1 means no inflation. Values above 5 to 10 are a common warning sign.
Interpretation
- The VIFs for Limit (229.2) and Rating (230.9) are enormous: the standard errors of these two coefficients are inflated by a factor of about √230 ≈ 15. This results from r = .997 between Limit and Rating, together with both variables’ strong correlation with Income (about .79).
- The other VIFs are small: Income 2.8, Cards 1.4, Age 1.0, Student 1.0.
- Decision: drop one of the pair. Rating is dropped and Limit is kept. Limit is the variable the business controls directly, and the two carry almost the same information.
final <- lm(Balance ~ Income + Limit + Cards + Age + Student, data = credit)
summary(final)
Call:
lm(formula = Balance ~ Income + Limit + Cards + Age + Student,
data = credit)
Residuals:
Min 1Q Median 3Q Max
-187.05 -79.57 -12.59 56.06 322.56
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.673e+02 2.199e+01 -21.250 < 2e-16 ***
Income -7.760e+00 2.341e-01 -33.149 < 2e-16 ***
Limit 2.661e-01 3.535e-03 75.296 < 2e-16 ***
Cards 2.355e+01 3.628e+00 6.492 2.55e-10 ***
Age -6.220e-01 2.933e-01 -2.120 0.0346 *
StudentYes 4.284e+02 1.655e+01 25.886 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 99.12 on 394 degrees of freedom
Multiple R-squared: 0.9541, Adjusted R-squared: 0.9535
F-statistic: 1638 on 5 and 394 DF, p-value: < 2.2e-16vif(final) Income Limit Cards Age Student
2.764201 2.703384 1.004962 1.039685 1.003462 Interpretation of the final Credit model
- Fit: R² = .954 and adjusted R² = .954. That is essentially the same as the full model (R² = .955), with half as many predictor variables.
- Collinearity is gone: every VIF is below 3 (Income 2.76, Limit 2.70, Cards 1.00, Age 1.04, Student 1.00).
- Coefficients, each holding the other predictor
variables constant:
- Limit: b = 0.266, p < .001. Each additional $1 of credit limit is associated with about $0.27 more in balance, so $1,000 more limit goes with about $266 more balance.
- Student: b = 428.38, p < .001. Students carry balances about $428 higher than non-students.
- Income: b = −7.76, p < .001. Each additional unit of income (thousands of dollars) is associated with about $7.76 less balance.
- Cards: b = 23.55, p < .001. Each additional card is associated with about $24 more balance.
- Age: b = −0.62, p = .035. Each additional year of age is associated with about $0.62 less balance: statistically significant but practically very small.
2.5.1 Residual checks for the final Credit model
The workflow requires residual checks as well as the VIF.
library(lmtest) # bptest()
par(mfrow = c(1, 3))
plot(final, which = 1) # linearity
plot(final, which = 3) # constant variance
hist(resid(final), breaks = 30, xlab = "Residual",
main = "Histogram of residuals") # normality
par(mfrow = c(1, 1))
bptest(final)
studentized Breusch-Pagan test
data: final
BP = 128.37, df = 5, p-value < 2.2e-16sum(credit$Balance == 0) # customers with a $0 balance[1] 90sum(fitted(final) < 0) # customers with a negative predicted balance[1] 62Interpretation. Unlike the Finance model in Case Study 2, the Credit model shows clear problems:
- Linearity: the residuals-versus-fitted plot shows a strong curved pattern rather than an even band.
- Constant variance: the height of the points in the scale-location plot changes across the fitted values, and the Breusch–Pagan test is significant: BP(5) = 128.37, p < .001.
- Normality: the residuals are right-skewed (skewness ≈ 0.81).
- The cause: 90 of the 400 customers have a balance of exactly $0, and a straight-line model cannot reproduce that floor. It predicts a negative balance for 62 customers.
What this means:
- The coefficient estimates: heteroscedasticity by itself does not bias the OLS coefficient estimates, but classical standard errors, tests and confidence intervals may be unreliable. The p-values above should be read with caution.
- What the model is still good for: it is a reasonable description of which variables are associated with balance.
- What it should not be used for: precise inference or prediction intervals, unless the problems are addressed, for example with robust (HC3) standard errors or a model designed for an outcome with many zeros.
2.6 Selection 4: best subsets (exhaustive)
bs_c <- regsubsets(Balance ~ Income + Limit + Rating + Cards + Age + Education +
Gender + Student + Married + Ethnicity,
data = credit, nvmax = 10)
sm <- summary(bs_c)
smSubset selection object
Call: regsubsets.formula(Balance ~ Income + Limit + Rating + Cards +
Age + Education + Gender + Student + Married + Ethnicity,
data = credit, nvmax = 10)
11 Variables (and intercept)
Forced in Forced out
Income FALSE FALSE
Limit FALSE FALSE
Rating FALSE FALSE
Cards FALSE FALSE
Age FALSE FALSE
Education FALSE FALSE
GenderMale FALSE FALSE
StudentYes FALSE FALSE
MarriedYes FALSE FALSE
EthnicityAsian FALSE FALSE
EthnicityCaucasian FALSE FALSE
1 subsets of each size up to 10
Selection Algorithm: exhaustive
Income Limit Rating Cards Age Education GenderMale StudentYes
1 ( 1 ) " " " " "*" " " " " " " " " " "
2 ( 1 ) "*" " " "*" " " " " " " " " " "
3 ( 1 ) "*" " " "*" " " " " " " " " "*"
4 ( 1 ) "*" "*" " " "*" " " " " " " "*"
5 ( 1 ) "*" "*" "*" "*" " " " " " " "*"
6 ( 1 ) "*" "*" "*" "*" "*" " " " " "*"
7 ( 1 ) "*" "*" "*" "*" "*" " " "*" "*"
8 ( 1 ) "*" "*" "*" "*" "*" " " "*" "*"
9 ( 1 ) "*" "*" "*" "*" "*" " " "*" "*"
10 ( 1 ) "*" "*" "*" "*" "*" " " "*" "*"
MarriedYes EthnicityAsian EthnicityCaucasian
1 ( 1 ) " " " " " "
2 ( 1 ) " " " " " "
3 ( 1 ) " " " " " "
4 ( 1 ) " " " " " "
5 ( 1 ) " " " " " "
6 ( 1 ) " " " " " "
7 ( 1 ) " " " " " "
8 ( 1 ) " " "*" " "
9 ( 1 ) "*" "*" " "
10 ( 1 ) "*" "*" "*" # The size that optimizes each criterion:
c(BIC_best = which.min(sm$bic), Cp_best = which.min(sm$cp),
adjR2_best = which.max(sm$adjr2)) BIC_best Cp_best adjR2_best
4 6 7 # Fit statistics for the best model of each size
round(data.frame(size = 1:10, adjr2 = sm$adjr2, cp = sm$cp, bic = sm$bic), 3) size adjr2 cp bic
1 1 0.745 1800.308 -535.947
2 2 0.874 685.197 -814.180
3 3 0.949 41.134 -1173.358
4 4 0.953 11.149 -1198.053
5 5 0.954 8.132 -1197.096
6 6 0.954 5.575 -1195.732
7 7 0.954 6.462 -1190.879
8 8 0.954 7.846 -1185.519
9 9 0.954 9.192 -1180.199
10 10 0.954 10.473 -1174.948# Coefficients of the BIC-best model
coef(bs_c, which.min(sm$bic)) (Intercept) Income Limit Cards StudentYes
-499.7272117 -7.8392288 0.2666445 23.1753794 429.6064203 What the code does
regsubsets()(from leaps) searches every combination of predictor variables and keeps the best model of each size, where “best” means the lowest RSS. It uses an efficient algorithm, so it finds the best subset without literally fitting every model.- Each indicator column counts separately: Ethnicity contributes two columns. So “size” means the number of predictor columns, not the number of original variables.
nvmax = 10asks for the best models with up to 10 predictor columns.summary(bs_c)prints a table of asterisks showing which predictor variables are in the best model of each size.which.min()andwhich.max()find the model size that optimizes each criterion.- As with the full model, the predictor variables are written out so that the 0/1 indicator columns are not entered twice.
Interpretation
- BIC is lowest at 4 predictor variables: Income, Limit, Cards and Student. This model keeps Limit and drops the redundant Rating automatically.
- Mallows’ Cp is lowest at 6 and adjusted R² is highest at 7. The best 6-predictor subset is Income, Limit, Rating, Cards, Age and Student, the same as the stepwise model, so it keeps Rating.
- So it is BIC, which penalizes extra predictor variables the most, that drops Rating. Cp and adjusted R² do not.
- Adjusted R² is flat: it reaches .953 at 4 predictor variables and .954 from 5 onward, so predictor variables beyond 4 or 5 add almost nothing.
- Comparing the candidates:
- The 4-predictor BIC model is the leanest.
- The 5-predictor model reported to management (with Age) has every predictor variable significant and no collinearity. It is not the best-subsets winner of size 5: that subset is Income, Limit, Rating, Cards and Student, which keeps Rating. The analyst chose the model with Age, rather than one that keeps the redundant Rating.
- Both are defensible. They differ only in Age, whose effect is small.
Note on the slide wording. The slide says best subsets keeps only Limit and drops Rating. That is true of the BIC-best model only. The Cp-best and adjusted-R²-best models keep Rating.
2.7 Explain it to your boss
Credit limit, income and student status are the strongest predictor variables of card balance. Rating was dropped because it duplicates Limit (r = .997). That leaves a simpler model that predicts just as well (R² = .954).
- Limit is the strongest predictor variable. Each additional $1 of credit limit is associated with about $0.27 higher balance, holding the other predictor variables constant.
- Income and student status matter too. Holding the others constant, each additional $1,000 of income is associated with about $7.76 lower balance, and students carry balances about $428 higher, which identifies a clear segment for targeted products or risk rules.
- Why this model. It matches the full model’s accuracy with five predictor variables and no collinearity, so it is a sound basis for describing these relationships.
These are associations from observational data, not causal effects. For example, issuers set credit limits partly from a customer’s history, so a higher limit does not necessarily cause a higher balance.
2.8 Backward, forward and stepwise selection with p-values (olsrr)
The olsrr functions select predictor variables using
p-value thresholds instead of AIC.
full <- lm(
Balance ~ Income + Limit + Rating + Cards + Age +
Education + Female + Student_Yes + Married_Yes,
data = credit)
# Backward: remove predictors with p > .05
backward <- ols_step_backward_p(full, p_val = .05)
# Forward: add predictors with p < .05
forward <- ols_step_forward_p(full, p_val = .05)
# Stepwise: add with p < .05; remove with p > .05
stepwise <- ols_step_both_p(
full,
p_enter = .05,
p_remove = .05
)
# View selection results
backward
Stepwise Summary
-----------------------------------------------------------------------------
Step Variable AIC SBC SBIC R2 Adj. R2
-----------------------------------------------------------------------------
0 Full Model 4820.855 4864.761 3686.216 0.95493 0.95389
1 Married_Yes 4819.333 4859.248 3684.621 0.95488 0.95396
2 Education 4817.901 4853.824 3683.116 0.95482 0.95401
3 Female 4817.039 4848.971 3682.168 0.95469 0.95400
-----------------------------------------------------------------------------
Final Model Output
------------------
Model Summary
-----------------------------------------------------------------
R 0.977 RMSE 97.745
R-Squared 0.955 MSE 9554.049
Adj. R-Squared 0.954 Coef. Var 18.963
Pred R-Squared 0.953 AIC 4817.039
MAE 78.121 SBC 4848.971
-----------------------------------------------------------------
RMSE: Root Mean Square Error
MSE: Mean Square Error
MAE: Mean Absolute Error
AIC: Akaike Information Criteria
SBC: Schwarz Bayesian Criteria
ANOVA
----------------------------------------------------------------------------
Sum of
Squares DF Mean Square F Sig.
----------------------------------------------------------------------------
Regression 80518292.240 6 13419715.373 1380.03 0.0000
Residual 3821619.670 393 9724.223
Total 84339911.910 399
----------------------------------------------------------------------------
Parameter Estimates
-----------------------------------------------------------------------------------------------
model Beta Std. Error Std. Beta t Sig lower upper
-----------------------------------------------------------------------------------------------
(Intercept) -493.734 24.825 -19.889 0.000 -542.540 -444.928
Income -7.795 0.233 -0.598 -33.395 0.000 -8.254 -7.336
Limit 0.194 0.032 0.972 5.981 0.000 0.130 0.257
Rating 1.091 0.485 0.367 2.251 0.025 0.138 2.044
Cards 18.212 4.319 0.054 4.217 0.000 9.721 26.702
Age -0.624 0.292 -0.023 -2.139 0.033 -1.198 -0.050
Student_Yes 425.610 16.510 0.278 25.780 0.000 393.152 458.068
-----------------------------------------------------------------------------------------------forward
Stepwise Summary
-----------------------------------------------------------------------------
Step Variable AIC SBC SBIC R2 Adj. R2
-----------------------------------------------------------------------------
0 Base Model 6042.711 6050.694 4903.834 0.00000 0.00000
1 Rating 5496.782 5508.756 4357.019 0.74585 0.74521
2 Income 5214.557 5230.523 4074.833 0.87512 0.87449
3 Student_Yes 4851.387 4871.344 3715.602 0.94988 0.94950
4 Limit 4834.524 4858.473 3699.038 0.95219 0.95170
5 Cards 4819.667 4847.607 3684.616 0.95416 0.95358
6 Age 4817.039 4848.971 3682.168 0.95469 0.95400
-----------------------------------------------------------------------------
Final Model Output
------------------
Model Summary
-----------------------------------------------------------------
R 0.977 RMSE 97.745
R-Squared 0.955 MSE 9554.049
Adj. R-Squared 0.954 Coef. Var 18.963
Pred R-Squared 0.953 AIC 4817.039
MAE 78.121 SBC 4848.971
-----------------------------------------------------------------
RMSE: Root Mean Square Error
MSE: Mean Square Error
MAE: Mean Absolute Error
AIC: Akaike Information Criteria
SBC: Schwarz Bayesian Criteria
ANOVA
----------------------------------------------------------------------------
Sum of
Squares DF Mean Square F Sig.
----------------------------------------------------------------------------
Regression 80518292.240 6 13419715.373 1380.03 0.0000
Residual 3821619.670 393 9724.223
Total 84339911.910 399
----------------------------------------------------------------------------
Parameter Estimates
-----------------------------------------------------------------------------------------------
model Beta Std. Error Std. Beta t Sig lower upper
-----------------------------------------------------------------------------------------------
(Intercept) -493.734 24.825 -19.889 0.000 -542.540 -444.928
Rating 1.091 0.485 0.367 2.251 0.025 0.138 2.044
Income -7.795 0.233 -0.598 -33.395 0.000 -8.254 -7.336
Student_Yes 425.610 16.510 0.278 25.780 0.000 393.152 458.068
Limit 0.194 0.032 0.972 5.981 0.000 0.130 0.257
Cards 18.212 4.319 0.054 4.217 0.000 9.721 26.702
Age -0.624 0.292 -0.023 -2.139 0.033 -1.198 -0.050
-----------------------------------------------------------------------------------------------stepwise
Stepwise Summary
---------------------------------------------------------------------------------
Step Variable AIC SBC SBIC R2 Adj. R2
---------------------------------------------------------------------------------
0 Base Model 6042.711 6050.694 4903.834 0.00000 0.00000
1 Rating (+) 5496.782 5508.756 4357.019 0.74585 0.74521
2 Income (+) 5214.557 5230.523 4074.833 0.87512 0.87449
3 Student_Yes (+) 4851.387 4871.344 3715.602 0.94988 0.94950
4 Limit (+) 4834.524 4858.473 3699.038 0.95219 0.95170
5 Cards (+) 4819.667 4847.607 3684.616 0.95416 0.95358
6 Age (+) 4817.039 4848.971 3682.168 0.95469 0.95400
---------------------------------------------------------------------------------
Final Model Output
------------------
Model Summary
-----------------------------------------------------------------
R 0.977 RMSE 97.745
R-Squared 0.955 MSE 9554.049
Adj. R-Squared 0.954 Coef. Var 18.963
Pred R-Squared 0.953 AIC 4817.039
MAE 78.121 SBC 4848.971
-----------------------------------------------------------------
RMSE: Root Mean Square Error
MSE: Mean Square Error
MAE: Mean Absolute Error
AIC: Akaike Information Criteria
SBC: Schwarz Bayesian Criteria
ANOVA
----------------------------------------------------------------------------
Sum of
Squares DF Mean Square F Sig.
----------------------------------------------------------------------------
Regression 80518292.240 6 13419715.373 1380.03 0.0000
Residual 3821619.670 393 9724.223
Total 84339911.910 399
----------------------------------------------------------------------------
Parameter Estimates
-----------------------------------------------------------------------------------------------
model Beta Std. Error Std. Beta t Sig lower upper
-----------------------------------------------------------------------------------------------
(Intercept) -493.734 24.825 -19.889 0.000 -542.540 -444.928
Rating 1.091 0.485 0.367 2.251 0.025 0.138 2.044
Income -7.795 0.233 -0.598 -33.395 0.000 -8.254 -7.336
Student_Yes 425.610 16.510 0.278 25.780 0.000 393.152 458.068
Limit 0.194 0.032 0.972 5.981 0.000 0.130 0.257
Cards 18.212 4.319 0.054 4.217 0.000 9.721 26.702
Age -0.624 0.292 -0.023 -2.139 0.033 -1.198 -0.050
-----------------------------------------------------------------------------------------------What the code does
fullis redefined here with the 0/1 indicators created earlier (Female,Student_Yes,Married_Yes) in place of the factors. Ethnicity is not included, so it is not a candidate.ols_step_backward_p()starts from the full model and repeatedly removes the predictor variable with the largest p-value, as long as that p-value is abovep_val(.05).ols_step_forward_p()starts from the intercept-only model and repeatedly adds the predictor variable with the smallest p-value, as long as it is belowp_val.ols_step_both_p()adds predictor variables with p <p_enterand, after each addition, removes any whose p-value has risen abovep_remove.- Why the thresholds are set explicitly: the olsrr
defaults are lenient (
p_val = 0.3;p_enter = 0.1,p_remove = 0.3). Without these arguments, the forward and backward searches (p_val = 0.3) also keep Female, giving seven predictor variables. The stepwise search with its defaults happens to give the same six as here. - Printing each object shows a step-by-step summary, then the final model’s fit statistics, ANOVA table and parameter estimates.
Interpretation
| Method | Path |
|---|---|
| Backward | Removed Married_Yes, then Education, then Female |
| Forward | Added Rating, Income, Student_Yes, Limit, Cards, then Age |
| Stepwise | The same additions as forward; nothing was removed |
- All three end at the same six predictor variables as the AIC searches: Income, Limit, Rating, Cards, Age and Student_Yes, with R² = .955 and adjusted R² = .954.
- Each keeps both Limit and Rating. Each is significant given the other (Limit p < .001, Rating p = .025), so a p-value search cannot detect the redundancy.
- The VIFs (about 230 each) are the tool that catches it.
- Student_Yes: b = 425.61. Holding the other predictor variables constant, students carry balances about $426 higher.
3 Case Study 2: Finance
3.1 Goal and the variables
We build several regression models to predict MR (market return) from 12 candidate predictor variables, and decide which model to report.
| Variable | Name | What it measures |
|---|---|---|
| MR | Market return | Return over the period. Outcome variable. |
| Beta | Market beta | Sensitivity of returns to overall market movements |
| P_E | Price-to-earnings ratio | Share price ÷ earnings per share |
| P_B | Price-to-book ratio | Market value ÷ book value of equity |
| DY | Dividend yield | Annual dividends ÷ share price |
| EG | Earnings growth | Rate at which earnings are growing |
| D_E | Debt-to-equity ratio | Total debt ÷ shareholders’ equity (leverage) |
| ROE | Return on equity | Net income ÷ shareholders’ equity (profitability) |
| Size | Firm size | Market capitalization |
| MOM | Momentum | Strength of the recent price trend |
| IR | Interest rate | Prevailing level of interest rates |
| VIX | Volatility index | Expected stock-market volatility |
| SR | Sharpe ratio | Return above the risk-free rate per unit of risk |
Notes on the data
- The definitions are assumed. Finance.csv has no codebook, so these are the standard finance meanings of the abbreviations. SR is taken to be the Sharpe ratio.
- The values are scaled. Every variable is centered near 0, so the numbers are not raw percentages or dollars, and 0 means roughly average. Coefficients are read as the change in MR for a one-unit change in the scaled predictor variable.
- No log transformation. Many values are negative, so a log transformation is not possible for these variables.
3.2 Step 0: Correlation matrix
finance <- read.csv("Finance.csv")
str(finance)'data.frame': 500 obs. of 13 variables:
$ MR : num 0.59 -1.01 -0.44 1.3 -0.59 -0.13 0.04 1.01 -0.94 -1.68 ...
$ Beta: num 0.49 -2.38 -2.5 1.54 1.33 -0.49 -0.38 0.53 -0.98 -1.95 ...
$ P_E : num -0.93 -2.1 0.48 3.16 1.03 -0.38 0.16 -0.35 -1.25 2.66 ...
$ P_B : num 0.7 -0.73 0.96 1.6 2.19 3.38 -1.41 -0.12 1.36 0.25 ...
$ DY : num 0.72 0.63 0.72 -0.34 0.04 0.87 0.53 -0.13 1.54 -0.16 ...
$ EG : num 1.02 -0.13 -0.05 1.06 0.53 -0.53 0.32 0.9 0.32 1.06 ...
$ D_E : num 0.12 0.14 -0.09 -0.6 -0.06 -0.02 -0.48 -0.1 -0.97 -0.05 ...
$ ROE : num -0.26 0.24 1.69 0.61 0.22 -0.77 0.38 2.51 -0.53 0.08 ...
$ Size: num -1.06 -2.36 -1.51 0.64 0.22 0.14 -0.3 0.46 0.43 1.24 ...
$ MOM : num 2.71 -1.54 1.36 2.57 0.53 -2.06 -0.67 1.22 -0.24 0.06 ...
$ IR : num 0.46 -0.46 -0.06 -0.66 -0.72 1.03 -0.72 0.32 -0.58 0 ...
$ VIX : num -0.78 2.19 -0.83 -1.36 -1.71 0.62 -2.41 -0.14 -2.6 -0.88 ...
$ SR : num 1.55 -2.16 8.7 3.18 0.61 -0.21 3.88 -3.43 -6.22 0.54 ...cor_matrix <- cor(finance) # all 13 variables are numeric
round(cor_matrix, 2) MR Beta P_E P_B DY EG D_E ROE Size MOM IR VIX
MR 1.00 0.63 0.18 0.20 0.13 0.16 -0.10 0.31 0.50 0.59 -0.27 -0.50
Beta 0.63 1.00 0.06 0.08 -0.06 0.13 0.06 0.18 0.36 0.53 -0.10 -0.39
P_E 0.18 0.06 1.00 0.51 -0.13 0.32 0.07 0.42 0.16 0.06 -0.12 -0.16
P_B 0.20 0.08 0.51 1.00 -0.03 0.27 0.02 0.36 0.24 0.10 -0.10 -0.14
DY 0.13 -0.06 -0.13 -0.03 1.00 -0.12 -0.26 -0.08 -0.12 -0.09 -0.16 -0.25
EG 0.16 0.13 0.32 0.27 -0.12 1.00 -0.04 0.31 0.12 0.20 -0.13 -0.27
D_E -0.10 0.06 0.07 0.02 -0.26 -0.04 1.00 -0.11 -0.10 -0.05 0.12 0.13
ROE 0.31 0.18 0.42 0.36 -0.08 0.31 -0.11 1.00 0.21 0.27 -0.10 -0.24
Size 0.50 0.36 0.16 0.24 -0.12 0.12 -0.10 0.21 1.00 0.36 -0.26 -0.39
MOM 0.59 0.53 0.06 0.10 -0.09 0.20 -0.05 0.27 0.36 1.00 -0.17 -0.28
IR -0.27 -0.10 -0.12 -0.10 -0.16 -0.13 0.12 -0.10 -0.26 -0.17 1.00 0.50
VIX -0.50 -0.39 -0.16 -0.14 -0.25 -0.27 0.13 -0.24 -0.39 -0.28 0.50 1.00
SR 0.26 0.32 0.45 0.21 0.03 0.24 -0.09 0.42 0.12 0.28 -0.05 -0.24
SR
MR 0.26
Beta 0.32
P_E 0.45
P_B 0.21
DY 0.03
EG 0.24
D_E -0.09
ROE 0.42
Size 0.12
MOM 0.28
IR -0.05
VIX -0.24
SR 1.00What the code does
read.csv()reads the 500 rows and 13 variables.str()confirms that every variable is numeric and that there are no text columns.cor()can therefore use the whole data frame, with no conversion to factors or indicators needed.- The matrix is symmetric (r for Beta and MOM equals r for MOM and Beta), and the diagonal is 1.00.
# Correlation of each predictor variable with MR, largest to smallest
# (the -1 drops MR's correlation with itself)
r_with_MR <- sort(cor_matrix["MR", -1], decreasing = TRUE)
round(r_with_MR, 2) Beta MOM Size ROE SR P_B P_E EG DY D_E IR VIX
0.63 0.59 0.50 0.31 0.26 0.20 0.18 0.16 0.13 -0.10 -0.27 -0.50 # Bar chart: blue = positive, orange = negative
barplot(rev(r_with_MR), horiz = TRUE, las = 1,
col = ifelse(rev(r_with_MR) > 0, "steelblue", "darkorange"),
xlim = c(-0.6, 0.8),
main = "Correlation with MR (Pearson r)", xlab = "r")
abline(v = 0)
# Largest correlations between pairs of predictor variables
pred_cor <- cor_matrix[-1, -1] # drop MR
pred_cor[upper.tri(pred_cor, diag = TRUE)] <- NA # keep each pair once
pairs <- na.omit(as.data.frame(as.table(pred_cor)))
names(pairs) <- c("Var1", "Var2", "r")
head(pairs[order(-abs(pairs$r)), ], 3) Var1 Var2 r
9 MOM Beta 0.5274981
15 P_B P_E 0.5131440
119 VIX IR 0.5017637What the code does
cor_matrix["MR", -1]takes the MR row of the matrix, without MR’s correlation with itself.sort(..., decreasing = TRUE)orders the correlations from largest to smallest.barplot(horiz = TRUE)draws horizontal bars.rev()puts the largest bar at the top,ifelse()colours bars by sign, andabline(v = 0)draws the zero line.- The last block finds the strongest predictor pairs:
- It removes MR from the matrix.
- It blanks out the upper triangle, so each pair appears only once.
- It reshapes the matrix into a list of pairs.
- It shows the three pairs with the largest |r|.
Interpretation
- Strongest associations with MR: Beta (r = .63), MOM (.59) and Size (.50) are positive; VIX (−.50) is negative.
- Weakest: D_E (−.10) and DY (.13).
- Largest pairs of predictor variables: MOM and Beta
(.53), P_E and P_B (.51), and VIX and IR (.50).
- These are moderate. None is anywhere near the .997 between Limit and Rating in the Credit data, so serious multicollinearity is unlikely here.
- Caution: r is a two-variable association. A predictor variable’s role can change once the others are controlled; we will see this for DY and SR in the final model.
3.3 Example 1: Hierarchical (blockwise) entry
The problem. An investment analyst enters predictor variables for MR in three blocks whose order comes from financial theory, before looking at the results:
- Block 1, risk: Beta, VIX.
- Block 2, market factors: MOM, Size.
- Block 3, growth, leverage and rates: EG, D_E, IR.
The question is not which single predictor variable matters most. It is whether each block as a whole adds explanatory power beyond the blocks already entered.
# Block 1 - risk baseline
m1 <- lm(MR ~ Beta + VIX, data = finance)
summary(m1)
Call:
lm(formula = MR ~ Beta + VIX, data = finance)
Residuals:
Min 1Q Median 3Q Max
-2.17774 -0.51482 0.02005 0.53286 2.31606
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.009364 0.036052 -0.26 0.795
Beta 0.286107 0.019457 14.71 < 2e-16 ***
VIX -0.165589 0.020047 -8.26 1.33e-15 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8014 on 497 degrees of freedom
Multiple R-squared: 0.4744, Adjusted R-squared: 0.4723
F-statistic: 224.3 on 2 and 497 DF, p-value: < 2.2e-16# Block 2 - add market factors
m2 <- lm(MR ~ Beta + VIX + MOM + Size, data = finance)
summary(m2)
Call:
lm(formula = MR ~ Beta + VIX + MOM + Size, data = finance)
Residuals:
Min 1Q Median 3Q Max
-1.92689 -0.45709 -0.00583 0.45594 2.22747
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.02350 0.03221 -0.730 0.466
Beta 0.17982 0.01989 9.040 < 2e-16 ***
VIX -0.11857 0.01867 -6.352 4.81e-10 ***
MOM 0.20577 0.02479 8.301 9.93e-16 ***
Size 0.16766 0.02788 6.013 3.54e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7152 on 495 degrees of freedom
Multiple R-squared: 0.5831, Adjusted R-squared: 0.5798
F-statistic: 173.1 on 4 and 495 DF, p-value: < 2.2e-16# Block 3 - add growth, leverage and interest rates
m3 <- lm(MR ~ Beta + VIX + MOM + Size +
EG + D_E + IR, data = finance)
summary(m3)
Call:
lm(formula = MR ~ Beta + VIX + MOM + Size + EG + D_E + IR, data = finance)
Residuals:
Min 1Q Median 3Q Max
-1.90055 -0.46439 0.01991 0.45931 2.30690
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.02342 0.03216 -0.728 0.4667
Beta 0.18910 0.02034 9.297 < 2e-16 ***
VIX -0.10407 0.02163 -4.811 2.00e-06 ***
MOM 0.20233 0.02509 8.063 5.69e-15 ***
Size 0.16067 0.02800 5.738 1.67e-08 ***
EG -0.02872 0.04681 -0.614 0.5398
D_E -0.10130 0.05529 -1.832 0.0675 .
IR -0.08777 0.07050 -1.245 0.2137
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7134 on 492 degrees of freedom
Multiple R-squared: 0.5877, Adjusted R-squared: 0.5819
F-statistic: 100.2 on 7 and 492 DF, p-value: < 2.2e-16What the code does
- The three models are nested: each contains every
predictor variable in the model before it.
- m2 repeats Beta and VIX from m1 and adds MOM and Size.
- m3 repeats all four and adds EG, D_E and IR.
- Leaving out an earlier predictor variable would break the nesting, and the comparison below would no longer be valid.
# Compare nested models
anova(m1, m2, m3)Analysis of Variance Table
Model 1: MR ~ Beta + VIX
Model 2: MR ~ Beta + VIX + MOM + Size
Model 3: MR ~ Beta + VIX + MOM + Size + EG + D_E + IR
Res.Df RSS Df Sum of Sq F Pr(>F)
1 497 319.18
2 495 253.17 2 66.011 64.8570 <2e-16 ***
3 492 250.38 3 2.795 1.8304 0.1407
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# R-squared change test for each block, one comparison at a time
anova(m1, m2)Analysis of Variance Table
Model 1: MR ~ Beta + VIX
Model 2: MR ~ Beta + VIX + MOM + Size
Res.Df RSS Df Sum of Sq F Pr(>F)
1 497 319.18
2 495 253.17 2 66.011 64.532 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m2, m3)Analysis of Variance Table
Model 1: MR ~ Beta + VIX + MOM + Size
Model 2: MR ~ Beta + VIX + MOM + Size + EG + D_E + IR
Res.Df RSS Df Sum of Sq F Pr(>F)
1 495 253.17
2 492 250.38 3 2.7945 1.8304 0.1407# Compare model fit metrics
list(m1 = m1, m2 = m2, m3 = m3) |>
map_dfr(glance, .id = "model") |>
dplyr::select(model, r.squared, adj.r.squared, AIC, BIC)# A tibble: 3 × 5
model r.squared adj.r.squared AIC BIC
<chr> <dbl> <dbl> <dbl> <dbl>
1 m1 0.474 0.472 1203. 1219.
2 m2 0.583 0.580 1091. 1116.
3 m3 0.588 0.582 1091. 1129.What the code does
anova(m1, m2, m3)runs an F test for each added block, comparing m2 with m1 and then m3 with m2.- The denominator differs between the two calls. When
anova()is given three models, it tests every row against the error term of the largest model, m3 (492 df). So its Block 2 row, F = 64.86, is F(2, 492). - The standard R² change test compares each block
with the model just before it.
anova(m1, m2)gives F(2, 495) = 64.53 for Block 2, andanova(m2, m3)gives F(3, 492) = 1.83 for Block 3. These are the values to report. The conclusions are the same either way. - This F test asks whether the new block explains significantly more variance (the R² change, ΔR²).
- The last command builds a table of fit statistics:
list()collects the three models.map_dfr(glance, .id = "model")runsglance()on each model and stacks the rows, labelling each row with the model’s name.dplyr::select()keeps the columns we want.|>is R’s built-in pipe (R 4.1 or later).
Interpretation
| Question | Answer |
|---|---|
| Variance explained by Block 1 alone | R² = .474. Beta and VIX together explain about 47% of the variance in MR. |
| Does Block 2 improve the model? | Yes. ΔR² = .583 − .474 = .109; F(2, 495) = 64.53, p < .001. Momentum and size explain about 11% more variance beyond risk. |
| Does Block 3 improve the model? | No. ΔR² = .588 − .583 = .005; F(3, 492) = 1.83, p = .141. |
| Which model to choose? | m2. Block 3 is not significant, and both AIC (1090.7 → 1091.1) and BIC (1116.0 → 1129.1) rise. |
- Degrees of freedom in the F tests: the numerator df is the number of predictor variables added (2, then 3). For the block-by-block tests, the denominator df is the residual df of the larger model in that comparison (495, then 492).
- Why the adjusted R² rise doesn’t count: adjusted R² rises slightly from .580 to .582, a gain of only .002. Adjusted R² penalizes added predictor variables less heavily than AIC or BIC do, so this tiny gain does not justify three more predictor variables.
Business conclusion. Risk measures explain about half of the variation in market returns. Adding momentum and firm size explains about 11% more. Earnings growth, leverage and interest rates add no significant explanatory power beyond these four, so the four-predictor model is preferred.
Why Example 3 ends with a different model: Example 1 considers only the seven predictor variables in its three blocks. Example 2 searches all 12 candidates, which is why Example 3 reports an eight-predictor model.
3.4 Example 2: Stepwise procedures and best subsets
The problem. A portfolio analyst wants a parsimonious model for MR: as few predictor variables as possible without losing explanatory power, chosen from all 12 candidates.
m_full <- lm(MR ~ ., data = finance)
m_null <- lm(MR ~ 1, data = finance)
# Forward - AIC-based
m_fwd <- step(m_null, direction = "forward",
scope = list(lower = m_null,
upper = m_full), trace = 0)
# Backward - AIC-based
m_bwd <- step(m_full, direction = "backward",
trace = 0)
# Stepwise (both directions) - AIC-based
m_step <- step(m_null, direction = "both",
scope = list(lower = m_null,
upper = m_full), trace = 0)
formula(m_fwd); formula(m_bwd); formula(m_step)MR ~ Beta + MOM + VIX + Size + DY + ROE + P_E + SRMR ~ Beta + P_E + DY + ROE + Size + MOM + VIX + SRMR ~ Beta + MOM + VIX + Size + DY + ROE + P_E + SRWhat the code does
MR ~ .uses every other column as a predictor variable. That is safe here because Finance contains only the outcome variable and the 12 candidates.step()is base R’s AIC search. It gives the same results asstepAIC().trace = 0hides the step-by-step tables.formula()prints the final model of each search.
Interpretation. All three AIC searches select the same eight predictor variables: Beta, P_E, DY, ROE, Size, MOM, VIX and SR. Forward selection adds them in the order Beta, MOM, VIX, Size, DY, ROE, P_E, SR. Beta enters first and SR last.
# p-value-based versions
ols_step_forward_p(m_full, p_val = 0.05)
Stepwise Summary
----------------------------------------------------------------------------
Step Variable AIC SBC SBIC R2 Adj. R2
----------------------------------------------------------------------------
0 Base Model 1520.162 1528.591 99.200 0.00000 0.00000
1 Beta 1264.840 1277.484 -155.850 0.40229 0.40109
2 MOM 1182.650 1199.508 -237.942 0.49491 0.49288
3 VIX 1123.921 1144.994 -296.329 0.55268 0.54997
4 Size 1090.671 1115.959 -329.247 0.58313 0.57976
5 DY 1060.187 1089.689 -359.180 0.60935 0.60540
6 ROE 1046.862 1080.579 -372.153 0.62114 0.61653
7 P_E 1042.409 1080.340 -376.411 0.62600 0.62068
8 SR 1037.158 1079.304 -381.383 0.63139 0.62538
----------------------------------------------------------------------------
Final Model Output
------------------
Model Summary
-----------------------------------------------------------------
R 0.795 RMSE 0.669
R-Squared 0.631 MSE 0.448
Adj. R-Squared 0.625 Coef. Var -1189.208
Pred R-Squared 0.618 AIC 1037.158
MAE 0.522 SBC 1079.304
-----------------------------------------------------------------
RMSE: Root Mean Square Error
MSE: Mean Square Error
MAE: Mean Absolute Error
AIC: Akaike Information Criteria
SBC: Schwarz Bayesian Criteria
ANOVA
---------------------------------------------------------------------
Sum of
Squares DF Mean Square F Sig.
---------------------------------------------------------------------
Regression 383.453 8 47.932 105.127 0.0000
Residual 223.866 491 0.456
Total 607.319 499
---------------------------------------------------------------------
Parameter Estimates
----------------------------------------------------------------------------------------
model Beta Std. Error Std. Beta t Sig lower upper
----------------------------------------------------------------------------------------
(Intercept) -0.036 0.030 -1.196 0.232 -0.096 0.023
Beta 0.206 0.019 0.374 10.607 0.000 0.168 0.244
MOM 0.211 0.024 0.297 8.757 0.000 0.163 0.258
VIX -0.063 0.019 -0.111 -3.292 0.001 -0.101 -0.025
Size 0.179 0.027 0.212 6.611 0.000 0.126 0.232
DY 0.295 0.043 0.205 6.840 0.000 0.210 0.379
ROE 0.064 0.020 0.104 3.184 0.002 0.024 0.103
P_E 0.075 0.022 0.111 3.357 0.001 0.031 0.118
SR -0.026 0.010 -0.091 -2.678 0.008 -0.044 -0.007
----------------------------------------------------------------------------------------ols_step_backward_p(m_full, p_val = 0.10)
Stepwise Summary
----------------------------------------------------------------------------
Step Variable AIC SBC SBIC R2 Adj. R2
----------------------------------------------------------------------------
0 Full Model 1042.935 1101.939 -375.311 0.63302 0.62398
1 P_B 1041.035 1095.825 -377.270 0.63295 0.62468
2 D_E 1039.592 1090.167 -378.791 0.63254 0.62503
3 IR 1038.353 1084.714 -380.110 0.63198 0.62522
4 EG 1037.158 1079.304 -381.383 0.63139 0.62538
----------------------------------------------------------------------------
Final Model Output
------------------
Model Summary
-----------------------------------------------------------------
R 0.795 RMSE 0.669
R-Squared 0.631 MSE 0.448
Adj. R-Squared 0.625 Coef. Var -1189.208
Pred R-Squared 0.618 AIC 1037.158
MAE 0.522 SBC 1079.304
-----------------------------------------------------------------
RMSE: Root Mean Square Error
MSE: Mean Square Error
MAE: Mean Absolute Error
AIC: Akaike Information Criteria
SBC: Schwarz Bayesian Criteria
ANOVA
---------------------------------------------------------------------
Sum of
Squares DF Mean Square F Sig.
---------------------------------------------------------------------
Regression 383.453 8 47.932 105.127 0.0000
Residual 223.866 491 0.456
Total 607.319 499
---------------------------------------------------------------------
Parameter Estimates
----------------------------------------------------------------------------------------
model Beta Std. Error Std. Beta t Sig lower upper
----------------------------------------------------------------------------------------
(Intercept) -0.036 0.030 -1.196 0.232 -0.096 0.023
Beta 0.206 0.019 0.374 10.607 0.000 0.168 0.244
P_E 0.075 0.022 0.111 3.357 0.001 0.031 0.118
DY 0.295 0.043 0.205 6.840 0.000 0.210 0.379
ROE 0.064 0.020 0.104 3.184 0.002 0.024 0.103
Size 0.179 0.027 0.212 6.611 0.000 0.126 0.232
MOM 0.211 0.024 0.297 8.757 0.000 0.163 0.258
VIX -0.063 0.019 -0.111 -3.292 0.001 -0.101 -0.025
SR -0.026 0.010 -0.091 -2.678 0.008 -0.044 -0.007
----------------------------------------------------------------------------------------ols_step_both_p(m_full,
p_enter = 0.05, p_remove = 0.10)
Stepwise Summary
----------------------------------------------------------------------------
Step Variable AIC SBC SBIC R2 Adj. R2
----------------------------------------------------------------------------
0 Base Model 1520.162 1528.591 99.200 0.00000 0.00000
1 Beta (+) 1264.840 1277.484 -155.850 0.40229 0.40109
2 MOM (+) 1182.650 1199.508 -237.942 0.49491 0.49288
3 VIX (+) 1123.921 1144.994 -296.329 0.55268 0.54997
4 Size (+) 1090.671 1115.959 -329.247 0.58313 0.57976
5 DY (+) 1060.187 1089.689 -359.180 0.60935 0.60540
6 ROE (+) 1046.862 1080.579 -372.153 0.62114 0.61653
7 P_E (+) 1042.409 1080.340 -376.411 0.62600 0.62068
8 SR (+) 1037.158 1079.304 -381.383 0.63139 0.62538
----------------------------------------------------------------------------
Final Model Output
------------------
Model Summary
-----------------------------------------------------------------
R 0.795 RMSE 0.669
R-Squared 0.631 MSE 0.448
Adj. R-Squared 0.625 Coef. Var -1189.208
Pred R-Squared 0.618 AIC 1037.158
MAE 0.522 SBC 1079.304
-----------------------------------------------------------------
RMSE: Root Mean Square Error
MSE: Mean Square Error
MAE: Mean Absolute Error
AIC: Akaike Information Criteria
SBC: Schwarz Bayesian Criteria
ANOVA
---------------------------------------------------------------------
Sum of
Squares DF Mean Square F Sig.
---------------------------------------------------------------------
Regression 383.453 8 47.932 105.127 0.0000
Residual 223.866 491 0.456
Total 607.319 499
---------------------------------------------------------------------
Parameter Estimates
----------------------------------------------------------------------------------------
model Beta Std. Error Std. Beta t Sig lower upper
----------------------------------------------------------------------------------------
(Intercept) -0.036 0.030 -1.196 0.232 -0.096 0.023
Beta 0.206 0.019 0.374 10.607 0.000 0.168 0.244
MOM 0.211 0.024 0.297 8.757 0.000 0.163 0.258
VIX -0.063 0.019 -0.111 -3.292 0.001 -0.101 -0.025
Size 0.179 0.027 0.212 6.611 0.000 0.126 0.232
DY 0.295 0.043 0.205 6.840 0.000 0.210 0.379
ROE 0.064 0.020 0.104 3.184 0.002 0.024 0.103
P_E 0.075 0.022 0.111 3.357 0.001 0.031 0.118
SR -0.026 0.010 -0.091 -2.678 0.008 -0.044 -0.007
----------------------------------------------------------------------------------------# AIC-based with reporting
ols_step_forward_aic(m_full)
Stepwise Summary
----------------------------------------------------------------------------
Step Variable AIC SBC SBIC R2 Adj. R2
----------------------------------------------------------------------------
0 Base Model 1520.162 1528.591 99.200 0.00000 0.00000
1 Beta 1264.840 1277.484 -155.850 0.40229 0.40109
2 MOM 1182.650 1199.508 -237.942 0.49491 0.49288
3 VIX 1123.921 1144.994 -296.329 0.55268 0.54997
4 Size 1090.671 1115.959 -329.247 0.58313 0.57976
5 DY 1060.187 1089.689 -359.180 0.60935 0.60540
6 ROE 1046.862 1080.579 -372.153 0.62114 0.61653
7 P_E 1042.409 1080.340 -376.411 0.62600 0.62068
8 SR 1037.158 1079.304 -381.383 0.63139 0.62538
----------------------------------------------------------------------------
Final Model Output
------------------
Model Summary
-----------------------------------------------------------------
R 0.795 RMSE 0.669
R-Squared 0.631 MSE 0.448
Adj. R-Squared 0.625 Coef. Var -1189.208
Pred R-Squared 0.618 AIC 1037.158
MAE 0.522 SBC 1079.304
-----------------------------------------------------------------
RMSE: Root Mean Square Error
MSE: Mean Square Error
MAE: Mean Absolute Error
AIC: Akaike Information Criteria
SBC: Schwarz Bayesian Criteria
ANOVA
---------------------------------------------------------------------
Sum of
Squares DF Mean Square F Sig.
---------------------------------------------------------------------
Regression 383.453 8 47.932 105.127 0.0000
Residual 223.866 491 0.456
Total 607.319 499
---------------------------------------------------------------------
Parameter Estimates
----------------------------------------------------------------------------------------
model Beta Std. Error Std. Beta t Sig lower upper
----------------------------------------------------------------------------------------
(Intercept) -0.036 0.030 -1.196 0.232 -0.096 0.023
Beta 0.206 0.019 0.374 10.607 0.000 0.168 0.244
MOM 0.211 0.024 0.297 8.757 0.000 0.163 0.258
VIX -0.063 0.019 -0.111 -3.292 0.001 -0.101 -0.025
Size 0.179 0.027 0.212 6.611 0.000 0.126 0.232
DY 0.295 0.043 0.205 6.840 0.000 0.210 0.379
ROE 0.064 0.020 0.104 3.184 0.002 0.024 0.103
P_E 0.075 0.022 0.111 3.357 0.001 0.031 0.118
SR -0.026 0.010 -0.091 -2.678 0.008 -0.044 -0.007
----------------------------------------------------------------------------------------ols_step_backward_aic(m_full)
Stepwise Summary
----------------------------------------------------------------------------
Step Variable AIC SBC SBIC R2 Adj. R2
----------------------------------------------------------------------------
0 Full Model 1042.935 1101.939 -375.311 0.63302 0.62398
1 P_B 1041.035 1095.825 -377.315 0.63295 0.62468
2 D_E 1039.592 1090.167 -378.853 0.63254 0.62503
3 IR 1038.353 1084.714 -380.178 0.63198 0.62522
4 EG 1037.158 1079.304 -381.451 0.63139 0.62538
----------------------------------------------------------------------------
Final Model Output
------------------
Model Summary
-----------------------------------------------------------------
R 0.795 RMSE 0.669
R-Squared 0.631 MSE 0.448
Adj. R-Squared 0.625 Coef. Var -1189.208
Pred R-Squared 0.618 AIC 1037.158
MAE 0.522 SBC 1079.304
-----------------------------------------------------------------
RMSE: Root Mean Square Error
MSE: Mean Square Error
MAE: Mean Absolute Error
AIC: Akaike Information Criteria
SBC: Schwarz Bayesian Criteria
ANOVA
---------------------------------------------------------------------
Sum of
Squares DF Mean Square F Sig.
---------------------------------------------------------------------
Regression 383.453 8 47.932 105.127 0.0000
Residual 223.866 491 0.456
Total 607.319 499
---------------------------------------------------------------------
Parameter Estimates
----------------------------------------------------------------------------------------
model Beta Std. Error Std. Beta t Sig lower upper
----------------------------------------------------------------------------------------
(Intercept) -0.036 0.030 -1.196 0.232 -0.096 0.023
Beta 0.206 0.019 0.374 10.607 0.000 0.168 0.244
P_E 0.075 0.022 0.111 3.357 0.001 0.031 0.118
DY 0.295 0.043 0.205 6.840 0.000 0.210 0.379
ROE 0.064 0.020 0.104 3.184 0.002 0.024 0.103
Size 0.179 0.027 0.212 6.611 0.000 0.126 0.232
MOM 0.211 0.024 0.297 8.757 0.000 0.163 0.258
VIX -0.063 0.019 -0.111 -3.292 0.001 -0.101 -0.025
SR -0.026 0.010 -0.091 -2.678 0.008 -0.044 -0.007
----------------------------------------------------------------------------------------ols_step_both_aic(m_full)
Stepwise Summary
----------------------------------------------------------------------------
Step Variable AIC SBC SBIC R2 Adj. R2
----------------------------------------------------------------------------
0 Base Model 1520.162 1528.591 99.200 0.00000 0.00000
1 Beta (+) 1264.840 1277.484 -155.850 0.40229 0.40109
2 MOM (+) 1182.650 1199.508 -237.942 0.49491 0.49288
3 VIX (+) 1123.921 1144.994 -296.329 0.55268 0.54997
4 Size (+) 1090.671 1115.959 -329.247 0.58313 0.57976
5 DY (+) 1060.187 1089.689 -359.180 0.60935 0.60540
6 ROE (+) 1046.862 1080.579 -372.153 0.62114 0.61653
7 P_E (+) 1042.409 1080.340 -376.411 0.62600 0.62068
8 SR (+) 1037.158 1079.304 -381.383 0.63139 0.62538
----------------------------------------------------------------------------
Final Model Output
------------------
Model Summary
-----------------------------------------------------------------
R 0.795 RMSE 0.669
R-Squared 0.631 MSE 0.448
Adj. R-Squared 0.625 Coef. Var -1189.208
Pred R-Squared 0.618 AIC 1037.158
MAE 0.522 SBC 1079.304
-----------------------------------------------------------------
RMSE: Root Mean Square Error
MSE: Mean Square Error
MAE: Mean Absolute Error
AIC: Akaike Information Criteria
SBC: Schwarz Bayesian Criteria
ANOVA
---------------------------------------------------------------------
Sum of
Squares DF Mean Square F Sig.
---------------------------------------------------------------------
Regression 383.453 8 47.932 105.127 0.0000
Residual 223.866 491 0.456
Total 607.319 499
---------------------------------------------------------------------
Parameter Estimates
----------------------------------------------------------------------------------------
model Beta Std. Error Std. Beta t Sig lower upper
----------------------------------------------------------------------------------------
(Intercept) -0.036 0.030 -1.196 0.232 -0.096 0.023
Beta 0.206 0.019 0.374 10.607 0.000 0.168 0.244
MOM 0.211 0.024 0.297 8.757 0.000 0.163 0.258
VIX -0.063 0.019 -0.111 -3.292 0.001 -0.101 -0.025
Size 0.179 0.027 0.212 6.611 0.000 0.126 0.232
DY 0.295 0.043 0.205 6.840 0.000 0.210 0.379
ROE 0.064 0.020 0.104 3.184 0.002 0.024 0.103
P_E 0.075 0.022 0.111 3.357 0.001 0.031 0.118
SR -0.026 0.010 -0.091 -2.678 0.008 -0.044 -0.007
----------------------------------------------------------------------------------------What the code does
- The p-value thresholds: .05 to enter and .10 to remove.
- Why the removal threshold is larger: a larger
p_removekeeps a borderline predictor variable from being dropped as soon as it enters. Ifp_removewere smaller thanp_enter, a variable could enter and then immediately qualify for removal, so the search could cycle. - The
_aicversions run the same AIC searches asstep(), with a fuller printed report.
Interpretation
- All six olsrr searches choose the same eight predictor variables, including the backward search with its more lenient p > .10 removal rule.
- Backward removal order: P_B, D_E, IR, then EG.
- On these data the choice of threshold does not change the result. In the Credit data, by contrast, looser thresholds (.30) also keep Gender (Female).
# Exhaustive best subsets
bs <- regsubsets(MR ~ ., data = finance,
nvmax = 12)
bs_sum <- summary(bs)
which.max(bs_sum$adjr2) # best adjusted R-squared[1] 8which.min(bs_sum$cp) # best Mallows' Cp[1] 8which.min(bs_sum$bic) # best BIC[1] 8round(data.frame(size = 1:12,
adjr2 = bs_sum$adjr2, cp = bs_sum$cp,
bic = bs_sum$bic), 3) size adjr2 cp bic
1 1 0.401 297.200 -244.892
2 2 0.493 176.279 -322.868
3 3 0.550 101.619 -377.382
4 4 0.592 46.617 -421.644
5 5 0.608 26.991 -435.996
6 6 0.617 16.763 -441.797
7 7 0.621 12.315 -442.036
8 8 0.625 7.170 -443.072
9 9 0.625 8.384 -437.662
10 10 0.625 9.640 -432.209
11 11 0.625 11.097 -426.552
12 12 0.624 13.000 -420.437coef(bs, which.min(bs_sum$bic))(Intercept) Beta P_E DY ROE Size
-0.03647434 0.20573600 0.07453615 0.29472713 0.06385921 0.17897854
MOM VIX SR
0.21068094 -0.06316503 -0.02551022 What the code does
regsubsets()fits the best model of every size from 1 to 12 predictor variables.nvmax = 12because there are 12 candidates.which.max()andwhich.min()find the best size under each criterion.data.frame()lays the three criteria side by side for every size.
Interpretation
- All three criteria choose 8 predictor variables, the same subset the stepwise searches chose.
- Adjusted R² is flat at .625 from 8 to 11 predictor variables, so the extra predictor variables add essentially nothing.
- Mallows’ Cp = 7.2 at size 8: a good model has Cp close to or below p + 1 (8 + 1 = 9), which indicates little bias.
- BIC is lowest at size 8 (−443.1) and worsens after that. BIC barely separates sizes 7 and 8 (−442.0 against −443.1).
- With equal fit, the simpler model is preferred.
Answers to the Example 2 questions
- Do the AIC-based searches agree? Yes: all three keep the same eight predictor variables.
- Order of entry: Beta enters first and SR last.
- Do the p-value searches agree? Yes. They choose the same model even with the more lenient removal threshold.
- Best size by each criterion: 8 by adjusted R², by Cp and by BIC.
- Why prefer 8 over 11? Adjusted R² does not improve, Cp is lowest, and BIC gets worse after 8.
- Excluded predictors: P_B, EG, D_E and IR. Recommendation: eight predictor variables explain about 63% of the variation in market returns; adding price-to-book, earnings growth, leverage and interest rates raises R² by only about .002.
Agreement across every method is evidence of a stable model, not proof that it is the best possible model.
3.5 Example 3: Final model and assumption checks
Selection procedures choose predictor variables. They do not check the regression assumptions, so the final model must be checked like any other regression model.
final <- lm(MR ~ Beta + P_E + DY + ROE + Size +
MOM + VIX + SR, data = finance)
summary(final)
Call:
lm(formula = MR ~ Beta + P_E + DY + ROE + Size + MOM + VIX +
SR, data = finance)
Residuals:
Min 1Q Median 3Q Max
-1.86481 -0.42115 0.00641 0.42340 2.10224
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.036474 0.030496 -1.196 0.232260
Beta 0.205736 0.019397 10.607 < 2e-16 ***
P_E 0.074536 0.022203 3.357 0.000849 ***
DY 0.294727 0.043092 6.840 2.37e-11 ***
ROE 0.063859 0.020055 3.184 0.001544 **
Size 0.178979 0.027073 6.611 9.98e-11 ***
MOM 0.210681 0.024058 8.757 < 2e-16 ***
VIX -0.063165 0.019185 -3.292 0.001065 **
SR -0.025510 0.009526 -2.678 0.007653 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.6752 on 491 degrees of freedom
Multiple R-squared: 0.6314, Adjusted R-squared: 0.6254
F-statistic: 105.1 on 8 and 491 DF, p-value: < 2.2e-16Interpretation
- Fit: R² = .631 and adjusted R² = .625. The overall model is significant: F(8, 491) = 105.1, p < .001.
- All eight predictor variables are significant (all p ≤ .008).
- Signs of the coefficients, holding the other
predictor variables constant:
- Positive: Beta, P_E, DY, ROE, Size and MOM. Higher values go with higher MR.
- Negative: VIX and SR.
- SR reverses sign. SR correlates positively with MR on its own (r = .26), but its coefficient here is negative (b = −0.026). This is why the zero-order correlations must not be read as the effects in the model.
3.5.1 Predictor importance
ols_regress(final) # Std. Beta column: standardized coefficients Model Summary
-----------------------------------------------------------------
R 0.795 RMSE 0.669
R-Squared 0.631 MSE 0.448
Adj. R-Squared 0.625 Coef. Var -1189.208
Pred R-Squared 0.618 AIC 1037.158
MAE 0.522 SBC 1079.304
-----------------------------------------------------------------
RMSE: Root Mean Square Error
MSE: Mean Square Error
MAE: Mean Absolute Error
AIC: Akaike Information Criteria
SBC: Schwarz Bayesian Criteria
ANOVA
---------------------------------------------------------------------
Sum of
Squares DF Mean Square F Sig.
---------------------------------------------------------------------
Regression 383.453 8 47.932 105.127 0.0000
Residual 223.866 491 0.456
Total 607.319 499
---------------------------------------------------------------------
Parameter Estimates
----------------------------------------------------------------------------------------
model Beta Std. Error Std. Beta t Sig lower upper
----------------------------------------------------------------------------------------
(Intercept) -0.036 0.030 -1.196 0.232 -0.096 0.023
Beta 0.206 0.019 0.374 10.607 0.000 0.168 0.244
P_E 0.075 0.022 0.111 3.357 0.001 0.031 0.118
DY 0.295 0.043 0.205 6.840 0.000 0.210 0.379
ROE 0.064 0.020 0.104 3.184 0.002 0.024 0.103
Size 0.179 0.027 0.212 6.611 0.000 0.126 0.232
MOM 0.211 0.024 0.297 8.757 0.000 0.163 0.258
VIX -0.063 0.019 -0.111 -3.292 0.001 -0.101 -0.025
SR -0.026 0.010 -0.091 -2.678 0.008 -0.044 -0.007
----------------------------------------------------------------------------------------ols_correlations(final) # zero-order, partial and part (semipartial) correlations Correlations
-------------------------------------------
Variable Zero Order Partial Part
-------------------------------------------
Beta 0.634 0.432 0.291
P_E 0.181 0.150 0.092
DY 0.134 0.295 0.187
ROE 0.314 0.142 0.087
Size 0.503 0.286 0.181
MOM 0.593 0.368 0.240
VIX -0.496 -0.147 -0.090
SR 0.265 -0.120 -0.073
-------------------------------------------What the code does
ols_regress()prints the model with a Std. Beta column of standardized coefficients.ols_correlations()prints, for each predictor variable:- the zero-order correlation, the ordinary r with MR;
- the partial correlation;
- the part (semipartial) correlation.
What each measure means
- Standardized coefficient (beta): the change in MR, in standard deviations, for a one-standard-deviation increase in the predictor variable, holding the others constant. Because every predictor variable is on the same scale, the betas can be compared.
- Partial correlation: the correlation between MR and the predictor variable after the other seven predictor variables have been removed from both.
- Part (semipartial) correlation, sr: the correlation between MR and the part of the predictor variable that is unique to it, with the other predictor variables removed from the predictor only. sr² is the drop in R² if that predictor variable were removed from the model.
Interpretation
| Predictor | Std. beta | Zero-order r | Partial | Part (sr) | sr² |
|---|---|---|---|---|---|
| Beta | .374 | .634 | .432 | .291 | .085 |
| MOM | .297 | .593 | .368 | .240 | .058 |
| Size | .212 | .503 | .286 | .181 | .033 |
| DY | .205 | .134 | .295 | .187 | .035 |
| P_E | .111 | .181 | .150 | .092 | .009 |
| VIX | −.111 | −.496 | −.147 | −.090 | .008 |
| ROE | .104 | .314 | .142 | .087 | .008 |
| SR | −.091 | .265 | −.120 | −.073 | .005 |
- The most important predictor variable: Beta is first and MOM second on every measure, after controlling for the other predictor variables in the model. Beta alone accounts for .085 of the variance in MR (8.5%) that no other predictor variable explains.
- Third place is a tie. Size and DY are essentially tied: Size has the larger beta (.212 against .205), DY the larger sr (.187 against .181).
- DY shows why controlling matters. It is only weakly correlated with MR on its own (r = .134), but it is one of the strongest predictor variables once the others are controlled.
- The unique pieces add up to only part of R². The sr² values sum to .240 (the rounded values in the table add to .241), while R² is .631. The remaining .392 is variance the predictor variables share, which cannot be credited to any single one of them.
3.5.2 Diagnostic plots and tests
par(mfrow = c(1, 3)) # three plots side by side
plot(final, which = 1) # linearity
plot(final, which = 3) # constant variance
hist(resid(final), breaks = 30, xlab = "Residual",
main = "Histogram of residuals") # normality
par(mfrow = c(1, 1)) # back to one plot per windowbptest(final) # constant variance
studentized Breusch-Pagan test
data: final
BP = 3.181, df = 8, p-value = 0.9225vif(final) # multicollinearity Beta P_E DY ROE Size MOM VIX SR
1.654137 1.465117 1.202264 1.410630 1.376216 1.534505 1.524353 1.554062 sum(abs(rstandard(final)) > 3) # count beyond ±3[1] 2max(abs(rstandard(final))) # largest[1] 3.14086What the code does
- The diagnostic plots:
plot(final, which = 1)is the residuals-versus-fitted plot, which checks linearity.plot(final, which = 3)is the scale-location plot: the square root of the absolute standardized residuals against the fitted values. It checks constant variance.hist(resid(final))checks whether the residuals are approximately normal.par(mfrow = c(1, 3))places the three plots side by side.
- The tests:
bptest()is the Breusch–Pagan test. Its null hypothesis is constant variance (homoscedasticity), so a small p-value is evidence of heteroscedasticity.vif()checks multicollinearity.rstandard()returns the standardized residuals.abs(...) > 3flags values beyond ±3,sum()counts them andmax()finds the largest.
- What the plots add automatically: the red lines are
trend lines added by
plot(), and the numbers mark the three largest residuals. In the residuals-versus-fitted plot, the red line helps reveal curvature. In both plots, do not over-read wiggles at the edges, where there are few points.
Interpretation
| Assumption | Check | Result | Conclusion |
|---|---|---|---|
| Linearity | Residuals vs fitted | Even band around 0; no curve | No evidence of a violation |
| Constant variance | Scale-location; bptest() |
Points at about the same height; BP(8) = 3.18, p = .922 | No evidence of a violation |
| Normality | Histogram of residuals | Symmetric and bell-shaped (skewness ≈ 0.03) | Approximately normal |
| Multicollinearity | vif() |
All VIFs between 1.20 and 1.65 | No concern |
| Independence | Study design | No time order or grouping is documented (the file has no codebook) | Assumed; it cannot be confirmed from these plots |
- Standardized residuals: 2 of the 500 lie beyond ±3, the largest at 3.14. With 500 observations, a couple of values slightly beyond ±3 are expected by chance.
- Word conclusions carefully: say there is no evidence of a violation, not that “the assumptions are satisfied”. Diagnostics can fail to detect a problem; they cannot prove an assumption holds.
- No transformation is needed. A transformation is a remedy for a detected problem, and none was detected. A log transformation is impossible anyway, because MR takes negative values (its minimum is −3.25). Adding a squared term would add complexity without fixing anything.
Recommendation. Report the eight-predictor model. Eight predictor variables explain about 63% of the variation in market returns (R² = .631). Adding price-to-book, earnings growth, leverage and interest rates would raise R² by only about .002.
4 Key takeaways
- Theory sets the order. Hierarchical entry tests whether each block adds explanatory power: report ΔR² and its F test.
- Searches propose models. Forward, backward, stepwise and best subsets compare models by AIC, BIC, Cp, adjusted R² or p-value thresholds.
- Agreement is evidence, not proof. In Finance, every procedure chose the same eight predictor variables. In Credit, every stepwise search kept both Limit and Rating, and only BIC best subsets and the VIF exposed the redundancy. The final Credit model also fails the residual checks.
- Check before you report. Use residual plots,
bptest(), a residual histogram andvif(), and transform only when a check shows a problem.