You were introduced to multiple regression in PSYC 122 and is an extension of simple regression. Now, we have more than one predictor variable for our outcome variable. We follow the same process as a simple regression, and test the variances explained by the model’s group of predictors with one \(F\) value in an “omnibus” test. This is unlike ANOVA where we can calculate an \(F\) value for each factor. For individual predictors and their significance on the outcome variable, we can calculate \(t\) values.
A multiple regression analysis finds the best possible set of values that predict values of the outcome variable given the observed values.
Remember the equation for simple regression:
\[ Y = b_0 + b_1*X+e \] You may notice that the intercept term here is \(b_0\) where in the simple regression equation it was \(a\). They are equivalent, but I have noted the differences in convention here so that if you come across equations, you are aware that different people use different notation conventions. There is usually an explanation of the notation when it is first introduced - so read carefully and you should be ok.
Now, the multiple regression equation simply adds further slope terms (extra \(b\)s) - one for each extra predictor variable
So a two predictor regression analysis equation becomes:
\[ Y_i = b_0 + b_1 * X_1 + b_2*X_2 + e_i \] As before:
All the relations of the terms in the equation remain the same as with simple regression, there are literally just more predictors. There is a little more to the interpretation of an individual beta weight, in the context of its neighbouring predictors, from that of simple regression (explained later).
You may be relieved to hear that the computation by hand becomes intractable at this point. We use the computer and functions to do the heavy lifting of model estimation and prediction for us.
This is a truncated set of values retrieved from Margriet’s first lecture last year. We used these to review correlation and simple regression. It makes sense to use them again, because we know the results of at least one of the predictors, so the addition of extra predictors may demonstrate that multiple regression is an extension of simple regression.
d <- read_csv("Covariance_Data.csv") # calling my dataframe d to save typing!Rows: 13 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (5): Participant_ID, Performance, IQ, Motivation, Social_Support
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.That is not to say that is all multiple regression is. There are practices that are needed with the presentation of more than one predictor but for now we will simple perform a regression with more than one predictor. The future lectures will visit these practices to fully inform the regression workflow.
head(d) # inspect the first 6 observations of `d`# A tibble: 6 × 5
Participant_ID Performance IQ Motivation Social_Support
<dbl> <dbl> <dbl> <dbl> <dbl>
1 1 85 109 89 73
2 2 84 106 84 80
3 3 87 125 59 67
4 4 69 84 60 58
5 5 81 89 60 67
6 6 71 109 62 75We know some information about this dataset from looking at the first few rows, however, without a supporting article or a codebook, there is little to tell us about the driving research question or the hypothesis for this data. So we cannot know for sure what the primary variables are and what the supporting /auxiliary / third variables may be.
We need more than two variables to complete an example. We know from the simple regression performed earlier the IQ has very little relationship with Performance, so we could just go for it and add Motivation and Social_Support.
Performance can still be \(Y\), with all other variables acting as \(X\) predictors in this model.
Here is a repeat of the simple regression linear model command:
summary(m1 <- lm(Performance ~ IQ, d))
Call:
lm(formula = Performance ~ IQ, data = d)
Residuals:
Min 1Q Median 3Q Max
-12.281 -3.755 -1.439 6.363 12.363
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 79.03336 17.87633 4.421 0.00103 **
IQ -0.01318 0.16492 -0.080 0.93775
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.727 on 11 degrees of freedom
Multiple R-squared: 0.00058, Adjusted R-squared: -0.09028
F-statistic: 0.006384 on 1 and 11 DF, p-value: 0.9378The multiple regression command is simply an addition of the other two predictor terms. I have called this m3 because it has three predictor terms (explained later).
m3 <- lm(Performance ~ IQ + Motivation + Social_Support, d)
summary(m3)
Call:
lm(formula = Performance ~ IQ + Motivation + Social_Support,
data = d)
Residuals:
Min 1Q Median 3Q Max
-9.7971 -6.9194 -0.5667 3.5231 14.4673
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 61.893215 27.601610 2.242 0.0516 .
IQ 0.004905 0.182974 0.027 0.9792
Motivation 0.062248 0.188307 0.331 0.7485
Social_Support 0.170331 0.281606 0.605 0.5602
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.163 on 9 degrees of freedom
Multiple R-squared: 0.08722, Adjusted R-squared: -0.217
F-statistic: 0.2867 on 3 and 9 DF, p-value: 0.834Let’s go through the summary output.
The Call of the model tells us the formula that we coded.
The Residuals give us some idea of whether we have a normal distribution of the error terms. Remember that a normal distribution is symmetric, which means that the min and max values should mirror each other, with the Min being negative and the max being positive. Here we have -9.79 for the Min values but 14.47 for the Max values. Maybe we have a problem with this assumption - we can check with some boxplots
The following two plots are identical in both visualisation and interpretation - except for the code used to call the residual values from the model object
boxplot(m3[['residuals']], # uses double brackets to identify the residuals from the model object
main='Boxplot: Residuals (index code method)', # main title label
ylab='residual value') # y axis label
boxplot(m3$residuals, # uses dollar sign to identify the residuals from the model object
main='Boxplot: Residuals (variable name code method)', # main title label
ylab='residual value') # y axis label
A shorter bottom whisker and a longer top whisker indicated a skewed distribution, which is indicative of a non-symettric distribution.
The Coefficients block: The columns:
Social_Support is a stronger predictor than the other predictors given its magnitude but remember that the units here for each predictor are on different scales (that we have no information about either). We need to standardise the predictors if we are to be able to make such comparisons (explained later).First, lets take a look at the difference of the intercept and IQ coefficient between the simple regression model and the multiple regression model. The intercept and IQ beta coefficients from the simple regression have changed in the context of a larger set of predictors:
IQ predictor is now positive (it was negative) but is of a smaller magnitude. The interpretation of this predictor, given its magnitude, is probably the same in m3 as m1. A one unit increase in IQ indicates a 0.005 increase in Performance Scores, holding Motivation and Social_Support predictors at a constant level. The t-value and associated p value remain non-significant.The added predictors of Motivation and Social_Support
Lets reprint the summary to save scrolling but also save the summary to an object so that we can call the values while interpreting and reporting the model
(m3_sum <- summary(m3))
Call:
lm(formula = Performance ~ IQ + Motivation + Social_Support,
data = d)
Residuals:
Min 1Q Median 3Q Max
-9.7971 -6.9194 -0.5667 3.5231 14.4673
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 61.893215 27.601610 2.242 0.0516 .
IQ 0.004905 0.182974 0.027 0.9792
Motivation 0.062248 0.188307 0.331 0.7485
Social_Support 0.170331 0.281606 0.605 0.5602
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.163 on 9 degrees of freedom
Multiple R-squared: 0.08722, Adjusted R-squared: -0.217
F-statistic: 0.2867 on 3 and 9 DF, p-value: 0.834When the \(F_b\) is not significant then \(b\) = 0 may be true. The data do not demonstrate any tendency for larger values of \(X\) to be associated with smaller or larger values of \(Y\).
None of the predictors show any significant ability to predict variance in Job Performance (IQ: t = 0.03, p = 0.98; Motivation: t = 0.33, p = 0.75 and Social Support: t = 0.6, p = 0.56). For this sample, the predictors and the job performance appear to be causally unrelated.
The model is not significant either (\(F{_(}{_3,_9}{_)}\) = 0.29, p = 0.834). \(F_{critical}\) for df = (3, 9) at the .05 level = 3.86. This means that \(F_{observed}\) is not larger than \(F_{critical}\) so we fail to reject the null hypothesis.
Plotting an individual effect (or lack of it) is very simple if we are not too worried about being pretty.
plot(predictorEffect("IQ", m3))
plot(predictorEffect("Motivation", m3))
plot(predictorEffect("Social_Support", m3))
Plotting it in this way, also allows us, with a simple additional argument to the plot call, to see the partial residuals - thereby checking that we have residual values that are normally distributed - I think the pink line probably tells us that we do not! We should be looking for an alternative way to analyse these two variables….or some transformations!
plot(predictorEffect("IQ", m3, partial.residuals = TRUE))
plot(predictorEffect("Motivation", m3, partial.residuals = TRUE))
plot(predictorEffect("Social_Support", m3, partial.residuals = TRUE))
Every regression model is built upon the following assumptions:
par(mfrow = c(2, 2)) # display plots in a 2 x 2 panel
plot(m3) # plot diagnostic plots for m3
Residuals vs Fitted: “fitted” here means the predicted values. We want the pink line to align pretty closely with the horizontal dashed line. Comparing this plot with that from the simple regression, this plot looks better. Take note of observation 13 - that was also labelled in the simple regression plots.
Normal Q-Q: If the residual points (open circles) follow the dashed line, you can assume the residuals are normally distributed
Scale-Location: This is checking for constant variance in the residuals - not much here. A good indication would be a horizontal pink line with equally spread points. Our graph is not good.
Residuals vs Leverage - are there any points that are having a large influence on the regression results. They will be numbered and you can then inspect them in your data file. Observations that show standardised residuals (see the table above) above 3 would be problematic. As would observations of a hat value above \(2(p+1)/n\) where \(p\) = is the number of predictors and \(n\) = is the number of observations
p <- 3
N <- nrow(d)
(m3_hat <- (2*(p+1))/N) # model hat value[1] 0.6153846And we can collect some diagnostic measures for each data point to help with assumption checking using the augment() function from the broom package:
m3_metrics <- augment(m3)
print.data.frame(m3_metrics) # this code makes all the columns in the tibble print out at once; it isn't pretty but it works for our purposes Performance IQ Motivation Social_Support .fitted .resid .hat
1 85 109 89 73 80.40217 4.59783441 0.2905074
2 84 106 84 80 81.26853 2.73147213 0.3098207
3 87 125 59 67 77.59123 9.40877367 0.3173281
4 69 84 60 58 75.91937 -6.91936802 0.4416438
5 81 89 60 67 77.47688 3.52312335 0.2418778
6 71 109 62 75 79.06213 -8.06213238 0.2457062
7 76 121 67 55 76.02561 -0.02561337 0.2303470
8 77 102 44 73 77.56667 -0.56666766 0.5160834
9 76 111 68 60 76.89046 -0.89046266 0.1118330
10 90 106 63 54 75.53271 14.46729219 0.2077073
11 74 107 93 75 80.98201 -6.98200905 0.3774818
12 74 97 52 58 75.48516 -1.48515549 0.2219124
13 65 133 60 50 74.79709 -9.79708713 0.4877509
.sigma .cooksd .std.resid
1 8.440766 4.576854e-02 0.668665424
2 8.580197 1.820451e-02 0.402759032
3 7.665622 2.261255e-01 1.394943074
4 8.015772 2.544376e-01 -1.134332135
5 8.539583 1.959622e-02 0.495664399
6 8.012472 1.053014e-01 -1.137127455
7 8.658576 9.570291e-07 -0.003576422
8 8.653791 2.654818e-03 -0.099786698
9 8.652136 4.217075e-04 -0.115743839
10 6.476835 2.598088e-01 1.991012125
11 8.073567 1.781360e-01 -1.084011442
12 8.638096 3.032965e-03 -0.206246770
13 7.179771 6.693121e-01 -1.676817738Columns 1 - 4 are our variables. fitted and residuals are the predicted values of Y and the error values. The .hat and .cooksd columns we use below with some explanation.
Looking down the .hat column, and comparing it to the model hat value of 0.62, none of the datapoints above have hat values larger than the model hat threshold value.
Checking for observations that are influential follows a similar pattern: observations that exceed the Cook’s distance value = \(4/(n-p-1)\) are likely to have high influence and the regression results may change if you exclude them. In the presence of such observations that exceed Cook’s distance, unless you know the observation are errors, you probably need to estimate the model without the observations and report both sets of results.
(m3_Cooks <- 4/(N-p-1)) # model Cook's distance value[1] 0.4444444Looking down the .cooksd column, and comparing it to the model Cooks Distance value of 0.44, observation no. 13 is looking influential. Remember the skewed distribution of the residuals with the long whisker… Time to inspect, possibly exclude 13 and rerun! As before however, the major problem here is small sample size.
Hopefully by now, you realise that writing your results section is quite formulaic. A results section comprises some descriptive statistics and inferential statistics. Descriptive statistics are number of observations and means and SDs of the variables - the same as any analysis method. The inferential statistics are reported at the model and predictor level. Because there are only three predictors in this model, you could present the descriptive statistics in the text with the predictor statistics.
In multiple regression, you may also want to present the relationships between the predictor variables as part of your results section. This is informative, as high correlations between predictor variables can suggest different modelling strategies, such as variable selection or factor analysis (see week 15).
# some preparation needed for a correlation matrix
# remove the participant ID column
d_cor <- d %>%
select(-Participant_ID)
res.cor <- correlate(d_cor, diagonal = 1)Correlation computed with
• Method: 'pearson'
• Missing treated using: 'pairwise.complete.obs'(res.cor2 <- res.cor %>% shave())# A tibble: 4 × 5
term Performance IQ Motivation Social_Support
<chr> <dbl> <dbl> <dbl> <dbl>
1 Performance 1 NA NA NA
2 IQ -0.0241 1 NA NA
3 Motivation 0.219 0.111 1 NA
4 Social_Support 0.274 -0.208 0.436 1This is a table of pairs of correlations. It shows how one predictor is associated with another predictor in the model. The first column is the row labels. From the second column onward, we see Pearson’s r correlations. The Performance column shows the correlation of each predictor with the outcome variable. If you recall, the simple correlation between IQ and Performance was -0.02. Here you see that when paired with Performance, Motivation and Social-Support seem to share much more variance with Performance - however when modelled together, neither of them show such strength, nor are they significant predictors.
In the final three columns, we see the correlational relationships that exist between pairs of predictors. Look at Motivation and Social Support - their relationship is positive, as one increases so does the other, and quite strong. Look at the relationship between IQ and Social Support: negative, as one increases the other decreases.
Above is the tibble; reporting standards are better approximated by going a little further.
kable(res.cor2, booktabs = T, digits = 2)| term | Performance | IQ | Motivation | Social_Support |
|---|---|---|---|---|
| Performance | 1.00 | NA | NA | NA |
| IQ | -0.02 | 1.00 | NA | NA |
| Motivation | 0.22 | 0.11 | 1.00 | NA |
| Social_Support | 0.27 | -0.21 | 0.44 | 1 |
We performed a multiple regression model with Job Performance Scores as our outcome variable and IQ, Motivation and Social Support as predictor variables. The bivariate correlation relationships for all variables in shown in table 1. The model was non-significant (F(3, 9) = 0.29, p = .834) and explained 0% of the variance on Job Performance Scores. IQ (mean = 107.6, SD = 13.5) was a non-significant predictor for Job Performance (b = 0.005, t = 0.03, p = .98), as was Motivation (mean = 66.2, SD = 14.3,b = 0.06, t = 0.33, p = .75), and Social Support (mean = 65, SD = 9.7,b = 0.17, t = 0.61, p = .56). Since neither the model nor any of the predictors were significant, we do not present any effects plots here.
There are lots of things to do as we run a multiple regression model - both before, while and after building the model, which will be the focus of the future lectures. We need to look at:
Your pre-lab activity is running through these three markdown documents for Correlation, Simple Regression and Multiple Regression review and testing your understanding with the WBA 1 Quiz. The Quiz runs through all the theory behind correlation and regression and is explicitly designed to make you agile with the r code. For some of you it will be intensely frustrating. Like all data analysis, start early and walk away from it when you are at your limit, returning to do a little more later on.