Crafting formulas

The basic structure of a formula:

y ~ x

Running a model with a formula is straightforward. Note that we don’t even have to put the formula in parentheses:

lm(life_expectancy ~ gdp, data = life_expec_dat)

## 
## Call:
## lm(formula = life_expectancy ~ gdp, data = life_expec_dat)
## 
## Coefficients:
## (Intercept)          gdp  
##    6.70e+01     3.12e-04

Note: Both sides of this equation go by many names. The left-hand side, often written as y, may be called the dependent variable, response variable, or target variable, depending on context.The right-hand side, often written as X, may be referred to as the independent variables, explanatory variables, or regressors. These terms have subtle differences in meaning depending on the context, but for our purposes today it’s enough to know that they all refer to one of the two sides of a regression equation.

Syntax

A quick recap of the syntax used within formulas. The ~ (tilde) sign generally differentiates between dependent and independent variables. The + sign is used to distinguish between independent variables. Here are a few more common formula inputs that are useful to know:

y ~ x1 + x2 # Standard formula

y ~ . # Including all other columns in data as x variables.

y ~ x1 - x2 # Ignore x2 in the analysis

y ~ x1 * x2 # Add the interaction term between x1 and x2
y ~ x1 + x2 + x1*x2 # The same as the above

y ~ x1:x2  # Only the interaction, no main effects (rare that you would want this)

y ~ x + I(x^2) # Adds the higher order term x^2. 

A More Explicit Way to Create Formulas

Formulas are generally straightforward to work with. To make things clearer and easier to automate we can split the process into two steps:

Step 1: Create a string containing the written formula.

formula_string <- paste("life_expectancy", "~ gdp")
formula_string

## [1] "life_expectancy ~ gdp"

Step 2: Transform the string into R’s formula class. Only then will your modeling function accept it as input.

formula <- as.formula(formula_string)
formula

## life_expectancy ~ gdp

Here we can see that the formula we created above is identical to doing it directly.

reg1 <- lm(formula, data = life_expec_dat) 
reg2 <- lm(life_expectancy ~ gdp, data = life_expec_dat)

reg1$coefficients == reg2$coefficients

## (Intercept)         gdp 
##        TRUE        TRUE

Iteration

Often, the modelling process requires you to run the same specification with multiple configurations of both dependent and independent variables. Model formulas make running many similar models much easier:

Step 1: Define a function that let’s you plug in different variables for x or y. The function should take the x and/or y variable(s) as a string and return a model object.

# function to include different y variables
lm_fun_iter_y <-  function(dep_var) {
  lm(as.formula(paste(dep_var, "~ gdp")), data = life_expec_dat)
}

# function to include different x variables 
lm_fun_iter_x <- function(indep_var) {
  lm(as.formula(paste("life_expectancy ~", paste(indep_var, collapse = "+"))), data = life_expec_dat)
}

Notice: Considering your model will likely include many independent variables, it’s unlikely that you’ll run a simple bivariate regression. Therefore, we use a nested paste function above that combines (with collapse()) all independent variables contained in a character vector with +.

Step 2: Use purrr::map() (refer to lab sessions 2 and 3) to iteratively apply the model to each variable contained in a character vector of variables:

# create vector of variables to iterate over
vars <- life_expec_dat |> 
  dplyr::select(-c("life_expectancy", "country")) |> 
  names()

vars

##  [1] "year"                            "status"                         
##  [3] "adult_mortality"                 "infant_deaths"                  
##  [5] "alcohol"                         "percentage_expenditure"         
##  [7] "hepatitis_b"                     "measles"                        
##  [9] "bmi"                             "under_five_deaths"              
## [11] "polio"                           "total_expenditure"              
## [13] "diphtheria"                      "hiv_aids"                       
## [15] "gdp"                             "population"                     
## [17] "thinness_1_19_years"             "thinness_5_9_years"             
## [19] "income_composition_of_resources" "schooling"                      
## [21] "large"

# run a bivariate model for each column
biv_model_out <-  vars |>
  purrr::map(lm_fun_iter_x)

biv_model_out

