crp.qmd

---
title: "Model selection"
subtitle: "C-reactive protein"
author:
  - name: "Nathan Constantine-Cooke" 
    url: https://scholar.google.com/citations?user=2emHWR0AAAAJ&hl=en&oi=ao
    corresponding: true
    affiliations:
      - ref: CGEM
      - ref: HGU
  - name: "Karla Monterrubio-Gómez"
    url: https://scholar.google.com/citations?user=YmyxSXAAAAAJ&hl=en
    affiliations:
      - ref: HGU
  - name: "Catalina A. Vallejos"
    url: https://scholar.google.com/citations?user=lkdrwm0AAAAJ&hl=en&oi=ao
    affiliations:
      - ref: HGU
#comments:
#  giscus: 
#    repo: quarto-dev/quarto-docs
---

## Introduction

```{R Setup}
#| message: false
set.seed(123)
if (file.exists("/.dockerenv")) { # Check if running in Docker
  # Assume igmm/cvallejo-predicct/libdr/ is passed to the data volume
  prefix <- "data/"
} else {
  # Assume running outside of a Docker container and the IGC(/IGMM) datastore is
  # mounted at /Volumes
  prefix <- "/Volumes/igmm/cvallejo-predicct/libdr/"
}


##########################
#--     Packages       --#
##########################

library(tidyverse)
# Support package (source found in libdr/)
library(libdr)
## Modelling ##
library(lcmm)
library(kml) # K-means
## Presentation ##
library(patchwork)
library(ggdist)
library(ggalluvial)
library(pander)
library(qqplotr)


##########################
#--     Data read      --#
#########################

dict <- readRDS(paste0(prefix, "processed/dict.RDS"))
crp_median <- readRDS(paste0(prefix, "processed/median-crp.RDS"))
crp <- readRDS(paste0(prefix, "processed/crp.RDS"))

dk.crp <- read.csv(paste0(prefix, "Denmark/2024-11-29/crp_8_models.csv"),
                    sep = ";")
```

## Model selection

```{R Load Processed CRP Models}
#| fig-width: 12
#| fig-height: 6.75
# set the number of groups
G.crp <- numeric()
models.crp.ma <- list()
G.cands <- seq(2, 10)
for (G.cand in G.cands) {
  file.name <- paste0(prefix, "/cache/crp-ma/crp-", G.cand, ".RDS")
  if (file.exists(file.name)) {
    G.crp <- c(G.crp, G.cand)
    models.crp.ma[[G.cand]] <- readRDS(file.name)
  }
}
rm(G.cand)
```

