Learning objectives

In this practical, we will use the built-in airquality data set to:

fit a multiple linear regression model;
identify and omit observations with missing values;
perform forward and backward variable selection;
use adjusted \(R^2\) and p-values as the main criteria for comparing models; and
examine a residual plot.

We will deliberately use only the numerical variables relevant to the regression exercise. The variables Month and Day are not used.

The data

The airquality data set is built into R. It contains daily air-quality measurements in New York, from May to September 1973.

# Load the built-in data
data(airquality)

# Look at the data
head(airquality)

str(airquality)

## 'data.frame':    153 obs. of  6 variables:
##  $ Ozone  : int  41 36 12 18 NA 28 23 19 8 NA ...
##  $ Solar.R: int  190 118 149 313 NA NA 299 99 19 194 ...
##  $ Wind   : num  7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
##  $ Temp   : int  67 72 74 62 56 66 65 59 61 69 ...
##  $ Month  : int  5 5 5 5 5 5 5 5 5 5 ...
##  $ Day    : int  1 2 3 4 5 6 7 8 9 10 ...

summary(airquality)

##      Ozone           Solar.R           Wind             Temp      
##  Min.   :  1.00   Min.   :  7.0   Min.   : 1.700   Min.   :56.00  
##  1st Qu.: 18.00   1st Qu.:115.8   1st Qu.: 7.400   1st Qu.:72.00  
##  Median : 31.50   Median :205.0   Median : 9.700   Median :79.00  
##  Mean   : 42.13   Mean   :185.9   Mean   : 9.958   Mean   :77.88  
##  3rd Qu.: 63.25   3rd Qu.:258.8   3rd Qu.:11.500   3rd Qu.:85.00  
##  Max.   :168.00   Max.   :334.0   Max.   :20.700   Max.   :97.00  
##  NAs    :37       NAs    :7                                       
##      Month            Day      
##  Min.   :5.000   Min.   : 1.0  
##  1st Qu.:6.000   1st Qu.: 8.0  
##  Median :7.000   Median :16.0  
##  Mean   :6.993   Mean   :15.8  
##  3rd Qu.:8.000   3rd Qu.:23.0  
##  Max.   :9.000   Max.   :31.0  
##

For this exercise, we will consider:

Response variable: Ozone
Candidate explanatory variables: Solar.R, Wind, and Temp

Thus, the general multiple regression model is

\[ Ozone = \beta_0 + \beta_1 Solar.R + \beta_2 Wind + \beta_3 Temp + \epsilon. \]

We will not use Month or Day in the regression model.

Missing observations

Before fitting the model, we should check for missing values.

colSums(is.na(airquality))

##   Ozone Solar.R    Wind    Temp   Month     Day 
##      37       7       0       0       0       0

There are missing observations in Ozone and Solar.R.

One simple approach for this practical is complete-case analysis: retain only observations that have no missing value in the variables used in the analysis.

aq <- airquality[, c("Ozone", "Solar.R", "Wind", "Temp")]

# Keep only complete observations
aq_complete <- na.omit(aq)

# Check the resulting data
dim(aq)

## [1] 153   4

dim(aq_complete)

## [1] 111   4

head(aq_complete)

The function na.omit() removes rows containing at least one missing value.

We can also verify that there are no missing values left:

colSums(is.na(aq_complete))

##   Ozone Solar.R    Wind    Temp 
##       0       0       0       0

Note: Missing-data handling should be decided before interpreting the regression results. Here we use complete-case analysis because it is simple and appropriate for demonstrating the regression and variable-selection procedures.

Exploratory look at the variables

It is useful to examine the variables before fitting the model.

pairs(aq_complete,
      main = "Scatterplot Matrix: Air Quality Variables")

We can also obtain correlations:

cor(aq_complete)

##              Ozone    Solar.R       Wind       Temp
## Ozone    1.0000000  0.3483417 -0.6124966  0.6985414
## Solar.R  0.3483417  1.0000000 -0.1271835  0.2940876
## Wind    -0.6124966 -0.1271835  1.0000000 -0.4971897
## Temp     0.6985414  0.2940876 -0.4971897  1.0000000

These plots and correlations are useful for understanding the data, but the variable-selection exercise below will be based on regression results.

Full multiple regression model

We first fit the model containing all three candidate explanatory variables.

full_model <- lm(Ozone ~ Solar.R + Wind + Temp,
                 data = aq_complete)

summary(full_model)

