Week 8: Mixture Models

Previously

• Multilevel/Hierarchical models
  • Impose structure on parameters by way of prior distributions
  • Information from grouping (partial pooling/“random effects”)
  • Strength of the prior determines the degree of pooling
• Model checking
  • Posterior predictive checks: can the model reproduce features of the data?
  • Comparing models: information criteria and cross-validation
• MrP: multilevel regression + post-stratification to adjust non-representative surveys

This week

• “Exploratory” modeling
  • Unsupervised “clustering” algorithms
  • K-means clustering
  • Finite mixture models
  • Obtaining MLEs or posterior modes via Expectation-Maximization
• Running example: Policy attitudes in the 2023 CES
  • Are there distinct “types” of voters in terms of their policy attitudes?
  • Next week we’ll ask the same question with a continuous latent variable (item response theory)

Finite Mixture Models

• In conventional regression models, we assume a known distribution for the outcome given a set of covariates

  • The covariates could be something like a discrete “group” label (e.g. \(X_i = x\))

• In this setting, conditional on \(X_i\), we would model each \(Y_i\) with a known distribution (e.g. Normal)

  \[Y_i | X_i = x \sim \text{Normal}(\mu_{x}, \sigma^2_{x})\]

  \(\mu_{x}\) and \(\sigma^2_{x}\) are the mean and variance parameters associated with the group \(X_i = x\).

• Now imagine that the group indicators are not observed but instead are latent variables

  • We can do inference on these conditional on the observed data and the model.

Finite Mixture Models

• In a finite mixture model, we assume that the outcome distribution for each observation \(i\) is governed by some discrete latent indicator \(z_i \in \{1, 2, \dotsc, K\}\)

  • \(z_1, z_2, \dotsc, z_N\) are i.i.d. \(\sim \text{Categorical}(\pi)\)

• Conditional on the latent variable \(z_i = k\), the outcome distribution is known.

  • For a continuous outcome, a gaussian mixture: \(Y_i | z_i = k \sim \text{Normal}(\mu_{k}, \sigma^2_{k})\)
  • For a binary outcome, a Bernoulli mixture: \(Y_i | z_i = k \sim \text{Bernoulli}(\beta_{k})\)

• But unconditionally, \(Y_i\) has a mixture distribution: its density is a weighted average of the component densities. In the gaussian case:

  \[Y_i \sim \sum_{k=1}^K \pi_k \times \text{Normal}(\mu_{k}, \sigma^2_{k})\]

  where \(\pi_k = Pr(z_i = k)\)

• Our goal is to estimate the parameters of each component and the mixing proportions \(\pi\) using likelihood inference.

  • Put priors on the component parameters and \(\pi\) to make it fully Bayesian.

Application: Policy attitudes

• The organization of voters’ policy views is at the heart of political science.
  • Are voters’ attitudes well-summarized by a single dimension (“left” vs. “right”)?
  • Or are many voters cross-pressured?
  • Is there any structure to voters’ attitudes (Converse 1964)?
• Recent methodological debate in the APSR
  • Fowler et al. (2023): most moderates are genuine centrists on a single dimension - few “Conversian” or Inattentive voters.
  • Broockman and Lauderdale (2025): this is a modelling artifact - moderates are actually much less constrained!
  • These next two weeks are about the statistical model(s) at the core of this debate.
• Surveys ask respondents batteries of yes/no policy questions
  • Can we recover groups of respondents with similar response patterns without telling the model anything about party or ideology?

Data: CES 2023

• The 2023 Cooperative Election Study (Schaffner et al. 2024): the full public release, as downloaded from Dataverse

library(haven)
ces_raw <- read_dta("data/CCES23_Common_OUTPUT.dta.gz") # haven reads the compressed Stata file directly
dim(ces_raw)

[1] 24500   320

• Stata files carry value labels: each policy question is coded 1 = Support, 2 = Oppose, 8 = skipped, 9 = not asked

count(ces_raw, CC23_324c) # "Make abortions illegal in all circumstances"

# A tibble: 3 × 2
  CC23_324c       n
  <dbl+lbl>   <int>
1 1 [Support]  4539
2 2 [Oppose]  19955
3 8 [skipped]     6

• The questionnaire (data/CES 2023 CC Questionnaire.pdf) has the full wording of every item

Data: CES 2023

• 24 “support/oppose” policy questions on guns and policing, immigration, abortion, the environment, and social welfare

items <- c(ban_assault_rifles = "CC23_321a", concealed_carry = "CC23_321b", background_checks = "CC23_321c",
           more_police = "CC23_321d", fewer_police = "CC23_321e", legal_status = "CC23_323a",
           border_patrols = "CC23_323b", border_wall = "CC23_323c", dreamers = "CC23_323d",
           abortion_choice = "CC23_324a", abortion_exceptions = "CC23_324b", abortion_ban = "CC23_324c",
           abortion_expand = "CC23_324d", epa_co2 = "CC23_326a", renewable_mandate = "CC23_326b",
           epa_enforcement = "CC23_326c", more_fossil_fuels = "CC23_326d", halt_leases = "CC23_326e",
           relax_zoning = "CC23_328a", housing_credits = "CC23_328b", medicaid_work_req = "CC23_328c",
           repeal_aca = "CC23_328d", expand_medicaid = "CC23_328e", forgive_loans = "CC23_328f")