```{R Processed CRP alluvial}
#| label: fig-crp-alluvial-ma
#| fig-cap: Alluvial plot of cluster membership across G for CRP

alluvial.df <- matrix(nrow = 0, ncol = 3)
colnames(alluvial.df) <- c("ids", "class", "G")

for (G in G.crp) {
  alluvial.df <- rbind(alluvial.df, cbind(models.crp.ma[[G]]$pprob[, 1:2],
    G = G
  ))
}
alluvial.df <- as.data.frame(alluvial.df)

alluvial.df$ids <- as.character(alluvial.df$ids)
alluvial.df$class <- as.factor(alluvial.df$class)

alluvial.df[alluvial.df[, "G"] == 3, "class"] <-
  plyr::mapvalues(alluvial.df[alluvial.df[, "G"] == 3, "class"],
    from = c(seq(1, 3)),
    to = c(1, 2, 3)
  )

alluvial.df[alluvial.df[, "G"] == 4, "class"] <-
  plyr::mapvalues(alluvial.df[alluvial.df[, "G"] == 4, "class"],
    from = c(seq(1, 4)),
    to = c(3, 1, 2, 4)
  )

alluvial.df[alluvial.df[, "G"] == 5, "class"] <-
  plyr::mapvalues(alluvial.df[alluvial.df[, "G"] == 5, "class"],
    from = c(seq(1, 5)),
    to = c(3, 5, 2, 1, 4)
  )

alluvial.df[alluvial.df[, "G"] == 6, "class"] <-
  plyr::mapvalues(alluvial.df[alluvial.df[, "G"] == 6, "class"],
    from = c(seq(1, 6)),
    to = c(2, 1, 3, 6, 4, 5)
  )

alluvial.df[alluvial.df[, "G"] == 7, "class"] <-
  plyr::mapvalues(alluvial.df[alluvial.df[, "G"] == 7, "class"],
    from = c(seq(1, 7)),
    to = c(2, 1, 7, 3, 5, 4, 6)
  )

alluvial.df[alluvial.df[, "G"] == 8, "class"] <-
  plyr::mapvalues(alluvial.df[alluvial.df[, "G"] == 8, "class"],
    from = c(seq(1, 8)),
    to = c(2, 6, 1, 7, 5, 4, 3, 8)
  )


alluvial.df[alluvial.df[, "G"] == 9, "class"] <-
  plyr::mapvalues(alluvial.df[alluvial.df[, "G"] == 9, "class"],
    from = c(seq(1, 9)),
    to = c(8, 2, 4, 1, 5, 9, 3, 6, 7)
  )

alluvial.df[alluvial.df[, "G"] == 10, "class"] <-
  plyr::mapvalues(alluvial.df[alluvial.df[, "G"] == 10, "class"],
    from = c(seq(1, 10)),
    to = c(9, 7, 1, 5, 3, 4, 2, 10, 6, 8)
  )

# eliminate label switching

p <- ggplot(
  alluvial.df,
  aes(
    x = G,
    stratum = class,
    alluvium = ids,
    fill = class,
    label = class
  )
) +
  scale_x_discrete(expand = c(.1, .1)) +
  geom_flow() +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 3) +
  theme_minimal() +
  theme(legend.position = "none") +
  scale_fill_manual(values = rainbow(10)) +
  xlab("Assumed number of clusters") +
  ylab("Frequency")
print(p)


for (format in c("png", "pdf")) {
  ggsave(paste0("plots/ma-nocor/alluvial.", format),
    p,
    width = 12,
    height = 6.75,
    units = "in"
  )
}
```

```{R Processed CRP alluvial no switching}
#| label: fig-crp-alluvial-ma-noswitch
#| fig-cap: Alluvial plot of cluster membership across G for CRP without label switching

alluvial.df <- matrix(nrow = 0, ncol = 3)
colnames(alluvial.df) <- c("ids", "class", "G")

for (G in G.crp) {
  alluvial.df <- rbind(alluvial.df, cbind(models.crp.ma[[G]]$pprob[, 1:2],
    G = G
  ))
}
alluvial.df <- as.data.frame(alluvial.df)

alluvial.df$ids <- as.character(alluvial.df$ids)
alluvial.df$class <- as.factor(alluvial.df$class)

p <- ggplot(
  alluvial.df,
  aes(
    x = G,
    stratum = class,
    alluvium = ids,
    fill = class,
    label = class
  )
) +
  scale_x_discrete(expand = c(.1, .1)) +
  geom_flow() +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 3) +
  theme_minimal() +
  theme(legend.position = "none") +
  scale_fill_manual(values = rainbow(10)) +
  xlab("Assumed number of clusters") +
  ylab("Frequency")
print(p)
```

## Posterior classifications

From the below data, we can see how posterior probabilities change as the number
of assumed clusters increase


:::{.column-page}
::: {.panel-tabset}

```{R CRP-ma postpro}
#| results: "asis"
postprob_quiet <- quietly(postprob)
for (G in G.crp) {
  cat(paste0("##### G = ", G, "\n"))
  output <- postprob_quiet(models.crp.ma[[G]])$result
  cat(pander::pander(output[[1]]))
  cat(pander::pander(output[[2]]))
  cat(pander::pander(output[[3]]))
}
rm(output)
```

