In SWotherspoon/MVNSeq: Normal Mixture and Hidden Markov Models

Some simple tests of the MVNSeq package.

set.seed(31)
library(MVNSeq)

Iris

Fit a simple mixture model to the iris data set and compare to the results from mclust.

y <- as.matrix(log(iris[,1:4]))

Use kmeans to get initial class memberships, then fit a mixture model

library(MVNSeq)
km <- kmeans(y,3)
fit <- mvnMix(y,km$cluster)

Fit equivalent mclust model

library(mclust)
fit.mclust <- Mclust(y,G=3,modelNames=c("VVV"))
summary(fit.mclust,parameters=T)

Compare classification

table(max.col(fit$pars$P),fit.mclust$classification)

Match labels and compare means and variances

j <- order(max.col(table(max.col(fit$pars$P),fit.mclust$classification)))
sapply(fit$pars$mvn[j],function(p) p$mu)
lapply(fit$pars$mvn[j],function(p) p$Sigma)

Simulation

Simulate m sequences of q dimensional responses from a K state hidden Markov or mixture model

## Number of components K, responses q and groups m
library(mvtnorm)
sim <- function(K,q,m,uscale=0,vscale=1,minl=100,meanl=500,hmm=T) {

  ## Vector of sequence lengths
  ns <- minl+rpois(m,max(meanl-minl,0))
  n <- sum(ns)

  ## Group structure
  gr <- rep(1:m,ns)
  J <- c(1,cumsum(ns))

  ## Responses, group means, population means, and covariances
  Y <- array(0,c(n,q,K))
  A <- array(0,c(m,q,K))
  Mu <- array(0,c(K,q))
  U <- vector(mode="list",K)
  V <- vector(mode="list",K)

  ## Simulate a response for each state and group
  for(k in 1:K) {
    ## Covariance of random effects
    U[[k]] <- crossprod(matrix(rnorm(q*q,0,uscale),q,q))
    ## Covariance of response about group means
    V[[k]] <- crossprod(matrix(rnorm(q*q,0,vscale),q,q))
    ## Population mean
    Mu[k,] <- runif(q,0,10)
    ## Group means
    A[,,k] <- rmvnorm(m,Mu[k,],U[[k]])
    ## Responses
    Y[,,k] <- A[gr,,k] + rmvnorm(n,rep(0,q),V[[k]])
  }

  ## state prior
  p <- runif(K)
  p <- p/sum(p)


  if(hmm) {
    ## Random transition matrix
    Q <- matrix(runif(K*K),K,K)
    Q <- Q/rowSums(Q)
  } else {
    Q <- matrix(p,K,K,byrow=T)
  }

  ## Simulate (Markovian) sequence of class membership and observations
  cl <- integer(n)
  y <- array(0,c(n,q))
  for(i in 1:m) {
    j <- J[i]
    cl[j] <- sample(K,1,prob=p)
    y[j,] <- Y[j,,cl[j]]
    for(j in (J[i]+1):J[i+1]) {
      cl[j] <- sample(K,1,prob=Q[cl[j-1],])
      y[j,] <- Y[j,,cl[j]]
    }
  }
  list(y=y,cl=cl,gr=gr,A=A,Mu=Mu,U=U,V=V,p=p,Q=Q)
}

Single sequence from a 3 component Mixture

Simulate a 3 component mixture.

s <- sim(K=3,q=3,m=1,hmm=F)

Use kmeans to get initial class memberships, then fit a mixture model

km <- kmeans(s$y,3)
fit <- mvnMix(s$y,km$cluster)

Fit equivalent mclust model

library(mclust)
fit.mclust <- Mclust(s$y,G=3,modelNames=c("VVV"))
summary(fit.mclust,parameters=T)

Compare to mclust classification

table(max.col(fit$pars$P),fit.mclust$classification)

Match labels and compare means and variances

j <- order(max.col(table(max.col(fit$pars$P),fit.mclust$classification)))
sapply(fit$pars$mvn[j],function(p) p$mu)
lapply(fit$pars$mvn[j],function(p) p$Sigma)