ces <- ces_raw %>%
  transmute(caseid, commonweight, pid3 = as.character(as_factor(pid3)), ideo5 = as.character(as_factor(ideo5)),
            across(all_of(items), ~ case_when(.x == 1 ~ 1, .x == 2 ~ 0))) %>% # skipped/not asked become NA
  drop_na(all_of(names(items)))
Y <- as.matrix(select(ces, all_of(names(items)))) # 1 = Support, 0 = Oppose

• 24,453 of the 24,500 respondents answered all 24 questions

Data: CES 2023

Clustering

• Each respondent is a vector of 24 ones and zeroes.

K-means

• A simple algorithm for finding clusters that has no underlying probabilistic model is the K-means algorithm
  • Straightforward to implement and still surprisingly popular and effective for a “first cut” at the data.
• K-means algorithm
  1. Initialize \(K\) distinct clusters by randomly assigning points to one of the \(K\) groups
  2. Calculate \(\mu_k\) as the mean of \(Y_i\) in observations in cluster \(k\).
  3. Reassign each point to the cluster \(k\) that has the smallest distance between \(Y_i\) and \(\mu_k\).
  4. Repeat 2-3 until no points change their assignments.
• With binary responses, \(\mu_k\) is the vector of proportions supporting each policy among respondents in cluster \(k\)

K-means example

• R’s kmeans() implements the algorithm; nstart re-runs it from multiple random starts and keeps the result with smallest within-cluster sum of squares.

set.seed(53703)
km2 <- kmeans(Y, centers = 2, nstart = 10)
km2$size

[1]  8803 15650

• Sorting our 150 respondents by their assigned cluster

K-means example

• The cluster centroids \(\mu_k\): share of each cluster supporting each policy

K-means example

• How do the clusters compare to respondents’ party identification?

table(cluster = km2$cluster, party = ces$pid3)

       party
cluster Democrat Independent Not sure Other Republican
      1      569        2605      348   450       4831
      2     8257        3812     1164   467       1950

• The clusters themselves don’t have any underlying interpretation or meaning
  • We assign the interpretation through inspection.

K-means

• Advantages
  • Easy to implement, fast updates, and quick convergence.
• Disadvantages
  • Multimodality (we can’t really solve this one!)
  • Choice of distance metric matters (e.g. should we standardize the items somehow?)
  • Euclidean distance is an awkward fit for binary/discrete outcomes
• With conventional K-means, we have no underlying probability model – each unit is assigned to a single cluster.
  • Sometimes called “hard” K-means.
• We can instead implement an approach sometimes referred to as a “soft” K-means algorithm
  • Our target of inference is the probability that a unit belongs to each cluster.

Bernoulli Mixture Models

• We observe \(N\) respondents each giving binary responses to \(J\) questions. \(Y_i\) is the vector of responses for unit \(i\) with components \(Y_{ij}\).

• Assume latent variables \(z = \{z_1, z_2, \dotsc, z_N\}\)

  \[z_i \underset{\text{i.i.d.}}{\sim} \text{Categorical}(\pi)\]

• Conditional on \(z_i = k\), each response is independently Bernoulli with a cluster-and-question-specific probability \(\beta_{kj}\)

  \[Y_{ij} | z_i = k \sim \text{Bernoulli}(\beta_{kj})\]

• Our targets of inference are the \(K \times J\) matrix of response probabilities \(\beta\) and the cluster prevalences \(\pi\)

  • And conditional on these estimates, we can obtain the cluster probabilities \(Pr(z_i = k | Y_i)\)

Local independence

• The key assumption: conditional on cluster membership, a respondent’s answers are independent

  \[f(Y_i | z_i = k) = \prod_{j=1}^J \beta_{kj}^{Y_{ij}}(1 - \beta_{kj})^{1 - Y_{ij}}\]

• Answers are certainly not independent marginally

  • Supporting an assault rifle ban is strongly correlated with supporting background checks

• The mixture model says all of that correlation comes from respondents belonging to different clusters

  • Within a cluster, knowing someone supports a ban tells you nothing about how they’ll answer on background checks
  • This is a strong assumption!

Expectation-Maximization

• Let’s write down the marginal likelihood, summing out the latent variables:

  \[\mathcal{L}(\beta, \pi; \mathbf{Y}) = \prod_{i=1}^N \sum_{k=1}^K f(Y_i| \beta_{k}, z_i = k) \times p(z_i = k | \pi)\]

• This is a tricky likelihood to maximize - if we take the log, we get a log of a sum (which doesn’t simplify as neatly as a log of a product)

  \[\ell(\beta, \pi ; \mathbf{Y}) = \sum_{i=1}^N \log\bigg(\sum_{k=1}^K \pi_k \prod_{j=1}^J \beta_{kj}^{Y_{ij}}(1 - \beta_{kj})^{1 - Y_{ij}} \bigg)\]

