Modeling heterogeneity

Weight-based update of the latent classes

The following weight-based updating scheme is analogue to Bauer, Büscher, and Batram (2019) and executed within the burn-in period:

• We remove class c, if s_c < ε_min, i.e. if the class weight s_c drops below some threshold ε_min. This case indicates that class c has a negligible impact on the mixing distribution.

• We split class c into two classes c₁ and c₂, if s_c > ε_max. This case indicates that class c has a high influence on the mixing distribution whose approximation can potentially be improved by increasing the resolution in directions of high variance. Therefore, the class means b_c1 and b_c2 of the new classes c₁ and c₂ are shifted in opposite directions from the class mean b_c of the old class c in the direction of the highest variance.

• We join two classes c₁ and c₂ to one class c, if ∥b_c1 − b_c2∥ < ε_distmin, i.e. if the euclidean distance between the class means b_c1 and b_c2 drops below some threshold ε_distmin. This case indicates location redundancy which should be repealed. The parameters of c are assigned by adding the values of s from c₁ and c₂ and averaging the values for b and Ω.

These rules contain choices on the values for ε_min, ε_max and ε_distmin. The adequate value for ε_distmin depends on the scale of the parameters. Per default, {RprobitB} sets

• epsmin = 0.01,

• epsmax = 0.99, and

• distmin = 0.1.

These values can be adapted through the latent_class list.

Dirichlet process-based update of the latent classes

As an alternative to the weight-based updating scheme to determine the correct number C of latent classes, {RprobitB} implements the Dirichlet process.⁴ The Dirichlet Process is a Bayesian nonparametric method, where nonparametric means that the number of model parameters can theoretically grow to infinity. The method allows to add more mixture components to the mixing distribution if needed for a better approximation, see Neal (2000) for a documentation of the general case. The literature offers many representations of the method, including the Chinese Restaurant Process (Aldous 1985), the stick-braking metaphor (Sethuraman 1994), and the Polya Urn model (Blackwell and MacQueen 1973).

In our case, we face the situation to find a distribution g that explains the decider-specific coefficients (β_n)_{n = 1, …, N}, where g is supposed to be a mixture of an unknown number C of Gaussian densities, i.e. g = ∑_{c = 1, …, C}s_cMVN(b_c, Ω_c).

Let z_n ∈ {1, …, C} denote the class membership of β_n. A priori, the mixture weights (s_c)_c are given a Dirichlet prior with concentration parameter δ/C, i.e. (s_c)_c ∣ δ ∼ D_C(δ/C, …, δ/C). Rasmussen (1999) shows that

$$ \Pr((z_n)_n\mid \delta) = \frac{\Gamma(\delta)}{\Gamma(N+\delta)} \prod_{c=1}^C \frac{\Gamma(m_c + \delta/C)}{\Gamma(\delta/C)}, $$ where Γ(⋅) denotes the gamma function and m_c = #{n : z_n = c} the number of elements that are currently allocated to class c. Crucially, the last equation is independent of the class weights (s_c)_c, yet it still depends on the finite number C of latent classes. However, Li, Schofield, and Gönen (2019) shows that

$$ \Pr(z_n = c \mid z_{-n}, \delta) = \frac{m_{c,-n} + \delta/C}{N-1+\delta},$$ where the notation −n means excluding the nth element. We can let C approach infinity to derive:

$$ \Pr(z_n = c \mid z_{-n}, \delta) \to \frac{m_{c,-n}}{N-1+\delta}. $$

Note that the allocation probabilities do not sum to 1, instead

$$ \sum_{c = 1}^C \frac{m_{c,-n}}{N-1+\delta} = \frac{N-1}{N-1+\delta}. $$

The difference to 1 equals

$$ \Pr(z_n \neq z_m ~ \forall ~ m \neq n \mid z_{-n}, \delta) = \frac{\delta}{N-1+\delta} $$

