Estimating unique users for genre combinations with latent class analysis

Introduction

Suppose you are a data scientist at a fictional video-on-demand service.

A media planner asks, “How many unique users (UU) watched either K-Drama or Romance in the past 90 days?” UU is an important metric because it helps determine whether a proposed targeting segment is large enough to support the planned advertising campaign. Because questions like this come up frequently, we want to answer them quickly and accurately without running an expensive query each time.

If all viewing logs are available, we can calculate UU by counting each qualifying viewer once. Scanning a large user-level dataset for every genre combination, however, can be expensive. One way to avoid these repeated scans is to approximate the joint distribution of genre indicators in advance and calculate UU from that distribution.

How might we model this situation?

Latent Class Analysis

I use latent class analysis (LCA) for this approximation. LCA is commonly used to identify and interpret latent subpopulations, but here I use it to approximate the joint distribution of multivariate binary data. Once the model has been fitted, its parameters are all we need to calculate UU for different genre combinations.

As an example, consider a panel of \(n = 50{,}000\) viewers. For each viewer, we observe \(J = 16\) binary indicators showing whether they watched each genre; viewers may have watched more than one genre. The 16 genres yield \(2^{16} - 1\) nonempty combinations in principle. For this example, I evaluate approximation accuracy for combinations of up to four genres, which covers a practical range of queries.

For a genre set \(S\) and a known population size \(M\), UU is

\[ \operatorname{UU}(S) = M \cdot P\!\left( \bigcup_{j \in S}\{X_j = 1\} \right) \]

Below, I show how to calculate this union probability directly from the fitted LCA and then evaluate the approximation error.

Synthetic data generation

To produce realistic co-occurrence, I introduce correlations among genres with similar viewing patterns.

Each viewer \(i\) has four latent factors, \(\mathbf{f}_i \sim N(\mathbf{0},\mathbf{I}_4)\). Conditional on these factors, the genre indicators are generated independently.

\[ P(X_{ij} = 1 \mid \mathbf{f}_i) = \sigma\!\left( \alpha_j + \boldsymbol{\lambda}_j^\top\mathbf{f}_i \right), \qquad \sigma(x) = \frac{1}{1 + e^{-x}} \]

Factor 1 affects crime and documentary viewing, Factor 2 romance, and Factor 3 anime, science fiction, fantasy, and horror. Factor 4 represents overall viewing volume and loads weakly on every genre. Heavy viewers tend to watch a broader range of genres, creating positive correlations across the panel. The model nevertheless generates binary genre indicators, not viewing time itself.

Genre	Serious \(\lambda_{j1}\)	Romance \(\lambda_{j2}\)	Fiction \(\lambda_{j3}\)	Volume \(\lambda_{j4}\)	\(\alpha_j\)	Target viewer share \(P(X_j = 1)\)
Crime Thriller	1.2	0.0	0.3	0.5	-0.27	0.45
True Crime	1.4	0.0	0.0	0.5	-1.18	0.30
Documentary	1.0	0.0	0.0	0.5	-0.77	0.35
K-Drama	0.0	1.4	0.0	0.5	-0.87	0.35
Romance	0.0	1.3	0.2	0.5	-0.56	0.40
Reality Dating	0.0	1.0	0.0	0.5	-1.71	0.20
Anime	0.0	0.2	1.3	0.5	-1.15	0.30
Sci-Fi	0.3	0.0	1.2	0.5	-1.45	0.25
Fantasy	0.0	0.3	1.1	0.5	-1.41	0.25
Horror	0.2	0.0	0.8	0.5	-2.31	0.12
Stand-up Comedy	0.0	0.0	0.0	0.6	-1.18	0.25
Sitcom	0.0	0.3	0.0	0.6	-0.93	0.30
Action	0.2	0.0	0.5	0.6	0.46	0.60
Kids & Family	0.0	0.0	0.2	0.4	-1.44	0.20
Cooking & Food	0.0	0.2	0.0	0.5	-0.90	0.30
Sports Doc	0.6	0.0	0.0	0.5	-1.55	0.20

The data-generating process is a continuous factor model rather than a finite mixture. The Bernoulli mixture fitted below therefore does not exactly match the process that generated the data. This reflects a realistic setting: observed data are unlikely to follow any particular finite mixture exactly.