## [[1]]
## 
## Call:
## lm(formula = as.formula(paste("life_expectancy ~", paste(indep_var, 
##     collapse = "+"))), data = life_expec_dat)
## 
## Coefficients:
## (Intercept)         year  
##    -635.872        0.351  
## 
## 
## [[2]]
## 
## Call:
## lm(formula = as.formula(paste("life_expectancy ~", paste(indep_var, 
##     collapse = "+"))), data = life_expec_dat)
## 
## Coefficients:
##      (Intercept)  statusDeveloping  
##             79.2             -12.1  
....

This returns output from all of the models that you’ve just run. Given that there are now quite a few models, it might be a good idea to keep the names of the independent variables that you feed into the model. You can do this with purrr::set_names().

# run a bivariate model for each column
biv_model_out_w_names <-  vars |>
  purrr::set_names() |> 
  purrr::map(lm_fun_iter_x)

biv_model_out_w_names$year # we can now easily index independent vars

## 
## Call:
## lm(formula = as.formula(paste("life_expectancy ~", paste(indep_var, 
##     collapse = "+"))), data = life_expec_dat)
## 
## Coefficients:
## (Intercept)         year  
##    -635.872        0.351

Say that you’re interested in seeing what effect gdp has on life_expectancy, but you suspect that other covariates such as alcohol consumption, health care expenditure, average body mass index and the propagation of aids also have an effect.

On top of this, it’s possible that these variables might have different effects on life_expectancy depending on how they interact with each other. We can determine which independent variables are robustly associated with the dependent variable by generating models for every possible combination of independent variables and then looking into the distribution of effects across these models.

Let’s review the code introduced in the lecture:

Step 1: Find out how many possible combinations of independent and dependent variables exist:

combinations <- 
  purrr::map(1:5, function(x) {combn(1:5, x)}) |> # create all possible combinations (assume we have 5 covariates)
  purrr::map(ncol) |> # extract number of combinations
  unlist() |> # pull out of list structure
  sum() # compute sum

combinations

## [1] 31

Step 2: Write a function that can run all possible combinations in one go:

combn_models <- function(depvar, covars, data) {
  
  # initialize empty list
  combn_list <- list()
  
  # generate list of covariate combinations
  for (i in seq_along(covars)) {
    combn_list[[i]] <- combn(covars, i, simplify = FALSE)
  } # what is combn?
  
  # unlist to make object easier to work with
  combn_list <-
    unlist(combn_list, recursive = FALSE) 
  
  # function to generate formulas
  gen_formula <- function(covars, depvar) {
    form <-
      as.formula(paste0(depvar, " ~ ", paste0(covars, collapse = "+")))
    form
  }
  
  # generate formulas
  formulas_list <-
    purrr::map(combn_list, gen_formula, depvar = depvar)
  
  # run models
  models_list <- purrr::map(formulas_list, lm, data = data)
  models_list
}

Step 3: Apply the function to our case:

depvar <- "life_expectancy"
covars <-  c("gdp", "percentage_expenditure", "alcohol", "bmi", "hiv_aids")
multiv_model_out <- combn_models(depvar = depvar, covars = covars, data = life_expec_dat)

# check whether we have the correct number of models
length(multiv_model_out)

## [1] 31

Tidying Model Objects

Luckily, we don’t have to deal with this. We can use the tidyverse’s broom package which is part of a collection of packages used for modeling that adhere to tidyverse principles referred to as tidymodels. Tidymodels offers a lot of cool packages that have their own tutorials, so check them out! 💪.

In today’s session, we will only focus on broom. Here are the three key broom functions that you should know:

1. broom::tidy(): Summarizes information about model components.

broom::tidy(multiv_model_out[[30]], conf.int = TRUE, conf.level = 0.95)

## # A tibble: 5 × 7
##   term                 estimate std.error statistic   p.value conf.low conf.high
##   <chr>                   <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
## 1 (Intercept)           6.15e+1 0.273         226.  0          6.10e+1  62.1    
## 2 percentage_expendit…  9.28e-4 0.0000602      15.4 1.94e- 51  8.10e-4   0.00105
## 3 alcohol               4.61e-1 0.0314         14.7 6.97e- 47  3.99e-1   0.522  
## 4 bmi                   1.66e-1 0.00638        26.0 3.46e-133  1.53e-1   0.178  
## 5 hiv_aids             -8.16e-1 0.0226        -36.1 2.43e-233 -8.60e-1  -0.772