:::
:::

```{R}
pprobs.2 <- c()
groups.2 <- models.crp.ma[[2]]
pprobs.3 <- c()
groups.3 <- models.crp.ma[[3]]
pprobs.4 <- c()
groups.4 <- models.crp.ma[[4]]
pprobs.5 <- c()
groups.5 <- models.crp.ma[[5]]
pprobs.6 <- c()
groups.6 <- models.crp.ma[[6]]
pprobs.7 <- c()
groups.7 <- models.crp.ma[[7]]
pprobs.8 <- c()
groups.8 <- models.crp.ma[[8]]
pprobs.9 <- c()
groups.9 <- models.crp.ma[[9]]
pprobs.10 <- c()
groups.10 <- models.crp.ma[[10]]


for (i in 1:nrow(models.crp.ma[[2]]$pprob)) {
  class.2 <- groups.2$pprob[i, 2]
  pprobs.2 <- c(pprobs.2, groups.2$pprob[i, class.2 + 2])
  class.3 <- groups.3$pprob[i, 2]
  pprobs.3 <- c(pprobs.3, groups.3$pprob[i, class.3 + 2])
  class.4 <- groups.4$pprob[i, 2]
  pprobs.4 <- c(pprobs.4, groups.4$pprob[i, class.4 + 2])
  class.5 <- groups.5$pprob[i, 2]
  pprobs.5 <- c(pprobs.5, groups.5$pprob[i, class.5 + 2])
  class.6 <- groups.6$pprob[i, 2]
  pprobs.6 <- c(pprobs.6, groups.6$pprob[i, class.6 + 2])
  class.7 <- groups.7$pprob[i, 2]
  pprobs.7 <- c(pprobs.7, groups.7$pprob[i, class.7 + 2])
  class.8 <- groups.8$pprob[i, 2]
  pprobs.8 <- c(pprobs.8, groups.8$pprob[i, class.8 + 2])
  class.9 <- groups.9$pprob[i, 2]
  pprobs.9 <- c(pprobs.9, groups.9$pprob[i, class.9 + 2])
  class.10 <- groups.10$pprob[i, 2]
  pprobs.10 <- c(pprobs.10, groups.10$pprob[i, class.10 + 2])
}
pprobs.2 <- tibble(prob = pprobs.2)
pprobs.3 <- tibble(prob = pprobs.3)
pprobs.4 <- tibble(prob = pprobs.4)
pprobs.5 <- tibble(prob = pprobs.5)
pprobs.6 <- tibble(prob = pprobs.6)
pprobs.7 <- tibble(prob = pprobs.7)
pprobs.8 <- tibble(prob = pprobs.8)
pprobs.9 <- tibble(prob = pprobs.9)
pprobs.10 <- tibble(prob = pprobs.10)

pprobs.2$Model <- as.factor(rep("Two clusters", nrow(pprobs.2)))
pprobs.3$Model <- as.factor(rep("Three clusters", nrow(pprobs.3)))
pprobs.4$Model <- as.factor(rep("Four clusters", nrow(pprobs.4)))
pprobs.5$Model <- as.factor(rep("Five clusters", nrow(pprobs.5)))
pprobs.6$Model <- as.factor(rep("Six clusters", nrow(pprobs.6)))
pprobs.7$Model <- as.factor(rep("Seven clusters", nrow(pprobs.7)))
pprobs.8$Model <- as.factor(rep("Eight clusters", nrow(pprobs.8)))
pprobs.9$Model <- as.factor(rep("Nine clusters", nrow(pprobs.9)))
pprobs.10$Model <- as.factor(rep("Ten clusters", nrow(pprobs.10)))
pprobs <- rbind(
  pprobs.2,
  pprobs.3,
  pprobs.4,
  pprobs.5,
  pprobs.6,
  pprobs.7,
  pprobs.8,
  pprobs.9,
  pprobs.10
)

p <- pprobs %>%
  ggplot(aes(x = prob, y = Model)) +
  # geom_histogram(bins = 40, fill = NA, position="identity")
  stat_slab(
    fill = "#5D5179", color = "gray",
    size = 0.8,
    alpha = 0.7
  ) +
  geom_dots(color = "#4F759B", dotsize = 1) +
  xlab("Posterior probability for cluster membership") +
  ylab("") +
  ggtitle(
    "Distribution of Posterior Probabilities Across Models",
    "Subject-specific posterior probabilities for assigned cluster"
  ) +
  theme_minimal() +
  scale_y_discrete(limits = rev)
ggsave("plots/pprob-distribution-crp-ma.png",
  p,
  width = 8.5,
  height = 8.5,
  units = "in"
)
ggsave("plots/pprob-distribution-crp-ma.pdf",
  p,
  width = 8.5,
  height = 8.5,
  units = "in"
)
print(p)
```