I calibrated each \(\alpha_j\) by simulation so that the share of viewers watching each genre ranges from 12% to 60%. The training and validation sets each contain 50,000 viewers. Final approximation accuracy is measured against an independent test set of 2,000,000 viewers generated by the same process.

genre_design <- tribble(
  ~genre            , ~serious , ~romance , ~fiction , ~volume , ~alpha ,
  "Crime Thriller"  , 1.2      , 0.0      , 0.3      , 0.5     , -0.27  ,
  "True Crime"      , 1.4      , 0.0      , 0.0      , 0.5     , -1.18  ,
  "Documentary"     , 1.0      , 0.0      , 0.0      , 0.5     , -0.77  ,
  "K-Drama"         , 0.0      , 1.4      , 0.0      , 0.5     , -0.87  ,
  "Romance"         , 0.0      , 1.3      , 0.2      , 0.5     , -0.56  ,
  "Reality Dating"  , 0.0      , 1.0      , 0.0      , 0.5     , -1.71  ,
  "Anime"           , 0.0      , 0.2      , 1.3      , 0.5     , -1.15  ,
  "Sci-Fi"          , 0.3      , 0.0      , 1.2      , 0.5     , -1.45  ,
  "Fantasy"         , 0.0      , 0.3      , 1.1      , 0.5     , -1.41  ,
  "Horror"          , 0.2      , 0.0      , 0.8      , 0.5     , -2.31  ,
  "Stand-up Comedy" , 0.0      , 0.0      , 0.0      , 0.6     , -1.18  ,
  "Sitcom"          , 0.0      , 0.3      , 0.0      , 0.6     , -0.93  ,
  "Action"          , 0.2      , 0.0      , 0.5      , 0.6     ,  0.46  ,
  "Kids & Family"   , 0.0      , 0.0      , 0.2      , 0.4     , -1.44  ,
  "Cooking & Food"  , 0.0      , 0.2      , 0.0      , 0.5     , -0.90  ,
  "Sports Doc"      , 0.6      , 0.0      , 0.0      , 0.5     , -1.55
)

genres <- genre_design$genre
s_pair <- c("K-Drama", "Romance")

gen_viewers <- function(n, design) {
  factor_names <- c("serious", "romance", "fiction", "volume")

  latent_factors <- matrix(
    rnorm(n * length(factor_names)),
    nrow = n,
    dimnames = list(NULL, factor_names)
  )

  factor_loadings <- design |>
    select(all_of(factor_names)) |>
    as.matrix()

  linear_predictor <- latent_factors %*%
    t(factor_loadings) |>
    sweep(
      MARGIN = 2,
      STATS = design$alpha,
      FUN = "+"
    )

  probabilities <- plogis(linear_predictor)

  matrix(
    rbinom(
      n = length(probabilities),
      size = 1,
      prob = probabilities
    ),
    nrow = n,
    dimnames = list(NULL, design$genre)
  )
}

empirical_union <- function(x, S) {
  mean(rowSums(x[, S, drop = FALSE]) > 0)
}

set.seed(2026)
X_train <- gen_viewers(50000, genre_design)

set.seed(2027)
X_valid <- gen_viewers(50000, genre_design)

The training data show clusters of correlated genres that are influenced by the same latent factors.

Figure 1: Correlations among genre indicators in the training data.

Overestimating the union under independence

The simplest way to estimate UU is to use the marginal viewer shares \(\hat p_j\) from the training data and assume that the genre indicators are mutually independent. Under this assumption, the share of users who watched at least one genre in \(S\) is

\[ \hat P_{\mathrm{indep}}\!\left( \bigcup_{j \in S}\{X_j = 1\} \right) = 1 - \prod_{j \in S}(1 - \hat p_j) \]

In the simulated data, however, genre preferences are not independent. Viewers who watch K-Drama are also more likely to watch Romance, so the overlap between the two audiences is larger than independence implies. The independence approximation consequently understates the overlap and overstates the share who watched at least one of the two genres.

Table 1 reports the marginal viewer shares from the training data, the union implied by independence, and the union observed in the validation data.

Table 1: Independence estimate and observed union for K-Drama and Romance.

K-Drama share	Romance share	Independence estimate	Observed union	Relative overestimate
34.7%	40.1%	60.9%	54.8%	11.1%

