Chapter 11
Unsupervised models

Lectures 21–23 · latent factors behind questionnaires, mixtures behind measurements, the posteriors that will not agree with themselves, and copulas that separate dependence from the margins

Every chapter so far has had an outcome: something the model was asked to predict or explain. This chapter drops it. We are left with a table of measurements and no label, and the question is whether the table itself contains structure: a few underlying dimensions along which people or objects vary, or a few groups into which they fall. Both kinds of structure are latent, never observed directly, and both are inferred by the same Bayesian machinery as before: write down a model in which the latent structure generates the data, and compute the posterior.

• Lecture 21 is about latent factors. A personality questionnaire asks 25 questions; psychologists claim the answers reflect five underlying traits. We build principal components, probabilistic PCA and the factor model, fit a Bayesian factor model in Stan, meet the rotation problem, and ask how many factors the data actually support.
• Lecture 22 is about mixtures. We hide the species labels of the Palmer penguins from Chapter 1 and ask whether four body measurements can find the species again. Along the way we meet label switching, a posterior that will not agree with itself, soft assignments, and the Dirichlet process, and then we look for breast-cancer subtypes in expression data.
• Lecture 23 is about copulas. Two genes can be strongly correlated in very different ways. We separate a joint distribution into its margins and its dependence with Sklar's theorem, meet the main copula families and their tails, fit them in Stan to co-amplified genes, use the Gaussian copula for questionnaires that mix ordinal, binary and continuous answers, and check whether home and away goals in football are independent.

Lecture 21 — Latent factors

11.1 Twenty-five questions

The Synthetic Aperture Personality Assessment project (SAPA) collected answers to personality items from thousands of volunteers online. William Revelle distributes 2,800 of them as the bfi data set (Revelle, 2024): 25 statements, five for each of the "Big Five" personality traits (Agreeableness, Conscientiousness, Extraversion, Neuroticism and Openness), each rated on a six-point scale from "very inaccurate" to "very accurate" as a description of oneself. "I am indifferent to the feelings of others" (A1), "I am exacting in my work" (C1), "I make friends easily" (E4), "I get angry easily" (N1), "I am full of ideas" (O1). Some statements are worded in reverse, so that agreeing indicates less of the trait; we reverse their scores so that every item points the same way.

The questions: do 25 answers really record 25 different things, or a handful of underlying dimensions? If a handful, how many, and what are they? Psychology's answer, the five-factor model, was itself found largely by factor analysis of questionnaires like this one. Here we ask whether a Bayesian factor model applied to these data agrees.

library(tidyverse)
library(rstan)
library(loo)
options(mc.cores = parallel::detectCores())

harbour <- c(ink = "#10191d", muted = "#4a5f68", grid = "#d5e3e6",
             sky = "#1d5187", amber = "#805907", rose = "#ad3b28",
             moss = "#1d6b54", lilac = "#5b3a94")
theme_harbour <- function(base_size = 11) {
  theme_minimal(base_size = base_size) +
    theme(text = element_text(colour = harbour[["ink"]]),
          axis.text = element_text(colour = harbour[["muted"]]),
          panel.grid.major = element_line(colour = harbour[["grid"]]),
          panel.grid.minor = element_blank(),
          strip.text = element_text(colour = harbour[["ink"]], hjust = 0),
          legend.position = "top")
}
theme_set(theme_harbour())

tmp <- tempfile()
download.file("https://raw.githubusercontent.com/cran/psychTools/master/data/bfi.rda", tmp, quiet = TRUE)
load(tmp)                                                   # creates `bfi`
reversed <- c("A1", "C4", "C5", "E1", "E2", "O2", "O5")     # negatively worded items
items <- bfi |> select(A1:O5) |> drop_na() |> mutate(across(all_of(reversed), ~ 7 - .x))
Y <- scale(as.matrix(items))                                # standardised: correlations, not raw scores
c(people = nrow(Y), items = ncol(Y))
round(cor(items)[c("N1", "N2", "E1", "E4", "C2"), c("N1", "N2", "E1", "E4", "C2")], 2)

people  items 
  2436     25 
      N1    N2    E1    E4    C2
N1  1.00  0.72 -0.01 -0.14 -0.02
N2  0.72  1.00 -0.02 -0.15 -0.01
E1 -0.01 -0.02  1.00  0.42 -0.02
E4 -0.14 -0.15  0.42  1.00  0.12
C2 -0.02 -0.01 -0.02  0.12  1.00

After dropping people with any missing answer, 2,436 of the 2,800 remain. The correlations show the pattern factor analysis looks for. Two Neuroticism items correlate 0.72 with each other and two Extraversion items 0.42, while items written for different traits correlate between −0.15 and 0.12. Items written for the same trait hang together; items written for different traits mostly do not.

11.2 Latent variables as causes

The factor model makes a causal claim, and it is worth stating it before fitting anything. A person has a level of neuroticism, which we cannot observe. That level causes her answers to "I get angry easily", "I get irritated easily" and "I panic easily". The five Neuroticism answers are correlated because they share a cause, and once the cause is known they are independent: this is local independence, the fork of Chapter 1 with an unobserved common cause.

This direction is an assumption, not a finding. Psychometricians call it a reflective model: the items reflect the trait. The alternative is a formative model, in which the items cause the construct: socioeconomic status is made up of income, education and occupation, and a rise in income raises status without any common cause behind it. A factor model fitted to formative indicators will happily find a "factor", because correlated variables always have a first principal axis, but that factor is a summary, not a cause, and interventions on it make no sense. A second warning applies to both: the latent variable is defined by the model. If we had asked different questions, we would have found different factors. The Big Five are a description of what people say about themselves in English-language questionnaires, and whether they are also causal entities in the brain is a question no factor analysis can settle (Borsboom, Mellenbergh and van Heerden, 2003).

11.3 From principal components to factors

Principal component analysis (PCA) finds the directions in the 25-dimensional space of answers along which people vary most. Mathematically, these are the eigenvectors of the correlation matrix, and each eigenvalue is the variance along its direction. No probability model is involved. If the five-trait story is right, we expect five eigenvalues to stand out and the other twenty to be small and similar.

ev <- eigen(cor(Y))$values
tibble(component = 1:10, eigenvalue = round(ev[1:10], 2),
       cumulative_share = round(cumsum(ev)[1:10] / sum(ev), 3))