#### Residual plots

::: {.panel-tabset}

```{R CRP MA residuals}
#| results: "asis"
for (G in G.crp) {
  cat(paste0("##### G = ", G, "\n"))
  plot(models.crp.ma[[G]], shades = TRUE)
}
```

:::

#### Model metrics

```{R CRP MA metrics plot}
#| label: fig-crp-ma-metrics
#| fig-cap: "Model metrics for CRP MA for G = 4-10"
summaryplot(models.crp.ma[[4]],
  models.crp.ma[[5]],
  models.crp.ma[[6]],
  models.crp.ma[[7]],
  models.crp.ma[[8]],
  models.crp.ma[[9]],
  models.crp.ma[[10]],
  which = c("loglik", "AIC", "BIC", "entropy", "ICL"),
  mfrow = c(2, 3),
  axis = "Class"
)
```

```{R Processed CRP metrics table}
crp.metrics <- makeMetrics(G.crp, models.crp.ma)
buildDT(crp.metrics)
```

#### Spaghetti plots per cluster



```{R}
for (G in G.crp) {
  # Data frame to hold processed data
  new.crp <- data.frame(
    ids = numeric(),
    crp_result = numeric(),
    crp_time = numeric(),
    class = numeric()
  )

  for (clust in 1:G) {
    ids.clust <- subset(models.crp.ma[[G]]$pprob, class == clust)$ids
    n.clust <- length(ids.clust)
    rand <- sample(n.clust, n.clust) # Randomise the order of the ids
    iters <- floor(n.clust / 6) # How many groups of six are there?

    # Matrix to hold the smoothed data
    crp.ma <- matrix(NA, nrow = iters, ncol = 7)
    for (i in 0:(iters - 1)) {
      # Find ids for group of five
      ids.select <- ids.clust[rand[((i * 6) + 1):((i * 6) + 6)]]
      crp.subset <- subset(crp_median, ids %in% ids.select)
      # Median process as per CRP preprocessing
      for (j in seq(0, 6)) {
        if (j == 6) {
          sub.obs <- subset(
            crp.subset,
            crp_time >= j - 0.5 & crp_time <= j + 1
          )
        } else {
          sub.obs <- subset(
            crp.subset,
            crp_time >= j - 0.5 & crp_time < j + 0.5
          )
        }
        if (nrow(sub.obs) > 0) {
          crp.ma[i + 1, j + 1] <- median(sub.obs$crp_result)
        }
      }
    }

    rownames(crp.ma) <- 1:iters
    crp.ma <- reshape2::melt(t(crp.ma),
      id.vars = row.names(crp.ma),
      na.rm = TRUE
    )
    colnames(crp.ma) <- c("crp_time", "ids", "crp_result")
    crp.ma <- crp.ma[, c(2, 3, 1)] # Make ids first column
    crp.ma$crp_time <- crp.ma$crp_time - 1
    # Take into account uneven spacing at start and end
    crp.ma$crp_time <- plyr::mapvalues(crp.ma$crp_time,
      from = c(0, 6),
      to = c(0.25, 6.25)
    )
    crp.ma$class <- clust # Identify cluster assignment
    new.crp <- rbind(new.crp, crp.ma)
  }

  cairo_pdf(paste0("plots/spaghetti/crp-", G, ".pdf"),
    width = 10,
    height = 17.5
  )
  grid::grid.newpage()
  spaghettiPlot(new.crp,
    models.crp.ma,
    G,
    clusters = TRUE,
    tmax = 6.25,
    sizes = TRUE,
    var.time = "crp_time",
    ylim = "data"
  )
  invisible(dev.off())

  png(paste0("plots/spaghetti/crp-", G, ".png"),
    width = 10,
    height = 17.5,
    units = "in",
    res = 300
  )
  grid::grid.newpage()
  spaghettiPlot(new.crp,
    models.crp.ma,
    G,
    clusters = TRUE,
    tmax = 6.25,
    sizes = TRUE,
    var.time = "crp_time",
    ylim = "data"
  )
  invisible(dev.off())

  grid::grid.newpage()
  spaghettiPlot(new.crp,
    models.crp.ma,
    G,
    clusters = TRUE,
    tmax = 6.25,
    sizes = TRUE,
    var.time = "crp_time",
    ylim = "data"
  )
}
```



