# Using diagnostic plots in order to decide variables to set up multiple linear regression - R

I'm working with a baseball dataset using R. The dataset contains the baseball season records for teams between the years 1871 and 2016. One of the columns has the number of wins for the season and all the other variables have other information like the number of stolen bases or number of home runs.

I'm trying to create a multiple linear regression. The dependent variable will be the number of wins for the season and the predictor variables will be some of the other columns in the dataset.

My questions are - How do I decide on which independent variables to use and also are diagnostic plots useful in deciding which columns to use as the independent variables in the multiple linear regression?

Here is a dput sample of my dataset:

dput(head(moneyball_training_data))

structure(list(INDEX = 1:6, TARGET_WINS = c(39L, 70L, 86L, 70L,
82L, 75L), TEAM_BATTING_H = c(1445L, 1339L, 1377L, 1387L, 1297L,
1279L), TEAM_BATTING_2B = c(194L, 219L, 232L, 209L, 186L, 200L
), TEAM_BATTING_3B = c(39L, 22L, 35L, 38L, 27L, 36L), TEAM_BATTING_HR = c(13L,
190L, 137L, 96L, 102L, 92L), TEAM_BATTING_BB = c(143L, 685L,
602L, 451L, 472L, 443L), TEAM_BATTING_SO = c(842L, 1075L, 917L,
922L, 920L, 973L), TEAM_BASERUN_SB = c(NA, 37L, 46L, 43L, 49L,
107L), TEAM_BASERUN_CS = c(NA, 28L, 27L, 30L, 39L, 59L), TEAM_BATTING_HBP = c(NA_integer_,
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_
), TEAM_PITCHING_H = c(9364L, 1347L, 1377L, 1396L, 1297L, 1279L
), TEAM_PITCHING_HR = c(84L, 191L, 137L, 97L, 102L, 92L), TEAM_PITCHING_BB = c(927L,
689L, 602L, 454L, 472L, 443L), TEAM_PITCHING_SO = c(5456L, 1082L,
917L, 928L, 920L, 973L), TEAM_FIELDING_E = c(1011L, 193L, 175L,
164L, 138L, 123L), TEAM_FIELDING_DP = c(NA, 155L, 153L, 156L,
168L, 149L)), row.names = c(NA, 6L), class = "data.frame")


For every single variable in the dataset I created a simple linear model - just one independent variable and one dependent variable.

Here you can also see the code I created to make diagnostic plots, specifically on one hand the fitted values vs. the residuals and on the other hand a histogram of the residuals.

Am I on the right track? If I go through all these diagnostic plots in order to check for homoscedasticity, etc, and only include the columns that have a good linear relationship to the dependent variable and combine them together, did I do it right in setting up my multiple linear regression?

For example, I see that when I create a linear model of number of stolen bases on the x and number of wins on the y, and I create diagnostic plots and see that the histogram of residuals for this linear model looks normally distributed... then I can include the variable 'number of stolen bases' in my multiple linear regression model that I'm creating... and then look to other variables to see what else I can include?

TEAM_BATTING_H_linear_model <- lm(TARGET_WINS ~ TEAM_BATTING_H, moneyball_training_data)

