This function computes the conditional probabilities t_{ik} that an observation i arises from the kth component for the current value of the mixture parameters.
1 
y 
(n x q) matrix of observed counts for n observations and q variables 
g 
Number of clusters 
conds 
Vector of length q defining the condition (treatment group) for each variable (column) in 
pi 
Vector of length 
s 
Vector of length q containing the estimates for the normalized library size parameters for each of the q variables in 
lambda 
(d x 
t 
(n x 
If all values of t_{ik} are 0 (or nearly zero), the observation is assigned with probability one to belong to the cluster with the closest mean (in terms of the Euclidean distance from the observation). To avoid calculation difficulties, extreme values of t_{ik} are smoothed, such that those smaller than 1e10 or larger than 11e10 are set equal to 1e10 and 11e10, respectively.
Andrea Rau <andrea.rau@jouy.inra.fr>
Rau, A., MaugisRabusseau, C., MartinMagniette, M.L., Celeux G. (2015). Coexpression analysis of highthroughput transcriptome sequencing data with Poisson mixture models. Bioinformatics, 31(9):14201427.
Rau, A., Celeux, G., MartinMagniette, M.L., MaugisRabusseau, C. (2011). Clustering highthroughput sequencing data with Poisson mixture models. Inria Research Report 7786. Available at http://hal.inria.fr/inria00638082.
PoisMixClus
for Poisson mixture model estimation and model selection;
PoisMixMean
to calculate the conditional percluster mean of each observation
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24  set.seed(12345)
## Simulate data as shown in Rau et al. (2011)
## Library size setting "A", high cluster separation
## n = 200 observations
simulate < PoisMixSim(n = 200, libsize = "A", separation = "high")
y < simulate$y
conds < simulate$conditions
s < colSums(y) / sum(y) ## TC estimate of lib size
## Run the PMMII model for g = 3
## "TC" library size estimate, EM algorithm
run < PoisMixClus(y, g = 3, norm = "TC",
conds = conds)
pi.est < run$pi
lambda.est < run$lambda
## Calculate the conditional probability of belonging to each cluster
proba < probaPost(y, g = 3, conds = conds, pi = pi.est, s = s,
lambda = lambda.est)
## head(round(proba,2))

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
All documentation is copyright its authors; we didn't write any of that.