### Chosen model

```{R}
new.crp <- data.frame(
  ids = numeric(),
  crp_result = numeric(),
  crp_time = numeric(),
  class = numeric()
)

for (clust in 1:8) {
  ids.clust <- subset(models.crp.ma[[8]]$pprob, class == clust)$ids
  n.clust <- length(ids.clust)
  rand <- sample(n.clust, n.clust) # Randomise the order of the ids
  iters <- floor(n.clust / 6) # How many groups of six are there?

  # Matrix to hold the smoothed data
  crp.ma <- matrix(NA, nrow = iters, ncol = 7)
  for (i in 0:(iters - 1)) {
    # Find ids for group of five
    ids.select <- ids.clust[rand[((i * 6) + 1):((i * 6) + 6)]]
    crp.subset <- subset(crp_median, ids %in% ids.select)
    # Median process as per CRP preprocessing
    for (j in seq(0, 6)) {
      if (j == 6) {
        sub.obs <- subset(
          crp.subset,
          crp_time >= j - 0.5 & crp_time <= j + 1
        )
      } else {
        sub.obs <- subset(
          crp.subset,
          crp_time >= j - 0.5 & crp_time < j + 0.5
        )
      }
      if (nrow(sub.obs) > 0) {
        crp.ma[i + 1, j + 1] <- median(sub.obs$crp_result)
      }
    }
  } 

  rownames(crp.ma) <- 1:iters
  crp.ma <- reshape2::melt(t(crp.ma),
    id.vars = row.names(crp.ma),
    na.rm = TRUE
  )
  colnames(crp.ma) <- c("crp_time", "ids", "crp_result")
  crp.ma <- crp.ma[, c(2, 3, 1)] # Make ids first column
  crp.ma$crp_time <- crp.ma$crp_time - 1
  # Take into account uneven spacing at start and end
  crp.ma$crp_time <- plyr::mapvalues(crp.ma$crp_time,
    from = c(0, 6),
    to = c(0.25, 6.25)
  )
  crp.ma$class <- clust # Identify cluster assignment
  new.crp <- rbind(new.crp, crp.ma)
}
```

#### Labelling

The labels/order of the clusters generated above are random. To improve the
readability and assist with interpretation of associations, it would be
beneficial to order clusters by some analogue for disease severity. This section
explores possible solutions for this task. 

##### Cumulative inflammation across follow-up

One possible solution is to order clusters by cumulative inflammation which is
given by the area under the curve for each mean cluster trajectory.

