Multiple Linear Regression and Shapiro-Wilk normality test

Question 1

#install.packages("ISLR")
library(ISLR)

Hitters_Fixed = na.omit(Hitters)

Hitters_Fixed

str(Hitters_Fixed)

## 'data.frame':    263 obs. of  20 variables:
##  $ AtBat    : int  315 479 496 321 594 185 298 323 401 574 ...
##  $ Hits     : int  81 130 141 87 169 37 73 81 92 159 ...
##  $ HmRun    : int  7 18 20 10 4 1 0 6 17 21 ...
##  $ Runs     : int  24 66 65 39 74 23 24 26 49 107 ...
##  $ RBI      : int  38 72 78 42 51 8 24 32 66 75 ...
##  $ Walks    : int  39 76 37 30 35 21 7 8 65 59 ...
##  $ Years    : int  14 3 11 2 11 2 3 2 13 10 ...
##  $ CAtBat   : int  3449 1624 5628 396 4408 214 509 341 5206 4631 ...
##  $ CHits    : int  835 457 1575 101 1133 42 108 86 1332 1300 ...
##  $ CHmRun   : int  69 63 225 12 19 1 0 6 253 90 ...
##  $ CRuns    : int  321 224 828 48 501 30 41 32 784 702 ...
##  $ CRBI     : int  414 266 838 46 336 9 37 34 890 504 ...
##  $ CWalks   : int  375 263 354 33 194 24 12 8 866 488 ...
##  $ League   : Factor w/ 2 levels "A","N": 2 1 2 2 1 2 1 2 1 1 ...
##  $ Division : Factor w/ 2 levels "E","W": 2 2 1 1 2 1 2 2 1 1 ...
##  $ PutOuts  : int  632 880 200 805 282 76 121 143 0 238 ...
##  $ Assists  : int  43 82 11 40 421 127 283 290 0 445 ...
##  $ Errors   : int  10 14 3 4 25 7 9 19 0 22 ...
##  $ Salary   : num  475 480 500 91.5 750 ...
##  $ NewLeague: Factor w/ 2 levels "A","N": 2 1 2 2 1 1 1 2 1 1 ...
##  - attr(*, "na.action")= 'omit' Named int [1:59] 1 16 19 23 31 33 37 39 40 42 ...
##   ..- attr(*, "names")= chr [1:59] "-Andy Allanson" "-Billy Beane" "-Bruce Bochte" "-Bob Boone" ...

1a) Run a regression analysis that uses “Salary” as the outcome and “Hits” and “Years” as predictors. Show all your work.

hitters_multiple_R = lm(Salary ~ Hits + Years, data=Hitters_Fixed)
summary(hitters_multiple_R)

## 
## Call:
## lm(formula = Salary ~ Hits + Years, data = Hitters_Fixed)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -838.06 -211.89  -43.40   76.07 2248.37 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -199.2510    67.4690  -2.953  0.00343 ** 
## Hits           4.3124     0.5013   8.603 7.46e-16 ***
## Years         36.9501     4.7187   7.831 1.24e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 366.1 on 260 degrees of freedom
## Multiple R-squared:  0.3465, Adjusted R-squared:  0.3415 
## F-statistic: 68.94 on 2 and 260 DF,  p-value: < 2.2e-16

Salary = -199.251 + 4.3124 * (Hits) + 36.9501 * (Years)

1a) Interpret the value of the Residual Standard Error that you obtained.

Note: Interpreting the value of RSE means commenting on what the value means. It is NOT whether the value is low or high and whether the equation is good. It is simply interpreting the meaning of its value.

summary(lm(Salary ~ Hits + Years, data=Hitters_Fixed))$sigma

## [1] 366.069

The Residual Standard Error(RSE) of 366.1 provides a measure of the average deviation of actual salaries from the salaries predicted by the model. On average, the values of actual Salary deviate from the regression equation by $366,069. So, we can expect the predictions of Salary based on the equation to be off by around $366,069.

coefficient of variation = RSE/sample mean

CV <- 366.069/mean(Hitters_Fixed$Salary)

