---
title: "Plotting and exporting"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Plotting and exporting}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5
)
has_ggplot2 <- requireNamespace("ggplot2", quietly = TRUE)
has_pheatmap <- requireNamespace("pheatmap", quietly = TRUE)
```

tidymatrix does not have its own plotting system. Instead it makes it easy to
get the data into the shape that plotting packages expect:

* the **metadata** (with PC scores, cluster labels, test results, ...) is a
  data frame, ready for ggplot2;
* `to_long()` gives **one row per matrix cell** with all metadata attached,
  for plots of the individual values;
* `plot_pheatmap()` draws a **heatmap** of the matrix with the metadata as
  annotations.

The same tools get the data out of R, which is covered at the end.

```{r load-packages}
library(tidymatrix)
library(dplyr, warn.conflicts = FALSE)
```

```{r load-ggplot2, eval = has_ggplot2}
library(ggplot2)
theme_set(theme_minimal())
```

## Data

We use the `big5` personality survey (see `?big5`): careless respondents
removed, reverse-keyed items re-scored, and a score per trait added to the
respondent metadata. The [Matrix operations](matrix-operations.html) and
[Row- and column-wise statistics](row-wise-statistics.html) vignettes explain
these steps.

```{r data}
tm <- tidymatrix(big5_responses, big5_respondents, big5_items) |>
  activate(rows) |>
  filter(completion_min > 3.5) |>
  transform_matrix(\(x, flip) ifelse(flip, 6L - x, x), flip = reversed) |>
  compute_across(
    \(answers, items) tapply(answers, items$trait, mean),
    add_to_data = TRUE,
    prefix = "score"
  ) |>
  mutate(age_group = cut(
    age,
    breaks = c(17, 29, 49, 64, Inf),
    labels = c("18-29", "30-49", "50-64", "65+")
  )) |>
  compute_prcomp(n_components = 2, store = FALSE) |>
  activate(columns) |>
  compute_prcomp(n_components = 2, store = FALSE)
```

The two `compute_prcomp()` calls add principal component scores to the
respondent and item metadata. `store = FALSE` skips storing the full `prcomp`
objects, which we do not need here.

## Plotting metadata

### PCA of respondents

After `compute_prcomp()`, the PC scores are columns of the row metadata. Plot
them with `as_tibble()` and ggplot2, and colour by any other variable:

```{r pca-rows, eval = has_ggplot2}
tm |>
  activate(rows) |>
  as_tibble() |>
  ggplot(aes(row_pca_PC1, row_pca_PC2, colour = score_Neuroticism)) +
  geom_point() +
  scale_colour_viridis_c() +
  labs(x = "PC1", y = "PC2", colour = "Neuroticism")
```

### PCA of items

With columns active, items become points. Labelling them with their IDs shows
that items of the same trait group together:

```{r pca-cols, eval = has_ggplot2}
tm |>
  activate(columns) |>
  as_tibble() |>
  ggplot(aes(column_pca_PC1, column_pca_PC2, colour = trait)) +
  geom_text(aes(label = item_id)) +
  labs(x = "PC1", y = "PC2", colour = NULL)
```

### Relationships between metadata variables

Trait scores computed from the matrix sit next to the demographics, so their
relationship is a regular scatterplot:

```{r scores-age, eval = has_ggplot2}
tm |>
  activate(rows) |>
  as_tibble() |>
  ggplot(aes(age, score_Conscientiousness)) +
  geom_jitter(height = 0.05, alpha = 0.5) +
  geom_smooth(method = "lm", formula = y ~ x) +
  labs(x = "Age", y = "Conscientiousness score")
```

### Test results

Results of `compute_ttest()` and friends are tibbles with the item metadata
included, so a volcano plot takes a few lines:

```{r volcano, eval = has_ggplot2}
tm |>
  activate(rows) |>
  filter(gender %in% c("Female", "Male")) |>
  activate(columns) |>
  compute_ttest(group_col = "gender", control = "Male", treatment = "Female") |>
  ggplot(aes(log2fc, -log10(p.value), colour = trait)) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", colour = "grey60") +
  geom_vline(xintercept = 0, colour = "grey60") +
  geom_text(aes(label = item_id)) +
  labs(
    x = "log2(mean women / mean men)",
    y = "-log10(p-value)",
    colour = NULL
  )
