## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
options(rmarkdown.html_vignette.check_title = FALSE)

## ----setup, message = FALSE, error=FALSE, warning=FALSE-----------------------
library(mda.biber)
library(corrplot)
library(tidyverse)
library(kableExtra)

## ----echo=FALSE---------------------------------------------------------------
knitr::kable(head(micusp_biber[,1:6]))

## ----echo=FALSE---------------------------------------------------------------
knitr::kable(head(micusp_biber[,63:68]))

## ----format_data, message = FALSE, error=FALSE, warning=FALSE-----------------
d <- micusp_biber |>
  mutate(doc_id = str_extract(doc_id, "^[A-Z]+")) |>
  mutate(doc_id = as.factor(doc_id)) |>
  select(where(~ any(. != 0))) # removes any columns containing all zeros as required for corrplot

## ----scree_plot, fig.height=4.5, fig.width=7----------------------------------
screeplot_mda(d)

## ----corr_plot, fig.height=5.5, fig.width=7-----------------------------------
# create a correlation matrix
cor_m <- cor(d[,-1], method = "pearson")

# plot the matrix
corrplot(cor_m, type = "upper", order = "hclust", tl.col = "black", tl.srt = 45, diag = F, tl.cex = 0.5)

## ----mda----------------------------------------------------------------------
m <- mda_loadings(d, n_factors = 2)

## ----results='hide'-----------------------------------------------------------
attributes(m)$group_means

## ----echo=FALSE---------------------------------------------------------------
knitr::kable(head(attributes(m)$group_means))

## ----results='hide'-----------------------------------------------------------
attributes(m)$loadings %>% arrange(-Factor1)

## ----echo=FALSE---------------------------------------------------------------
knitr::kable(head(attributes(m)$loadings %>% arrange(-Factor1)))

## ----echo=FALSE---------------------------------------------------------------
knitr::kable(tail(attributes(m)$loadings %>% arrange(-Factor1)))

## ----stickplot_mda, fig.height=4.5, fig.width=2.5-----------------------------
stickplot_mda(m, n_factor = 1)

## ----heatmap_mda, fig.height=5.5, fig.width=7---------------------------------
heatmap_mda(m, n_factor = 1)

## ----boxplot_mda, fig.height=6.5, fig.width=7---------------------------------
boxplot_mda(m, n_factor = 1)

## -----------------------------------------------------------------------------

# Carry out regression
f1_lm <- lm(Factor1 ~ group, data = m)
f2_lm <- lm(Factor2 ~ group, data = m)

# Convert ANOVA results into data.frames allows for easier name manipulation
f1_aov <- data.frame(anova(f1_lm), r.squared = c(summary(f1_lm)$r.squared*100, NA))
f2_aov <- data.frame(anova(f2_lm), r.squared = c(summary(f2_lm)$r.squared*100, NA))

# Putting all into one data.frame/table
anova_results <- data.frame(rbind(c("DF", "Sum Sq", "Mean Sq", "F value", "Pr(>F)", "*R*^2^",
                                    "DF", "Sum Sq", "Mean Sq", "F value", "Pr(>F)", "*R*^2^"),
                                   cbind(round(f1_aov, 2), round(f2_aov, 2))))
colnames(anova_results) <- c("", "", "", "", "", "", "", "", "", "", "", "")
row.names(anova_results)[1] <- ""
anova_results[is.na(anova_results)] <- "--"

## ----results='asis'-----------------------------------------------------------
anova_results %>% knitr::kable("html") %>%
  kableExtra::kable_styling(bootstrap_options = "striped", full_width = F) %>%
  kableExtra::add_header_above(c("", "Dimension 1" = 6, "Dimension 2" = 6))