• Suppose we knew the \(z_i\) as well (they were like “data”), then the “complete” log-likelihood would be easier to maximize!

  \[\ell(\beta, \pi ; \mathbf{Y}, \mathbf{Z}) = \sum_{i=1}^N \sum_{k=1}^K I(z_i = k) \bigg[\log \pi_k + \sum_{j=1}^J \Big(Y_{ij} \log \beta_{kj} + (1 - Y_{ij}) \log (1 - \beta_{kj})\Big)\bigg]\]

Expectation-Maximization

• The Expectation-Maximization algorithm (Dempster, Laird, and Rubin 1977) is an iterative method for obtaining a maximum likelihood estimate (or maximum a posteriori (MAP) in a Bayesian setting) when the likelihood consists of a sum/integral over the distribution of some latent variables

  • It uses the idea of iteratively optimizing a lower bound on the likelihood until convergence is reached.
  • The key trick is Jensen’s inequality. For a concave function like the logarithm:

  \[\log(E[X]) \ge E[\log(X)]\]

• Consider the general setting where we have a parameter \(\theta\), data \(\mathbf{X}\) and discrete latent variables \(\mathbf{Z}\)

• We want to optimize the log-likelihood:

  \[\ell(\theta; \mathbf{X}) = \log \sum_{\mathbf{Z}} f(\mathbf{X}, \mathbf{Z} | \theta) = \log \sum_{\mathbf{Z}} f(\mathbf{X} | \mathbf{Z}, \theta) p(\mathbf{Z} | \theta)\]

Expectation-Maximization

• Instead of optimizing the log-likelihood directly, we’ll optimize its lower bound. Define some arbitrary distribution on the latent variables \(q(\mathbf{Z})\).

• Multiply by 1:

  \[\ell(\theta; \mathbf{X}) = \log \sum_{\mathbf{Z}} f(\mathbf{X} | \mathbf{Z}, \theta) p(\mathbf{Z} | \theta) \frac{q(\mathbf{Z})}{q(\mathbf{Z})}\]

• Rearranging terms, we can see that this can be written as an expectation over the distribution \(q\)

  \[\ell(\theta; \mathbf{X}) = \log E_{q} \bigg[\frac{ f(\mathbf{X} | \mathbf{Z}, \theta) p(\mathbf{Z} | \theta)}{q(\mathbf{Z})}\bigg]\]

• By Jensen’s inequality, we have a lower bound

  \[\ell(\theta; \mathbf{X}) \ge E_{q} \bigg[ \log \bigg(\frac{ f(\mathbf{X} | \mathbf{Z}, \theta) p(\mathbf{Z} | \theta)}{q(\mathbf{Z})}\bigg)\bigg]\]

Expectation-Maximization

• And by properties of logs, our lower-bound is:

  \[\ell(\theta; \mathbf{X}) \ge E_{q}[\log f(\mathbf{X}|\mathbf{Z}, \theta)] + E_{q}[\log p(\mathbf{Z}| \theta)] - E_{q}[\log q(\mathbf{Z})]\]

• We can optimize this iteratively by switching between finding an optimal distribution \(q^{(t+1)}\) given parameter values \(\theta^{(t)}\) and finding parameter values \(\theta^{(t+1)}\) given an existing choice of \(q^{(t+1)}\).

  1. Expectation step: Find \(q^{(t+1)}\) (what distribution are we taking the expectation over)
  2. Maximization step: Find \(\theta^{(t+1)}\) (given our \(q\) distribution, what is the value of the parameter values that maximizes the lower bound).

E-step

• Can we find a closed-form “optimal” solution for our update \(q^{(t+1)}\) given \(\theta^{(t)}\)?

  • Yes, find the \(q\) that makes the inequality an equality!

  \[\ell(\theta^{(t)}; \mathbf{X}) \ge E_{q} \bigg[ \log \bigg(\frac{ f(\mathbf{X} | \mathbf{Z}, \theta^{(t)}) p(\mathbf{Z} | \theta^{(t)})}{q(\mathbf{Z})}\bigg)\bigg]\]

  \[\ell(\theta^{(t)}; \mathbf{X}) \ge \sum_{\mathbf{Z}}\log \bigg(\frac{ f(\mathbf{X}, \mathbf{Z} | \theta^{(t)}) }{q(\mathbf{Z})}\bigg) q(\mathbf{Z})\]

E-step

• Let’s try to get the \(\mathbf{Z}\) out of the log. Start by factoring the joint distribution in the numerator conditioning on \(\mathbf{X}\) instead of \(\mathbf{Z}\)

  \[\ell(\theta^{(t)}; \mathbf{X}) \ge \sum_{\mathbf{Z}}\log \bigg(\frac{ f(\mathbf{Z} | \mathbf{X}, \theta^{(t)}) p(\mathbf{X} | \theta^{(t)})}{q(\mathbf{Z})}\bigg) q(\mathbf{Z})\]