Compare to actual classification

table(max.col(fit$pars$P),s$cl)

Match labels and compare means

j <- order(max.col(table(max.col(fit$pars$P),s$cl)))
t(s$Mu)
sapply(fit$pars$mvn[j],function(p) p$mu)

and variances

s$V
lapply(fit$pars$mvn[j],function(p) p$Sigma)

Single HHM sequence

Simulate one long sequence

s <- sim(K=3,q=3,m=1,minl=10000,hmm=T)
pairs(s$y,pch=".",col=s$cl)

Mixture model

Use kmeans to get initial class memberships, then fit a mixture model

km <- kmeans(s$y,3)
fit <- mvnMix(s$y,km$cluster)

Compare to actual classification

table(max.col(fit$pars$P),s$cl)

Match labels and compare means

j <- order(max.col(table(max.col(fit$pars$P),s$cl)))
t(s$Mu)
sapply(fit$pars$mvn[j],function(p) p$mu)

and variances

s$V
lapply(fit$pars$mvn[j],function(p) p$Sigma)

HMM

Use kmeans to get initial class memberships, then fit an hidden Markov model

km <- kmeans(s$y,3)
fit <- mvnHMM(s$y,km$cluster)

Compare to actual classification

table(max.col(fit$pars$P),s$cl)

Match labels and compare means

j <- order(max.col(table(max.col(fit$pars$P),s$cl)))
t(s$Mu)
sapply(fit$pars$mvn[j],function(p) p$mu)

and variances

s$V
lapply(fit$pars$mvn[j],function(p) p$Sigma)

and transition probabilities

fit$Q[j,j]
s$Q

Multiple HHM Sequence

Simulate many shorter sequences all with varying means

s <- sim(K=3,q=3,m=100,minl=100,meanl=500,hmm=T,uscale = 0.8)
pairs(s$y,pch=".",col=s$cl)

Fixed Effects Mixture Model

Use kmeans to get initial class memberships, but use actual means to avoid label permutation issues, then fit a mixture model to each sequence forcing common mixing fractions across sequences

km <- kmeans(s$y,s$Mu)
fit <- gmvnMix(s$y,km$cluster,s$gr,common.fractions = TRUE)

Compare classifications across all groups

table(unlist(lapply(fit$pars,function(p) max.col(p$P))),s$cl)

Compare group means

## Check parameter estimates match simulation
opar <- par(mfrow=c(2,2),mar=c(1,1,1,1)+0.1)
matplot(s$A[,,1],t(sapply(fit$pars,function(p) p$mvn[[1]]$mu)),pch=16,cex=0.5,xlab="",ylab="")
matplot(s$A[,,2],t(sapply(fit$pars,function(p) p$mvn[[2]]$mu)),pch=16,cex=0.5,xlab="",ylab="")
matplot(s$A[,,3],t(sapply(fit$pars,function(p) p$mvn[[3]]$mu)),pch=16,cex=0.5,xlab="",ylab="")
par(opar)

Compare covariances

fit$pars[[1]]$mvn[[1]]$Sigma
s$V[[1]]
fit$pars[[1]]$mvn[[2]]$Sigma
s$V[[2]]
fit$pars[[1]]$mvn[[3]]$Sigma
s$V[[3]]

Fixed Effects Hidden Markov Model

Use kmeans to get initial class memberships, but use actual means to avoid label permutation issues, then fit a hidden Markov model to each sequence forcing common transition probabilities across sequences

km <- kmeans(s$y,s$Mu)
fit <- gmvnHMM(s$y,km$cluster,s$gr,common.transition = TRUE)

Compare classifications across all groups

table(unlist(lapply(fit$pars,function(p) max.col(p$P))),s$cl)

Compare group means

