```#' Compute the Maximization step calculation for features still active.
#'
#' Maximization step is solved by weighted least squares.  The function also
#' computes counts residuals.
#'
#' 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)*f_0(y_ij) +(1-pi_j (S_j)) *
#' 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 (s_j))\$. The responsibilities are defined
#' as \$z_ij = pr(delta_ij=1 | data)\$.
#'
#' @param z Matrix (m x n) of estimate responsibilities (probabilities that a
#' count comes from a spike distribution at 0).
#' @param y Matrix (m x n) of count observations.
#' @param mmCount Model matrix for the count distribution.
#' @param stillActive Boolean vector of size M, indicating whether a feature
#' converged or not.
#' @param fit2 Previous fit of the count model.
#' @param dfMethod Either 'default' or 'modified' (by responsibilities)
#' @return Update matrix (m x n) of estimate responsibilities (probabilities
#' that a count comes from a spike distribution at 0).
doCountMStep <-
function(z, y, mmCount, stillActive,fit2=NULL,dfMethod="modified"){

if (is.null(fit2)){
fit=limma::lmFit(y[stillActive,],mmCount,weights = (1-z[stillActive,]))
if(dfMethod=="modified"){
df = rowSums(1-z[stillActive,,drop=FALSE]) - ncol(mmCount)
fit\$df[stillActive] = df
fit\$df.residual[stillActive] = df
}
countCoef = fit\$coefficients
countMu=tcrossprod(countCoef, mmCount)
residuals=sweep((y[stillActive,,drop=FALSE]-countMu),1,fit\$sigma,"/")
dat = list(fit = fit, residuals = residuals)
return(dat)
} else {

residuals = fit2\$residuals
fit2 = fit2\$fit

fit=limma::lmFit(y[stillActive,,drop=FALSE],mmCount,weights = (1-z[stillActive,,drop=FALSE]))

fit2\$coefficients[stillActive,] = fit\$coefficients
fit2\$stdev.unscaled[stillActive,]=fit\$stdev.unscaled
fit2\$sigma[stillActive] = fit\$sigma
fit2\$Amean[stillActive] = fit\$Amean

if(dfMethod=="modified"){
df = rowSums(1-z[stillActive,,drop=FALSE]) - ncol(mmCount)
fit\$df = df
fit\$df.residual = df
}
fit2\$df[stillActive]    = fit\$df
fit2\$df.residual[stillActive]    = fit\$df.residual

countCoef = fit\$coefficients
countMu=tcrossprod(countCoef, mmCount)
r=sweep((y[stillActive,,drop=FALSE]-countMu),1,fit\$sigma,"/")
residuals[stillActive,]=r

dat = list(fit = fit2, residuals=residuals)

return(dat)
}
}
```