2. broom::glance(): Reports information about the entire model.

broom::glance(multiv_model_out[[30]])

## # A tibble: 1 × 12
##   r.squared adj.r.squared sigma statistic p.value    df logLik    AIC    BIC
##       <dbl>         <dbl> <dbl>     <dbl>   <dbl> <dbl>  <dbl>  <dbl>  <dbl>
## 1     0.608         0.607  5.99     1052.       0     4 -8727. 17465. 17501.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

3. augment(): Adds information about observations to a dataset.

broom::augment(multiv_model_out[[30]], se_fit = TRUE) 

## # A tibble: 2,720 × 13
##    .rownames life_expectancy percentage_expenditure alcohol   bmi hiv_aids
##    <chr>               <dbl>                  <dbl>   <dbl> <dbl>    <dbl>
##  1 1                    65                    71.3     0.01  19.1      0.1
##  2 2                    59.9                  73.5     0.01  18.6      0.1
##  3 3                    59.9                  73.2     0.01  18.1      0.1
##  4 4                    59.5                  78.2     0.01  17.6      0.1
##  5 5                    59.2                   7.10    0.01  17.2      0.1
##  6 6                    58.8                  79.7     0.01  16.7      0.1
##  7 7                    58.6                  56.8     0.01  16.2      0.1
##  8 8                    58.1                  25.9     0.03  15.7      0.1
##  9 9                    57.5                  10.9     0.02  15.2      0.1
## 10 10                   57.3                  17.2     0.03  14.7      0.1
## # ℹ 2,710 more rows
## # ℹ 7 more variables: .fitted <dbl>, .se.fit <dbl>, .resid <dbl>, .hat <dbl>,
## #   .sigma <dbl>, .cooksd <dbl>, .std.resid <dbl>

Doesn’t this look so much tidier?

Broom with Many Models

Where broom really shines is when it comes to dealing with multiple models like we have in our multiv_model_out. It allows us to easily move away from hard-to-deal-with lists while retaining each of the three important model outputs. Since our example replicates several variables in different models, we include the .id argument in purrr::map(). This adds a column to our tibble which will contain a model identifier. In our example, the identifier of each model corresponds with the index number of the model output.

multiv_model_out_broom <- purrr::map_dfr(multiv_model_out[29:31], # the model objects
                                         broom::tidy, # the function to iterate with
                                         .id = "model_type") # id column

multiv_model_out_broom

## # A tibble: 16 × 6
##    model_type term                     estimate  std.error statistic   p.value
##    <chr>      <chr>                       <dbl>      <dbl>     <dbl>     <dbl>
##  1 1          (Intercept)            61.5       0.292        210.    0        
##  2 1          gdp                     0.000160  0.00000927    17.2   1.13e- 62
##  3 1          alcohol                 0.392     0.0339        11.6   4.20e- 30
##  4 1          bmi                     0.168     0.00693       24.2   2.72e-115
##  5 1          hiv_aids               -0.779     0.0227       -34.4   1.69e-209
##  6 2          (Intercept)            61.5       0.273        226.    0        
##  7 2          percentage_expenditure  0.000928  0.0000602     15.4   1.94e- 51
##  8 2          alcohol                 0.461     0.0314        14.7   6.97e- 47
##  9 2          bmi                     0.166     0.00638       26.0   3.46e-133
## 10 2          hiv_aids               -0.816     0.0226       -36.1   2.43e-233
## 11 3          (Intercept)            61.5       0.293        210.    0        
## 12 3          gdp                     0.000166  0.0000219      7.59  4.69e- 14
## 13 3          percentage_expenditure -0.0000464 0.000145      -0.321 7.48e-  1
## 14 3          alcohol                 0.394     0.0344        11.5   1.27e- 29
## 15 3          bmi                     0.167     0.00695       24.1   1.26e-114
## 16 3          hiv_aids               -0.779     0.0227       -34.4   2.18e-209

Now that we have multiple model outputs in a format that we can work with, it is time to start thinking about visualizing our results. You can plot the residuals, the goodness of fit statistics etc. Similarly, you could also use broom::glance() to provide some robustness statistics or descriptive tables.

The next section will focus on how to present your results!

Using Regression Tables