• Now we can see the optimal choice for \(q^{(t+1)}\) revealed to us. Suppose we plug in \(q = f(\mathbf{Z} | \mathbf{X}, \theta^{(t)})\), the conditional distribution of \(\mathbf{Z}\) given \(\mathbf{X}\) and \(\theta^{(t)}\):

  \[\ell(\theta^{(t)}; \mathbf{X}) \ge \sum_{\mathbf{Z}}\log p(\mathbf{X} | \theta^{(t)}) \times f(\mathbf{Z} | \mathbf{X}, \theta^{(t)})\]

• The logged term no longer depends on \(\mathbf{Z}\), so pull it out of the sum

  \[\ell(\theta^{(t)}; \mathbf{X}) \ge \log p(\mathbf{X} | \theta^{(t)}) \times \sum_{\mathbf{Z}} f(\mathbf{Z} | \mathbf{X}, \theta^{(t)})\]

E-step

• The sum is equal to \(1\) (sum over a PMF/PDF) and we’re left with the definition of the marginal log-likelihood, so this expression holds with equality

  \[\ell(\theta^{(t)}; \mathbf{X}) = \log p(\mathbf{X} | \theta^{(t)})\]

• So our choice of distribution over which to take the expectation of the joint likelihood is \(q^{(t+1)} = f(\mathbf{Z} | \mathbf{X}, \theta^{(t)})\)

M-step

• Let’s go back to our lower bound - for a given value of \(q^{(t+1)}\), we want to find the \(\theta^{(t+1)}\) that maximize the “complete data” log-likelihood

  \[\ell(\theta; \mathbf{X}) \ge E_{q^{(t+1)}}[\log f(\mathbf{X}|\mathbf{Z}, \theta)] + E_{q^{(t+1)}}[\log p(\mathbf{Z}| \theta)] - E_{q^{(t+1)}}[\log q^{(t+1)}(\mathbf{Z})]\]

• Since the third term doesn’t depend on \(\theta\), we only need to worry about the first two.

  • This expectation is sometimes called the “Q-Function” \(Q(\theta | \theta^{(t)})\) or the expected “complete data” log-likelihood

  \[Q(\theta | \theta^{(t)}) = E_{\mathbf{Z}|\mathbf{X}, \theta^{(t)}}[\log f(\mathbf{X}, \mathbf{Z}| \theta)]\]

• The M-step sets \(\theta^{(t+1)}\) to the value of \(\theta\) that maximizes this Q-function

  \[\theta^{(t+1)} = \argmax_{\theta} Q(\theta | \theta^{(t)})\]

• Each E-step makes the bound touch the log-likelihood at \(\theta^{(t)}\), and each M-step pushes the bound up, so the log-likelihood never decreases from one iteration to the next

Deriving EM for the Bernoulli mixture

• For current values \(\beta^{(t)}\), \(\pi^{(t)}\), let’s derive the conditional distribution \(z_i | Y_i, \beta^{(t)}, \pi^{(t)}\) using Bayes’ rule

  \[p(z_i = k | Y_i, \beta^{(t)}, \pi^{(t)}) = \frac{p(z_i = k | \pi^{(t)}) f(Y_i | \beta^{(t)}, z_i = k)}{f(Y_i | \beta^{(t)}, \pi^{(t)})}\]

• Given \(z_i = k\), the responses are independent Bernoulli draws. And the denominator is just the sum of the numerator over \(k\)

  \[\gamma_{ik} = p(z_i = k | Y_i, \beta^{(t)}, \pi^{(t)}) = \frac{\overbrace{\pi_k^{(t)}}^{\text{prior}} \overbrace{\prod_{j=1}^J \big(\beta_{kj}^{(t)}\big)^{Y_{ij}}\big(1 - \beta_{kj}^{(t)}\big)^{1 - Y_{ij}}}^{\text{likelihood of } Y_i \text{ under cluster } k}}{\underbrace{\sum_{l=1}^K \pi_l^{(t)} \prod_{j=1}^J \big(\beta_{lj}^{(t)}\big)^{Y_{ij}}\big(1 - \beta_{lj}^{(t)}\big)^{1 - Y_{ij}}}_{\text{prior} \times \text{likelihood, summed over all clusters}}}\]

• These weights, \(\gamma_{ik}\) are sometimes called the “responsibility” parameters as they denote the extent to which each cluster is “responsible” for an observation.

Deriving EM for the Bernoulli mixture

• Recall the “complete” log-likelihood

  \[\ell(\beta, \pi ; \mathbf{Y}, \mathbf{Z}) = \sum_{i=1}^N \sum_{k=1}^K I(z_i = k) \bigg[\log \pi_k + \sum_{j=1}^J \Big(Y_{ij} \log \beta_{kj} + (1 - Y_{ij}) \log (1 - \beta_{kj})\Big)\bigg]\]

• Now, taking the expectation over \(z_i\), the only component that is not a constant is \(I(z_i = k)\) and \(E[I(z_i = k)] = p(z_i = k)\) (fundamental bridge).

  • And we got that (conditional) expectation in the previous slide: \(\gamma_{ik}\)

