1.1 Check the Dependent and Independent Variables for Linear Correlation

First, we can use the cor() function from the "stats" package in the base version of R to check for linear correlation between the dependent variable ("rating") and the independent variables, and among the independent variables. The function produces a matrix of Pearson’s correlation coefficients between every possible pair.

cor(attitude)

              rating complaints privileges  learning    raises  critical
rating     1.0000000  0.8254176  0.4261169 0.6236782 0.5901390 0.1564392
complaints 0.8254176  1.0000000  0.5582882 0.5967358 0.6691975 0.1877143
privileges 0.4261169  0.5582882  1.0000000 0.4933310 0.4454779 0.1472331
learning   0.6236782  0.5967358  0.4933310 1.0000000 0.6403144 0.1159652
raises     0.5901390  0.6691975  0.4454779 0.6403144 1.0000000 0.3768830
critical   0.1564392  0.1877143  0.1472331 0.1159652 0.3768830 1.0000000
advance    0.1550863  0.2245796  0.3432934 0.5316198 0.5741862 0.2833432
             advance
rating     0.1550863
complaints 0.2245796
privileges 0.3432934
learning   0.5316198
raises     0.5741862
critical   0.2833432
advance    1.0000000

The first row shows the results for the linear relationship between the dependent variable ("rating") and the independent variables. The results suggest a strong linear relationship between "rating" and some of the independent variables, particularly, "complaints" (0.83), "learning" (0.62), and "raises" (0.59).

We can also use the plot() function to visually check for linear correlations based on the scatter plots of all the variable pairs.

plot(attitude)

Scatter Plots of All Variable Pairs in R

The first row shows the scatter plots of the dependent variable ("rating") versus the independent variables. The results also suggest a strong linear relationship between "rating" and the variables "complaints", "learning", and "raises" which are the first, third, and fourth plots on the first row respectively. Hence, these three variables are candidates in the final model.

1.2 Perform Variable Selection

We could include some of the most correlated variables with "rating" in the model and remove or add variables based on whether they are significant. Alternatively, we could automate the variable selection process by using stepwise regression.

The step() function from the base "stats" package can be used to perform various stepwise regression based on AIC.

Forward Stepwise Regression

Without any special preference for some variables, we can use the forward stepwise regression starting with no variable.

First, we set the intercept model and the full model. The intercept model has no variable in the model, while the full model contains all the variables available.

# Specify the intercept model
intercept_model = lm(rating ~ 1, data = attitude)

# Specify the full model
full_model = lm(rating ~ ., data = attitude)
#or
full_model = lm(rating ~ complaints + privileges + learning +
                raises + critical + advance, data = attitude)

Our variable selection then starts from the intercept (or null) model. Then from the intercept model, based on AIC, we add one variable repeatedly, until no other variable based on the AIC criterion significantly improves model performance when added. Potentially, we can pick variables up to the full set of variables available (or scope) if they all significantly add to the predictive power of the model. Set trace = 1 to see the steps of variable additions, it is set to 0 here.

# Perform forward stepwise regression
step(intercept_model, direction = 'forward', 
     scope = formula(full_model), trace = 0)

Call:
lm(formula = rating ~ complaints + learning, data = attitude)

Coefficients:
(Intercept)   complaints     learning  
     9.8709       0.6435       0.2112  

As can be seen above, the final selection includes only "complaints" and "learning". This means that other variables such as "raises" are not significantly increasing the predictive power of the model when added to the two, and these two are most likely the most predictive pair.

1.3 Run the Regression Model with the Chosen Variables

Run the multiple linear regression model using the lm() function, and print the results using the summary() function.

model = lm(rating ~ complaints + learning, data = attitude)
summary(model)

Call:
lm(formula = rating ~ complaints + learning, data = attitude)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.5568  -5.7331   0.6701   6.5341  10.3610 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   9.8709     7.0612   1.398    0.174    
complaints    0.6435     0.1185   5.432 9.57e-06 ***
learning      0.2112     0.1344   1.571    0.128    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.817 on 27 degrees of freedom
Multiple R-squared:  0.708, Adjusted R-squared:  0.6864 
F-statistic: 32.74 on 2 and 27 DF,  p-value: 6.058e-08

Interpretation of the Results

• Coefficients:
  • The estimated intercept (\(\widehat \beta_0\)) is \(\text{summary(model)\$coefficients[1, 1]}\) \(= 9.871\).
  • The estimated coefficient for "complaints" or \(x_1\) (\(\widehat \beta_1\)) is \(\text{summary(model)\$coefficients[2, 1]}\) \(= 0.644\).
  • The estimated coefficient for "learning" or \(x_2\) (\(\widehat \beta_2\)) is \(\text{summary(model)\$coefficients[3, 1]}\) \(= 0.211\).
• P-values:
  For level of significance, \(\alpha = 0.05\).
  • The p-value for "complaints" is \(\text{summary(model)\$coefficients[2, 4]}\) \(= 9.57\times 10^{-6}\).
    Since the p-value is less than \(0.05\), we say that "complaints" is a statistically significant predictor of "rating".
  • The p-value for "learning" is \(\text{summary(model)\$coefficients[3, 4]}\) \(= 0.128\).
    Since the p-value is greater than \(0.05\), we say that "learning" is NOT a statistically significant predictor of "rating" in the presence of "complaints". However, this variable was included in the model because it increases the predictive power of the model when added to the "complaints" variable and passes the AIC criterion.
  • Remember that "learning" is highly correlated with "rating", with correlation of \(0.6236782\). Hence, when it is the only variable in the model the p-value is as low as \(0.000231\). Also, the r-squared is as high as \(0.389\) (or \(0.6236782^2\)). This means while "learning" is not statistically significant in the presence of "complaints", it is a useful predictor when added as a second variable as supported by the AIC criterion.