(366.069/mean(Hitters_Fixed$Salary))*100

## [1] 68.3059

The coefficient of variation is 68.3%. An RSE of $366,069 is showing measures of standard deviation of residuals (or errors) in the model predictions that are widely spread from the observed data, indicating a poor fit. Therefore, the RSE is suggesting the model does not fit the actual data very well. A coefficient of variation of 68.3% is way higher than a 10% or 20% coefficient, which would be the ideal magnitude of error compared to a typical value of Y.

1b) Does the equation that you obtained in part a allow you to eliminate at least 50% of the variation existing in the values of Salary? Justify.

summary(hitters_multiple_R)$adj.r.squared

## [1] 0.3415173

Adjusted R squared is 0.3415 (34.15%)

If we use the linear equation instead of the sample mean to make predictions of Y, we can reduce only 34.15% of the variability in the values of Y. Since we reduce less than 50% of the variability compared to using the sample mean, this indicates that the model does not explain half or more of the variation in Salary.

1c) For the equation you obtained in 1 a), assume that Assumption 2 regarding the validity of the linear regression analysis is satisfied. Run an analysis to check for the validity of assumptions 3 and 4 (DO NOT check assumption 1). Show all your work AND comment on your findings.

Assumptions 3

hypothesis test: the Shapiro test to check for normality

Ho: The residuals follow a Normal distribution Ha: The residuals do NOT follow a Normal distribution

shapiro.test(residuals(hitters_multiple_R))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(hitters_multiple_R)
## W = 0.84691, p-value = 1.975e-15

PV= 1.975e-15, which is a lot smaller than alpha (0.05). Thus, we reject Ho and support Ha. The data is giving us evidence that the residuals do NOT follow a Normal distribution. Therefore, Assumption 3 is not satisfied.

Assumptions 4

plot(predict(hitters_multiple_R), residuals(hitters_multiple_R), xlab = "Predict", ylab="Residuals")

abline(h=0, col = "red")

The plot is showing some outliers and the residuals getting spread as the predicted value increase. Therefore, the variability is not the same for all values. Assumption 4 is not satisfied

1d) Add the variable “RBI” as a third predictor to the equation that you obtained in part a. Show all your work.

hitters_multiple_R1 = lm(Salary ~ Hits + Years + RBI, data=Hitters_Fixed)
summary(hitters_multiple_R1)

## 
## Call:
## lm(formula = Salary ~ Hits + Years + RBI, data = Hitters_Fixed)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -724.40 -190.04  -45.45   85.01 2248.38 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -194.5418    67.1578  -2.897 0.004092 ** 
## Hits           3.0545     0.8183   3.733 0.000233 ***
## Years         35.2177     4.7782   7.371 2.28e-12 ***
## RBI            2.7890     1.4385   1.939 0.053618 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 364.1 on 259 degrees of freedom
## Multiple R-squared:  0.3559, Adjusted R-squared:  0.3484 
## F-statistic:  47.7 on 3 and 259 DF,  p-value: < 2.2e-16

1d1) Write the equation.

Salary = -194.5418 + 3.0545 * (Hits) + 35.2177 * (Years) + 2.789 * (RBI)

1d2) Use Adjusted R Squared to decide whether this equation is better than the one you obtained in part a.

summary(lm(Salary ~ Hits + Years, data=Hitters_Fixed))$adj.r.squared

## [1] 0.3415173

summary(lm(Salary ~ Hits + Years + RBI, data=Hitters_Fixed))$adj.r.squared

## [1] 0.3484308

(summary(lm(Salary ~ Hits + Years + RBI, data=Hitters_Fixed))$adj.r.squared -  summary(lm(Salary ~ Hits + Years, data=Hitters_Fixed))$adj.r.squared) / summary(lm(Salary ~ Hits + Years, data=Hitters_Fixed))$adj.r.squared

## [1] 0.0202436

No, the new equation with “RBI” added produces almost the same results in Adjusted R Squared when compared to the original equation that only includes “Hits” and “Years” as predictors. The adjusted R squared for the second equation is only 2% larger than the first equation and therefore does not meet the at least 5% increase minimum to be considered better.