```{R}
rank.full <- rankCumulative(models.crp.ma[[8]],
  tmax = 6.25,
  var.time = "crp_time"
)
rank.full %>%
  ggplot(aes(x = paste0("CRP", New), y = Area)) +
  geom_bar(stat = "identity", fill = "#745C97", color = "#39375B") +
  ylab("Cumulative inflammation (Area)") +
  xlab("Cluster") +
  theme_minimal()
```

```{R}
png(paste0("plots/spaghetti/crp-ncs-reordered.png"),
  width = 10,
  height = 16,
  units = "in",
  res = 300
)
grid::grid.newpage()
spaghettiPlot(new.crp,
  models.crp.ma,
  G = 8,
  log = TRUE,
  tmax = 6.25,
  sizes = TRUE,
  knots = FALSE,
  var.time = "crp_time",
  clusters = TRUE,
  ylim = "data",
  mapping = rank.full$Original
)
invisible(dev.off())

cairo_pdf(paste0("plots/spaghetti/crp-ncs-reordered.pdf"),
  width = 10,
  height = 16
)
grid::grid.newpage()
spaghettiPlot(new.crp,
  models.crp.ma,
  G = 8,
  log = TRUE,
  tmax = 6.25,
  sizes = TRUE,
  knots = FALSE,
  var.time = "crp_time",
  clusters = TRUE,
  ylim = "data",
  mapping = rank.full$Original
)
invisible(dev.off())

spaghettiPlot(new.crp,
  models.crp.ma,
  G = 8,
  log = TRUE,
  tmax = 6.25,
  sizes = TRUE,
  knots = FALSE,
  var.time = "crp_time",
  clusters = TRUE,
  ylim = "data",
  mapping = rank.full$Original
)
```

With Danish predictions

```{R}
dk.crp <- dk.crp %>%
  mutate(cluster = str_replace_all(cluster, "Ypred_class", "")) %>%
  filter(time <= 4.25)
```

```{R}
dk.crp$cluster <- plyr::mapvalues(dk.crp$cluster,
                                   from = seq(1, 8),
                                   to = c(1, 2, 8, 3, 5, 7, 6, 4))
```

```{R}
cairo_pdf(paste0("plots/spaghetti/crp-scot-dk.pdf"),
  width = 10,
  height = 16
)
grid::grid.newpage()
spaghettiPlot(new.crp,
  models.crp.ma,
  G = 8,
  log = TRUE,
  tmax = 6.25,
  sizes = TRUE,
  knots = FALSE,
  var.time = "crp_time",
  clusters = TRUE,
  ylim = "data",
  external = dk.crp,
  mapping = rank.full$Original
)
invisible(dev.off())

png(paste0("plots/spaghetti/crp-scot-dk.png"),
  width = 10,
  height = 16,
  units = "in",
  res = 300
)
grid::grid.newpage()
spaghettiPlot(new.crp,
  models.crp.ma,
  G = 8,
  log = TRUE,
  tmax = 6.25,
  sizes = TRUE,
  knots = FALSE,
  var.time = "crp_time",
  clusters = TRUE,
  ylim = "data",
  external = dk.crp,
  mapping = rank.full$Original
)
invisible(dev.off())

grid::grid.newpage()
spaghettiPlot(new.crp,
  models.crp.ma,
  G = 8,
  log = TRUE,
  tmax = 6.25,
  sizes = TRUE,
  knots = FALSE,
  var.time = "crp_time",
  clusters = TRUE,
  ylim = "data",
  external = dk.crp,
  mapping = rank.full$Original
)
```