# A tibble: 10 x 3
   component eigenvalue cumulative_share
       <int>      <dbl>            <dbl>
 1         1       5.13            0.205
 2         2       2.75            0.315
 3         3       2.14            0.401
 4         4       1.85            0.475
 5         5       1.55            0.537
 6         6       1.07            0.58 
 7         7       0.84            0.614
 8         8       0.8             0.646
 9         9       0.72            0.674
10        10       0.69            0.702

Five eigenvalues stand clearly above 1 (5.13, 2.75, 2.14, 1.85 and 1.55), the sixth is 1.07, and from the seventh on they decline slowly from 0.84. Together the first five components account for 54% of the variance in the 25 items. That supports a handful of dimensions, and five is a reasonable reading, though the cut-off is not sharp: the classical rule of keeping eigenvalues above 1 would say six. Notice too that the first eigenvalue is much larger than the rest; the rotation widget below shows why.

Probabilistic PCA (Tipping and Bishop, 1999) turns this into a generative model. Each person $i$ has a vector of $K$ latent scores $f_i\sim \text{Normal}(0,I_K)$, and her standardised answers are a noisy linear function of them:

The $25\times K$ matrix $\Lambda$ holds the loadings: ${\Lambda }_{jk}$ says how strongly factor $k$ drives item $j$. In probabilistic PCA the noise has the same variance for every item, $\Psi ={\sigma }^{2}I$; in the factor analysis model each item has its own unique variance ${\psi }_j^2$, the part of the answer that the factors do not explain. Integrating out the unobserved $f_i$ gives the distribution of the data:

All the information about the factors is in the covariance matrix: the model says that 300 correlations between 25 items are generated by $25K$ loadings and 25 unique variances. It also says something uncomfortable. If $R$ is any $K\times K$ rotation matrix ($R{R}^{\top }=I$), then $(\Lambda R)(\Lambda R{)}^{\top }=\Lambda {\Lambda }^{\top }$: rotated loadings fit the data exactly as well. The likelihood cannot tell a factor model from any rotation of itself. The widget shows this with two factors.

At 0° the plot shows the solution as fitted, with the loadings constrained to lower-triangular form and E3 and N1 as anchors (Section 11.4 explains the constraint). Because the anchors are strong items of two different traits, this solution is already close to simple structure, and the varimax button moves it by only about a degree: the Neuroticism items lie along factor 2, and the Extraversion and Agreeableness items, with the weaker Conscientiousness and Openness items, along factor 1. That so many items load positively on one direction is why the first eigenvalue is so large. Now rotate by 45°: every item loads on both factors and the picture is hard to read, yet the readout shows that nothing about the fit has changed. Rotation is a choice made for interpretation, and the varimax criterion (Kaiser, 1958), which maximises the variance of the squared loadings within each factor, is one convention for making it.

11.4 A Bayesian factor model

For a Bayesian fit, the rotation problem is more than cosmetic. If the likelihood is unchanged along a continuum of rotations, so is the posterior (with rotation-invariant priors), and the sampler wanders around a ridge of equally good solutions; the chains never agree and the loadings have no meaningful posterior mean. The standard fix (Geweke and Zhou, 1996) is to choose one representative from each set of rotations: constrain $\Lambda$ to be lower triangular with a positive diagonal. Item 1 loads only on factor 1, item 2 only on factors 1 and 2, and so on up to item $K$; from item $K+1$ on, all loadings are free. Any $\Lambda$ can be rotated into this form, so the constraint restricts the parameterisation, not the fit. The first $K$ items are called anchors.

Because the model depends on the data only through the sample covariance matrix $S$, the likelihood can be computed from $S$ in one step rather than person by person: for $N$ people, $\log p(Y)=-\frac{N}{2}[\log |\Sigma |+tr({\Sigma }^{-1}S)]$ plus a constant, with $\Sigma =\Lambda {\Lambda }^{\top }+\Psi$. That makes the model fast even for thousands of respondents. We still compute each person's log-likelihood in generated quantities, for leave-one-out comparison.

