Save multiple tab_model html objects into 1 file (e.g., PDF, Word, or html file)

I'd like to 3 sjPlot::tab_model object into 1 file (e.g., PDF, Word, or html file). I'm not sure how to do this. Please help.

I used the top 300 rows of the diamonds dataset.

diamonds_top300 <- data.frame(dplyr::top_n(ggplot2::diamonds, 300, y))

I repeated the same glmer on 3 DVs, price, x, and y.

## DV price for diamonds_top300, using poisson family
(
  mlm_top300_poisson_price <- 
    lme4::glmer(
      price ~ cut + color + carat + (1 | clarity) + (1 | depth),
      data = diamonds_top300,
      family = poisson()
    )
)

Here is the code that creates the 3 tab_model files I used.

## creates sjplot::tab_model objects of models of interest
# ---- NOTE: tab_model of mlm_top300_poisson_price
tab_model_mlm_top300_poisson_price <- 
  sjPlot::tab_model(mlm_top300_poisson_price)
# ---- NOTE: tab_model of mlm_top300_poisson_x
tab_model_mlm_top300_poisson_x <- 
  sjPlot::tab_model(mlm_top300_poisson_x)
# ---- NOTE: tab_model of mlm_top300_poisson_y
tab_model_mlm_top300_poisson_y <- 
  sjPlot::tab_model(mlm_top300_poisson_y)

Please let me know how I can do this.

FYI, I use a 2013 Macbook Pro with a 2.4 GHz dual-core intel chip, 8 GB of ram, macOS big sur 11.2.2, RStudio Version 1.4.1106, and the R Base Package 4.04.

Thanks.

Here is the script I used to create this analysis:

#### follow up tests for continuous variables associated with general linear mixed models analyses ####



# Loads packages
# ---- NOTE: making plots and diamonds dataset
if(!require(ggplot2)){install.packages("ggplot2")}
# ---- NOTE: run mixed effects models
if(!require(lme4)){install.packages("lme4")}
# ---- NOTE: for data wrangling
if(!require(dplyr)){install.packages("dplyr")}
# ---- NOTE: for emmeans command
if(!require(emmeans)){install.packages("emmeans")}
# ---- NOTE: for sjPlot command
if(!require(sjPlot)){install.packages("sjPlot")}




# dataset creation

## for dataset with top 300 rows
# ---- NOTE: selects only the top 300 rows of the dataset
diamonds_top300 <- data.frame(dplyr::top_n(ggplot2::diamonds, 300, y))
# ---- NOTE: gives dataset info
head(diamonds_top300)
str(diamonds_top300)
colnames(diamonds_top300)
nrow(diamonds_top300)
# ---- NOTE: gives unique values of Fixed and Random effects, and dvs
unique(diamonds_top300$price)
unique(diamonds_top300$x)
unique(diamonds_top300$y)
unique(diamonds_top300$cut)
unique(diamonds_top300$color)
unique(diamonds_top300$carat)
unique(diamonds_top300$clarity)
unique(diamonds_top300$x)
unique(diamonds_top300$table)




# model creation

## DV price for diamonds_top300, using poisson family
(
  mlm_top300_poisson_price <- 
    lme4::glmer(
      price ~ cut + color + carat + (1 | clarity) + (1 | depth),
      data = diamonds_top300,
      family = poisson()
    )
)

### gives summary of model
summary(mlm_top300_poisson_price)

### gives Anova of model
Anova(mlm_top300_poisson_price)

### does follow up tests, based on significant fixed effects in Anova(mlm_top300_poisson)

#### when fixed effect is cut
emmeans::emmeans(
  mlm_top300_poisson_price, 
  "cut")

#### when fixed effect is color
emmeans::emmeans(
  mlm_top300_poisson_price, 
  "color")

#### when fixed effect is carat
emmeans::emmeans(
  mlm_top300_poisson_price, 
  "carat")

## DV x for diamonds_top300, using poisson family
(
  mlm_top300_poisson_x <- 
    lme4::glmer(
      x ~ cut + color + carat + (1 | clarity) + (1 | depth),
      data = diamonds_top300,
      family = poisson()
    )
)

### gives summary of model
summary(mlm_top300_poisson_x)

### gives Anova of model
Anova(mlm_top300_poisson_x)

### does follow up tests, based on significant fixed effects in Anova(mlm_top300_poisson)

#### when fixed effect is cut
emmeans::emmeans(
  mlm_top300_poisson_x, 
  "cut")

#### when fixed effect is color
emmeans::emmeans(
  mlm_top300_poisson_x, 
  "color")

#### when fixed effect is carat
emmeans::emmeans(
  mlm_top300_poisson_x, 
  "carat")

## DV y for diamonds_top300, using poisson family
(
  mlm_top300_poisson_y <- 
    lme4::glmer(
      y ~ cut + color + carat + (1 | clarity) + (1 | depth),
      data = diamonds_top300,
      family = poisson()
    )
)

### gives summary of model
summary(mlm_top300_poisson_y)

### gives Anova of model
Anova(mlm_top300_poisson_y)

### does follow up tests, based on significant fixed effects in Anova(mlm_top300_poisson)

#### when fixed effect is cut
emmeans::emmeans(
  mlm_top300_poisson_y, 
  "cut")

#### when fixed effect is color
emmeans::emmeans(
  mlm_top300_poisson_y, 
  "color")

#### when fixed effect is carat
emmeans::emmeans(
  mlm_top300_poisson_y, 
  "carat")

## creates sjplot::tab_model objects of models of interest
# ---- NOTE: tab_model of mlm_top300_poisson_price
tab_model_mlm_top300_poisson_price <- 
  sjPlot::tab_model(mlm_top300_poisson_price)
# ---- NOTE: tab_model of mlm_top300_poisson_x
tab_model_mlm_top300_poisson_x <- 
  sjPlot::tab_model(mlm_top300_poisson_x)
# ---- NOTE: tab_model of mlm_top300_poisson_y
tab_model_mlm_top300_poisson_y <- 
  sjPlot::tab_model(mlm_top300_poisson_y)