Formulas and model families

suppressPackageStartupMessages(library("tidyverse"))

package 㤼㸱tidyverse㤼㸲 was built under R version 3.6.3

suppressPackageStartupMessages(library("modelr"))

package 㤼㸱modelr㤼㸲 was built under R version 3.6.3

1. What happens if you repeat the analysis of sim2 using a model without an intercept. What happens to the model equation? What happens to the predictions?

To run a model without an intercept, add - 1 or + 0 to the right-hand-side of the formula:

mod2a <- lm(y ~ x - 1, data = sim2)
mod2 <- lm(y ~ x, data = sim2)

The predictions are exactly the same in the models with and without an intercept:

grid <- sim2 %>%
  data_grid(x) %>%
  spread_predictions(mod2, mod2a)
grid

2. Use model_matrix() to explore the equations generated for the models I fit to sim3 and sim4. Why is * a good shorthand for interaction?

For x1 * x2 when x2 is a categorical variable produces indicator variables x2b, x2c, x2d and variables x1:x2b, x1:x2c, and x1:x2d which are the products of x1 and x2* variables:

x3 <- model_matrix(y ~ x1 * x2, data = sim3)
x3

We can confirm that the variables x1:x2b is the product of x1 and x2b,

all(x3[["x1:x2b"]] == (x3[["x1"]] * x3[["x2b"]]))

[1] TRUE

and similarly for x1:x2c and x2c, and x1:x2d and x2d:

all(x3[["x1:x2c"]] == (x3[["x1"]] * x3[["x2c"]]))

[1] TRUE

all(x3[["x1:x2d"]] == (x3[["x1"]] * x3[["x2d"]]))

[1] TRUE

For x1 * x2 where both x1 and x2 are continuous variables, model_matrix() creates variables x1, x2, and x1:x2:

x4 <- model_matrix(y ~ x1 * x2, data = sim4)
x4

Confirm that x1:x2 is the product of the x1 and x2,

all(x4[["x1"]] * x4[["x2"]] == x4[["x1:x2"]])

[1] TRUE

The asterisk * is good shorthand for an interaction since an interaction between x1 and x2 includes terms for x1, x2, and the product of x1 and x2.

3. Using the basic principles, convert the formulas in the following two models into functions. (Hint: start by converting the categorical variable into 0-1 variables.)

mod1 <- lm(y ~ x1 + x2, data = sim3)
mod2 <- lm(y ~ x1 * x2, data = sim3)

The problem is to convert the formulas in the models into functions. I will assume that the function is only handling the conversion of the right hand side of the formula into a model matrix. The functions will take one argument, a data frame with x1 and x2 columns, and it will return a data frame. In other words, the functions will be special cases of the model_matrix() function.

Consider the right hand side of the first formula, ~ x1 + x2. In the sim3 data frame, the column x1 is an integer, and the variable x2 is a factor with four levels.

levels(sim3$x2)

[1] "a" "b" "c" "d"

Since x1 is numeric it is unchanged. Since x2 is a factor it is replaced with columns of indicator variables for all but one of its levels. I will first consider the special case in which x2 only takes the levels of x2 in sim3. In this case, “a” is considered the reference level and omitted, and new columns are made for “b”, “c”, and “d”.

model_matrix_mod1 <- function(.data) {
  mutate(.data,
    x2b = as.numeric(x2 == "b"),
    x2c = as.numeric(x2 == "c"),
    x2d = as.numeric(x2 == "d"),
    `(Intercept)` = 1
  ) %>%
    select(`(Intercept)`, x1, x2b, x2c, x2d)
}
model_matrix_mod1(sim3)

A more general function for ~ x1 + x2 would not hard-code the specific levels in x2.

model_matrix_mod1b <- function(.data) {
  # the levels of x2
  lvls <- levels(.data$x2)
  # drop the first level
  # this assumes that there are at least two levels
  lvls <- lvls[2:length(lvls)]
  # create an indicator variable for each level of x2
  for (lvl in lvls) {
    # new column name x2 + level name
    varname <- str_c("x2", lvl)
    # add indicator variable for lvl
    .data[[varname]] <- as.numeric(.data$x2 == lvl)
  }
  # generate the list of variables to keep
  x2_variables <- str_c("x2", lvls)
  # Add an intercept
  .data[["(Intercept)"]] <- 1
  # keep x1 and x2 indicator variables
  select(.data, `(Intercept)`, x1, one_of(x2_variables))
}
model_matrix_mod1b(sim3)

Consider the right hand side of the first formula, ~ x1 * x2. The output data frame will consist of x1, columns with indicator variables for each level (except the reference level) of x2, and columns with the x2 indicator variables multiplied by x1.

As with the previous formula, first I’ll write a function that hard-codes the levels of x2.

model_matrix_mod2 <- function(.data) {
  mutate(.data,
    `(Intercept)` = 1,
    x2b = as.numeric(x2 == "b"),
    x2c = as.numeric(x2 == "c"),
    x2d = as.numeric(x2 == "d"),
    `x1:x2b` = x1 * x2b,
    `x1:x2c` = x1 * x2c,
    `x1:x2d` = x1 * x2d
  ) %>%
    select(`(Intercept)`, x1, x2b, x2c, x2d, `x1:x2b`, `x1:x2c`, `x1:x2d`)
}
model_matrix_mod2(sim3)

For a more general function which will handle arbitrary levels in x2, I will extend the model_matrix_mod1b() function that I wrote earlier.

model_matrix_mod2b <- function(.data) {
  # get dataset with x1 and x2 indicator variables
  out <- model_matrix_mod1b(.data)
  # get names of the x2 indicator columns
  x2cols <- str_subset(colnames(out), "^x2")
  # create interactions between x1 and the x2 indicator columns
  for (varname in x2cols) {
    # name of the interaction variable
    newvar <- str_c("x1:", varname)
    out[[newvar]] <- out$x1 * out[[varname]]
  }
  out
}
model_matrix_mod2b(sim3)

These functions could be further generalized to allow for x1 and x2 to be either numeric or factors. However, generalizing much more than that and we will soon start reimplementing all of the matrix_model() function.

4. For sim4, which of mod1 and mod2 is better? I think mod2 does a slightly better job at removing patterns, but it’s pretty subtle. Can you come up with a plot to support my claim?

Estimate models mod1 and mod2 on sim4,

mod1 <- lm(y ~ x1 + x2, data = sim4)
mod2 <- lm(y ~ x1 * x2, data = sim4)

and add the residuals from these models to the sim4 data,

sim4_mods <- gather_residuals(sim4, mod1, mod2)

Frequency plots of both the residuals,

ggplot(sim4_mods, aes(x = resid, colour = model)) +
  geom_freqpoly(binwidth = 0.5) +
  geom_rug()

and the absolute values of the residuals,

ggplot(sim4_mods, aes(x = abs(resid), colour = model)) +
  geom_freqpoly(binwidth = 0.5) +
  geom_rug()

does not show much difference in the residuals between the models. However, mod2 appears to have fewer residuals in the tails of the distribution between 2.5 and 5 (although the most extreme residuals are from mod2.

This is confirmed by checking the standard deviation of the residuals of these models,

sim4_mods %>%
  group_by(model) %>%
  summarise(resid = sd(resid))

The standard deviation of the residuals of mod2 is smaller than that of mod1.