## 
## Call:
## lm(formula = Ozone ~ Solar.R + Wind + Temp, data = aq_complete)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -40.485 -14.219  -3.551  10.097  95.619 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -64.34208   23.05472  -2.791  0.00623 ** 
## Solar.R       0.05982    0.02319   2.580  0.01124 *  
## Wind         -3.33359    0.65441  -5.094 1.52e-06 ***
## Temp          1.65209    0.25353   6.516 2.42e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.18 on 107 degrees of freedom
## Multiple R-squared:  0.6059, Adjusted R-squared:  0.5948 
## F-statistic: 54.83 on 3 and 107 DF,  p-value: < 2.2e-16

The regression output contains, among other quantities:

estimated regression coefficients;
standard errors;
t-statistics;
p-values;
\(R^2\);
adjusted \(R^2\).

For model selection in this practical, we will focus on adjusted \(R^2\) and p-values.

Adjusted \(R^2\)

The ordinary \(R^2\) never decreases when an additional explanatory variable is added.

Adjusted \(R^2\) includes a penalty for adding unnecessary variables. Therefore, it is more useful than ordinary \(R^2\) when comparing models with different numbers of predictors.

We can extract adjusted \(R^2\) using:

summary(full_model)$adj.r.squared

## [1] 0.5948449

Forward variable selection

Idea

Forward selection starts with a model containing only the intercept.

At each step, we consider adding one of the available explanatory variables.

For this practical, we will use the following teaching rule:

add a variable if its addition gives a meaningful improvement in adjusted \(R^2\) and its p-value is acceptable;
among candidate additions, prefer the variable giving the larger adjusted \(R^2\);
continue until no remaining variable provides a useful improvement.

We will do this manually, so that students can see what happens at every step.

Step 0: Intercept-only model

m0 <- lm(Ozone ~ 1, data = aq_complete)

summary(m0)

## 
## Call:
## lm(formula = Ozone ~ 1, data = aq_complete)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -41.1  -24.1  -11.1   19.9  125.9 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   42.099      3.158   13.33   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.28 on 110 degrees of freedom

summary(m0)$adj.r.squared

## [1] 0

The intercept-only model contains no explanatory variable.

Step 1: Consider each variable separately

Now fit three one-predictor models.

m_solar <- lm(Ozone ~ Solar.R, data = aq_complete)
m_wind  <- lm(Ozone ~ Wind, data = aq_complete)
m_temp  <- lm(Ozone ~ Temp, data = aq_complete)

summary(m_solar)

## 
## Call:
## lm(formula = Ozone ~ Solar.R, data = aq_complete)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -48.292 -21.361  -8.864  16.373 119.136 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 18.59873    6.74790   2.756 0.006856 ** 
## Solar.R      0.12717    0.03278   3.880 0.000179 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 31.33 on 109 degrees of freedom
## Multiple R-squared:  0.1213, Adjusted R-squared:  0.1133 
## F-statistic: 15.05 on 1 and 109 DF,  p-value: 0.0001793

summary(m_wind)

## 
## Call:
## lm(formula = Ozone ~ Wind, data = aq_complete)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -43.513 -18.597  -5.035  15.814  88.437 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  99.0413     7.4724   13.25  < 2e-16 ***
## Wind         -5.7288     0.7082   -8.09 9.09e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 26.42 on 109 degrees of freedom
## Multiple R-squared:  0.3752, Adjusted R-squared:  0.3694 
## F-statistic: 65.44 on 1 and 109 DF,  p-value: 9.089e-13

summary(m_temp)

## 
## Call:
## lm(formula = Ozone ~ Temp, data = aq_complete)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -40.922 -17.459  -0.874  10.444 118.078 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -147.6461    18.7553  -7.872 2.76e-12 ***
## Temp           2.4391     0.2393  10.192  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.92 on 109 degrees of freedom
## Multiple R-squared:  0.488,  Adjusted R-squared:  0.4833 
## F-statistic: 103.9 on 1 and 109 DF,  p-value: < 2.2e-16

Compare their adjusted \(R^2\) values:

c(
  Solar.R = summary(m_solar)$adj.r.squared,
  Wind    = summary(m_wind)$adj.r.squared,
  Temp    = summary(m_temp)$adj.r.squared
)

##   Solar.R      Wind      Temp 
## 0.1132809 0.3694195 0.4832625

We can also look at the p-value for the predictor in each model:

