Predictions based on multiple variables

What if we wanted to know the average amount of working hours per week for females with 18 years of education? In the same way as before, when we predicted the average hours of work per week given a specific value of years of education, we could simply add the specification of the gender variable to our code. Importantly: you must provide values for all variables in the model for the calculations to make sense, so predict() will throw an error if you have not.

# Creating a new dataframe with the correct specifications
woman_educ <- data.frame(yeareduc = 18,
                         gender = "Female")

# Predicting the mean outcome value given the specifications above
predict(m2, newdata = woman_educ)

##        1 
## 37.06204

As we can see from the output, the expected mean amount of working hours per week for a woman with 18 years of completed education is 37.06204. You can also see that all we did to obtain this information, was simply to add another argument inside the data.frame() function where we specified which gender we wanted to predict at.

Checking linearity

The linearity assumption can be checked either by plotting a scatterplot and checking the correlation between the explanatory variable and the outcome variable (continuous variables), or you can use a residuals vs fitted plot. This plot shows whether there are non-linear patterns of the residuals in your model. If the plot shows equally spread residuals around a horizontal line without any distinct patterns, it is a good indication that you don’t have any non-linear relationship.

Let’s illustrate this using our multiple regression model from before, m2. To make a quick-and-dirty plot, we can use the function plot().

For linear models, there are six automatic plots, but in this case we are only interested in the first of these. Thus, we specify which = 1 as well:

# Plotting a residuals vs fitted plot
plot(m2, which = 1)

The plot shows no clear non-linear patterns, e.g. there is no curvature on the line, and the residuals are pretty equally spread along a horizontal line. In this case, this is a good indication that we do not have any issues with non-linearity in our model.

ggplot(m2, aes(x = .fitted, y = .resid)) +
  geom_point()+
  geom_smooth(method = "lm", se = F)

## `geom_smooth()` using formula 'y ~ x'

Checking normality of residuals

The residuals should also be normally distributed. We can check this assumption by plotting a histogram of the residuals and/or creating a quantile-quantile-plot.

# Plotting a histogram of residuals 
ggplot(m2, aes(x = .resid)) +
  geom_histogram(bins = 20, fill = "steelblue", color = "black") +
  labs(title = "Histograms of residuals", x = "Residuals", y = "Frequency")

# Plotting a qq-plot
plot(m2, which = 2)

ggplot(m2, aes(sample = .resid)) +
  geom_qq() +
  geom_qq_line()

Checking homoscedasticity

Another assumption about the residuals is the assumption of homoscedasticity, i.e that the residuals should have constant variance at all values of the explanatory variables. We can check whether or not this assumption is met by plotting a scale-location plot (also called spread-location plot). It is a good indication that homoscedasticity is present if the plot shows that the residuals are equally (randomly) spread along a horizontal line and there’s no distinct patterns. We can obtain this plot by using the plot() function and specifying the which = argument to number 3:

# Plotting a scale-location plot 
plot(m2, which = 3)

What do you think? Does the plot show equal variance?

To do this and the next plot in ggplot() requires a bit more work, so we leave that out for now.

Checking for influential outliers

Finally, we could check for influential outliers, i.e outliers that have an impact on the regression line, by plotting a residuals vs leverage plot. Not all outliers may be influential, meaning that the results would stay the same whether we included or excluded them from our analysis. However, if the results change depending on whether we keep or drop an observation, this observation would be defined as an influential observation.

We can interpret the plot by looking out for outlying values at the lower or upper right corner. These are the spots where observations may be influential to the regression line. We should be looking for observations that fall outside of a red, dashed line. When observations are located outside the dashed lines, they are usually influential to the regression results. This essentially means that the results would be different if we dropped those observations from the analysis.

By using the plot() function, and specifying the which = argument to number 5, we can plot a residuals vs leverage plot:

# Plotting a residuals vs leverage plot
plot(m2, which = 5)

What do you think? Does it appear to be any influential outliers? Since we cannot see any of the red, dashed lines, this means that there are no influential outliers. Even though we see that observations number 1216, 83 and 17 are marked as outliers, they are not on the outside of any red, dashed lines and are therefore not influential to the regression results.

What now?

The keen observer may have noticed that some observations are marked in each plot, observations number 1216 and 17. When this happens it may be useful to take a further look at these cases individually. Are there anything remarkably different about these observations? May it be caused by data entry errors?

If you do notice patterns in residuals in your data, the answer to what to do next is not straight-forward. You should go back to your theory, e.g. is it reasonable to assume a linear relationship between the predictors and your outcome? Maybe there are important variables that are missing from your model?