fa_code <- "
data { int<lower=1> N; int<lower=1> P; int<lower=1> K; matrix[N, P] Y; }
transformed data { matrix[P, P] S = crossprod(Y) / N; }   // sample covariance (Y is centred)
parameters {
  vector<lower=0>[K] ldiag;                         // positive diagonal of the loadings
  vector[K * (K - 1) / 2 + (P - K) * K] lfree;      // loadings below the diagonal
  vector<lower=0>[P] psi;                           // unique standard deviations
}
transformed parameters {
  matrix[P, K] L = rep_matrix(0, P, K);             // lower triangular: identifies the rotation
  { int idx = 1;
    for (j in 1:K) { L[j, j] = ldiag[j];
      for (i in (j + 1):P) { L[i, j] = lfree[idx]; idx += 1; } } }
}
model {
  matrix[P, P] Sigma = add_diag(tcrossprod(L), square(psi));
  ldiag ~ normal(0, 1); lfree ~ normal(0, 1); psi ~ exponential(1);
  target += -0.5 * N * (log_determinant(Sigma) + trace(mdivide_left_spd(Sigma, S)));
}
generated quantities {
  vector[N] log_lik;
  { matrix[P, P] Lc = cholesky_decompose(add_diag(tcrossprod(L), square(psi)));
    for (n in 1:N) log_lik[n] = multi_normal_cholesky_lpdf(Y[n]' | rep_vector(0, P), Lc); }
}"
fa_mod <- stan_model(model_code = fa_code)

Which items should be the anchors? This matters more than the theory suggests. The constraint forces anchor $k$ to load positively on factor $k$, so an anchor that barely relates to the factor we hope it defines leaves the posterior with two nearly equal modes, one for each sign, and chains can settle in different ones. In a first attempt (before reverse-scoring the negatively worded items, and with anchors chosen by hand), one chain found a solution in which the "Agreeableness" anchor had a loading of 0.01 while the other found 0.41, and $\widehat{R}$ reached 37. The cure is to choose, for each factor, the item that loads most strongly on it in a quick classical solution (maximum-likelihood factor analysis with factanal()), and to start every chain at that classical solution rotated into lower-triangular form.

fit_fa <- function(K) {
  path <- sprintf("fits/ch11_fa%d.rds", K); if (file.exists(path)) return(readRDS(path))
  ml <- factanal(Y, K, rotation = if (K > 1) "varimax" else "none")    # classical solution
  Lml <- unclass(ml$loadings); anchors <- integer(0)
  for (j in 1:K) { a <- abs(Lml[, j]); a[anchors] <- -1; anchors <- c(anchors, which.max(a)) }
  ord <- c(anchors, setdiff(seq_len(ncol(Y)), anchors))              # anchors first
  Lm <- Lml[ord, , drop = FALSE]
  Q <- qr.Q(qr(t(Lm[1:K, , drop = FALSE])))                          # rotate to lower-triangular form
  Li <- Lm %*% Q; Li <- sweep(Li, 2, sign(diag(Li[1:K, , drop = FALSE])), "*")
  init <- function() list(ldiag = as.array(diag(Li[1:K, , drop = FALSE])),
                          lfree = as.array(unlist(lapply(1:K, function(j) Li[(j + 1):ncol(Y), j]))),
                          psi = sqrt(ml$uniquenesses[ord]))
  fit <- sampling(fa_mod, data = list(N = nrow(Y), P = ncol(Y), K = K, Y = Y[, ord]),
                  chains = 4, iter = 1000, seed = 11, refresh = 0, init = init)
  out <- list(fit = fit, items = colnames(Y)[ord]); saveRDS(out, path); out
}
fa5 <- fit_fa(5)
fa5$items[1:5]                                           # the anchors
s5 <- summary(fa5$fit, pars = c("ldiag", "lfree", "psi"))$summary
c(max_rhat = round(max(s5[, "Rhat"]), 3), min_ess = round(min(s5[, "n_eff"])),
  divergences = get_num_divergent(fa5$fit))

[1] "N1" "E2" "C4" "A3" "O3"
   max_rhat     min_ess divergences 
       1.01      953.00        0.00 

The anchors chosen from the classical solution are N1 ("Get angry easily"), E2 ("Find it difficult to approach others", reverse-scored), C4 ("Do things in a half-way manner", reverse-scored), A3 ("Know how to comfort others") and O3 ("Carry the conversation to a higher level"). Started from the classical solution, the four chains agree: the largest $\widehat{R}$ is 1.01, the smallest effective sample size is 953, and there are no divergent transitions. The fit takes under a minute.

The lower-triangular loadings are identified but not interpretable: they depend on which items we chose as anchors. To read the factors, we rotate. We take the varimax rotation of the posterior mean loadings as a reference, and then rotate every posterior draw to match the reference as closely as possible (an orthogonal Procrustes rotation, one singular value decomposition per draw). Because rotation does not change the fit, every rotated draw is an equally valid draw from the posterior; we have only chosen how to look at it.

rotate_draws <- function(res) {
  L <- rstan::extract(res$fit, "L")$L                        # draws x items x factors
  ref <- unclass(varimax(apply(L, c(2, 3), mean))$loadings)
  key <- substr(res$items[apply(abs(ref), 2, which.max)], 1, 1)  # trait letter of each factor's top item
  ref <- sweep(ref, 2, sign(ref[cbind(apply(abs(ref), 2, which.max), seq_len(ncol(ref)))]), "*")
  rot <- apply(L, 1, function(M) { sv <- svd(crossprod(M, ref)); M %*% sv$u %*% t(sv$v) })
  rot <- array(t(rot), dim(L)); dimnames(rot) <- list(NULL, res$items, key); rot
}
lam <- rotate_draws(fa5)
lam_tab <- as.data.frame.table(apply(lam, c(2, 3), mean), responseName = "mean") |>
  bind_cols(lo = as.vector(apply(lam, c(2, 3), quantile, 0.05)),
            hi = as.vector(apply(lam, c(2, 3), quantile, 0.95))) |>
  rename(item = Var1, factor = Var2) |> mutate(trait = substr(item, 1, 1))
lam_tab |> group_by(factor) |> slice_max(abs(mean), n = 3) |> mutate(across(c(mean, lo, hi), ~ round(.x, 2)))

# A tibble: 15 x 6
# Groups:   factor [5]
   item  factor  mean    lo    hi trait
   <fct> <fct>  <dbl> <dbl> <dbl> <chr>
 1 N1    N       0.82  0.79  0.85 N    
 2 N2    N       0.79  0.76  0.82 N    
 3 N3    N       0.72  0.69  0.75 N    
 4 E2    E       0.68  0.64  0.71 E    
 5 E4    E       0.62  0.59  0.65 E    
 6 E1    E       0.59  0.55  0.63 E    
 7 C4    C       0.66  0.62  0.69 C    
 8 C2    C       0.63  0.59  0.66 C    
 9 C5    C       0.58  0.54  0.61 C    
10 A3    A       0.66  0.63  0.7  A    
11 A2    A       0.6   0.57  0.64 A    
12 A5    A       0.58  0.55  0.61 A    
13 O3    O       0.62  0.58  0.65 O    
14 O1    O       0.53  0.49  0.56 O    
15 O5    O       0.51  0.47  0.55 O    

After rotation, each factor is led by the items of one trait, so we can name the factors N, E, C, A and O. The strongest loadings are those of N1 and N2 on Neuroticism (0.82 and 0.79), E2 on Extraversion (0.68), C4 on Conscientiousness (0.66), A3 on Agreeableness (0.66) and O3 on Openness (0.62). The 90% intervals are narrow, about ±0.03: with 2,436 people, the loadings are well determined once the rotation is fixed. The names are ours, not the model's. The model found five directions, and we named them after the items that load on them.

lam_tab |> mutate(item = fct_rev(factor(item)), big = abs(mean) > 0.3) |>
  ggplot(aes(y = item)) +
  geom_vline(xintercept = 0, colour = harbour[["grid"]]) +
  geom_errorbarh(aes(xmin = lo, xmax = hi, colour = big), height = 0) +
  geom_point(aes(x = mean, colour = big), size = 1.3) +
  facet_grid(trait ~ factor, scales = "free_y", space = "free_y",
             labeller = labeller(factor = function(x) paste("factor", x))) +
  scale_colour_manual(values = c(`FALSE` = harbour[["muted"]], `TRUE` = harbour[["sky"]]), guide = "none") +
  scale_x_continuous(breaks = c(0, 0.4, 0.8)) + labs(x = "loading", y = NULL)

Every item loads most strongly on the factor of its own trait, so the five-factor structure is recovered. But it is not perfectly simple. The largest cross-loadings link Extraversion and Agreeableness: E4 ("Make friends easily") loads 0.36 on Agreeableness, and A5 ("Make people feel at ease") 0.36 on Extraversion. N4 ("Often feel blue") loads −0.37 on Extraversion, so it measures low mood as much as emotional reactivity. O4 ("Spend time reflecting on things") is the weakest Openness item (0.37), with small loadings on Neuroticism (0.21) and, negatively, Extraversion (−0.22). Departures like these are what a sixth or seventh factor can absorb.

11.5 How many factors?

Five factors was the theory's choice. The data can be asked, too: fit models with $K=1,\dots ,8$ and compare their expected log predictive density by PSIS-LOO, which here measures how well each model predicts a new person's 25 answers.

fits_k <- map(1:8, fit_fa)
loos_k <- map(fits_k, ~ loo(extract_log_lik(.x$fit)))
cmp <- map_dfr(1:8, function(k) {
  d <- loos_k[[k]]$pointwise[, "elpd_loo"] - loos_k[[5]]$pointwise[, "elpd_loo"]
  tibble(K = k, diff = sum(d), se = sqrt(length(d)) * sd(d),
         max_rhat = max(summary(fits_k[[k]]$fit, pars = c("ldiag", "psi"))$summary[, "Rhat"]))
})
ggplot(cmp, aes(K, diff)) + geom_hline(yintercept = 0, colour = harbour[["muted"]], linetype = 2) +
  geom_pointrange(aes(ymin = diff - 2 * se, ymax = diff + 2 * se), colour = harbour[["sky"]], size = 0.3) +
  geom_line(colour = harbour[["sky"]]) + scale_x_continuous(breaks = 1:8) +
  labs(x = "number of factors K", y = "elpd difference from K = 5")
cmp |> mutate(across(c(diff, se), ~ round(.x, 1)), max_rhat = round(max_rhat, 3))
L2 <- apply(rstan::extract(fits_k[[2]]$fit, "L")$L, c(2, 3), mean); items2 <- fits_k[[2]]$items
rm(fits_k, loos_k); invisible(gc())                      # free memory before compiling the next model

# A tibble: 8 x 4
      K   diff    se max_rhat
  <int>  <dbl> <dbl>    <dbl>
1     1 -4486  138.      1   
2     2 -2482. 102       1   
3     3 -1461.  75.2     1.00
4     4  -728.  57       1.00
5     5     0    0       1.00
6     6   277.  32.4     1.00
7     7   395.  39       1.00
8     8   463.  42.4     1.01

LOO does not stop at five. Each extra factor up to five improves the expected log predictive density by hundreds or thousands of units, and beyond five the gains shrink but stay much larger than their standard errors: +277 for six factors, +395 for seven and +463 for eight, with standard errors from 32 to 42. All eight fits mixed well (largest $\widehat{R}$ 1.01). With 2,436 people, small departures from the five-factor model become detectable: pairs of items with similar wording, items that belong partly to two traits, and six-point answers that are not normally distributed. A sixth or seventh factor may describe these better without being a new personality trait. LOO answers the question "which model predicts a new person's answers best?", and its answer is "more than five factors". Whether the extra factors are psychology or artefact is a substantive question, and Exercise 1 looks at what five factors miss.

Lecture 22 — Mixtures

11.6 Penguins without labels

In Chapter 1 we used the Palmer penguins (Gorman, Williams and Fraser, 2014; Horst, Hill and Gorman, 2020) to predict body mass. Now we hide the one variable a biologist would never be without, the species, and ask whether four measurements (bill length, bill depth, flipper length and body mass) can find it again. The question is a fair test of clustering, because we know the answer: 151 Adélie, 68 Chinstrap and 123 Gentoo penguins, measured on three islands near Palmer Station, Antarctica, between 2007 and 2009.

library(palmerpenguins)
birds <- penguins |> drop_na(bill_length_mm, bill_depth_mm, flipper_length_mm, body_mass_g)
Yp <- scale(as.matrix(select(birds, bill_length_mm, bill_depth_mm, flipper_length_mm, body_mass_g)))
colnames(Yp) <- c("bill_length", "bill_depth", "flipper", "mass")
c(penguins = nrow(Yp), measurements = ncol(Yp))

    penguins measurements 
         342            4 

A mixture model says that each penguin belongs to one of $K$ unobserved groups, $z_i\in \{1,\dots ,K\}$, with probabilities $w_1,\dots ,w_K$, and that its measurements are drawn from a distribution specific to its group:

The group labels $z_i$ are discrete, and Hamiltonian Monte Carlo cannot sample discrete parameters. As with the missing data of Chapter 7, the remedy is to sum them out: the likelihood of one penguin is

which depends only on continuous parameters. Each penguin's group is then recovered after fitting by Bayes' rule, draw by draw:

Averaging over the posterior gives each penguin a soft assignment: a probability for each group, not a single label. Before HMC, mixtures were fitted by two older algorithms that are still worth knowing. The EM algorithm (Dempster, Laird and Rubin, 1977) alternates between computing these membership probabilities (the E step) and re-estimating each group's weight, mean and variance from the penguins weighted by their memberships (the M step); it climbs to a local maximum of the likelihood. Gibbs sampling (Diebolt and Robert, 1994) instead draws each $z_i$ from its membership probabilities and then draws the parameters given the labels, and produces a sample from the posterior. The widget runs both on flipper length alone.

With two components, both algorithms quickly find a small-flipper group (Adélie and Chinstrap) and a large-flipper group (Gentoo). Try a few new starts with three components: EM sometimes stops at a poor solution, for instance splitting the Gentoo into two and leaving Adélie and Chinstrap merged, because it only ever climbs. That is a local mode, and the Bayesian posterior has it too, as the next section shows.

11.7 Chains that disagree

We write the mixture in Stan with full covariance matrices for each group (an LKJ prior on each correlation matrix, as in Chapter 6), a Dirichlet(2) prior on the weights, and the per-penguin membership probabilities computed in generated quantities. Each chain starts from a k-means solution.

mix_code <- "
data { int<lower=1> N; int<lower=1> D; int<lower=1> K; array[N] vector[D] y; }
parameters {
  simplex[K] w;                                   // group weights
  array[K] vector[D] mu;                          // group means
  array[K] cholesky_factor_corr[D] Omega;         // group correlation matrices
  array[K] vector<lower=0>[D] sigma;              // group standard deviations
}
model {
  array[K] matrix[D, D] Lk;
  w ~ dirichlet(rep_vector(2, K));
  for (k in 1:K) {
    mu[k] ~ normal(0, 1.5); sigma[k] ~ exponential(2); Omega[k] ~ lkj_corr_cholesky(2);
    Lk[k] = diag_pre_multiply(sigma[k], Omega[k]);
  }
  for (n in 1:N) {                                // the group label is summed out
    vector[K] lp;
    for (k in 1:K) lp[k] = log(w[k]) + multi_normal_cholesky_lpdf(y[n] | mu[k], Lk[k]);
    target += log_sum_exp(lp);
  }
}
generated quantities {
  vector[N] log_lik; array[N] vector[K] prob;     // membership probabilities by Bayes' rule
  { array[K] matrix[D, D] Lk;
    for (k in 1:K) Lk[k] = diag_pre_multiply(sigma[k], Omega[k]);
    for (n in 1:N) { vector[K] lp;
      for (k in 1:K) lp[k] = log(w[k]) + multi_normal_cholesky_lpdf(y[n] | mu[k], Lk[k]);
      log_lik[n] = log_sum_exp(lp); prob[n] = softmax(lp); } }
}"
mix_path <- "fits/ch11_mix_model.rds"                  # cache the compiled model between sessions
if (file.exists(mix_path)) mix_mod <- readRDS(mix_path) else { mix_mod <- stan_model(model_code = mix_code); saveRDS(mix_mod, mix_path) }
fit_mix <- function(Y, K, tag) {
  path <- sprintf("fits/ch11_%s%d.rds", tag, K); if (file.exists(path)) return(readRDS(path))
  set.seed(11); km <- kmeans(Y, K, nstart = 50)            # starting values from k-means
  init <- function() list(w = as.array(as.numeric(table(km$cluster)) / nrow(Y)),
    mu = lapply(1:K, function(k) as.array(km$centers[k, ])),
    sigma = lapply(1:K, function(k) as.array(pmax(apply(Y[km$cluster == k, , drop = FALSE], 2, sd), 0.1))),
    Omega = lapply(1:K, function(k) diag(ncol(Y))))
  fit <- sampling(mix_mod, data = list(N = nrow(Y), D = ncol(Y), K = K, y = lapply(seq_len(nrow(Y)), function(i) Y[i, ])),
                  chains = 4, iter = 1000, seed = 11, refresh = 0, init = init)
  saveRDS(fit, path); fit
}
m3 <- fit_mix(Yp, 3, "penguin")
c(max_rhat = round(max(summary(m3, pars = c("w", "mu"))$summary[, "Rhat"]), 2), divergences = get_num_divergent(m3))

   max_rhat divergences 
      23.48        0.00 

The largest $\widehat{R}$ is 23.5, although there are no divergent transitions. By every rule used in this course, these chains have not converged.

The diagnostics say the four chains disagree. Before trying to fix anything, look at what each chain found: its average log posterior density, its group weights, and how the penguins it puts in each group line up with the species we hid.

per_chain <- function(fit, K) {                  # posterior-mean memberships, one N x K matrix per chain
  arr <- as.array(fit); pn <- grep("^prob\\[", dimnames(arr)[[3]], value = TRUE)
  map(seq_len(dim(arr)[2]), ~ matrix(colMeans(arr[, .x, pn]), ncol = K))
}
lp3 <- colMeans(as.array(m3, pars = "lp__")[, , 1]); P3 <- per_chain(m3, 3)
round(rbind(lp = lp3, sapply(1:4, function(ch) colMeans(as.array(m3, pars = "w")[, ch, ]))), 2)
for (chain in 1:4) { cat("chain", chain, "\n"); print(table(birds$species, group = max.col(P3[[chain]]))) }

      chain:1  chain:2  chain:3  chain:4
lp   -1210.43 -1210.78 -1210.59 -1210.94
w[1]     0.45     0.44     0.36     0.44
w[2]     0.20     0.36     0.45     0.20
w[3]     0.36     0.20     0.20     0.36
chain 1 
           group
              1   2   3
  Adelie    149   2   0
  Chinstrap   4  64   0
  Gentoo      0   0 123
chain 2 
           group
              1   2   3
  Adelie    149   0   2
  Chinstrap   4   0  64
  Gentoo      0 123   0
chain 3 
           group
              1   2   3
  Adelie      0 149   2
  Chinstrap   0   4  64
  Gentoo    123   0   0
chain 4 
           group
              1   2   3
  Adelie    149   2   0
  Chinstrap   4  64   0
  Gentoo      0   0 123

All four chains reached the same log posterior density (−1210.4 to −1210.9) and put the same penguins together: one group holds 149 Adélie and 4 Chinstrap, another 64 Chinstrap and 2 Adélie, and the third all 123 Gentoo. Only the numbering differs. Gentoo is group 3 in chains 1 and 4, group 2 in chain 2 and group 1 in chain 3. The chains agree completely about the penguins and disagree only about the names.

Two different problems are hiding in that output. The first is label switching. The likelihood is unchanged if we swap the names of two groups: nothing in the data says which group is "1". So a posterior with $K$ well-separated groups has $K!$ identical copies of every mode, one for each labelling, and different chains can land in different copies. Averaging their draws then blends Gentoo with Chinstrap into nonsense, and $\widehat{R}$ rightly complains. This is not a problem with the fit, only with its bookkeeping. A common remedy is to constrain one coordinate of the means to be ordered (${\mu }_{1,1}<{\mu }_{2,1}<{\mu }_{3,1}$), but it works only if the groups are well separated in that coordinate; Chinstrap and Gentoo penguins have almost the same average bill length, so ordering by bill length would not have helped. The widget shows label switching within a single chain, which happens when groups overlap.

Relabelling fixes the first problem. We take the chain with the highest log posterior as the reference and, for each other chain, choose the permutation of its labels that best matches the reference's membership probabilities (Stephens, 2000, does this draw by draw; chain by chain suffices when labels do not switch within a chain, which the clean per-chain tables above indicate).

The second problem did not occur in this run, but it does in others, and it is not bookkeeping. In an earlier run with a different random seed, one of the four chains settled in a genuinely different solution, a local mode like the one EM gets stuck in: Adélie and Chinstrap merged into one group and a few birds in a small group of their own, with a log posterior density about 50 units below the others. No relabelling turns that into the same answer. Such a mode carries essentially no posterior probability, since $e^{-50}$ is negligible, but a chain that falls into it cannot climb out. Our align_chains() function therefore keeps only the chains whose average log density is within 10 units of the best, and reports how many it kept. With a minor mode closer in log density, the principled alternative is to weight the chains by their predictive performance (stacking; Yao, Vehtari and Gelman, 2022). We will need the rule for four groups.

perms <- function(K) if (K == 1) list(1) else do.call("c", lapply(1:K, function(i) lapply(perms(K - 1), function(p) c(i, ifelse(p >= i, p + 1, p)))))
align_chains <- function(fit, K, gap = 10) {
  lp <- colMeans(as.array(fit, pars = "lp__")[, , 1]); P <- per_chain(fit, K)
  ref <- which.max(lp); keep <- which(lp > max(lp) - gap)            # chains in the dominant mode
  pm <- map(seq_along(P), function(ch) { pp <- perms(K); pp[[which.max(map_dbl(pp, ~ sum(P[[ref]] * P[[ch]][, .x])))]] })
  prob <- Reduce(`+`, map(keep, ~ P[[.x]][, pm[[.x]]])) / length(keep)
  ll <- as.array(fit, pars = "log_lik")[, keep, , drop = FALSE]
  list(prob = prob, keep = keep, lp = lp, loo = suppressWarnings(loo(ll, r_eff = relative_eff(exp(ll)))))
}
a3 <- align_chains(m3, 3)
c(chains_kept = length(a3$keep))
table(birds$species, group = max.col(a3$prob))
c(certain_0.9 = sum(apply(a3$prob, 1, max) > 0.9), uncertain = sum(apply(a3$prob, 1, max) <= 0.9))

chains_kept 
          4 
           group
              1   2   3
  Adelie    149   2   0
  Chinstrap   4  64   0
  Gentoo      0   0 123
certain_0.9   uncertain 
        329          13 

After relabelling, all four chains are kept and the pooled assignments recover the species almost exactly: 149 of 151 Adélie, 64 of 68 Chinstrap and all 123 Gentoo fall in groups of their own, so only 6 of the 342 penguins are grouped with another species, all of them between Adélie and Chinstrap. 329 penguins are assigned with probability above 0.9; 13 are borderline.

11.8 Soft assignments and choosing K

The figure shows every penguin in the two measurements that separate the species best, coloured by its most probable group, with paler points for less certain assignments, and with the hidden species as the marker shape.

birds |> mutate(group = factor(max.col(a3$prob)), certainty = apply(a3$prob, 1, max)) |>
  ggplot(aes(bill_length_mm, flipper_length_mm, colour = group, shape = species, alpha = certainty)) +
  geom_point(size = 1.9) +
  scale_colour_manual(values = unname(harbour[c("sky", "amber", "rose")])) +
  scale_alpha_continuous(range = c(0.25, 1), limits = c(0.5, 1)) +
  guides(alpha = "none") + theme(legend.position = "right") +
  labs(x = "bill length (mm)", y = "flipper length (mm)", colour = "group", shape = "species")

The Gentoo group is unmistakable, with its long flippers (and, in the measurements not shown, shallower bills and greater weight). Adélie and Chinstrap differ mainly in bill length, and the pale points, the uncertain assignments, sit where their clouds meet.

unc <- order(apply(a3$prob, 1, max))[1:15]
as_tibble(a3$prob[unc, ], .name_repair = ~ paste("group", 1:3)) |>
  mutate(bird = factor(paste0("#", unc, " (", birds$species[unc], ", ", birds$sex[unc], ")"), levels = rev(paste0("#", unc, " (", birds$species[unc], ", ", birds$sex[unc], ")")))) |>
  pivot_longer(-bird, names_to = "group", values_to = "p") |>
  ggplot(aes(p, bird, fill = group)) + geom_col(width = 0.75) +
  scale_fill_manual(values = unname(harbour[c("sky", "amber", "rose")])) +
  labs(x = "posterior probability", y = NULL, fill = NULL)

The borderline penguins are all on the Adélie–Chinstrap boundary: 11 of the 15 are female Chinstrap and 4 are male Adélie. Small Chinstrap females and large Adélie males overlap in every measurement, so much of the uncertainty is explained by sex, a variable the model never sees. A biologist would recognise these birds as hard cases, and the model says so instead of hiding them in a confident label.

How many groups would the data have chosen without our hint? We fit $K=1,\dots ,4$ and compare by PSIS-LOO (on one core, the four-group model alone takes about 20 minutes, which is itself a warning sign), using the chains of the dominant mode for each.

mk <- map(1:4, ~ fit_mix(Yp, .x, "penguin"))
ak <- map2(mk, 1:4, align_chains)
tibble(K = 1:4, chains_kept = map_int(ak, ~ length(.x$keep)),
       lp_range = map_dbl(ak, ~ round(max(.x$lp) - min(.x$lp), 1)),
       elpd = map_dbl(ak, ~ .x$loo$estimates["elpd_loo", "Estimate"]),
       se = map_dbl(ak, ~ .x$loo$estimates["elpd_loo", "SE"]),
       k_above_0.7 = map_int(ak, ~ sum(.x$loo$diagnostics$pareto_k > 0.7))) |>
  mutate(elpd_vs_3 = elpd - elpd[3], across(c(elpd, se, elpd_vs_3), ~ round(.x, 1)))
table(paste(birds$species, birds$sex), group = max.col(ak[[4]]$prob))

# A tibble: 4 x 7
      K chains_kept lp_range   elpd    se k_above_0.7 elpd_vs_3
  <int>       <int>    <dbl>  <dbl> <dbl>       <int>     <dbl>
1     1           4      0.6 -1530   23.2           0    -336  
2     2           4      0.1 -1236.  27.2           0     -42.6
3     3           4      0.5 -1194.  28.9           1       0  
4     4           1     70.7 -1073.  40.1           1     121. 
                  group
                    1  2  3  4
  Adelie NA         4  0  1  0
  Adelie female    68  0  5  0
  Adelie male      70  2  1  0
  Chinstrap female  6 25  3  0
  Chinstrap male    0 33  1  0
  Gentoo NA         0  0  0  4
  Gentoo female     0  0  0 58
  Gentoo male       0  0  0 61

The comparison does not stop at three. Two groups are much worse than three (−42.6), but four groups beat three by 121 units of expected log predictive density, about three times the standard errors. That is not the end of the story. With four groups the chains disagreed: their average log densities spanned 70.7 units, only one chain was in the best mode, and so the four-group estimate rests on 500 draws from a single chain, with no guarantee that an even better mode was not missed. Nor is the fourth group a new species, or one sex of a species. Gentoo stays one group; most Chinstrap (with two Adélie males) form another; most Adélie (with six female Chinstrap) a third; and the fourth holds 11 birds, mostly females of the two smaller species. A component like this improves the description of the density by catching birds that sit away from their species' main cloud, but it is not a natural kind. The lesson is the factor model's lesson again: LOO says how many components describe the data best, which is not the same as how many kinds of thing there are.

Choosing $K$ by comparing finite mixtures has a Bayesian alternative: put a prior on partitions with no fixed number of groups. The Dirichlet process (Ferguson, 1973) does this, and its behaviour is easiest to see through the Chinese restaurant process: penguins arrive one at a time; each joins an existing group with probability proportional to its size, or starts a new group with probability proportional to a concentration parameter $\alpha$. The expected number of groups after $n$ arrivals is ${\sum }_{i=1}^n\alpha /(\alpha +i-1)$, which grows like $\alpha \log n$. In a Dirichlet process mixture the data then decide how many of these groups are occupied.

Two features of the prior matter in practice. Group sizes are very unequal: a few large groups and a tail of small ones, a "rich get richer" pattern that suits some data (species abundances, word frequencies) and not others. And the number of groups keeps growing with $n$, so a Dirichlet process mixture fitted to a large data set from three true groups tends to report several extra small ones (Miller and Harrison, 2013). That is not a bug in the method, but it means the posterior number of groups is not an estimate of "the" number of clusters; it is the number of components needed to describe the density.

11.9 Subtypes from expression

The cancer gene-expression project asks the penguin question of tumours: do expression profiles fall into groups, and are they the known molecular subtypes? We return to the TCGA breast-cancer data of Chapter 10 (Cancer Genome Atlas Network, 2012; Goldman et al., 2020), now taking every primary tumour with a PAM50 subtype call, 844 in all. PAM50 is a 50-gene classifier that assigns each tumour to one of five intrinsic subtypes, Luminal A, Luminal B, HER2-enriched, Basal-like and Normal-like (Parker et al., 2009); the calls play the part of the hidden penguin species. No labels are used in fitting, so there is no training and test split. With 1,000 genes and 844 tumours a full-covariance mixture is out of reach, so we combine this chapter's two tools: reduce the 1,000 most variable genes to their first three principal components, then fit a three-group mixture in that space with the code of Section 11.7.

xena <- "https://tcga-xena-hub.s3.us-east-1.amazonaws.com/download/TCGA.BRCA.sampleMap%2F"
dir.create("data/tcga", showWarnings = FALSE, recursive = TRUE)
for (f in c("HiSeqV2.gz", "BRCA_clinicalMatrix"))
  if (!file.exists(file.path("data/tcga", f))) download.file(paste0(xena, f), file.path("data/tcga", f))
expr <- read_tsv("data/tcga/HiSeqV2.gz", show_col_types = FALSE)
clin <- read_tsv("data/tcga/BRCA_clinicalMatrix", show_col_types = FALSE, guess_max = 5000)
X <- t(as.matrix(expr[, -1])); colnames(X) <- expr$sample
pam <- clin |> filter(substr(sampleID, 14, 15) == "01", !is.na(PAM50Call_RNAseq))    # primary tumours
X <- X[intersect(rownames(X), pam$sampleID), ]
subtype <- pam$PAM50Call_RNAseq[match(rownames(X), pam$sampleID)]
her2 <- pam$HER2_Final_Status_nature2012[match(rownames(X), pam$sampleID)]
table(subtype)
G <- X[, order(apply(X, 2, var), decreasing = TRUE)[1:1000]]                       # no labels used
pc <- prcomp(G, scale. = TRUE)
round(summary(pc)$importance[2:3, 1:5], 3)
Z <- scale(pc$x[, 1:3])
mt <- fit_mix(Z, 3, "xena"); at <- align_chains(mt, 3)
c(chains_kept = length(at$keep), lp_range = round(max(at$lp) - min(at$lp), 1))
table(subtype, group = max.col(at$prob))
c(uncertain = sum(apply(at$prob, 1, max) <= 0.9))

subtype
 Basal   Her2   LumA   LumB Normal 
   141     67    421    192     23 
                         PC1   PC2   PC3   PC4   PC5
Proportion of Variance 0.161 0.074 0.043 0.032 0.026
Cumulative Proportion  0.161 0.235 0.278 0.310 0.336
chains_kept    lp_range 
        4.0         0.6 
        group
subtype    1   2   3
  Basal  135   0   6
  Her2     0   0  67
  LumA     0 314 107
  LumB     0  39 153
  Normal   5  12   6
uncertain 
      279 

The first three principal components explain 16%, 7% and 4% of the variance in the 1,000 genes, 28% together. The four chains agree after relabelling (log densities within 0.6 of each other). The three groups follow the best-known divisions in breast cancer. Group 1 is Basal-like: it holds 135 of the 141 Basal tumours and nothing else apart from five Normal-like ones. Basal tumours differ from the rest across thousands of genes, so they dominate the first component and separate cleanly, just as the Gentoo penguins did. Group 2 is mostly Luminal A (314 tumours, with 39 Luminal B). Group 3 collects all 67 HER2-enriched tumours and 153 of the 192 Luminal B, but also 107 Luminal A. What separates groups 2 and 3 is roughly the difference between slowly and rapidly dividing tumours, and the PAM50 classes on that axis overlap rather than forming islands. And 279 of the 844 tumours have no group with probability above 0.9: a third of the tumours sit between groups, against 13 of 342 penguins.

mt4 <- fit_mix(Z, 4, "xena"); at4 <- align_chains(mt4, 4)
c(chains_kept = length(at4$keep), lp_range = round(max(at4$lp) - min(at4$lp), 1))
loo_compare(list(K3 = at$loo, K4 = at4$loo))
table(subtype, group = max.col(at4$prob))
c(uncertain = sum(apply(at4$prob, 1, max) <= 0.9))
table(clinical_her2 = her2, group_K3 = max.col(at$prob), useNA = "ifany")

chains_kept    lp_range 
        4.0         0.3 
   elpd_diff se_diff
K4   0.0       0.0  
K3 -25.1       8.2  
        group
subtype    1   2   3   4
  Basal    1   0 136   4
  Her2    26   0   0  41
  LumA   139 280   0   2
  LumB   164  27   0   1
  Normal   2  12   4   5
uncertain 
      394 
             group_K3
clinical_her2   1   2   3
    Equivocal   1   0   8
    Negative  125 299 210
    Positive    3  21  89
    <NA>       11  45  32

PSIS-LOO prefers four groups to three by 25.1 units of expected log predictive density, about three standard errors. The fourth group is mostly HER2-enriched: 41 of its 53 tumours, and 41 of the 67 HER2-enriched tumours overall, with the other 26 joining a group of Luminal B and Luminal A. The Basal group is unchanged, but the extra group buys a better density at the price of more ambiguity: 394 tumours are now uncertain. The last table compares the three groups with clinical HER2 status, the outcome of Chapter 10: group 3 holds 89 of the 113 clinically positive tumours, but they are a minority of it (210 are negative). That is an honest, unglamorous result. HER2 status is decided by the amplification of one small region of chromosome 17, a strong signal in a handful of genes, and a variance-based summary of a thousand genes does not emphasise it; the supervised horseshoe of Chapter 10 found it in one gene. Unsupervised groups are groups in whatever the reduction emphasises.

Lecture 23 — Copulas

11.10 Two genes on one amplicon

Chapter 10 left a puzzle that this lecture takes up. In its Exercise 6 the horseshoe could not decide between ERBB2 and GRB7, two genes a few tens of kilobases apart on chromosome 17 that are amplified together in HER2-positive tumours. Their correlation across tumours is high, but a correlation coefficient is a poor description of how they move together. Look at them across the 1,097 primary tumours of the TCGA data, and at a second pair, ESR1, the gene for the oestrogen receptor, and GATA3, a transcription factor that works with it.

library(copula)
xena <- "https://tcga-xena-hub.s3.us-east-1.amazonaws.com/download/TCGA.BRCA.sampleMap%2F"
dir.create("data/tcga", showWarnings = FALSE, recursive = TRUE)
if (!file.exists("data/tcga/HiSeqV2.gz")) download.file(paste0(xena, "HiSeqV2.gz"), "data/tcga/HiSeqV2.gz")
expr <- read_tsv("data/tcga/HiSeqV2.gz", show_col_types = FALSE)          # log2(RSEM + 1)
G <- expr |> filter(sample %in% c("ERBB2", "GRB7", "ESR1", "GATA3")) |>
  pivot_longer(-sample, names_to = "tumour") |> pivot_wider(names_from = sample) |>
  filter(substr(tumour, 14, 15) == "01")                                  # primary tumours
nrow(G)
dep <- function(x, y) c(pearson = cor(x, y), spearman = cor(x, y, method = "spearman"),
                        kendall = cor(x, y, method = "kendall"))
rbind(`ERBB2-GRB7, log scale` = dep(G$ERBB2, G$GRB7), `ERBB2-GRB7, raw scale` = dep(2^G$ERBB2, 2^G$GRB7),
      `ESR1-GATA3, log scale` = dep(G$ESR1, G$GATA3), `ESR1-GATA3, raw scale` = dep(2^G$ESR1, 2^G$GATA3)) |> round(3)

[1] 1097
                      pearson spearman kendall
ERBB2-GRB7, log scale   0.843    0.678   0.510
ERBB2-GRB7, raw scale   0.871    0.678   0.510
ESR1-GATA3, log scale   0.835    0.717   0.525
ESR1-GATA3, raw scale   0.560    0.717   0.525

The two pairs have almost the same rank correlations: Kendall's $\tau$ is 0.510 for ERBB2 and GRB7 and 0.525 for ESR1 and GATA3, and Spearman's 0.678 and 0.717. Pearson's correlation is higher, 0.843 and 0.835, on the log scale on which expression is usually analysed; on the original scale of normalised counts it is 0.871 for the first pair but only 0.560 for the second. The rank correlations do not change at all. A dependence measure that changes when the data are re-expressed is describing the margins as much as the dependence.

pobs <- function(x) rank(x) / (length(x) + 1)                 # pseudo-observations
pairs_long <- bind_rows(
  tibble(pair = "ERBB2 and GRB7", x = G$ERBB2, y = G$GRB7),
  tibble(pair = "ESR1 and GATA3", x = G$ESR1, y = G$GATA3)) |>
  mutate(u = pobs(x), v = pobs(y), .by = pair)
p_log <- ggplot(pairs_long, aes(x, y)) + geom_point(size = 0.5, alpha = 0.45, colour = harbour[["sky"]]) +
  facet_wrap(~ pair, ncol = 1, scales = "free") + labs(x = "first gene, log2 expression", y = "second gene")
p_rank <- ggplot(pairs_long, aes(u, v)) + geom_point(size = 0.5, alpha = 0.45, colour = harbour[["sky"]]) +
  annotate("rect", xmin = 0.95, xmax = 1, ymin = 0.95, ymax = 1, fill = NA, colour = harbour[["rose"]], linetype = 2) +
  annotate("rect", xmin = 0, xmax = 0.05, ymin = 0, ymax = 0.05, fill = NA, colour = harbour[["rose"]], linetype = 2) +
  facet_wrap(~ pair, ncol = 1) + coord_equal() + labs(x = "rank of first gene", y = "rank of second gene")
patchwork::wrap_plots(p_log, p_rank, widths = c(1.15, 1))