c(
  Solar.R = coef(summary(m_solar))["Solar.R", "Pr(>|t|)"],
  Wind    = coef(summary(m_wind))["Wind", "Pr(>|t|)"],
  Temp    = coef(summary(m_temp))["Temp", "Pr(>|t|)"]
)

##      Solar.R         Wind         Temp 
## 1.793109e-04 9.089415e-13 1.552677e-17

At this stage, select the variable that provides the strongest improvement according to the chosen criteria.

For the airquality data, students should inspect the output rather than simply accepting a pre-written answer.

Step 2: Add a second variable

Suppose the variable selected in Step 1 is Temp. We now compare the models obtained by adding each of the two remaining variables.

m_temp_solar <- lm(Ozone ~ Temp + Solar.R, data = aq_complete)
m_temp_wind  <- lm(Ozone ~ Temp + Wind, data = aq_complete)

summary(m_temp_solar)

## 
## Call:
## lm(formula = Ozone ~ Temp + Solar.R, data = aq_complete)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -36.610 -15.976  -2.928  12.371 115.555 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -145.70316   18.44672  -7.899 2.53e-12 ***
## Temp           2.27847    0.24600   9.262 2.22e-15 ***
## Solar.R        0.05711    0.02572   2.221   0.0285 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.5 on 108 degrees of freedom
## Multiple R-squared:  0.5103, Adjusted R-squared:  0.5012 
## F-statistic: 56.28 on 2 and 108 DF,  p-value: < 2.2e-16

summary(m_temp_wind)

## 
## Call:
## lm(formula = Ozone ~ Temp + Wind, data = aq_complete)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -42.156 -13.216  -3.123  10.598  98.492 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -67.3220    23.6210  -2.850  0.00524 ** 
## Temp          1.8276     0.2506   7.294 5.29e-11 ***
## Wind         -3.2948     0.6711  -4.909 3.26e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.73 on 108 degrees of freedom
## Multiple R-squared:  0.5814, Adjusted R-squared:  0.5736 
## F-statistic: 74.99 on 2 and 108 DF,  p-value: < 2.2e-16

Compare adjusted \(R^2\):

c(
  adjRsq_Temp_Solar.R = summary(m_temp_solar)$adj.r.squared,
  adjRsq_Temp_Wind    = summary(m_temp_wind)$adj.r.squared
)

## adjRsq_Temp_Solar.R    adjRsq_Temp_Wind 
##           0.5012485           0.5736257

And examine the p-values:

coef(summary(m_temp_solar))

##                  Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept) -145.70315510 18.44671758 -7.898595 2.529334e-12
## Temp           2.27846684  0.24599582  9.262218 2.215559e-15
## Solar.R        0.05710959  0.02571885  2.220534 2.847063e-02

coef(summary(m_temp_wind))

##               Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) -67.321953 23.6210451 -2.850084 5.236066e-03
## Temp          1.827554  0.2505520  7.294112 5.294630e-11
## Wind         -3.294839  0.6711482 -4.909257 3.261728e-06

The variable that improves the model according to our criteria can be retained.

Step 3: Consider the remaining variable

After selecting two variables, consider whether the remaining variable should be added.

For example, if the two-variable model is:

m_two <- lm(Ozone ~ Temp + Wind, data = aq_complete)

summary(m_two)

## 
## Call:
## lm(formula = Ozone ~ Temp + Wind, data = aq_complete)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -42.156 -13.216  -3.123  10.598  98.492 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -67.3220    23.6210  -2.850  0.00524 ** 
## Temp          1.8276     0.2506   7.294 5.29e-11 ***
## Wind         -3.2948     0.6711  -4.909 3.26e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.73 on 108 degrees of freedom
## Multiple R-squared:  0.5814, Adjusted R-squared:  0.5736 
## F-statistic: 74.99 on 2 and 108 DF,  p-value: < 2.2e-16

summary(m_two)$adj.r.squared

## [1] 0.5736257

Then compare it with the full model:

summary(full_model)$adj.r.squared

## [1] 0.5948449

Also inspect the p-value of the additional variable:

coef(summary(full_model))

##                 Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept) -64.34207893 23.05472435 -2.790841 6.226638e-03
## Solar.R       0.05982059  0.02318647  2.579979 1.123664e-02
## Wind         -3.33359131  0.65440710 -5.094063 1.515934e-06
## Temp          1.65209291  0.25352979  6.516366 2.423506e-09

If the additional variable does not provide a useful improvement in adjusted \(R^2\) and has a large p-value, we have a reason to stop the forward-selection procedure.

