fit_mvn_PG_smooth: Multivariate normal mixture model clustering - PG multinom...

In spruce: Spatial Random Effects Clustering of Single Cell Data

Multivariate normal mixture model clustering - PG multinom regression Spatial smoothing

Description

Implement Gibbs sampling for MVN model with spatial smoothing prior. Includes fixed effects multinomial regression on cluster indicators using Polya-Gamma data augmentation.

Usage

fit_mvn_PG_smooth(
  Y,
  W,
  coords_df,
  K,
  r = 3,
  nsim = 2000,
  burn = 1000,
  z_init = NULL,
  verbose = FALSE
)

Arguments

Y	An n x g matrix of gene expression values. n is the number of cell spots and g is the number of features.
W	An n x v matrix of covariates to predict cluster membership. Should include an intercept (i.e., first column is 1)
coords_df	An n x 2 data frame or matrix of 2d spot coordinates.
K	The number of mixture components to fit.
r	Empirical spatial smoothing
nsim	Number of total MCMC iterations to run.
burn	Number of MCMC iterations to discard as burn in. The number of saved samples is nsim - burn.
z_init	Optional initialized allocation vector. Initialized with hierarchical clustering if NULL.
verbose	Logical for printing cluster allocations at each iteration.

Value

a list of posterior samples

Examples

# parameters
data(coords_df_sim)
coords_df <- coords_df_sim[,1:2]
z <- remap_canonical2(coords_df_sim$z)
                                 
n <- nrow(coords_df) # number of observations
g <- 3 # number of features
K <- length(unique(coords_df_sim$z)) # number of clusters (mixture components)
pi <- table(z)/length(z) # cluster membership probability

W <- matrix(0, nrow = n, ncol = 2)
W[,1] <- 1
W[,2] <- sample(c(0,1),size = n, replace = TRUE, prob = c(0.5,0.5))

# Cluster Specific Parameters
Mu <- list(
  Mu1 = rnorm(g,-5,1),
  Mu2 = rnorm(g,0,1),
  Mu3 = rnorm(g,5,1),
  Mu4 = rnorm(g,-2,3)
)
# cluster specific variance-covariance
S <- matrix(1,nrow = g,ncol = g) # y covariance matrix
diag(S) <- 1.5
Sig <- list(
  Sig1 = S,
  Sig2 = S, 
  Sig3 = S,
  Sig4 = S
)

Y <- matrix(0, nrow = n, ncol = g)
for(i in 1:n)
{
  Y[i,] <- mvtnorm::rmvnorm(1,mean = Mu[[z[i]]],sigma = Sig[[z[i]]])
}

# fit model
# in practice use more mcmc iterations
fit <- fit_mvn_PG_smooth(Y = Y, coords_df = coords_df, W = W, K = K, nsim = 10, burn = 0)

spruce documentation built on March 18, 2022, 7:01 p.m.