In closing this section on regression diagnostics, we would like to emphasize that the methods described in this part, i.e regression diagnostic plots of residuals, do not give the full picture of whether the assumptions of linear regression are met. As stated in BPS: “(…) regression diagnostics is a subject that could fill an entire book” (chapter 29, p. 50). Obviously, we are not able to cover the full topic here. However, the methods described in this section is a good place to start and sufficient for the purposes of this course.

Including a interaction in regression model

Let’s say we wanted to include an interaction of yeareduc*gender in our model to check whether the effect of years of education on working hours depends on whether the respondents are males or females. In Norway, women tend to outdo men in higher education, and maybe the differences in working hours may be explained by the fact that women tend to have more years of education than men? All we have to do is simply to add an interaction term to our formula = argument.

To estimate a model of the form \(y = \alpha + \beta_1 yeareduc + \beta_2gender + beta_3(yeareduc * gender)\) we add each term to the formula =.

m3 <- lm(formula = hours_work ~ yeareduc + gender + yeareduc*gender,
         data = issp)
summary(m3)

## 
## Call:
## lm(formula = hours_work ~ yeareduc + gender + yeareduc * gender, 
##     data = issp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -38.690  -1.690   1.084   4.803  52.893 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            35.9767     2.2483  16.001   <2e-16 ***
## yeareduc                0.2261     0.1506   1.501   0.1337    
## genderFemale           -6.6673     3.2684  -2.040   0.0417 *  
## yeareduc:genderFemale   0.2209     0.2125   1.039   0.2990    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.49 on 711 degrees of freedom
## Multiple R-squared:  0.03228,    Adjusted R-squared:  0.02819 
## F-statistic: 7.905 on 3 and 711 DF,  p-value: 3.423e-05

Notice that we simply add yeareduc*gender to the formula used in the previous model. By combining two variables using the * symbol, we create an interaction term.

Again, we can obtain the result of our new model using the summary() function:

summary(m3)

## 
## Call:
## lm(formula = hours_work ~ yeareduc + gender + yeareduc * gender, 
##     data = issp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -38.690  -1.690   1.084   4.803  52.893 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            35.9767     2.2483  16.001   <2e-16 ***
## yeareduc                0.2261     0.1506   1.501   0.1337    
## genderFemale           -6.6673     3.2684  -2.040   0.0417 *  
## yeareduc:genderFemale   0.2209     0.2125   1.039   0.2990    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.49 on 711 degrees of freedom
## Multiple R-squared:  0.03228,    Adjusted R-squared:  0.02819 
## F-statistic: 7.905 on 3 and 711 DF,  p-value: 3.423e-05

As we can see from the output, we still get a slope coefficient for both yeareduc and genderFemale, but in addition we also get a coefficient for the interaction term (yeareduc:genderFemale). The coefficients for yeareduc and genderFemale are interpreted the same way as before:

• yeareduc: for each additional year of education, the amount of working hours increases by 0.2261.
• genderFemale: On average, females work 6.6673 hours less than males

The interpretation of the interaction term is a bit different: Simply put, the coefficient for the interaction term tells us that the effect of years of education on working hours per week is 0.2209 higher for females than for males (reference group). However, the interaction term is not significant, meaning that we cannot say that there is an actual interaction effect in the population - it might be due to some random variation in this particular sample.

Visualizing interactions with ggplot2

It is often useful to visualize an interaction. An interaction in regression results in two regression lines, and in cases where there is an interaction effect, the two regression lines will be different from each other. Maybe one of them is much steeper than the other etc.

We can do this by using the ggplot2 package, combined with the geom geom_smooth(), which you are already familiar with:

# Plot with two regression lines 
ggplot(issp, aes(x=yeareduc, y=hours_work, color=gender)) +
  geom_smooth(method = "lm", se = FALSE) + 
  labs(x = "Years of education",
       y = "Working hours per week",
       color = "Moderator") +
  theme_bw()

## `geom_smooth()` using formula 'y ~ x'

However, we usually want to clearly communicate as much information as possible, and combining a scatterplot with regression lines for both genders may be more useful:

# Scatterplot combined with two regression lines
ggplot(issp, aes(x=yeareduc, y=hours_work, color=gender)) +
  geom_jitter(alpha = 0.5) +
  geom_smooth(method = "lm", se = FALSE) + 
  labs(x = "Years of education",
       y = "Working hours per week",
       color = "Moderator") +
  theme_bw()

## `geom_smooth()` using formula 'y ~ x'

This plot communicates more clearly that gender does not substantially moderate the relationship between years of education and working hours per week than the previous one where we only included the regression lines.

Modifying apperance of tables

To be able to merge together several regression models into a single table, we first have to generate each of the models as gtsummary objects and do some modifications to each of them.