```{R}
for (G in c(7, 9)) {
  new.crp <- data.frame(
    ids = numeric(),
    crp_result = numeric(),
    crp_time = numeric(),
    class = numeric()
  )

  for (clust in 1:G) {
    ids.clust <- subset(models.crp.ma[[G]]$pprob, class == clust)$ids
    n.clust <- length(ids.clust)
    rand <- sample(n.clust, n.clust) # Randomise the order of the ids
    iters <- floor(n.clust / 6) # How many groups of six are there?

    # Matrix to hold the smoothed data
    crp.ma <- matrix(NA, nrow = iters, ncol = 7)
    for (i in 0:(iters - 1)) {
      # Find ids for group of five
      ids.select <- ids.clust[rand[((i * 6) + 1):((i * 6) + 6)]]
      crp.subset <- subset(crp_median, ids %in% ids.select)
      # Median process as per CRP preprocessing
      for (j in seq(0, 6)) {
        if (j == 6) {
          sub.obs <- subset(
            crp.subset,
            crp_time >= j - 0.5 & crp_time <= j + 1
          )
        } else {
          sub.obs <- subset(
            crp.subset,
            crp_time >= j - 0.5 & crp_time < j + 0.5
          )
        }
        if (nrow(sub.obs) > 0) {
          crp.ma[i + 1, j + 1] <- median(sub.obs$crp_result)
        }
      }
    }

    rownames(crp.ma) <- 1:iters
    crp.ma <- reshape2::melt(t(crp.ma),
      id.vars = row.names(crp.ma),
      na.rm = TRUE
    )
    colnames(crp.ma) <- c("crp_time", "ids", "crp_result")
    crp.ma <- crp.ma[, c(2, 3, 1)] # Make ids first column
    crp.ma$crp_time <- crp.ma$crp_time - 1
    # Take into account uneven spacing at start and end
    crp.ma$crp_time <- plyr::mapvalues(crp.ma$crp_time,
      from = c(0, 6),
      to = c(0.25, 6.25)
    )
    crp.ma$class <- clust # Identify cluster assignment
    new.crp <- rbind(new.crp, crp.ma)
  }


  rank.full <- rankCumulative(models.crp.ma[[G]],
    tmax = 6.25,
    var.time = "crp_time"
  )
  png(paste0("plots/spaghetti/crp-ncs-reordered-", G, ".png"),
    width = 10,
    height = 16,
    units = "in",
    res = 300
  )

  grid::grid.newpage()
  spaghettiPlot(
    new.crp,
    models.crp.ma,
    G = G,
    log = TRUE,
    tmax = 6.25,
    sizes = TRUE,
    knots = FALSE,
    var.time = "crp_time",
    clusters = TRUE,
    ylim = "data",
    mapping = rank.full$Original
  )
  invisible(dev.off())

  cairo_pdf(paste0("plots/spaghetti/crp-ncs-reordered-", G, ".pdf"),
    width = 10,
    height = 16
  )
  grid::grid.newpage()
  spaghettiPlot(
    new.crp,
    models.crp.ma,
    G = G,
    log = TRUE,
    tmax = 6.25,
    sizes = TRUE,
    knots = FALSE,
    var.time = "crp_time",
    clusters = TRUE,
    ylim = "data",
    mapping = rank.full$Original
  )
  invisible(dev.off())
}
```

### Alluvial anchored plot