Question 2

Apply Linear Regression to predict “Salary” based on these predictors: Hits, League, and Division. Show all your work.

hitters_multiple_R2 = lm(Salary ~ Hits + League + Division, data=Hitters_Fixed)
summary(hitters_multiple_R2)

## 
## Call:
## lm(formula = Salary ~ Hits + League + Division, data = Hitters_Fixed)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -844.52 -264.61  -59.75  169.95 1958.07 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  118.1472    78.3232   1.508  0.13266    
## Hits           4.3350     0.5574   7.778 1.75e-13 ***
## LeagueN       46.7805    50.1133   0.933  0.35143    
## DivisionW   -140.7468    49.6199  -2.837  0.00492 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 400.8 on 259 degrees of freedom
## Multiple R-squared:  0.2196, Adjusted R-squared:  0.2105 
## F-statistic: 24.29 on 3 and 259 DF,  p-value: 6.985e-14

Salary = 118.1472 + 4.335 * (Hits) + 46.7805 * (League) - 140.7468 * (Division)

Based on the results, write the equation ONLY with the predictors that show a significant partial t-test. Justify your selections/decisions.

Salary = 118.1472 + 4.335 * (Hits) - 140.7468 * (Division)

The p-values associated to the t value for Hits and Division are 1.75e-13 and 0.00492 respectively. These p-values are smaller than alpha (0.05). We reject Ho, thus, we conclude that the relationship between Hits, Division and Salary are statistically significant. Adding Hits and Division to the equation reduces an amount of variability in Salary that is statistically significant and therefore, the amount of variability reduced by adding Hits and Division is not random, it is a real amount of variability. League on the other hand has a large p-value of 0.35143 that does NOT show a significant partial t-test and therefore, it is excluded from the equation.

Use the equation with the predictors you chose in part a) to predict the Salary for a player who connected 120 hits and played in the National League in the West division in 1986. Show your work.

prediction3 = 118.1472 + 4.335 * (120) - 140.7468 * (as.numeric(Hitters_Fixed$Division[1]))
prediction3

## [1] 356.8536

A player who connected 120 hits and played in National League in the West division in 1986 is estimated to have a salary of $356,853.60

Does the equation with the predictors you chose in part a) improve when you add the “NewLeague” predictor? Justify your answer using RSE. Show all your work.

hitters_multiple_R3 = lm(Salary ~ Hits + Division + NewLeague, data=Hitters_Fixed)
summary(hitters_multiple_R3)

## 
## Call:
## lm(formula = Salary ~ Hits + Division + NewLeague, data = Hitters_Fixed)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -833.27 -252.26  -64.79  174.46 1959.69 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  129.4426    76.6772   1.688  0.09259 .  
## Hits           4.2907     0.5536   7.750  2.1e-13 ***
## DivisionW   -141.1239    49.6555  -2.842  0.00484 ** 
## NewLeagueN    33.9060    49.8290   0.680  0.49683    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 401.1 on 259 degrees of freedom
## Multiple R-squared:  0.2183, Adjusted R-squared:  0.2093 
## F-statistic: 24.11 on 3 and 259 DF,  p-value: 8.539e-14

(summary(lm(Salary ~ Hits + Division + NewLeague, data= Hitters_Fixed))$sigma / mean(Hitters_Fixed$Salary))*100

## [1] 74.85093

We want a coefficient of variance of at most <=20%. In this case, it is quite above that (74.85%). Therefore, we consider that the RSE is not as low as we want. We want an equation with a lower RSE.

(summary(lm(Salary ~ Hits + Division, data=Hitters_Fixed))$sigma -  summary(lm(Salary ~ Hits + Division + NewLeague, data=Hitters_Fixed))$sigma) / summary(lm(Salary ~ Hits + Division, data=Hitters_Fixed))$sigma

## [1] -0.001034281

In addition, the equation with the third predictor NewLeague added produces an increase in RSE of 0.1% compared to the previous equation. This is the opposite of what we want and well below the 5% decrease required for it to be useful (> 5%) when compared to the equation with two predictors.