# R/doEStep.R In metagenomeSeq: Statistical analysis for sparse high-throughput sequencing

#### Documented in doEStep

```#' Compute the Expectation step.
#'
#' Estimates the responsibilities \$z_ij = fracpi_j cdot I_0(y_ijpi_j cdot
#' I_0(y_ij + (1-pi_j) cdot f_count(y_ij
#'
#' Maximum-likelihood estimates are approximated using the EM algorithm where
#' we treat mixture membership \$delta_ij\$ = 1 if \$y_ij\$ is generated from the
#' zero point mass as latent indicator variables. The density is defined as
#' \$f_zig(y_ij = pi_j(S_j) cdot f_0(y_ij) +(1-pi_j (S_j))cdot
#' f_count(y_ij;mu_i,sigma_i^2)\$. The log-likelihood in this extended model is
#' \$(1-delta_ij) log f_count(y;mu_i,sigma_i^2 )+delta_ij log
#' pi_j(s_j)+(1-delta_ij)log (1-pi_j (sj))\$. The responsibilities are defined
#' as \$z_ij = pr(delta_ij=1 | data)\$.
#'
#' @param countResiduals Residuals from the count model.
#' @param zeroResiduals Residuals from the zero model.
#' @param zeroIndices Index (matrix m x n) of counts that are zero/non-zero.
#' @return Updated matrix (m x n) of estimate responsibilities (probabilities
#' that a count comes from a spike distribution at 0).
doEStep <-
function(countResiduals,  zeroResiduals, zeroIndices)
{
pi_prop=getPi(zeroResiduals)
w1=sweep(zeroIndices, 2, pi_prop, FUN="*")

countDensity=getCountDensity(countResiduals)
w2=sweep(countDensity, 2, 1-pi_prop, FUN="*")
z=w1/(w1+w2)
z[z>1-1e-6]=1-1e-6
z[!zeroIndices]=0
z
}
```

## Try the metagenomeSeq package in your browser

Any scripts or data that you put into this service are public.

metagenomeSeq documentation built on Nov. 1, 2018, 4:36 a.m.