COMIX: Coarsened mixtures of hierarchical skew normal kernels

S. Gorsky, C. Chan and L. Ma

knitr::opts_chunk$set(
  comment = '', fig.width = 6, fig.height = 6
)

library(COMIX)
# Number of observations for each sample (row) and cluster (column):
njk <- 
  matrix(
    c(
      150, 300,
      250, 200
    ),
    nrow = 2,
    byrow = TRUE
  )

# Dimension of data:
p <- 3

# Scale and skew parameters for first cluster:
Sigma1 <- matrix(0.5, nrow = p, ncol = p) + diag(0.5, nrow = p)
alpha1 <- rep(0, p)
alpha1[1] <- -5
# location parameter for first cluster in first sample:
xi11 <- rep(0, p)
# location parameter for first cluster in second sample (aligned with first):
xi21 <- rep(0, p)


# Scale and skew parameters for second cluster:
Sigma2 <- matrix(-1/3, nrow = p, ncol = p) + diag(1 + 1/3, nrow = p)
alpha2 <- rep(0, p)
alpha2[2] <- 5
# location parameter for second cluster in first sample:
xi12 <- rep(3, p)
# location parameter for second cluster in second sample (misaligned with first):
xi22 <- rep(4, p)

# Sample data:
set.seed(1)
Y <- 
  rbind(
    sn::rmsn(njk[1, 1], xi = xi11, Omega = Sigma1, alpha = alpha1),
    sn::rmsn(njk[1, 2], xi = xi12, Omega = Sigma2, alpha = alpha2),
    sn::rmsn(njk[2, 1], xi = xi21, Omega = Sigma1, alpha = alpha1),
    sn::rmsn(njk[2, 2], xi = xi22, Omega = Sigma2, alpha = alpha2)
  )

C <- c(rep(1, rowSums(njk)[1]), rep(2, rowSums(njk)[2]))

prior <- list(zeta = 1, K = 10)
pmc <- list(naprt = 5, nburn = 200, nsave = 200)
# Fit the model:
res <- comix(Y, C, pmc = pmc, prior = prior)

initializing all particles...
Done
Merged clusters (iteration 17)
Merged clusters (iteration 48)
Merged clusters (iteration 76)
Merged clusters (iteration 77)
Merged clusters (iteration 94)
Iteration: 100 of 400
Iteration: 200 of 400
Iteration: 300 of 400
Iteration: 400 of 400

# Relabel to resolve potential label switching issues:
res_relab <- relabelChain(res)

K=10
T=200
N=900
J=2
p=3

# Generate calibrated data:
cal <- calibrateNoDist(res_relab)

# Compare raw and calibrated data:
oldpar <- par(mfrow = c(1, 2)) 
plot(Y, col = C, xlim = range(Y[,1]), ylim = range(Y[,2]), main = "Raw Data",
     xlab = expression(Y[1]), ylab = expression(Y[2]))
plot(cal$Y_cal, col = C, main = "Calibrated Data",
     xlab = expression(Y[1]), ylab = expression(Y[2]))

par(oldpar)

# Get posterior estimates for the model parameters:
res_summary <- summarizeChain(res_relab)
# Check for instance, the cluster assignment labels:
table(res_summary$t)

  2   3 
400 500 

# Indeed the same as 
colSums(njk)

[1] 400 500

# Or examine the skewness parameter for the non-trivial clusters:
res_summary$alpha[ , unique(res_summary$t)]

           [,1]        [,2]
[1,] -4.6033765  0.09181518
[2,] -0.5486982  4.92952169
[3,]  0.3524381 -0.24937654

# And compare those to
cbind(alpha1, alpha2)

     alpha1 alpha2
[1,]     -5      0
[2,]      0      5
[3,]      0      0