fit_mvn_MCAR: Multivariate normal spatial mixture model clustering

In spruce: Spatial Random Effects Clustering of Single Cell Data

Multivariate normal spatial mixture model clustering

Description

Implement Gibbs sampling for MVN model with MCAR spatial random effects

Usage

fit_mvn_MCAR(Y, coords_df, K, nsim = 2000, burn = 1000, z_init = NULL)

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.
coords_df	An n x 2 data frame or matrix of 2d spot coordinates.
K	The number of mixture components to fit.
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. Randomly initialized if NULL.

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)
A <- build_knn_graph(as.matrix(coords_df),k = 4)
                                 
n <- nrow(coords_df) # number of observations
g <- 2 # number of features
K <- length(unique(coords_df_sim$z)) # number of clusters (mixture components)
pi <- table(z)/length(z) # cluster membership probability

# 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
)

# generate phi - not cluster specific
# conditional covariance of phi_i given phi_noti
m <- colSums(A)
M <- diag(m)
V <- matrix(0.4,nrow = g, ncol = g) # CAR covariance
diag(V) <- 0.6
V_true <- V
rho <- 0.999999 # Spatial dependence parameter ~ 1 for intrinsic CAR
Q <- diag(m) - rho*A	# m is number of neighbors for each spot
covphi <- solve(Q) %x% V # gn x gn covariance of phis
phi <- mvtnorm::rmvnorm(1, sigma=covphi)		# gn vector of spatial effects
PHI <- matrix(phi, ncol=g, byrow=TRUE) 	# n x g matrix of spatial effects
PHI <- t(scale(t(PHI)))

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

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

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