Important: There is no universal rule that says adjusted \(R^2\) and p-values must always give the same decision. They measure different aspects of model fit. In this practical, students should report both rather than relying on only one number.

A compact forward-selection table

The following code can be used to compare the models considered during forward selection.

forward_results <- data.frame(
  Model = c(
    "Ozone ~ 1",
    "Ozone ~ Solar.R",
    "Ozone ~ Wind",
    "Ozone ~ Temp",
    "Ozone ~ Temp + Solar.R",
    "Ozone ~ Temp + Wind",
    "Ozone ~ Solar.R + Wind + Temp"
  ),
  Adjusted_R2 = c(
    summary(m0)$adj.r.squared,
    summary(m_solar)$adj.r.squared,
    summary(m_wind)$adj.r.squared,
    summary(m_temp)$adj.r.squared,
    summary(m_temp_solar)$adj.r.squared,
    summary(m_temp_wind)$adj.r.squared,
    summary(full_model)$adj.r.squared
  )
)

forward_results

Use the regression summaries to add the relevant p-values to their interpretation.

Backward variable selection

Idea

Backward selection starts with the full model.

At each step:

fit the current model;
examine the p-values of the explanatory variables;
identify a variable that is not contributing sufficiently;
remove it;
compare the resulting model using adjusted \(R^2\) and p-values;
continue until no further removal is justified by the selected criteria.

Start with the full model:

summary(full_model)

## 
## Call:
## lm(formula = Ozone ~ Solar.R + Wind + Temp, data = aq_complete)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -40.485 -14.219  -3.551  10.097  95.619 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -64.34208   23.05472  -2.791  0.00623 ** 
## Solar.R       0.05982    0.02319   2.580  0.01124 *  
## Wind         -3.33359    0.65441  -5.094 1.52e-06 ***
## Temp          1.65209    0.25353   6.516 2.42e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.18 on 107 degrees of freedom
## Multiple R-squared:  0.6059, Adjusted R-squared:  0.5948 
## F-statistic: 54.83 on 3 and 107 DF,  p-value: < 2.2e-16

Look at the coefficient table:

coef(summary(full_model))

##                 Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept) -64.34207893 23.05472435 -2.790841 6.226638e-03
## Solar.R       0.05982059  0.02318647  2.579979 1.123664e-02
## Wind         -3.33359131  0.65440710 -5.094063 1.515934e-06
## Temp          1.65209291  0.25352979  6.516366 2.423506e-09

Suppose one variable has the largest p-value and its removal does not reduce adjusted \(R^2\) substantially. We can fit the reduced model.

For example:

model_reduced <- lm(Ozone ~ Temp + Wind,
                    data = aq_complete)

summary(model_reduced)

## 
## Call:
## lm(formula = Ozone ~ Temp + Wind, data = aq_complete)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -42.156 -13.216  -3.123  10.598  98.492 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -67.3220    23.6210  -2.850  0.00524 ** 
## Temp          1.8276     0.2506   7.294 5.29e-11 ***
## Wind         -3.2948     0.6711  -4.909 3.26e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.73 on 108 degrees of freedom
## Multiple R-squared:  0.5814, Adjusted R-squared:  0.5736 
## F-statistic: 74.99 on 2 and 108 DF,  p-value: < 2.2e-16

Compare adjusted \(R^2\):

c(
  Full_model = summary(full_model)$adj.r.squared,
  Reduced_model = summary(model_reduced)$adj.r.squared
)

##    Full_model Reduced_model 
##     0.5948449     0.5736257

Then compare the p-values in the reduced model:

coef(summary(model_reduced))

##               Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) -67.321953 23.6210451 -2.850084 5.236066e-03
## Temp          1.827554  0.2505520  7.294112 5.294630e-11
## Wind         -3.294839  0.6711482 -4.909257 3.261728e-06

If all remaining predictors have useful p-values and the reduced model has a suitable adjusted \(R^2\), we can stop.

Note: Backward selection starts with all candidate variables, whereas forward selection starts with no explanatory variables.

If both procedures produce the same final model, that is a useful observation to report.

If they produce different models, students should report the variables in each model and compare the adjusted \(R^2\) and p-values.

Final regression model

After variable selection, suppose the final model is stored as final_model.

For example:

final_model <- lm(Ozone ~ Temp + Wind,
                  data = aq_complete)

summary(final_model)