The workhorse function in modelsummary is unsurprisingly modelsummary(). 🤔

Here we only need to input the models list as the function does all of the tidying for us.

modelsummary::modelsummary(multiv_model_out[29:31], 
                           output = "kableExtra",
                           escape = FALSE) # can be set to multiple different output formats (see the documentation for more information)

It looks pretty solid already, but can definitely be improved upon.

For example, we can set acceptable number of digits:

modelsummary::modelsummary(multiv_model_out[[31]], 
                           output = "kableExtra",
                           fmt = "%.2f", # 2-digits and trailing 0
                           escape = FALSE) # to avoid escaping special characters in kableExtra

We can also report only the C.I. by getting rid of the p-values and intercept term:

model_table <- modelsummary::modelsummary(multiv_model_out[29:31], 
                                          output = "kableExtra",
                                          fmt = "%.2f",  # 2-digits and trailing 0  
                                          estimate  = "{estimate}", 
                                          statistic = "conf.int",
                                          coef_omit = "Intercept",
                                          escape = FALSE) 
model_table

Wouldn’t it look more professional if our coefficients began with capital letter? And how about we get rid of some of the bloated goodness of fit statistics?

mod_table <- modelsummary::modelsummary(multiv_model_out[29:31], 
                                        output = "kableExtra",
                                        fmt = "%.2f",  # 2-digits and trailing 0  
                                        estimate  = "{estimate}",
                                        statistic = "conf.int",
                                        coef_omit = "Intercept",
                                        coef_rename=c("gdp"="Gdp", 
                                                      "bmi"="Avg. BMI", 
                                                      "alcohol" = "Alcohol Consum.",
                                                      "hiv_aids"= "HIV cases", 
                                                      "percentage_expenditure" = "Health Expenditure (% of GDP)",
                                                      "life_expectancy" = "Life Expectancy"),
                                        gof_omit = 'DF|Deviance|Log.Lik|AIC|BIC',
                                        title = 'An Improved Regression Table',
                                        escape = FALSE)

mod_table

When using the kableExtra package, you can even post-process your table:

mod_table |> 
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = F, fixed_thead = T) |> 
  row_spec(3, color = 'red') |>
  row_spec(5, background = 'lightblue')

Using Plots 🎨

If you are looking for an alternative way to present your results, you might also consider a nice coefficient plot.

Again, modelsummary provides a very accessible function for this: modelplot(). As with the regression table, you input the Model List and not the tidied output. The tidying happens under the hood.

# plot 
modelsummary::modelplot(multiv_model_out[29:31], coef_omit = "Intercept") 

Now we just need to customize the plot to our liking:

modelsummary::modelplot(multiv_model_out[29:31], coef_omit = "Intercept") +
  labs(x = 'Coefficients', 
       y = 'Term names',
       title = 'Linear regression models of "Life expectancy"',
       caption = "Comparing multiple models of life expectancy at the country level, data originated at the WHO") 

We will cover data visualization with R in much greater detail in the next session, so hopefully the section on coefficient plot will be made clearer then!

More Formula info

Create a formula for a model with a large number of variables:

xnam <- paste0("x", 1:25)
(fmla <- as.formula(paste("y ~ ", paste(xnam, collapse= "+"))))

## y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11 + 
##     x12 + x13 + x14 + x15 + x16 + x17 + x18 + x19 + x20 + x21 + 
##     x22 + x23 + x24 + x25

The update function

The update function lets us incrementally modify/update a formula:

#?update

Working with model outputs in base R

Inspecting models in R is straightforward with the summary() function.

lm(life_expectancy ~ gdp, data = life_expec_dat) |> summary()

## 
## Call:
## lm(formula = life_expectancy ~ gdp, data = life_expec_dat)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -30.94  -4.97   2.01   5.82  21.84 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 6.70e+01   1.94e-01   345.7   <2e-16 ***
## gdp         3.12e-04   1.20e-05    25.9   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.6 on 2483 degrees of freedom
##   (453 observations deleted due to missingness)
## Multiple R-squared:  0.213,  Adjusted R-squared:  0.213 
## F-statistic:  672 on 1 and 2483 DF,  p-value: <2e-16

But often you want to post-process estimation results. So let’s examine the output of a model function more closely.