• So our Q function is

  \[Q(\beta, \pi | \beta^{(t)}, \pi^{(t)}) = \sum_{i=1}^N \sum_{k=1}^K \gamma_{ik} \bigg[\log \pi_k + \sum_{j=1}^J \Big(Y_{ij} \log \beta_{kj} + (1 - Y_{ij}) \log (1 - \beta_{kj})\Big)\bigg]\]

• \(\beta^{(t)}\) and \(\pi^{(t)}\) only enter through the responsibilities \(\gamma_{ik}\)

Deriving EM for the Bernoulli mixture

• Closed form solutions for the M-step: each \(\beta_{kj}\) is a weighted Bernoulli MLE, and \(\pi\) comes from maximizing \(\sum_k n_k \log \pi_k\) subject to \(\sum_k \pi_k = 1\)

  \[n_k = \sum_{i=1}^N \gamma_{ik}\] \[\pi_k = \frac{n_k}{N}\] \[\beta_{kj} = \frac{1}{n_k}\sum_{i=1}^N \gamma_{ik} Y_{ij}\]

• Compare to K-means: \(\beta_{kj}\) is the share supporting policy \(j\) in cluster \(k\), but every respondent counts toward every cluster in proportion to \(\gamma_{ik}\)

  • K-means is the special case where each \(\gamma_{ik}\) is forced to be \(0\) or \(1\)

MAP estimation with EM

• EM works just as well for finding the posterior mode – add the log prior to the Q-function

  \[Q(\beta, \pi | \beta^{(t)}, \pi^{(t)}) + \log p(\beta) + \log p(\pi)\]

• With conjugate priors, the M-step stays in closed form

  • \(\beta_{kj} \sim \text{Beta}(a, b)\) and \(\pi \sim \text{Dirichlet}(\alpha, \dotsc, \alpha)\)

  \[\beta_{kj} = \frac{\sum_{i=1}^N \gamma_{ik} Y_{ij} + a - 1}{n_k + a + b - 2} \qquad \pi_k = \frac{n_k + \alpha - 1}{N + K(\alpha - 1)}\]

• The prior acts like pseudo-observations added to each cluster

  • Keeps \(\beta_{kj}\) away from exactly \(0\) or \(1\) (where the log-likelihood blows up) when a cluster is small

• Below we use \(a = b = 2\) and a flat prior on \(\pi\) (\(\alpha = 1\))

  • With 24,453 respondents, this has little effect on the estimates, but it keeps EM from crashing when a cluster has no units.

Implementing EM

• For the complete data log-likelihood, entry \((i, k)\) is \(\log \pi_k + \sum_j Y_{ij}\log \beta_{kj} + (1 - Y_{ij})\log(1 - \beta_{kj})\)

bmm_complete_loglik <- function(Y, beta, pi){ # N x K matrix: complete-data log-lik of unit i if z_i = k
  Y %*% t(log(beta)) + (1 - Y) %*% t(log(1 - beta)) + matrix(log(pi), nrow(Y), length(pi), byrow = T)
}

• The observed-data log-likelihood sums \(\pi_k f(Y_i | z_i = k)\) over \(k\), then takes the log

bmm_loglik <- function(Y, beta, pi){
  sum(log(rowSums(exp(bmm_complete_loglik(Y, beta, pi)))))
}

Implementing EM

• The E-step: normalize each row of \(\pi_k f(Y_i | z_i = k)\) to sum to \(1\)

bmm_estep <- function(Y, beta, pi){
  joint <- exp(bmm_complete_loglik(Y, beta, pi))
  joint / rowSums(joint) # N x K matrix of responsibilities
}

• The M-step: weighted shares supporting each policy, plus the Beta(\(a\), \(b\)) pseudo-observations

bmm_mstep <- function(Y, gamma, a = 2, b = 2){
  n_k <- colSums(gamma)
  list(pi = n_k / nrow(Y),
       beta = (t(gamma) %*% Y + a - 1) / (n_k + a + b - 2)) # K x J matrix
}

Implementing EM

bmm_em <- function(Y, K, a = 2, b = 2, maxit = 5000, tol = 1e-6){
  ## Initialize with random responsibilities
  gamma <- matrix(runif(nrow(Y) * K), ncol = K)
  gamma <- gamma / rowSums(gamma)
  curr_post <- -Inf
  for (iter in 1:maxit){
    # M-step
    params <- bmm_mstep(Y, gamma, a, b)
    # E-step
    gamma <- bmm_estep(Y, params$beta, params$pi)
    # Check convergence: change in the per-observation log posterior (up to a constant)
    log_lik <- bmm_loglik(Y, params$beta, params$pi)
    new_post <- log_lik + (a - 1) * sum(log(params$beta)) + (b - 1) * sum(log(1 - params$beta))
    if (abs(new_post - curr_post) / nrow(Y) < tol) break
    curr_post <- new_post
  }
  if (iter == maxit) warning("Maxit reached")
  list(beta = params$beta, pi = params$pi, gamma = gamma, log_lik = log_lik, log_post = new_post, iter = iter)
}