```

## Plotting individual values with `to_long()`

`to_long()` converts the whole tidymatrix into a long table: one row per cell,
with the row metadata, the column metadata and the `value`.

```{r to-long}
long <- to_long(tm)
dim(long)
long |> select(respondent_id, gender, age_group, item_id, trait, value)
```

If a column name occurs in both the row and the column metadata, all columns
are prefixed with `row.` and `col.` to keep them apart.

From here, any ggplot is possible. The distribution of answers to the
neuroticism items, by gender:

```{r long-bar, eval = has_ggplot2, fig.height = 5}
long |>
  filter(trait == "Neuroticism", gender != "Non-binary") |>
  count(item_id, gender, value) |>
  group_by(item_id, gender) |>
  mutate(share = n / sum(n)) |>
  ggplot(aes(factor(value), share, fill = gender)) +
  geom_col(position = "dodge") +
  facet_wrap(~item_id) +
  labs(x = "Answer (re-scored)", y = "Share of respondents", fill = NULL)
```

Average answers per item and age group:

```{r long-line, eval = has_ggplot2}
long |>
  group_by(trait, item_id, age_group) |>
  summarise(mean = mean(value), .groups = "drop") |>
  ggplot(aes(age_group, mean, group = item_id, colour = trait)) +
  geom_line() +
  facet_wrap(~trait, nrow = 1) +
  labs(x = "Age group", y = "Mean answer") +
  theme(legend.position = "none", axis.text.x = element_text(angle = 45))
```

`to_long()` always converts the full matrix. Filter the tidymatrix first to
keep the long table small.

## Heatmaps

`plot_pheatmap()` passes the matrix to `pheatmap::pheatmap()` and fills in
the annotations and clusterings from the tidymatrix:

* `row_names` and `col_names` choose the metadata columns used as labels;
* `row_annotation` and `col_annotation` choose the metadata columns shown as
  colour bars (by default all of them, `FALSE` for none);
* if a stored `compute_hclust()` result exists for rows or columns, it is used
  to order them. `TRUE` lets pheatmap cluster, `FALSE` turns clustering off;
* all other arguments go to `pheatmap::pheatmap()`.

### Respondents × items

Standardising the items and clipping extreme values gives a readable colour
scale. Clustering both dimensions with `compute_hclust()` makes the trait
structure visible:

```{r heatmap, eval = has_pheatmap, fig.height = 7}
tm |>
  activate(columns) |>
  scale() |>
  activate(matrix) |>
  clip_values(min = -2, max = 2) |>
  activate(columns) |>
  compute_hclust(method = "ward.D2") |>
  activate(rows) |>
  compute_hclust(method = "ward.D2") |>
  plot_pheatmap(
    col_names = "item_id",
    row_annotation = c("age_group", "gender"),
    col_annotation = "trait",
    show_rownames = FALSE,
    treeheight_row = 20
  )
```

### Aggregated heatmap

A heatmap does not need to show the raw data. Summarising the respondents by
age group gives a compact overview, here in questionnaire order and without
clustering:

```{r heatmap-agg, eval = has_pheatmap, fig.height = 3}
tm |>
  activate(rows) |>
  group_by(age_group) |>
  summarise(n = n()) |>
  activate(columns) |>
  arrange(trait, item_id) |>
  plot_pheatmap(
    row_names = "age_group",
    col_names = "item_id",
    row_annotation = FALSE,
    col_annotation = "trait",
    row_cluster = FALSE,
    col_cluster = FALSE,
    display_numbers = TRUE,
    number_format = "%.1f",
    fontsize_number = 6
  )
```

## Exporting

### The matrix

`pull_active()` returns the active component. With the matrix active, that is
a plain matrix:

```{r export-matrix}
m <- tm |>
  activate(matrix) |>
  pull_active()

m[1:3, 1:6]
```

### The metadata

With rows or columns active, `pull_active()`, `as.data.frame()` and
`as_tibble()` return the metadata, including anything added by analyses:

```{r export-metadata}
tm |>
  activate(rows) |>
  as_tibble() |>
  select(respondent_id, starts_with("score_"), starts_with("row_pca"))
```

### Everything at once

`to_long()` is the most complete export: every value together with all its
metadata, in a format that spreadsheets, databases and other tools understand.

```{r export-file}
out <- tempfile(fileext = ".csv")

tm |>
  activate(columns) |>
  select(item_id, trait, item_text) |>
  activate(rows) |>
  select(respondent_id, age, gender, country) |>
  to_long() |>
  write.csv(out, row.names = FALSE)

read.csv(out, nrows = 3)
```

The first two steps keep only the metadata columns that are needed; the matrix
is not affected.

To save the tidymatrix itself, including stored analyses, use `saveRDS()` and
`readRDS()`.

## See also

* [PCA and clustering](statistical-analysis.html)
* [MDS, t-SNE and UMAP](dimensionality-reduction.html)
* [Row- and column-wise statistics](row-wise-statistics.html)