and constitutes the probability that a new cluster for observation n is created. Neal (2000) points out that this probability is proportional to the prior parameter δ: A greater value for δ encourages the creation of new clusters, a smaller value for δ increases the probability of an allocation to an already existing class.

In summary, the Dirichlet process updates the allocation of each β coefficient vector one at a time, dependent on the other allocations. The number of clusters can theoretically rise to infinity, however, as we delete unoccupied clusters, C is bounded by N. As a final step after the allocation update, we update the class means b_c and covariance matrices Ω_c by means of their posterior predictive distribution. The mean and covariance matrix for a new generated cluster is drawn from the prior predictive distribution. The corresponding formulas are given in Li, Schofield, and Gönen (2019).

The Dirichlet process directly integrates into our existing Gibbs sampler. Given β values, it updated the class means b_c and class covariance matrices Ω_c. The Dirichlet process updating scheme is implemented in the function update_classes_dp(). In the following, we give a small example in the bivariate case P_r = 2. We sample true class means b_true and class covariance matrices Omega_true for C_true = 3 true latent classes. The true (unbalanced) class sizes are given by the vector N, and z_true denotes the true allocations.

set.seed(1)
P_r <- 2
C_true <- 3
N <- c(100, 70, 30)
(b_true <- matrix(replicate(C_true, rnorm(P_r)), nrow = P_r, ncol = C_true))
#>            [,1]       [,2]       [,3]
#> [1,] -0.6264538 -0.8356286  0.3295078
#> [2,]  0.1836433  1.5952808 -0.8204684
(Omega_true <- matrix(replicate(C_true, rwishart(P_r + 1, 0.1 * diag(P_r))$W, simplify = TRUE),
  nrow = P_r * P_r, ncol = C_true
))
#>           [,1]        [,2]       [,3]
#> [1,] 0.3093652  0.14358543  0.2734617
#> [2,] 0.1012729 -0.07444148 -0.1474941
#> [3,] 0.1012729 -0.07444148 -0.1474941
#> [4,] 0.2648235  0.05751780  0.2184029
beta <- c()
for (c in 1:C_true) {
  for (n in 1:N[c]) {
    beta <- cbind(beta, rmvnorm(mu = b_true[, c, drop = F], Sigma = matrix(Omega_true[, c, drop = F], ncol = P_r)))
  }
}
z_true <- rep(1:3, times = N)

We specify the following prior parameters (for their definition see the vignette on model fitting):

delta <- 0.1
xi <- numeric(P_r)
D <- diag(P_r)
nu <- P_r + 2
Theta <- diag(P_r)

Initially, we start with C = 1 latent classes. The class mean b is set to zero, the covariance matrix Omega to the identity matrix:

z <- rep(1, ncol(beta))
C <- 1
b <- matrix(0, nrow = P_r, ncol = C)
Omega <- matrix(rep(diag(P_r), C), nrow = P_r * P_r, ncol = C)

The following call to update_classes_dp() updates the latent classes in 100 iterations. Note that we specify the arguments Cmax and s_desc. The former denotes the maximum number of latent classes. This specification is not a requirement for the Dirichlet process per se, but rather for its implementation. Knowing the maximum possible class number, we can allocate the required memory space, which leads to a speed improvement. We later can verify that we won’t exceed the number of Cmax = 10 latent classes at any point of the Dirichlet process. Setting s_desc = TRUE ensures that the classes are ordered by their weights in a descending order to ensure identifiability.

for (r in 1:100) {
  dp <- RprobitB:::update_classes_dp(
    Cmax = 10, beta, z, b, Omega, delta, xi, D, nu, Theta, s_desc = TRUE
  )
  z <- dp$z
  b <- dp$b
  Omega <- dp$Omega
}

The Dirichlet process was able to infer the true number C_true = 3 of latent classes:

par(mfrow = c(1, 2))
plot(t(beta), xlab = bquote(beta[1]), ylab = bquote(beta[2]), pch = 19)
RprobitB:::plot_class_allocation(beta, z, b, Omega, r = 100, perc = 0.95)