Implementing EM

• EM only finds a local maximum, so run it from several random starts and keep the run with the highest log posterior

bmm_em_runs <- function(Y, K, runs = 10){
  fits <- map(1:runs, function(r) bmm_em(Y, K))
  fits[[which.max(map_dbl(fits, "log_post"))]]
}

• Cluster labels are arbitrary, so it’s helpful to put them in a consistent order. We’ll summarize each by the average share giving the liberal answer

conservative_items <- c("concealed_carry", "more_police", "border_patrols", "border_wall", "abortion_exceptions",
                        "abortion_ban", "more_fossil_fuels", "medicaid_work_req", "repeal_aca")
liberal_answer <- 1 - colnames(Y) %in% conservative_items # 1 = "Support" is the liberal answer
order_clusters <- function(fit){
  liberal_share <- colMeans(liberal_answer * t(fit$beta) + (1 - liberal_answer) * (1 - t(fit$beta)))
  o <- order(liberal_share, decreasing = T)
  list(beta = fit$beta[o, , drop = F], pi = fit$pi[o], gamma = fit$gamma[, o, drop = F], log_lik = fit$log_lik, log_post = fit$log_post)
}

Two clusters

set.seed(60637)
bmm2 <- order_clusters(bmm_em_runs(Y, K = 2))
bmm2$pi

[1] 0.666 0.334

Two clusters

• Cluster 2 (33%): conservative answers on almost everything. Cluster 1 (67%): liberal on most questions
  • Some questions split the clusters much less than others: 82% of cluster 2 support background checks; 44% of cluster 1 support more police
• Almost the same split as K-means

table(kmeans = km2$cluster, mixture = apply(bmm2$gamma, 1, which.max))

      mixture
kmeans     1     2
     1   698  8105
     2 15650     0

• It lines up mostly with Party ID – 33% of Republicans land in “liberal” cluster 1

round(prop.table(table(party = ces$pid3, cluster = apply(bmm2$gamma, 1, which.max)), 1), 2)

             cluster
party            1    2
  Democrat    0.95 0.05
  Independent 0.62 0.38
  Not sure    0.81 0.19
  Other       0.53 0.47
  Republican  0.33 0.67

Two clusters

• Each respondent has a probability of belonging to each cluster

• 86% of respondents are assigned to a cluster with more than 99% probability
  • Under local independence, 24 answers are 24 independent pieces of evidence about \(z_i\), so the posterior piles up near \(0\) and \(1\)

Five clusters

• Let’s add more clusters. With \(K = 5\), the log posterior (up to a constant) across 10 random starts:

set.seed(60637)
fits5 <- map(1:10, function(r) bmm_em(Y, K = 5))
map_dbl(fits5, "log_post")

 [1] -296879 -296879 -296879 -296879 -296978 -296879 -296879 -296879 -296879
[10] -296879

bmm5 <- order_clusters(fits5[[which.max(map_dbl(fits5, "log_post"))]])

• Different starts get stuck at different local maxima - so be sure to use multiple runs!
• Label switching: runs that reach the same maximum may number the clusters differently. Permuting the labels doesn’t change the likelihood.

Five clusters

Five clusters

• Cluster 1: Consistent liberals
• Cluster 5: Consistent conservatives
• Cluster 2: Liberal, except on police and the border
• Cluster 4: Probabilities between 30% and 70% on most questions
  • Seems like a cluster that has ideologically unconstrained preferences (Converse 1964)
• Cluster 3: A majority supports 22 of the 24 proposals – including both sides of the same question
  • 65% support more police and 46% support fewer; 93% support expanding Medicaid and 62% support repealing the ACA
  • Looks like acquiescence bias - a tendency to say “support” regardless of the content

Five clusters

• Who is in each cluster?

hard5 <- apply(bmm5$gamma, 1, which.max)
round(prop.table(table(cluster = hard5, ideology = ces$ideo5), 1), 2)[, c("Very liberal", "Liberal", "Moderate",
                                                                         "Conservative", "Very conservative", "Not sure")]

       ideology
cluster Very liberal Liberal Moderate Conservative Very conservative Not sure
      1         0.39    0.36     0.18         0.01              0.00     0.05
      2         0.10    0.29     0.46         0.06              0.01     0.07
      3         0.07    0.14     0.44         0.16              0.08     0.11
      4         0.03    0.09     0.40         0.23              0.09     0.16
      5         0.00    0.01     0.17         0.44              0.36     0.02

• The middle clusters (2-4) have the most self-described moderates; cluster 4 has the largest share of “Not sure”

Five clusters

• The CES has survey weights. We can estimate each cluster’s prevalence in the population by weighting the responsibilities with these weights (commonweight)

rbind(unweighted = colMeans(bmm5$gamma),
      weighted = colSums(ces$commonweight * bmm5$gamma) / sum(ces$commonweight))

            [,1]  [,2]  [,3]  [,4]  [,5]
unweighted 0.189 0.159 0.193 0.246 0.213
weighted   0.186 0.152 0.194 0.264 0.204