```{R}
alluvial.df <- matrix(nrow = 0, ncol = 3)
colnames(alluvial.df) <- c("ids", "class", "G")

for (G in G.crp) {
  alluvial.df <- rbind(alluvial.df, cbind(models.crp.ma[[G]]$pprob[, 1:2],
    G = G
  ))
}
alluvial.df <- as.data.frame(alluvial.df)

alluvial.df$ids <- as.character(alluvial.df$ids)
alluvial.df$class <- as.factor(alluvial.df$class)


alluvial.df[alluvial.df[, "G"] == 2, "class"] <-
  plyr::mapvalues(alluvial.df[alluvial.df[, "G"] == 2, "class"],
    from = c(1, 2),
    to = c(1, 2)
  )

alluvial.df[alluvial.df[, "G"] == 3, "class"] <-
  plyr::mapvalues(alluvial.df[alluvial.df[, "G"] == 3, "class"],
    from = c(seq(1, 3)),
    to = c(1, 2, 6)
  )

alluvial.df[alluvial.df[, "G"] == 4, "class"] <-
  plyr::mapvalues(alluvial.df[alluvial.df[, "G"] == 4, "class"],
    from = c(seq(1, 4)),
    to = c(6, 1, 2, 7)
  )

alluvial.df[alluvial.df[, "G"] == 5, "class"] <-
  plyr::mapvalues(alluvial.df[alluvial.df[, "G"] == 5, "class"],
    from = c(seq(1, 5)),
    to = c(6, 5, 2, 1, 7)
  )

alluvial.df[alluvial.df[, "G"] == 6, "class"] <-
  plyr::mapvalues(alluvial.df[alluvial.df[, "G"] == 6, "class"],
    from = c(seq(1, 6)),
    to = c(2, 1, 6, 3, 7, 5)
  )

alluvial.df[alluvial.df[, "G"] == 7, "class"] <-
  plyr::mapvalues(alluvial.df[alluvial.df[, "G"] == 7, "class"],
    from = c(seq(1, 7)),
    to = c(2, 1, 4, 6, 5, 7, 3)
  )

alluvial.df[alluvial.df[, "G"] == 8, "class"] <-
  plyr::mapvalues(alluvial.df[alluvial.df[, "G"] == 8, "class"],
    from = rank.full$Original,
    to = rank.full$New
  )


alluvial.df[alluvial.df[, "G"] == 9, "class"] <-
  plyr::mapvalues(alluvial.df[alluvial.df[, "G"] == 9, "class"],
    from = c(seq(1, 9)),
    to = c(8, 2, 7, 1, 5, 9, 6, 3, 4)
  )


alluvial.df[alluvial.df[, "G"] == 10, "class"] <-
  plyr::mapvalues(alluvial.df[alluvial.df[, "G"] == 10, "class"],
    from = c(seq(1, 10)),
    to = c(9, 4, 1, 5, 6, 7, 2, 10, 3, 8)
  )


# eliminate label switching

p <- ggplot(
  alluvial.df,
  aes(
    x = G,
    stratum = class,
    alluvium = ids,
    fill = class,
    label = class
  )
) +
  scale_x_discrete(expand = c(.1, .1)) +
  geom_flow() +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 3) +
  theme_minimal() +
  theme(legend.position = "none") +
  scale_fill_manual(values = c(viridis::viridis(8), "#FFF282", "#F8F3C3")) +
  xlab("Assumed number of clusters") +
  ylab("Frequency")
print(p)


for (format in c("png", "pdf")) {
  ggsave(paste0("plots/ma-nocor/anchored-alluvial.", format),
    p,
    width = 12,
    height = 6.75,
    units = "in"
  )
}

saveRDS(p, paste0(prefix, "processed/plots/crp-alluvial.RDS"))
```

## Acknowledgments {.appendix}

This work is funded by the Medical Research Council & University of Edinburgh
via a Precision Medicine PhD studentship (MR/N013166/1, to **NC-C**).

This work has made use of the resources provided by the [Edinburgh Compute and
Data Facility (ECDF)](http://www.ecdf.ed.ac.uk/).

## Author contributions {.appendix}



## References {.appendix}

::: {#refs}
:::

## Reuse {.appendix}

Licensed by 
<a href="https://creativecommons.org/licenses/by/4.0/">CC BY</a>
 unless otherwise stated.
 
## {.appendix}

<div class = "center">
  <img class = "center" src="../../images/MRC_HGU_Edinburgh RGB.png" alt="MRC Human Genetics Unit logo" height = 50px>
  <img src="../../images/cgem-logo.png" alt="Centre for Genomic & Experimental Medicine logo" height = 50px>
</div>

---

## Session info {.appendix}

```{R Session info}
pander(sessionInfo())
```