ggplot(data = TEAM_BATTING_H_linear_model, aes(x = .fitted, y = .resid)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed") +
xlab("Fitted values") +
ylab("Residuals") +
labs(title = "", subtitle = "Scatterplot - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

ggplot(data = TEAM_BATTING_H_linear_model, aes(x = .resid)) +
geom_histogram() +
xlab("Residuals") +
labs(title = "", subtitle = "Histogram - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

TEAM_BATTING_2B_linear_model <- lm(TARGET_WINS ~ TEAM_BATTING_2B, moneyball_training_data)

ggplot(data = TEAM_BATTING_2B_linear_model, aes(x = .fitted, y = .resid)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed") +
xlab("Fitted values") +
ylab("Residuals") +
labs(title = "", subtitle = "Scatterplot - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

ggplot(data = TEAM_BATTING_2B_linear_model, aes(x = .resid)) +
geom_histogram() +
xlab("Residuals") +
labs(title = "", subtitle = "Histogram - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

TEAM_BATTING_3B_linear_model <- lm(TARGET_WINS ~ TEAM_BATTING_3B, moneyball_training_data)

ggplot(data = TEAM_BATTING_3B_linear_model, aes(x = .fitted, y = .resid)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed") +
xlab("Fitted values") +
ylab("Residuals") +
labs(title = "", subtitle = "Scatterplot - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

ggplot(data = TEAM_BATTING_3B_linear_model, aes(x = .resid)) +
geom_histogram() +
xlab("Residuals") +
labs(title = "", subtitle = "Histogram - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

TEAM_BATTING_HR_linear_model <- lm(TARGET_WINS ~ TEAM_BATTING_HR, moneyball_training_data)

ggplot(data = TEAM_BATTING_HR_linear_model, aes(x = .fitted, y = .resid)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed") +
xlab("Fitted values") +
ylab("Residuals") +
labs(title = "", subtitle = "Scatterplot - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

ggplot(data = TEAM_BATTING_HR_linear_model, aes(x = .resid)) +
geom_histogram() +
xlab("Residuals") +
labs(title = "", subtitle = "Histogram - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

TEAM_BATTING_BB_linear_model <- lm(TARGET_WINS ~ TEAM_BATTING_BB, moneyball_training_data)

ggplot(data = TEAM_BATTING_BB_linear_model, aes(x = .fitted, y = .resid)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed") +
xlab("Fitted values") +
ylab("Residuals") +
labs(title = "", subtitle = "Scatterplot - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

ggplot(data = TEAM_BATTING_BB_linear_model, aes(x = .resid)) +
geom_histogram() +
xlab("Residuals") +
labs(title = "", subtitle = "Histogram - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

TEAM_BATTING_SO_linear_model <- lm(TARGET_WINS ~ TEAM_BATTING_SO, moneyball_training_data)

ggplot(data = TEAM_BATTING_SO_linear_model, aes(x = .fitted, y = .resid)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed") +
xlab("Fitted values") +
ylab("Residuals") +
labs(title = "", subtitle = "Scatterplot - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

ggplot(data = TEAM_BATTING_SO_linear_model, aes(x = .resid)) +
geom_histogram() +
xlab("Residuals") +
labs(title = "", subtitle = "Histogram - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

TEAM_BASERUN_SB_linear_model <- lm(TARGET_WINS ~ TEAM_BASERUN_SB, moneyball_training_data)

ggplot(data = TEAM_BASERUN_SB_linear_model, aes(x = .fitted, y = .resid)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed") +
xlab("Fitted values") +
ylab("Residuals") +
labs(title = "", subtitle = "Scatterplot - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

ggplot(data = TEAM_BASERUN_SB_linear_model, aes(x = .resid)) +
geom_histogram() +
xlab("Residuals") +
labs(title = "", subtitle = "Histogram - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

TEAM_BASERUN_CS_linear_model <- lm(TARGET_WINS ~ TEAM_BASERUN_CS, moneyball_training_data)

ggplot(data = TEAM_BASERUN_CS_linear_model, aes(x = .fitted, y = .resid)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed") +
xlab("Fitted values") +
ylab("Residuals") +
labs(title = "", subtitle = "Scatterplot - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

ggplot(data = TEAM_BASERUN_CS_linear_model, aes(x = .resid)) +
geom_histogram() +
xlab("Residuals") +
labs(title = "", subtitle = "Histogram - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

TEAM_BATTING_HBP_linear_model <- lm(TARGET_WINS ~ TEAM_BATTING_HBP, moneyball_training_data)

ggplot(data = TEAM_BATTING_HBP_linear_model, aes(x = .fitted, y = .resid)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed") +
xlab("Fitted values") +
ylab("Residuals") +
labs(title = "", subtitle = "Scatterplot - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

ggplot(data = TEAM_BATTING_HBP_linear_model, aes(x = .resid)) +
geom_histogram() +
xlab("Residuals") +
labs(title = "", subtitle = "Histogram - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

TEAM_PITCHING_H_linear_model <- lm(TARGET_WINS ~ TEAM_PITCHING_H, moneyball_training_data)

ggplot(data = TEAM_PITCHING_H_linear_model, aes(x = .fitted, y = .resid)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed") +
xlab("Fitted values") +
ylab("Residuals") +
labs(title = "", subtitle = "Scatterplot - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

ggplot(data = TEAM_PITCHING_H_linear_model, aes(x = .resid)) +
geom_histogram() +
xlab("Residuals") +
labs(title = "", subtitle = "Histogram - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

TEAM_PITCHING_HR_linear_model <- lm(TARGET_WINS ~ TEAM_PITCHING_HR, moneyball_training_data)

ggplot(data = TEAM_PITCHING_HR_linear_model, aes(x = .fitted, y = .resid)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed") +
xlab("Fitted values") +
ylab("Residuals") +
labs(title = "", subtitle = "Scatterplot - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

ggplot(data = TEAM_PITCHING_HR_linear_model, aes(x = .resid)) +
geom_histogram() +
xlab("Residuals") +
labs(title = "", subtitle = "Histogram - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

TEAM_PITCHING_BB_linear_model <- lm(TARGET_WINS ~ TEAM_PITCHING_BB, moneyball_training_data)

ggplot(data = TEAM_PITCHING_BB_linear_model, aes(x = .fitted, y = .resid)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed") +
xlab("Fitted values") +
ylab("Residuals") +
labs(title = "", subtitle = "Scatterplot - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

ggplot(data = TEAM_PITCHING_BB_linear_model, aes(x = .resid)) +
geom_histogram() +
xlab("Residuals") +
labs(title = "", subtitle = "Histogram - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

TEAM_PITCHING_SO_linear_model <- lm(TARGET_WINS ~ TEAM_PITCHING_SO, moneyball_training_data)

ggplot(data = TEAM_PITCHING_SO_linear_model, aes(x = .fitted, y = .resid)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed") +
xlab("Fitted values") +
ylab("Residuals") +
labs(title = "", subtitle = "Scatterplot - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

ggplot(data = TEAM_PITCHING_SO_linear_model, aes(x = .resid)) +
geom_histogram() +
xlab("Residuals") +
labs(title = "", subtitle = "Histogram - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

TEAM_FIELDING_E_linear_model <- lm(TARGET_WINS ~ TEAM_FIELDING_E, moneyball_training_data)

ggplot(data = TEAM_FIELDING_E_linear_model, aes(x = .fitted, y = .resid)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed") +
xlab("Fitted values") +
ylab("Residuals") +
labs(title = "", subtitle = "Scatterplot - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

ggplot(data = TEAM_FIELDING_E_linear_model, aes(x = .resid)) +
geom_histogram() +
xlab("Residuals") +
labs(title = "", subtitle = "Histogram - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

TEAM_FIELDING_DP_linear_model <- lm(TARGET_WINS ~ TEAM_FIELDING_DP, moneyball_training_data)

ggplot(data = TEAM_FIELDING_DP_linear_model, aes(x = .fitted, y = .resid)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed") +
xlab("Fitted values") +
ylab("Residuals") +
labs(title = "", subtitle = "Scatterplot - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")

ggplot(data = TEAM_FIELDING_DP_linear_model, aes(x = .resid)) +
geom_histogram() +
xlab("Residuals") +
labs(title = "", subtitle = "Histogram - Linear Model") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")


Here are two sample images:

Your question pertains to the "selection of variables" problem in regression, for which there is a large statistical literature that is quite a lot to summarise. There are a few issues to consider here, so I will go through them slowly.

A preliminary note on training and testing: One of the difficulties that one encounters in conducting selection of variables is that many of the methods for doing this "optimise" over sets of available variables. One problem with this is that most testing methods constructed for regression analysis (e.g., T-tests, F-tests, etc.) assume that no optimisation of selected variables has occurred, which means that if you apply these tests then they will be biased towards showing relationships between the selected explanatory variables and the response variable.

In order to avoid this, you can use one of two general approaches. One approach is to use a model form like LASSO regression that has a built-in selection mechanism as part of the model. These types of models build in a selection mechanism by imposing a penalty for non-zero regression coefficients in the model. The subsequent tests for the model take account of this in-built selection mechanism and so they are not biased by it. Alternatively, one can use regular regression methods (that do not have in-built selection mechanisms) but split your data into a training set used to fit the models and perform selection of variables, and then a testing set used to formally test the relationship between the selected explanatory variables and the response variable.

In any case, whatever method you use, you will need to consider whether you are performing your selection of variables "within the model" and if not, you will need to consider using a separate training and testing phase, which generally involves some random partition of your data for this purpose.

Use one multiple regression, not lots of simple regressions: It is not useful in these cases to construct simple regression models for the individual variables. That approach ignores the colinearity between the explanatory variables, and it does not assist in understanding their contribution to the multiple regression model. The ultimate goal in regression analysis is to model the response as a function of all relevant explanatory variables, so you should be using a multiple regression model for this.

From what you have written, it appears that your goal in using the simple regression models was to construct a plot that will show the relationship between each explanatory variable and the response variable. Fortunately, this can be done in a multiple regression model by using added variable plots. This plot shows the relationship between these variables conditional on having removed the effect of all other explanatory variables in the model. Consequently, it shows you the conditional relationship of interest. (Note that the line-of-best-fit in the added variable plot uses the estimated slope coefficient from the multiple regression model, which is not true when you plot the marginal relationships using simplie regressions.)

You can construct added variable plots using the avPlots function in the car package in R. This function takes an input model and automatically prints all its added variable plots. It is also possible to extract the information for these plots and plot them using ggplot2 if you prefer the latter graphical system (see instructions here).

Some methods for selection of variables: There are a few broad methods still in use for selection of variables. The simplest method is to avoid it all together, using the "kitchen sink" approach where you fit a regression model will all the available variables in it. Under this approach you eschew the selection of variables analysis and just use all the explanatory variables that are plausible predictors a priori. If you do this then you may end up including some explanatory variables that are not related to the response variable. For models constructed for predictive analysis, the only real cost of this is loss of a few degrees-of-freedom, which is usually not a big problem. For models constructed for expalantory or causal reasons, this method can be problematic. Some analysts do not like this approach, but I think it is quite reasonable in many cases.

Another way to do selection of variables in your regression analysis is to use a model form with an in-built selection mechanism, such as LASSO regression. This is equivalent to a Bayesian analysis where you give a prior to your regression coefficients, with some non-zero probability that they are zero. Estimation in LASSO regression will estimate coefficients to be zero if there is sufficiently weak evidence of a relationship between the corresponding explanatory variable and the response variable. The model has an adjustable "penalty" function that allows you to set the appropriate bar for how much evidence you require for includion of a variable in the model.

Finally, another common way to do selection of variables is to use standard regression models (i.e., those without an in-built selection mechanism) and then select variables to include based on optimisation of some kind of "goodness-of-fit" statistic. Commonly used statistics for this purpose include the Akaike information criterion, the Bayesian information criterion or the adjusted coefficient of variation. If your model space is not too large (i.e., you do not have too many available explanatory variables) it is usually possible to search over all possible models. Since this involves model selection using optimisation over a large number of models, it is important to separate training of the model from formal testing of relationships in the selected model.

• Hi Ben, thanks for the thorough response here. The part about simple regressions ignoring the collinearity and therefore not being useful was especially to the point of my question. I'm wondering about a couple of other clarifying points. 1) You mention AIC as a goodness-of-fit statistic. When I look at the summary of a linear model in R, I see p-values listed for different independent variables. Would p-value here also be a goodness of fit statistic? Mar 1, 2021 at 7:32
• and 2) if I use a kitchen sink approach in deciding variables to include in the model or I go through different goodness-of-fit statistics or an a-priori thinking / eyeballing it -- is there a standard workflow method to show in an r markdown document "my thinking"? I could present "here's my multiple linear regression model" but how would I show "this is the best possible model based on permutations of independent variables" Mar 1, 2021 at 7:34

Someone else with more expertise may speak more clearly to this, but I thought I'd note a few considerations (not really hazarding an answer). I'm not sure how useful bivariable diagnostic plots are for determining whether to include a variable in a multivariable model. I'd be interested to hear opinions on that. But I do think that there are general limitations in being guided by bivariable relationships for choosing which variables to include in a larger model. Two key problems with this are:

1. Masked relationships. Predictors can have a relationship with the outcome which is not visible in a bivariable analysis. Say your 'team batting H' variable is weakly negatively correlated with the outcome, while number of stolen bases variable is weakly positively correlated with the outcome. Both of these relationships might appear weak when inspected in a plot, and thus you may choose to not include them in the multivariable model. If there is a positive correlation between these two predictors, however, those bivariable plots will underestimate the true effect of each predictor on the outcome. So by not including these predictors in the multivariable model, you have potentially missed two variables strongly related to the outcome.

2. Spurious associations. The inverse of (1). Predictors which appear to be strongly related to the outcome in bivariable plots, may in fact have little relationship with the outcome in the multivariable model.

p.s. if you're concerned about non-linearity, which is a reasonable expectation, why not include spline terms for your continuous predictors?

• Hi Lachlan, thanks for the note here. I see that if I analyze these bivariate relationships separately, I may totally miss the idea that these variable can be collinear, as I understand it, and there could be an interaction taking place. Based on both of your responses, I think I should only use diagnostic plots after creating a linear model to assess the model as a whole. Your suggestion about spline terms seems like a good one but I don't know much about them :/ Mar 1, 2021 at 9:19