HW: Find a data set of which you can fit multiple linear regression and interpret your results

A good dataset for Multiple Linear Regression (MLR) should have:

.One continuous dependent variable (response variable).

.Two or more independent variables (predictors).

.Enough observations to identify relationships.

Let’s use Zee Entertainment Enterprises as in the last work we did

Dependent Variable (Y)

.Close

Independent Variables (X)

.Prev Close
.Open
.High
.Low
.Volume
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.1     ✔ readr     2.1.6
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.2     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
ZEEL <- read.csv("ZEEL.csv")

model <- lm(
  Close ~ `Prev.Close` + Open + High + Low + Volume,
  data = ZEEL
)

summary(model)
## 
## Call:
## lm(formula = Close ~ Prev.Close + Open + High + Low + Volume, 
##     data = ZEEL)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -92.735  -1.512  -0.285   1.253  67.805 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.500e-01  1.329e-01   3.387 0.000712 ***
## Prev.Close  -6.050e-02  1.145e-02  -5.285 1.31e-07 ***
## Open        -4.981e-01  1.390e-02 -35.830  < 2e-16 ***
## High         7.368e-01  9.363e-03  78.689  < 2e-16 ***
## Low          8.209e-01  7.854e-03 104.523  < 2e-16 ***
## Volume       2.585e-08  8.287e-09   3.119 0.001823 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.843 on 5300 degrees of freedom
## Multiple R-squared:  0.9992, Adjusted R-squared:  0.9992 
## F-statistic: 1.401e+06 on 5 and 5300 DF,  p-value: < 2.2e-16

A multiple linear regression model was fitted to examine the relationship between the closing stock price (Close) of Zee Entertainment Enterprises Ltd (ZEEL) and the predictor variables Prev.Close, Open, High, Low, and Volume. The results indicate that all predictor variables are statistically significant since their p-values are less than 0.05. The variables High and Low have strong positive effects on the closing price, while Prev.Close and Open have negative coefficients after controlling for the effects of the other variables.

Trading volume also has a positive and statistically significant effect, although its coefficient is very small due to the large scale of volume values. The model achieved an R-squared value of 0.9992, meaning that approximately 99.92% of the variation in the closing price is explained by the independent variables included in the model. Furthermore, the overall model is highly significant (F-statistic = 1.401 × 10⁶, p-value < 2.2 × 10⁻¹⁶), indicating that the predictors collectively have a significant relationship with the closing price. These results suggest that the model has excellent predictive performance for ZEEL stock prices.

However, because the stock price variables (Prev.Close, Open, High, Low, and Close) are likely to be highly correlated, multicollinearity may be present, which can make individual coefficient estimates less stable and more difficult to interpret even though the overall predictive accuracy remains very high.

Introduction

Variable selection identifies the most relevant predictors in a model, improving accuracy and interpretability while reducing overfitting.

Methods

1. Filter Methods

Rank variables by statistical scores before modeling.

# Correlation matrix on mtcars
round(cor(mtcars), 2)
##        mpg   cyl  disp    hp  drat    wt  qsec    vs    am  gear  carb
## mpg   1.00 -0.85 -0.85 -0.78  0.68 -0.87  0.42  0.66  0.60  0.48 -0.55
## cyl  -0.85  1.00  0.90  0.83 -0.70  0.78 -0.59 -0.81 -0.52 -0.49  0.53
## disp -0.85  0.90  1.00  0.79 -0.71  0.89 -0.43 -0.71 -0.59 -0.56  0.39
## hp   -0.78  0.83  0.79  1.00 -0.45  0.66 -0.71 -0.72 -0.24 -0.13  0.75
## drat  0.68 -0.70 -0.71 -0.45  1.00 -0.71  0.09  0.44  0.71  0.70 -0.09
## wt   -0.87  0.78  0.89  0.66 -0.71  1.00 -0.17 -0.55 -0.69 -0.58  0.43
## qsec  0.42 -0.59 -0.43 -0.71  0.09 -0.17  1.00  0.74 -0.23 -0.21 -0.66
## vs    0.66 -0.81 -0.71 -0.72  0.44 -0.55  0.74  1.00  0.17  0.21 -0.57
## am    0.60 -0.52 -0.59 -0.24  0.71 -0.69 -0.23  0.17  1.00  0.79  0.06
## gear  0.48 -0.49 -0.56 -0.13  0.70 -0.58 -0.21  0.21  0.79  1.00  0.27
## carb -0.55  0.53  0.39  0.75 -0.09  0.43 -0.66 -0.57  0.06  0.27  1.00
library(corrplot)
corrplot(cor(mtcars), method = "color", tl.cex = 0.8)

# Variance of each predictor (drop near-zero variance)
sapply(mtcars, var)
##          mpg          cyl         disp           hp         drat           wt 
## 3.632410e+01 3.189516e+00 1.536080e+04 4.700867e+03 2.858814e-01 9.573790e-01 
##         qsec           vs           am         gear         carb 
## 3.193166e+00 2.540323e-01 2.489919e-01 5.443548e-01 2.608871e+00

2. Wrapper Methods

Use model performance to evaluate subsets.

  • Forward Selection – start empty, add one variable at a time
  • Backward Elimination – start full, remove one variable at a time
  • Stepwise – combination of both