model_out <- lm(life_expectancy ~ gdp, data = life_expec_dat)
class(model_out)

## [1] "lm"

str(model_out)

## List of 13
##  $ coefficients : Named num [1:2] 6.70e+01 3.12e-04
##   ..- attr(*, "names")= chr [1:2] "(Intercept)" "gdp"
##  $ residuals    : Named num [1:2485] -2.22 -7.33 -7.33 -7.74 -7.85 ...
##   ..- attr(*, "names")= chr [1:2485] "1" "2" "3" "4" ...
##  $ effects      : Named num [1:2485] -3458.11 221.84 -7.22 -7.63 -7.74 ...
##   ..- attr(*, "names")= chr [1:2485] "(Intercept)" "gdp" "" "" ...
##  $ rank         : int 2
##  $ fitted.values: Named num [1:2485] 67.2 67.2 67.2 67.2 67.1 ...
##   ..- attr(*, "names")= chr [1:2485] "1" "2" "3" "4" ...
##  $ assign       : int [1:2] 0 1
##  $ qr           :List of 5
##   ..$ qr   : num [1:2485, 1:2] -49.8498 0.0201 0.0201 0.0201 0.0201 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2485] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:2] "(Intercept)" "gdp"
##   .. ..- attr(*, "assign")= int [1:2] 0 1
##   ..$ qraux: num [1:2] 1.02 1.01
##   ..$ pivot: int [1:2] 1 2
##   ..$ tol  : num 1e-07
##   ..$ rank : int 2
##   ..- attr(*, "class")= chr "qr"
##  $ df.residual  : int 2483
##  $ na.action    : 'omit' Named int [1:453] 161 162 163 164 165 166 167 168 169 170 ...
##   ..- attr(*, "names")= chr [1:453] "161" "162" "163" "164" ...
##  $ xlevels      : Named list()
##  $ call         : language lm(formula = life_expectancy ~ gdp, data = life_expec_dat)
##  $ terms        :Classes 'terms', 'formula'  language life_expectancy ~ gdp
##   .. ..- attr(*, "variables")= language list(life_expectancy, gdp)
##   .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:2] "life_expectancy" "gdp"
##   .. .. .. ..$ : chr "gdp"
##   .. ..- attr(*, "term.labels")= chr "gdp"
##   .. ..- attr(*, "order")= int 1
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "predvars")= language list(life_expectancy, gdp)
##   .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
##   .. .. ..- attr(*, "names")= chr [1:2] "life_expectancy" "gdp"
##  $ model        :'data.frame':   2485 obs. of  2 variables:
##   ..$ life_expectancy: num [1:2485] 65 59.9 59.9 59.5 59.2 58.8 58.6 58.1 57.5 57.3 ...
##   ..$ gdp            : num [1:2485] 584.3 612.7 631.7 670 63.5 ...
##   ..- attr(*, "terms")=Classes 'terms', 'formula'  language life_expectancy ~ gdp
##   .. .. ..- attr(*, "variables")= language list(life_expectancy, gdp)
##   .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. ..$ : chr [1:2] "life_expectancy" "gdp"
##   .. .. .. .. ..$ : chr "gdp"
##   .. .. ..- attr(*, "term.labels")= chr "gdp"
##   .. .. ..- attr(*, "order")= int 1
##   .. .. ..- attr(*, "intercept")= int 1
##   .. .. ..- attr(*, "response")= int 1
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. .. ..- attr(*, "predvars")= language list(life_expectancy, gdp)
##   .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
##   .. .. .. ..- attr(*, "names")= chr [1:2] "life_expectancy" "gdp"
##   ..- attr(*, "na.action")= 'omit' Named int [1:453] 161 162 163 164 165 166 167 168 169 170 ...
##   .. ..- attr(*, "names")= chr [1:453] "161" "162" "163" "164" ...
##  - attr(*, "class")= chr "lm"

Ugh, that’s a lot of information in a list structure. Let’s use some helper functions to unpack these:

coef(model_out)
fitted.values(model_out)
residuals(model_out)
model.matrix(model_out)

We can also dig in the summary of the model.

#str(summary(model_out))
summary(model_out)$coefficients

##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  6.7e+01    1.9e-01     346  0.0e+00
## gdp          3.1e-04    1.2e-05      26 2.7e-131