A basic example for a multiple linear regression

Let’s look at another basic example, again using the {palmerpenguins} data:

library(equatiomatic)

# fit a basic multiple linear regression model
m <- lm(bill_length_mm ~ bill_depth_mm + flipper_length_mm, penguins)

Now we can pull the TeX code with extract_eq

extract_eq(m)

\[ \operatorname{bill\_length\_mm} = \alpha + \beta_{1}(\operatorname{bill\_depth\_mm}) + \beta_{2}(\operatorname{flipper\_length\_mm}) + \epsilon \]

You can of course set echo = FALSE as well, and then you’ll get just the equation, which will look like the below.

\[ \operatorname{bill\_length\_mm} = \alpha + \beta_{1}(\operatorname{bill\_depth\_mm}) + \beta_{2}(\operatorname{flipper\_length\_mm}) + \epsilon \]

Let’s dive into a bit more complex example.

An example with a slightly more complicated model

The above was super simple, but—often—you have categorical variables with lots of levels and interactions, which can make the output of your models (and writing a model equation for your models) a bit more complicated. {equatiomatic} is intended to “smooth out” some of these issues, requiring the same process as in the above example to extract (and present) the model’s equation:

For example, here is an example using a categorical variable (island) in an interaction with the bill_depth_mm variable we used in the example above:

m2 <- lm(bill_length_mm ~ bill_depth_mm * island, penguins)
extract_eq(m2)

\[ \operatorname{bill\_length\_mm} = \alpha + \beta_{1}(\operatorname{bill\_depth\_mm}) + \beta_{2}(\operatorname{island}_{\operatorname{Dream}}) + \beta_{3}(\operatorname{island}_{\operatorname{Torgersen}}) + \beta_{4}(\operatorname{bill\_depth\_mm} \times \operatorname{island}_{\operatorname{Dream}}) + \beta_{5}(\operatorname{bill\_depth\_mm} \times \operatorname{island}_{\operatorname{Torgersen}}) + \epsilon \]

Sometimes, for such models, the equations can get overly long.

That’s where the wrap and terms_per_line arguments come in.

extract_eq(m2, wrap = TRUE) # default terms_per_line = 4

\[ \begin{aligned} \operatorname{bill\_length\_mm} &= \alpha + \beta_{1}(\operatorname{bill\_depth\_mm}) + \beta_{2}(\operatorname{island}_{\operatorname{Dream}}) + \beta_{3}(\operatorname{island}_{\operatorname{Torgersen}})\ + \\ &\quad \beta_{4}(\operatorname{bill\_depth\_mm} \times \operatorname{island}_{\operatorname{Dream}}) + \beta_{5}(\operatorname{bill\_depth\_mm} \times \operatorname{island}_{\operatorname{Torgersen}}) + \epsilon \end{aligned} \]

extract_eq(m2, wrap = TRUE, terms_per_line = 2)

\[ \begin{aligned} \operatorname{bill\_length\_mm} &= \alpha + \beta_{1}(\operatorname{bill\_depth\_mm})\ + \\ &\quad \beta_{2}(\operatorname{island}_{\operatorname{Dream}}) + \beta_{3}(\operatorname{island}_{\operatorname{Torgersen}})\ + \\ &\quad \beta_{4}(\operatorname{bill\_depth\_mm} \times \operatorname{island}_{\operatorname{Dream}}) + \beta_{5}(\operatorname{bill\_depth\_mm} \times \operatorname{island}_{\operatorname{Torgersen}})\ + \\ &\quad \epsilon \end{aligned} \]

Maybe you want different intercept notation, such as \(\beta_0\)? Simply pass “beta” to the intercept argument, as follows.

extract_eq(m2, wrap = TRUE, intercept = "beta")

\[ \begin{aligned} \operatorname{bill\_length\_mm} &= \beta_{0} + \beta_{1}(\operatorname{bill\_depth\_mm}) + \beta_{2}(\operatorname{island}_{\operatorname{Dream}}) + \beta_{3}(\operatorname{island}_{\operatorname{Torgersen}})\ + \\ &\quad \beta_{4}(\operatorname{bill\_depth\_mm} \times \operatorname{island}_{\operatorname{Dream}}) + \beta_{5}(\operatorname{bill\_depth\_mm} \times \operatorname{island}_{\operatorname{Torgersen}}) + \epsilon \end{aligned} \]