## 
## Call:
## lm(formula = Ozone ~ Temp + Wind, data = aq_complete)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -42.156 -13.216  -3.123  10.598  98.492 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -67.3220    23.6210  -2.850  0.00524 ** 
## Temp          1.8276     0.2506   7.294 5.29e-11 ***
## Wind         -3.2948     0.6711  -4.909 3.26e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.73 on 108 degrees of freedom
## Multiple R-squared:  0.5814, Adjusted R-squared:  0.5736 
## F-statistic: 74.99 on 2 and 108 DF,  p-value: < 2.2e-16

The fitted regression equation can be written using the estimated coefficients:

coef(final_model)

## (Intercept)        Temp        Wind 
##  -67.321953    1.827554   -3.294839

Students should write the fitted equation in the form

\[ \widehat{Ozone} = b_0 + b_1 Temp + b_2 Wind. \]

The exact numerical values should be taken from the R output.

Residuals

A residual is

\[ e_i = y_i - \hat{y}_i. \]

We can obtain the fitted values and residuals using:

fitted_values <- fitted(final_model)
residual_values <- residuals(final_model)

head(
  data.frame(
    Fitted = fitted_values,
    Residual = residual_values
  )
)

Residual plot

A basic residual plot is a plot of residuals against fitted values.

plot(fitted(final_model),
     residuals(final_model),
     xlab = "Fitted values",
     ylab = "Residuals",
     main = "Residuals vs Fitted Values")

abline(h = 0, lty = 2)

How to read the residual plot

For a reasonable linear regression model, we generally look for:

residuals scattered around zero;
no obvious curved pattern;
roughly constant spread of residuals across fitted values;
no single observation that appears extremely unusual.

A visible pattern in the residual plot may suggest that the simple linear model is not adequately describing the data.

Note: The residual plot is a diagnostic tool. It should not be used as a variable-selection criterion in this exercise because the syllabus criteria specified here are adjusted \(R^2\) and p-values.

A useful diagnostic version

R also provides standard diagnostic plots.

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

par(mfrow = c(1, 1))

For this practical, however, the Residuals vs Fitted plot is the main plot students need to interpret.

Suggested reporting format

Students can report their analysis in the following order:

Data: State that the built-in airquality data set was used.
Variables: State that Ozone is the response and Solar.R, Wind, and Temp are candidate explanatory variables.
Missing values: State that incomplete observations were removed using na.omit().
Full model: Report the fitted multiple regression model and its adjusted \(R^2\).
Forward selection: Show the variables considered at each step and explain the decision using adjusted \(R^2\) and p-values.
Backward selection: Start from the full model and show which variable(s) were considered for removal.
Final model: Give the final regression equation and interpret the coefficients.
Residual plot: Show the residual-versus-fitted plot and briefly comment on its pattern.

Important distinction: selection criterion versus diagnostic

In this exercise, keep the following distinction clear:

Purpose	Quantity
Compare model fit while accounting for model size	Adjusted \(R^2\)
Assess evidence against a zero coefficient	p-value
Examine model assumptions/patterns	Residual plot

The residual plot is therefore not being used to decide which variable enters or leaves the model. It is used after fitting the model to assess whether the fitted linear model shows obvious problems.

Practice questions

How many observations are available before removing missing values?
How many observations remain after na.omit()?
Which single predictor gives the largest adjusted \(R^2\)?
Which variable is selected first in forward selection?
What happens to adjusted \(R^2\) when a second variable is added?
In the full model, which predictor has the largest p-value?
What happens to adjusted \(R^2\) when that predictor is removed?
Do forward and backward selection lead to the same final model?
Write the final fitted regression equation.
What does the residual plot suggest about the adequacy of the linear model?

Reproducibility

The analysis uses only the built-in airquality data set, so no external data file is required.

Multiple Linear Regression and Variable Selection in R: Guide line

Dr Kuntal Bakuli

Learning objectives

The data

Missing observations

Note: Missing-data handling should be decided before interpreting the regression results. Here we use complete-case analysis because it is simple and appropriate for demonstrating the regression and variable-selection procedures.

Exploratory look at the variables

Full multiple regression model

Adjusted \(R^2\)

Forward variable selection

Idea

Step 0: Intercept-only model

Step 1: Consider each variable separately

Step 2: Add a second variable

Step 3: Consider the remaining variable

A compact forward-selection table

Backward variable selection

Idea

Final regression model

Residuals

Residual plot

How to read the residual plot

A useful diagnostic version

Suggested reporting format

Important distinction: selection criterion versus diagnostic

Practice questions

Reproducibility