In a population of 20 million, an error of just a few percentage points in the union probability can lead to an overestimate of roughly a million users.

Approximating the joint distribution with LCA

The independence approximation makes union probabilities easy to calculate, but it cannot capture co-occurrence across genres. An LCA with \(K\) classes—equivalently, a Bernoulli mixture model—keeps the calculation simple while allowing genres to be dependent at the population level.

\[ P(\mathbf{x}) = \sum_{k=1}^K\pi_k\prod_{j=1}^J \theta_{jk}^{x_j}(1 - \theta_{jk})^{1-x_j}, \qquad \mathbf{x} \in \{0,1\}^J \]

Here, \(\pi_k\) is the mixing proportion for class \(k\), and \(\theta_{jk}\) is the probability that a viewer in class \(k\) watches genre \(j\). Genre indicators are conditionally independent within each class, but mixing multiple classes induces dependence in the population. I estimate the model with the EM algorithm implemented in {flexmix}.

An LCA with \(K\) classes has \(K(J + 1) - 1\) parameters, far fewer than the \(2^{16} - 1 = 65{,}535\) required by the saturated model. When \(K = 1\), the model has the same 16 parameters as the independence model; larger values of \(K\) use more parameters to capture richer dependence structures.

I do not try to interpret the latent classes themselves. Instead, I use the mixture to approximate \(P(\mathbf{x})\) and judge it by the accuracy of its UU estimates across genre combinations.

Computing UU for arbitrary conditions

Conditional on class \(k\), the genre indicators are independent. Let \(A\) be the set of genres required to be present and \(B\) the set required to be absent, with no genre belonging to both sets. Then

\[ P\!\left( \bigcap_{j \in A}\{X_j = 1\}\; \cap \;\bigcap_{j \in B}\{X_j = 0\} \right) = \sum_{k=1}^K\pi_k \prod_{j \in A}\theta_{jk}\prod_{j \in B}(1 - \theta_{jk}) \]

The probability of watching at least one genre in \(S\) is the complement of watching none of them:

\[ P\!\left( \bigcup_{j \in S}\{X_j = 1\} \right) = 1 - \sum_{k=1}^K\pi_k\prod_{j \in S}(1 - \theta_{jk}) \]

The same model can also handle intersections and exclusions, such as “watched Anime but not Horror,” by using the corresponding products of \(\theta_{jk}\) and \(1 - \theta_{jk}\).

Selecting the number of classes using validation data

I compare \(K \in \{1,2,3,4,6,8,12,16\}\). Each model is fitted using only the 50,000-viewer training set. The evaluation covers all 120 two-way, 560 three-way, and 1,820 four-way genre unions. For each query, I compare the estimated union probability with its empirical counterpart in the validation data. I exclude one-way queries because their marginal probabilities can be estimated directly from the training data.

library(flexmix)

Ks <- c(2, 3, 4, 6, 8, 12, 16)

set.seed(2029)
fits <- stepFlexmix(
  X_train ~ 1,
  k = Ks,
  nrep = 3,
  model = FLXMCmvbinary(),
  verbose = FALSE
)

models <- map(Ks, \(k) getModel(fits, which = as.character(k)))

params <- map(
  models,
  \(model) {
    list(
      pi = prior(model),
      theta = parameters(model)
    )
  }
)

params <- c(
  list(list(
    pi = 1,
    theta = matrix(
      colMeans(X_train),
      ncol = 1,
      dimnames = list(genres, NULL)
    )
  )),
  params
)

K_grid <- c(1, Ks)

queries <- map(2:4, \(n_way) {
  combn(seq_along(genres), n_way, simplify = FALSE)
}) |>
  list_flatten()

query_size <- lengths(queries)

uu_union <- function(pi_hat, theta_hat, S) {
  1 - sum(pi_hat * apply(1 - theta_hat[S, , drop = FALSE], 2, prod))
}

valid_observed <- map_dbl(queries, \(S) empirical_union(X_valid, S))

valid_error <- map2(K_grid, params, \(k, par) {
  estimate <- map_dbl(
    queries,
    \(S) uu_union(par$pi, par$theta, S)
  )

  tibble(
    query_id = seq_along(queries),
    k = k,
    n_way = query_size,
    observed = valid_observed,
    estimate = estimate,
    error = estimate - observed
  )
}) |>
  list_rbind()