• One interesting feature is that the low-structure cluster 4 is under-represented in the raw sample

Challenges with mixture models

• Highly model dependent - the clusters we find depend on the distribution we assume for \(Y_i\) given \(z_i\)
  • Local independence is a substantial assumption: all of the correlation between answers has to be explained by the clusters
  • But some questions are logically dependent on others even conditional on ideology (e.g. abortion, policing).
• Multi-modality - Mixture model likelihoods often have multiple modes and EM is only guaranteed to converge to a local optimum
  • Can get stuck in a “bad” EM run (Solution: Run EM from multiple starting values and pick the one with the best log-likelihood (or log posterior))
  • “Label-switching” problem: Permuting the labels doesn’t change the likelihood
• Challenges with Bayes - Sampling-based inference can be tricky with mixtures due to the label-switching problem.
  • Different chains (or the same chain) can visit different permutations of the labels, so averaging \(\beta_{1j}\) across draws will mix clusters together
  • Common to use an approximation to the posterior around the mode obtained via EM.

Mixture models in Stan

• Stan can’t sample discrete parameters like \(z_i\) (no gradient)
• But we can marginalize out the \(z_i\). Same model and priors as our EM: Beta(2, 2) on each \(\beta_{kj}\), flat on \(\pi\)

bmm_stan_code <- "
data {
  int N; int J; int K;
  matrix[N, J] Y; // 1 = Support, 0 = Oppose
}
parameters {
  simplex[K] pi; // cluster prevalences (flat prior by default)
  matrix<lower=0, upper=1>[J, K] beta; // response probabilities
}
model {
  matrix[N, K] lp = Y * log(beta) + (1 - Y) * log1m(beta); // log f(Y_i | z_i = k), as in bmm_complete_loglik
  to_vector(beta) ~ beta(2, 2);
  for (i in 1:N) target += log_sum_exp(log(pi) + lp[i]'); // log of the sum over clusters
}
"

Mixture models in Stan

• Fit \(K = 5\) (this takes a few minutes – the likelihood sums over all 24,453 respondents at every step)

bmm_stan_fit <- stan(
  model_code = bmm_stan_code,  # Stan code
  data = list(N = nrow(Y), J = ncol(Y), K = 5, Y = Y), # named list of data
  chains = 4,             # number of Markov chains
  warmup = 500,           # number of warmup iterations per chain
  iter = 1000,            # total number of iterations per chain
  cores = 4,              # number of cores
  refresh = 0,            # no progress shown
  seed = 53706            # random seed
  )

• Label switching: each chain’s posterior mean of \(\pi\). The chains found the same clusters, but numbered them differently

round(get_posterior_mean(bmm_stan_fit, pars = "pi")[, 1:4], 2)

      mean-chain:1 mean-chain:2 mean-chain:3 mean-chain:4
pi[1]         0.28         0.17         0.17         0.15
pi[2]         0.18         0.22         0.18         0.17
pi[3]         0.17         0.18         0.22         0.18
pi[4]         0.15         0.15         0.28         0.22
pi[5]         0.23         0.28         0.15         0.28

Stan vs. EM

• Relabel every draw the same way we ordered the EM clusters (by the share giving the liberal answer)

relabel <- function(beta, pi){ # beta is J x K in Stan
  liberal_share <- colMeans(liberal_answer * beta + (1 - liberal_answer) * (1 - beta))
  o <- order(liberal_share, decreasing = T)
  list(beta = t(beta[, o]), pi = pi[o]) # K x J, like the EM estimates
}
relabel_draws <- function(fit){
  draws <- rstan::extract(fit)
  relabeled <- map(1:nrow(draws$pi), function(s) relabel(draws$beta[s, , ], draws$pi[s, ]))
  list(beta = simplify2array(map(relabeled, "beta")), # K x J x draws
       pi = simplify2array(map(relabeled, "pi")))     # K x draws
}
stan_random <- relabel_draws(bmm_stan_fit)

bmm5_worst <- order_clusters(fits5[[which.min(map_dbl(fits5, "log_post"))]]) # the EM start that found a worse local maximum
rbind(`EM (best start)` = bmm5$pi, `EM (worst start)` = bmm5_worst$pi, `Stan (posterior mean)` = rowMeans(stan_random$pi))

                       [,1]  [,2]  [,3]  [,4]  [,5]
EM (best start)       0.189 0.159 0.193 0.246 0.213
EM (worst start)      0.276 0.179 0.149 0.172 0.224
Stan (posterior mean) 0.276 0.177 0.149 0.173 0.225

Stan vs. EM

• All four chains align on the worse local maximum (log-likelihood -296,726 vs. -296,609)
• MCMC can’t jump between well-separated modes, so like EM, there’s no guarantee we’ve found the best mode.

Starting Stan at the EM mode

• Solution: Find the mode with EM, then run MCMC for uncertainty

em_init <- list(pi = bmm5$pi, beta = t(bmm5$beta)) # beta is J x K in Stan
bmm_stan_fit_em <- stan(
  model_code = bmm_stan_code, data = list(N = nrow(Y), J = ncol(Y), K = 5, Y = Y),
  chains = 4, warmup = 500, iter = 1000, cores = 4, refresh = 0, seed = 53706,
  init = rep(list(em_init), 4) # the same starting values for all four chains
  )
stan_em <- relabel_draws(bmm_stan_fit_em)

rbind(`EM (best start)` = bmm5$pi, `Stan (started at EM)` = rowMeans(stan_em$pi))

                      [,1]  [,2]  [,3]  [,4]  [,5]
EM (best start)      0.189 0.159 0.193 0.246 0.213
Stan (started at EM) 0.188 0.159 0.194 0.246 0.213

• No label switching (every chain starts with the same labels) and we’ve guaranteed a good mode.

Starting Stan at the EM mode

• Largest 95% interval width across the 120 \(\beta_{kj}\): 0.048

Choosing K

• Choosing \(K\) is a problem of model comparison and selection.
  • Just like before, more parameters = better in-sample fit.
  • But we want to avoid over-fitting
• Just as last week, two approaches
  • Information criteria - Evaluate the in-sample log-likelihood at the maximum, penalized by some factor for the number of parameters
    • Here, \(KJ + (K - 1)\) parameters
  • Cross-validation - Estimate the model on a training set. Compute the log-likelihood on a held-out test set.

Choosing K

set.seed(53703)
K_choice <- 1:10
fits_K <- map(K_choice, function(K) order_clusters(bmm_em_runs(Y, K, runs = 5)))
k_fit <- tibble(K = K_choice, log_lik = map_dbl(fits_K, "log_lik"), n_par = K * ncol(Y) + K - 1) %>%
  mutate(AIC = 2 * n_par - 2 * log_lik, BIC = log(nrow(Y)) * n_par - 2 * log_lik)

• Big gains up to \(K = 3\), then smaller ones – but even BIC keeps improving through \(K = 10\)

Choosing K

• 5-fold cross-validation: average held-out log-likelihood per respondent

set.seed(60637)
n_folds <- 5
cv_fold <- sample(rep(1:n_folds, length.out = nrow(Y)))
heldout <- expand_grid(K = K_choice, fold = 1:n_folds) %>%
  mutate(ll = map2_dbl(K, fold, function(K, f){
    fold_fit <- bmm_em_runs(Y[cv_fold != f, ], K, runs = 3)
    bmm_loglik(Y[cv_fold == f, ], fold_fit$beta, fold_fit$pi) / sum(cv_fold == f)
  }))

• With 24,453 respondents, more clusters keep predicting better out of sample. No obvious “right” \(K\)

Checking local independence

• Under the mixture model, the implied covariance between questions is \[\text{Cov}(Y_{ij}, Y_{il}) = \sum_{k=1}^K \pi_k \beta_{kj}\beta_{kl} - p_j p_l\]

• Like a posterior predictive check, but plugging in the EM estimates: does the model reproduce the observed correlations among the various pairs of questions?

Checking local independence

• With \(K = 2\), the model under-predicts many of the positive correlations: answers are still correlated within clusters
  • Worst of all: always allowing abortion and expanding abortion access (observed correlation 0.7, implied 0.29)
  • Local independence is an implausible assumption - lots of overconfidence in the classification probabilities.
• Adding clusters soaks up some of that correlation.

Logically exclusive questions

• A coherent respondent can’t support both 10% more police and 10% fewer police – these answers should be perfectly dependent

table(more_police = ces$more_police, fewer_police = ces$fewer_police)

           fewer_police
more_police     0     1
          0  4869  5572
          1 11970  2042

• Yet 8% of respondents support both. Within each \(K = 5\) cluster, the share supporting both:

• Local independence over-predicts “support both”, as it should – but interestingly the observed shares are nowhere near the \(0\%\) that perfect dependence implies.

Classifying the clusters

• Count the liberal answers each respondent gives (0-24) based on their “hard” cluster at \(K = 5\)

Latent classes vs. latent traits

• Mixture model (“latent class”): a discrete \(z_i\) for each respondent

  \[Pr(Y_{ij} = 1 | z_i = k) = \beta_{kj}\]

• Item response theory (“latent trait”): a continuous \(\theta_i\) for each respondent, and a slope and intercept for each question

  \[Pr(Y_{ij} = 1 | \theta_i) = F(\beta_j\theta_i - \tau_j)\]

• Both assume local independence conditional on the underlying latent variable

• But IRT uses a single continuous trait rather than a discrete classification.

• Maybe a single dimension works better here?

  • But it probably won’t do so well for cluster 4 (no structure) or cluster 3 (supporting both sides)

• Next week: fit IRT to the same 24 questions and see what we learn.

Next week

• Item response theory and ideal point models
  • Latent trait models for binary responses
  • “Ideal point” models: inferring latent continuous preferences from discrete responses
  • Adding temporal structure through a random-walk prior.
• Readings: Clinton, Jackman, and Rivers (2004) and Imai, Lo, and Olmsted (2016)