## Check parameter estimates match simulation
opar <- par(mfrow=c(2,2),mar=c(1,1,1,1)+0.1)
matplot(s$A[,,1],t(sapply(fit$pars,function(p) p$mvn[[1]]$mu)),pch=16,cex=0.5,xlab="",ylab="")
matplot(s$A[,,2],t(sapply(fit$pars,function(p) p$mvn[[2]]$mu)),pch=16,cex=0.5,xlab="",ylab="")
matplot(s$A[,,3],t(sapply(fit$pars,function(p) p$mvn[[3]]$mu)),pch=16,cex=0.5,xlab="",ylab="")
par(opar)

Compare covariances

fit$pars[[1]]$mvn[[1]]$Sigma
s$V[[1]]
fit$pars[[1]]$mvn[[2]]$Sigma
s$V[[2]]
fit$pars[[1]]$mvn[[3]]$Sigma
s$V[[3]]

Compare transition probabilities

fit$pars[[1]]$Q
s$Q

Random Effects Mixture Model

Use kmeans to get initial class memberships, but use actual means to avoid label permutation issues, then fit a mixture model to each sequence forcing common mixing fractions across sequences

km <- kmeans(s$y,s$Mu)
fit <- grmvnMix(s$y,km$cluster,s$gr,common.fractions = TRUE)

Compare classifications across all groups

table(unlist(lapply(fit$pars,function(p) max.col(p$P))),s$cl)

Compare group means

## Check parameter estimates match simulation
opar <- par(mfrow=c(2,2),mar=c(1,1,1,1)+0.1)
matplot(s$A[,,1],t(sapply(fit$pars,function(p) p$mvn[[1]]$mu)),pch=16,cex=0.5,xlab="",ylab="")
matplot(s$A[,,2],t(sapply(fit$pars,function(p) p$mvn[[2]]$mu)),pch=16,cex=0.5,xlab="",ylab="")
matplot(s$A[,,3],t(sapply(fit$pars,function(p) p$mvn[[3]]$mu)),pch=16,cex=0.5,xlab="",ylab="")
par(opar)

Compare covariances

fit$pars[[1]]$mvn[[1]]$Sigma
s$V[[1]]
fit$pars[[1]]$mvn[[2]]$Sigma
s$V[[2]]
fit$pars[[1]]$mvn[[3]]$Sigma
s$V[[3]]

Compare covariances of group means

fit$muv[[1]]$U
s$U[[1]]
fit$muv[[2]]$U
s$U[[2]]
fit$muv[[3]]$U
s$U[[3]]

Random Effects Hidden Markov Model

Use kmeans to get initial class memberships, but use actual means to avoid label permutation issues, then fit a hidden Markov model to each sequence forcing common transition probabilities across sequences

km <- kmeans(s$y,s$Mu)
fit <- grmvnHMM(s$y,km$cluster,s$gr,common.transition = TRUE)

Compare classifications across all groups

table(unlist(lapply(fit$pars,function(p) max.col(p$P))),s$cl)

Compare group means

## Check parameter estimates match simulation
opar <- par(mfrow=c(2,2),mar=c(1,1,1,1)+0.1)
matplot(s$A[,,1],t(sapply(fit$pars,function(p) p$mvn[[1]]$mu)),pch=16,cex=0.5,xlab="",ylab="")
matplot(s$A[,,2],t(sapply(fit$pars,function(p) p$mvn[[2]]$mu)),pch=16,cex=0.5,xlab="",ylab="")
matplot(s$A[,,3],t(sapply(fit$pars,function(p) p$mvn[[3]]$mu)),pch=16,cex=0.5,xlab="",ylab="")
par(opar)

Compare covariances

fit$pars[[1]]$mvn[[1]]$Sigma
s$V[[1]]
fit$pars[[1]]$mvn[[2]]$Sigma
s$V[[2]]
fit$pars[[1]]$mvn[[3]]$Sigma
s$V[[3]]

Compare covariances of group means

fit$muv[[1]]$U
s$U[[1]]
fit$muv[[2]]$U
s$U[[2]]
fit$muv[[3]]$U
s$U[[3]]

Compare transition probabilities

fit$pars[[1]]$Q
s$Q

SWotherspoon/MVNSeq documentation built on June 1, 2022, 10:49 p.m.