Model mieszanin wzrostu (GMM) i LCGA – jak wyodrębnić klasy trajektorii?

suppressPackageStartupMessages(library(mclust))
set.seed(3561)
N <- 300; t <- 0:4                                  # 300 osób, 5 pomiarów
klasa <- sample(1:3, N, replace = TRUE, prob = c(.60, .25, .15))
b0 <- c(10, 12, 25)[klasa] + rnorm(N, 0, 1.5)      # zmienność wewnątrz klas (jak w GMM)
b1 <- c(0, 3, -3)[klasa]   + rnorm(N, 0, 0.5)
Y  <- outer(b0, rep(1, 5)) + outer(b1, t) + matrix(rnorm(N * 5, 0, 2), N)
# LCGA: mieszanina trajektorii liniowych, algorytm EM, bez efektów losowych
lcga <- function(Y, K, starty = 20) {
  X <- cbind(1, t); naj <- NULL
  for (s in 1:starty) {
    w <- matrix(runif(nrow(Y) * K), ncol = K); w <- w / rowSums(w)
    ll_old <- -Inf
    for (it in 1:500) {
      B <- sapply(1:K, function(k) coef(lm.wfit(X[rep(1:5, nrow(Y)), ], c(t(Y)), rep(w[, k], each = 5))))
      s2 <- sapply(1:K, function(k) sum(w[, k] * rowSums((Y - outer(rep(1, nrow(Y)), c(X %*% B[, k])))^2)) / (5 * sum(w[, k])))
      p  <- colMeans(w)
      lf <- sapply(1:K, function(k) log(p[k]) + rowSums(dnorm(Y, outer(rep(1, nrow(Y)), c(X %*% B[, k])), sqrt(s2[k]), log = TRUE)))
      lf <- matrix(lf, ncol = K); m <- apply(lf, 1, max)
      ll <- sum(m + log(rowSums(exp(lf - m)))); w <- exp(lf - m); w <- w / rowSums(w)
      if (abs(ll - ll_old) < 1e-8) break; ll_old <- ll
    }
    if (is.null(naj) || ll > naj$ll) naj <- list(ll = ll, B = B, p = p, w = w)
  }
  npar <- 4 * K - 1                                 # 2 współczynniki + wariancja na klasę, K - 1 proporcji
  E <- if (K == 1) NA else 1 + sum(naj$w * log(pmax(naj$w, 1e-300))) / (nrow(Y) * log(K))
  c(naj, BIC = -2 * naj$ll + npar * log(nrow(Y)), entropia = E)
}
wyniki <- lapply(1:5, function(K) lcga(Y, K))
round(sapply(wyniki, function(r) c(BIC = r$BIC, entropia = r$entropia)), 2)
r3 <- wyniki[[3]]
round(rbind(B = r3$B, p = r3$p), 2)                # trajektorie i udziały klas (K = 3)
table(prawdziwa = klasa, LCGA = max.col(r3$w))
# podejście dwuetapowe: wyraz wolny i nachylenie każdej osoby, potem mclust
wsp <- t(apply(Y, 1, function(y) coef(lm(y ~ t))))
mc  <- Mclust(wsp, G = 1:5, verbose = FALSE)
c(G = mc$G, BIC_mclust = round(mc$bic, 1))
table(prawdziwa = klasa, mclust = mc$classification)
# GMM z efektami losowymi wewnątrz klas – wymaga pakietu lcmm, niewykonany:
# m1 <- lcmm::hlme(y ~ t, random = ~ t, subject = "id", ng = 1, data = dane_dlugie)
# m3 <- lcmm::hlme(y ~ t, mixture = ~ t, random = ~ t, subject = "id", ng = 3,
#                  data = dane_dlugie, B = m1)

            [,1]    [,2]    [,3]    [,4]   [,5]
BIC      9404.74 8369.04 7911.38 7743.00 7691.9
entropia      NA    0.96    0.99    0.91    0.9
  [,1]  [,2]  [,3]
  9.92 25.15 11.94
t 0.00 -2.97  2.92
p 0.62  0.13  0.24
         LCGA
prawdziwa   1   2   3
        1 186   0   1
        2   0   0  73
        3   0  40   0
         G BIC_mclust 
       3.0    -2563.1 
         mclust
prawdziwa   1   2   3
        1 187   0   0
        2   1   0  72
        3   0  40   0