We also wrap all the variable names in \operatorname by default so they show up as plain text, but if you’d like your variable names to be italicized just set ital_vars = TRUE.

extract_eq(m2, wrap = TRUE, ital_vars = TRUE)

\[ \begin{aligned} bill\_length\_mm &= \alpha + \beta_{1}(bill\_depth\_mm) + \beta_{2}(island_{Dream}) + \beta_{3}(island_{Torgersen})\ + \\ &\quad \beta_{4}(bill\_depth\_mm \times island_{Dream}) + \beta_{5}(bill\_depth\_mm \times island_{Torgersen}) + \epsilon \end{aligned} \]

Raw TeX code

Currently, the intercept argument defaults to "alpha" and only takes one additional argument, "beta". However, the raw_tex and greek arguments allows you to specify whatever syntax you would like both for the intercept and for the coefficients. For example

extract_eq(m2,
  wrap = TRUE,
  intercept = "\\hat{\\phi}",
  greek = "\\hat{\\gamma}",
  raw_tex = TRUE
)

\[ \begin{aligned} \operatorname{bill\_length\_mm} &= \hat{\phi} + \hat{\gamma}_{1}(\operatorname{bill\_depth\_mm}) + \hat{\gamma}_{2}(\operatorname{island}_{\operatorname{Dream}}) + \hat{\gamma}_{3}(\operatorname{island}_{\operatorname{Torgersen}})\ + \\ &\quad \hat{\gamma}_{4}(\operatorname{bill\_depth\_mm} \times \operatorname{island}_{\operatorname{Dream}}) + \hat{\gamma}_{5}(\operatorname{bill\_depth\_mm} \times \operatorname{island}_{\operatorname{Torgersen}}) + \epsilon \end{aligned} \]

Coefficients

In many cases, you’re interested more in including the coefficient estimates (e.g., 3.04), than the Greek symbols (e.g., \(\\beta\_1\)). This may be helpful for communicating the results of a model (and, possibly, for teaching about the statistical model). To do this, simply change the use_coefs argument:

extract_eq(m2, wrap = TRUE, use_coefs = TRUE)

\[ \begin{aligned} \operatorname{\widehat{bill\_length\_mm}} &= 63.72 - 1.16(\operatorname{bill\_depth\_mm}) - 54.12(\operatorname{island}_{\operatorname{Dream}}) - 35.13(\operatorname{island}_{\operatorname{Torgersen}})\ + \\ &\quad 3.05(\operatorname{bill\_depth\_mm} \times \operatorname{island}_{\operatorname{Dream}}) + 1.73(\operatorname{bill\_depth\_mm} \times \operatorname{island}_{\operatorname{Torgersen}}) \end{aligned} \]

Logistic Regression

lr <- glm(sex ~ species * bill_length_mm,
  data = penguins,
  family = binomial(link = "logit")
)

extract_eq(lr, wrap = TRUE)

\[ \begin{aligned} \log\left[ \frac { P( \operatorname{sex} = \operatorname{male} ) }{ 1 - P( \operatorname{sex} = \operatorname{male} ) } \right] &= \alpha + \beta_{1}(\operatorname{species}_{\operatorname{Chinstrap}}) + \beta_{2}(\operatorname{species}_{\operatorname{Gentoo}}) + \beta_{3}(\operatorname{bill\_length\_mm})\ + \\ &\quad \beta_{4}(\operatorname{species}_{\operatorname{Chinstrap}} \times \operatorname{bill\_length\_mm}) + \beta_{5}(\operatorname{species}_{\operatorname{Gentoo}} \times \operatorname{bill\_length\_mm}) \end{aligned} \]

You can also (optionally) show the how the data are assumed to be distributed.

extract_eq(lr, wrap = TRUE, show_distribution = TRUE)