valid_by_size <- valid_error |>
  summarise(rmse = sqrt(mean(error^2)), .by = c(k, n_way))

valid_overall <- valid_error |>
  summarise(rmse = sqrt(mean(error^2)), .by = k)

selected_k <- 8

Let \(u_q\) denote the union probability observed in the validation data and \(\hat u_{q,K}\) the estimate from the model with \(K\) classes. I select \(K\) using RMSE on the probability scale across all two-way, three-way, and four-way queries. Every combination receives equal weight, so query sizes with more combinations contribute more to the overall RMSE. With \(Q\) total queries, the criterion is

\[ \operatorname{RMSE}(K) = \sqrt{\frac{1}{Q}\sum_{q=1}^{Q}(\hat u_{q,K} - u_q)^2} \]

Figure 2: Validation RMSE by number of latent classes.

RMSE declines through \(K = 8\), with only small improvements beyond that point. I therefore choose \(K = 8\) as a balance between approximation accuracy and model complexity.

Approximation accuracy on the test set

After fixing \(K\), I evaluate approximation accuracy on an independently generated test set of two million viewers.

Even with synthetic data, a union probability calculated from a finite sample varies slightly from one draw to another. Using two million test viewers reduces this variation and brings the observed union probabilities close to their population values. Most of the remaining error can therefore be attributed to the model approximation rather than random variation in the test data. The test set is not used to fit the model or select \(K\).

rm(X_valid)

set.seed(2028)
X_test <- gen_viewers(2000000, genre_design)

test_observed <- map_dbl(queries, \(S) empirical_union(X_test, S))
selected_params <- params[[match(selected_k, K_grid)]]

test_estimate <- map_dbl(
  queries,
  \(S) {
    uu_union(
      selected_params$pi,
      selected_params$theta,
      S
    )
  }
)

test_error <- tibble(
  n_way = query_size,
  observed = test_observed,
  estimate = test_estimate,
  error = estimate - observed,
  rel_error = abs(error) / observed
)

test_summary <- test_error |>
  summarise(
    rmse = sqrt(mean(error^2)),
    mean_relative_error = mean(rel_error),
    p95_relative_error = quantile(rel_error, 0.95),
    max_relative_error = max(rel_error),
    .by = n_way
  )

The results below compare estimates from the selected LCA with the union probabilities observed in the test data. The estimates lie close to the 45-degree line, indicating small approximation errors.

Figure 3: Observed and estimated union probabilities in the test data.

Table 2: Test error by number of genres in the query.

Query size	RMSE	Mean relative error	95th percentile relative error	Maximum relative error
2-way	0.19%	0.33%	0.80%	1.08%
3-way	0.19%	0.25%	0.58%	1.01%
4-way	0.19%	0.21%	0.49%	0.96%

Summary

When targeting conditions overlap, simply adding their individual audience sizes will double-count some users. Modeling the joint distribution of genre indicators accounts for this overlap without repeatedly scanning user-level data.

In this example, I used latent class analysis not to discover and interpret audience segments, but to represent a complex joint distribution with a relatively small number of parameters. I therefore selected the number of classes using validation error on the UU queries of interest rather than an information criterion.

Once fitted, the model can answer arbitrary genre-union queries using only the estimated parameters. Its practical value lies in providing a reusable way to estimate the sizes of many overlapping audiences.

Implementation

Once the model is fitted, the only quantities required for aggregation are the class proportions \(\hat\pi_k\) and the class-specific genre probabilities \(\hat\theta_{jk}\).

library(tidyverse)
library(flexmix)

set.seed(2029)
sf <- stepFlexmix(
  X_train ~ 1,
  k = selected_k,
  nrep = 3,
  model = FLXMCmvbinary(),
  verbose = FALSE
)
fit <- getModel(sf, which = as.character(selected_k))

pi_hat <- prior(fit)
theta_hat <- parameters(fit)

uu_union <- function(pi_hat, theta_hat, S) {
  1 - sum(pi_hat * apply(1 - theta_hat[S, , drop = FALSE], 2, prod))
}

s_pair <- c("K-Drama", "Romance")
S <- match(s_pair, colnames(X_train))
uu_union(pi_hat, theta_hat, S)