library(MASS)

full  <- lm(mpg ~ ., data = mtcars)
empty <- lm(mpg ~ 1, data = mtcars)

# Forward selection
fwd <- stepAIC(empty, scope = list(upper = full), direction = "forward", trace = FALSE)
summary(fwd)
## 
## Call:
## lm(formula = mpg ~ wt + cyl + hp, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9290 -1.5598 -0.5311  1.1850  5.8986 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 38.75179    1.78686  21.687  < 2e-16 ***
## wt          -3.16697    0.74058  -4.276 0.000199 ***
## cyl         -0.94162    0.55092  -1.709 0.098480 .  
## hp          -0.01804    0.01188  -1.519 0.140015    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.512 on 28 degrees of freedom
## Multiple R-squared:  0.8431, Adjusted R-squared:  0.8263 
## F-statistic: 50.17 on 3 and 28 DF,  p-value: 2.184e-11
# Backward elimination
bwd <- stepAIC(full, direction = "backward", trace = FALSE)
summary(bwd)
## 
## Call:
## lm(formula = mpg ~ wt + qsec + am, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4811 -1.5555 -0.7257  1.4110  4.6610 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   9.6178     6.9596   1.382 0.177915    
## wt           -3.9165     0.7112  -5.507 6.95e-06 ***
## qsec          1.2259     0.2887   4.247 0.000216 ***
## am            2.9358     1.4109   2.081 0.046716 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.459 on 28 degrees of freedom
## Multiple R-squared:  0.8497, Adjusted R-squared:  0.8336 
## F-statistic: 52.75 on 3 and 28 DF,  p-value: 1.21e-11

3. Embedded Methods

Selection happens during model training.

  • LASSO – shrinks some coefficients exactly to zero (L1 penalty)
  • Ridge – shrinks but retains all variables (L2 penalty)
  • Elastic Net – mix of LASSO and Ridge
library(glmnet)

x <- model.matrix(mpg ~ ., data = mtcars)[, -1]
y <- mtcars$mpg

# --- LASSO (alpha = 1) ---
lasso_cv <- cv.glmnet(x, y, alpha = 1)
plot(lasso_cv)                          # CV error vs log(lambda)

coef(lasso_cv, s = "lambda.min")        # non-zero = selected variables
## 11 x 1 sparse Matrix of class "dgCMatrix"
##              lambda.min
## (Intercept) 36.31068553
## cyl         -0.87115711
## disp         .         
## hp          -0.01394790
## drat         .         
## wt          -2.72254841
## qsec         .         
## vs           .         
## am           0.33946305
## gear         .         
## carb        -0.05780686
# --- Ridge (alpha = 0) ---
ridge_cv <- cv.glmnet(x, y, alpha = 0)
coef(ridge_cv, s = "lambda.min")
## 11 x 1 sparse Matrix of class "dgCMatrix"
##               lambda.min
## (Intercept) 21.051283516
## cyl         -0.374112703
## disp        -0.005318127
## hp          -0.011506803
## drat         1.055629523
## wt          -1.204585685
## qsec         0.160391657
## vs           0.787069385
## am           1.591536197
## gear         0.541785546
## carb        -0.534626533
# --- Elastic Net (alpha = 0.5) ---
enet_cv <- cv.glmnet(x, y, alpha = 0.5)
coef(enet_cv, s = "lambda.min")
## 11 x 1 sparse Matrix of class "dgCMatrix"
##               lambda.min
## (Intercept) 31.507946409
## cyl         -0.656456607
## disp        -0.001865773
## hp          -0.013916343
## drat         0.563405569
## wt          -2.058202045
## qsec         .          
## vs           0.231158096
## am           1.053549722
## gear         .          
## carb        -0.290622942

4. Random Forest Importance

Tree-based methods rank variables by their contribution to reducing impurity.

library(randomForest)

set.seed(42)
rf_model <- randomForest(mpg ~ ., data = mtcars, importance = TRUE)
varImpPlot(rf_model, main = "Variable Importance (Random Forest)")

# Numeric importance scores
importance(rf_model)
##        %IncMSE IncNodePurity
## cyl  12.620655     178.16498
## disp 13.790918     253.10985
## hp   11.548795     169.12588
## drat  4.385795      65.69089
## wt   12.737424     259.44080
## qsec  2.985945      33.48435
## vs    4.408521      34.18138
## am    2.127998      15.04090
## gear  3.340548      20.09017
## carb  7.184614      32.49537

5. Boruta (All-Relevant Selection)

Boruta compares each variable against random “shadow” features to decide importance.

library(Boruta)

set.seed(42)
boruta_out <- Boruta(mpg ~ ., data = mtcars, doTrace = 0)
print(boruta_out)
## Boruta performed 24 iterations in 0.23522 secs.
##  10 attributes confirmed important: am, carb, cyl, disp, drat and 5
## more;
##  No attributes deemed unimportant.
plot(boruta_out, las = 2, cex.axis = 0.7, main = "Boruta Feature Selection")

# Final confirmed variables
getSelectedAttributes(boruta_out, withTentative = FALSE)
##  [1] "cyl"  "disp" "hp"   "drat" "wt"   "qsec" "vs"   "am"   "gear" "carb"