• R-squared:
  The r-squared value is \(\text{summary(model)\$r.squared}\) \(= 0.708\).
  This means that the model, using the independent variables, "complaints" and "learning", explains \(70.8\%\) of the variations in the dependent variable \(y\) (rating) from the sample mean of \(y\), \(\bar y\).

1.4 Run the Regression Model Diagnostics

R allows you to plot standard regression diagnostics plots as below, and our focus here will be on the first two plots:

par(mfrow= c(2,2))
plot(model)

Standard Linear Regression Diagnostics Plots in R

The multiple linear regression framework has the theoretical assumption that the error terms (or residuals) are 1) independent, 2) normal distributed, 3) have constant variance, and 4) have mean zero.

By construction of the ordinary least squares model, the error terms have mean zero. Hence, we only need to check for 1) independence, 2) normal distribution, and 3) constant variance of the residuals.

Check for Independence and Constant Variance of Residuals

To check for independence and constant variance of the residuals, we can:

• Look at the scatter plot of the residuals versus the fitted values and confirm that the residuals do not oscillate around zero in any noticeable or clear pattern. The plot of the residuals versus the fitted values below (or above) shows no pattern, hence, there is no evidence auto-correlation and they are independent.

• Look at the scatter plot of the residuals versus the fitted values and confirm that the residuals do not depend on the predicted values. That is, the residuals do not increase (widen) or decrease (narrow) as the predicted value increases or decreases. The plot of the residuals versus the fitted values below (or above) does not show any widening or narrowing of the residuals, hence, there is no heteroskedasticity, and the constant variance assumption of residuals is not violated.

• Look at the autocorrelation plot from the acf() function in the base "stats" package. If the vertical lines, excluding the first one, are between the two blue horizontal lines, we can conclude that the residuals are not auto-correlated. All the lines are between the blue lines, this is another evidence of independence of residuals.

plot(model$fitted.values, model$residuals,
     main = "Residuals versus Fitted Values")
abline(h=0)

acf(model$residuals)

Regression Residual versus Fitted Values and Auto-correlation Plot in R

Check for Normality of the Residuals

Here we check for the normality of the residuals.

Normal Q-Q Plot

To check for the normality of the residuals, we can look at the quantile-quantile plot of the residuals versus the theoretical quantiles of from the standard normal distribution, and confirm that there is an alignment between the two.

If the values of the residuals are close to the theoretical values from normal distribution, the residuals will be very close to the quantile-quantile line on the plot. The closer they are to the line, the more the residuals fit to a normal distribution.

qqnorm(model$residuals)
qqline(model$residuals)

Quantile versus Quantile Plot of Regression Residuals in R

The result suggests that the residuals, given the sample size, fairly fit to a normal distribution as the dots are close to the diagonal except for about two values on the left tail and two values on the right tail.

Shapiro-Wilk Test

We can also use the Shapiro-Wilk test via the shapiro.test() function from the base "stats" package. At \(\alpha = 0.05\), that is, \(5\%\) level of significance, we reject the null hypothesis that the residuals are normal distributed, if the p-value is less than \(0.05\). Else, we fail to reject that they are normal distributed.

shapiro.test(model$residuals)

    Shapiro-Wilk normality test

data:  model$residuals
W = 0.94486, p-value = 0.123

With the p-value \(0.123\) being greater than \(0.05\), we fail to reject the normal distribution assumption. Hence, our residuals are in line with the normal distribution assumption.

1.5 Use the Model for Prediction

To perform predictions, we can use the predict.lm() function from the base "stats" package, or calculate directly.

Calculate Predictions Directly

For complaints = 70 (\(x_1 = 70\)), and learning = 60 (\(x_2 = 60\)).

#Set the coefficient and intercept
b0 = 9.8709; b1 = 0.6435; b2 = 0.2112

#or

b0 = model$coefficients[1]
b1 = model$coefficients[2]
b2 = model$coefficients[3]

prediction = b0 + b1*70 + b2*60
prediction

(Intercept) 
   67.58862 

Use the predict.lm() Function for Prediction

For complaints = 70 (\(x_1 = 70\)), and learning = 60 (\(x_2 = 60\)).

predict.lm(model, data.frame(complaints = 70, learning = 60))

       1 
67.58862 

For multiple values of both complaints and learning, for example, complaints = \(\{70, 60, 50\}\), and learning = \(\{60, 50, 45\}\).

complaints = c(70, 60, 50); learning = c(60, 50, 45)
pred_data = data.frame(complaints, learning)
pred_data

  complaints learning
1         70       60
2         60       50
3         50       45

predict.lm(model, pred_data)

       1        2        3 
67.58862 59.04153 51.55039 

We can also see the fitted values from the "attitude" data itself:

predict.lm(model, attitude)

       1        2        3        4        5        6        7        8 
50.92676 62.46037 69.48935 60.33851 74.00392 54.55679 64.81330 69.75025 
       9       10       11       12       13       14       15       16 
76.78918 59.05147 56.22644 56.71842 58.63903 72.78648 74.62755 82.99328 
      17       18       19       20       21       22       23       24 
79.14211 64.32132 66.95505 58.59926 42.79211 62.21935 62.90263 45.93016 
      25       26       27       28       29       30 
54.75804 72.72682 73.76290 56.05502 79.56450 75.09964