\[ \begin{aligned} \operatorname{sex} &\sim Bernoulli\left(\operatorname{prob}_{\operatorname{sex} = \operatorname{male}}= \hat{P}\right) \\ \log\left[ \frac { \hat{P} }{ 1 - \hat{P} } \right] &= \alpha + \beta_{1}(\operatorname{species}_{\operatorname{Chinstrap}}) + \beta_{2}(\operatorname{species}_{\operatorname{Gentoo}}) + \beta_{3}(\operatorname{bill\_length\_mm})\ + \\ &\quad \beta_{4}(\operatorname{species}_{\operatorname{Chinstrap}} \times \operatorname{bill\_length\_mm}) + \beta_{5}(\operatorname{species}_{\operatorname{Gentoo}} \times \operatorname{bill\_length\_mm}) \end{aligned} \]

Probit Regression

Probit regression works similarly to logistic regression:

pr <- glm(sex ~ species * bill_length_mm,
  data = penguins,
  family = binomial(link = "probit")
)

extract_eq(pr, wrap = TRUE)

\[ \begin{aligned} P( \operatorname{sex} = \operatorname{male} ) &= \Phi[\alpha + \beta_{1}(\operatorname{species}_{\operatorname{Chinstrap}}) + \beta_{2}(\operatorname{species}_{\operatorname{Gentoo}}) + \beta_{3}(\operatorname{bill\_length\_mm})\ + \\ &\qquad\ \beta_{4}(\operatorname{species}_{\operatorname{Chinstrap}} \times \operatorname{bill\_length\_mm}) + \beta_{5}(\operatorname{species}_{\operatorname{Gentoo}} \times \operatorname{bill\_length\_mm})] \end{aligned} \]

Again, you can (optionally) show the how the data are assumed distributed.

extract_eq(pr, wrap = TRUE, show_distribution = TRUE)

\[ \begin{aligned} \operatorname{sex} &\sim Bernoulli\left(\operatorname{prob}_{\operatorname{sex} = \operatorname{male}}= \hat{P}\right) \\ \hat{P} &= \Phi[\alpha + \beta_{1}(\operatorname{species}_{\operatorname{Chinstrap}}) + \beta_{2}(\operatorname{species}_{\operatorname{Gentoo}}) + \beta_{3}(\operatorname{bill\_length\_mm})\ + \\ &\qquad\ \beta_{4}(\operatorname{species}_{\operatorname{Chinstrap}} \times \operatorname{bill\_length\_mm}) + \beta_{5}(\operatorname{species}_{\operatorname{Gentoo}} \times \operatorname{bill\_length\_mm})] \end{aligned} \]

Ordinal Regression

For these examples, we’ll use wine tasting data from the {ordinal} package. If you don’t already have this package, you can download it with

install.packages("ordinal")

The data look like this

We’ll fit a model from the documentation, with the ordinal rating response predicted by an interaction between the temperature and and contact.

You can fit ordinal response models with either MASS::polr or ordinal::clm, and {equatiomatic} works with either, and with logistic or probit link functions.

mass_ologit <- MASS::polr(rating ~ temp * contact,
  data = ordinal::wine
)

extract_eq(mass_ologit, wrap = TRUE, terms_per_line = 2)
#> 
#> Re-fitting to get Hessian
#> 
#> 
#> Re-fitting to get Hessian

mass_oprobit <- MASS::polr(rating ~ temp * contact,
  data = ordinal::wine,
  method = "probit"
)

extract_eq(mass_oprobit, wrap = TRUE, terms_per_line = 2)
#> 
#> Re-fitting to get Hessian
#> 
#> 
#> Re-fitting to get Hessian

ordinal_ologit <- ordinal::clm(rating ~ temp * contact,
  data = ordinal::wine,
  link = "logit"
)

extract_eq(ordinal_ologit, wrap = TRUE, terms_per_line = 2)

ordinal_probit <- ordinal::clm(rating ~ temp * contact,
  data = ordinal::wine,
  link = "probit"
)

extract_eq(ordinal_probit, wrap = TRUE, terms_per_line = 2)