Converting a model into latex equation

equation_lm.R

equation_lm <- function(model) {
  rhs  <- colnames(model.matrix(model))[-1]
  lhs  <- all.vars(formula(model))[1]
  lhs_eq <- paste(lhs, "= ")

  betas <- paste0("\\beta_{", seq_along(rhs), "}(")
	rhs_eq <- paste0(betas, rhs, ")")
 	rhs_eq <- paste("\\alpha +", paste(rhs_eq, collapse = " + "))
 	
 	error <- "+ \\epsilon" 

 	cat(
 	    paste("$$\n", 
 	      paste0(lhs_eq, rhs_eq), 
 	      error, 
 	      "\n$$")
 	    )
}

# fit a couple simple models
mod1 <- lm(mpg ~ cyl + disp, mtcars)
mod2 <- lm(mpg ~ ., mtcars)

# test it out
equation_lm(mod1)

## $$
##  mpg = \alpha + \beta_{1}(cyl) + \beta_{2}(disp) + \epsilon 
## $$

equation_lm(mod2)

## $$
##  mpg = \alpha + \beta_{1}(cyl) + \beta_{2}(disp) + \beta_{3}(hp) + \beta_{4}(drat) + \beta_{5}(wt) + \beta_{6}(qsec) + \beta_{7}(vs) + \beta_{8}(am) + \beta_{9}(gear) + \beta_{10}(carb) + \epsilon 
## $$

# Test with a categorical predictor
mod3 <- lm(Sepal.Length ~ Sepal.Width + Species, iris)
equation_lm(mod3)

## $$
##  Sepal.Length = \alpha + \beta_{1}(Sepal.Width) + \beta_{2}(Speciesversicolor) + \beta_{3}(Speciesvirginica) + \epsilon 
## $$