As you can see from the help documentation on tbl_regression() (scroll down to the section: “help documentation”), there are a lot of different functions you can use to modify the appearance of the tables. In this session, we’ll focus on the simple things and for those of you who are interested, check out the help documentation and try it out yourself in R.

Let’s begin with our first model, m1, to show you the default output of tbl_regression():

# Default output from tbl_regression()
m1 %>%
  tbl_regression()

In the code above, we send the model m1 through tbl_regression() by using the pipe operator %>% from the magrittr package (included in tidyverse).

From the output we can see the coefficient for the explanatory variable, yeareduc, and it’s corresponding 95% confidence interval and p-value. However, the intercept is not shown in the default output. We can include the intercept in the table by setting the intercept = argument to TRUE:

# Including the intercept in the table
m1 %>%
  tbl_regression(intercept = TRUE)

There we go! Now the intercept is also shown in the table. We would also like to modify the labels of the variables so that they do not just show the variable names. We can do this by using the labels = argument:

# Modifying label names of variables
m1 %>%
  tbl_regression(intercept = TRUE, label = list(yeareduc = "Years of education"))

If we wanted to highlight the variables, we could either bold the labels of the variables by adding the bold_labels() function, or we could use the italicize_labels() function to italicize the variable labels:

# Bold labels of variables
m1 %>%
  tbl_regression(intercept = TRUE, label = list(yeareduc = "Years of education")) %>%
  bold_labels()

# Italicize labels of variables
m1 %>%
  tbl_regression(intercept = TRUE, label = list(yeareduc = "Years of education")) %>%
  italicize_labels()

Preparing each of the models

Now we have shown you some of the basic aesthetic functions to use with tbl_regression(), and we can move on to preparing each of our three models as gtsummary objects. All we have to do is simply to set up the models using the lm() function and then send it through the tbl_regression() function and store it in an object (below we call them t1, t2 and t3). In the code below we do this for all three models, and we do some modifications to each of them:

# Model 1: simple linear regression
t1 <-
   lm(formula = hours_work ~ yeareduc, data = issp) %>%
  tbl_regression(intercept = TRUE, label = list(yeareduc = "Years of education")) %>%
  bold_labels()

t1

# Model 2: Multiple regression, including gender
t2 <-
    lm(formula = hours_work ~ yeareduc + gender, data = issp) %>%
  tbl_regression(intercept = TRUE, label = list(yeareduc = "Years of education", gender = "Gender")) %>%
  bold_labels() %>%
  italicize_levels()

t2

# Model 3: Interaction term, yeareduc*gender
t3 <-
    lm(formula = hours_work ~ yeareduc*gender, data = issp) %>%
  tbl_regression(intercept = TRUE, label = list(yeareduc = "Years of education", gender = "Gender")) %>%
  bold_labels() %>%
  italicize_levels()

t3

As you can see from the output of t2 and t3, we use the function italicize_levels() to italicize the levels of categorical variables, in our case, the gender variable. You can also see in the output from t3 that when we specify the variable labels, the new labels also automatically appears in the interaction term!

Merging the models together

Now that we have stored each of the models as gtsummary objects, we can merge them together in a single table output using the tbl_merge() function:

merged <- tbl_merge(tbls = list(t1, t2, t3), tab_spanner = c("**Model 1**", "**Model 2**", "**Model 3**"))
 
merged

In the code above, we specify which gtsummary objects to merge by setting the tbls = argument, and then we use the tab_spanner = argument to give labels to each of the models.

As you can see, the results for all the models are now merged together side-by-side into a single table!

Help documentation

To see the help-pages for the functions tbl_regression() and tbl_merge(), just run the following codes and the help-pages will open in your browser.

help(tbl_regression, package = "gtsummary")

help(tbl_merge, package = "gtsummary")

Saving tables with models

After we’ve created the table outputs for our models in tbl_regression, we also have to save them to our working folder. We can do this by using the functions as_flex_table() from the gtsummary package in combination with the function save_as_docx() from the flextable package. We want to save the table which merged all of our models, merged. The following code saves the merged table as a word document to our working folder. so first off we have to store the merged table as an R object:

# Saving the table as a word file to our working folder 
merged %>%
  as_flex_table() %>%
  save_as_docx(path = "output.docx")

In the code above we use the object merged, then we send it through the function as_flex_table() to convert the object into a flextable object, and finally, we then send that output through the function save_as_docx() to save the table to our working folder. We specify the path = argument with what we want to name the file (output) and the correct file extension (.docx) to save it as a word file.

Help documentation

To see the help-pages for the functions as_flex_table() and save_as_docx(), just run the following codes and the help-pages will open in your browser.

help(as_flex_table, package = "gtsummary")

help(save_as_docx, package = "flextable")