R/EMmixt.R

In isglobal-brge/CNVassoc: Association analysis of CNVs and imputed SNPs incorporating uncertainty

EMmixt <- function(x, mu.surrog, sd.surrog, 
        pi, tol = 10^-5, max.iter = 1000, verbose = FALSE, 
        var.equal = TRUE) {
        rn <- rownames(x)
        if (length(dim(x)) > 1) 
            x <- x[, 1]
        mean.weight <- function(x, w) sum(x * 
            w)/sum(w)
        sd.weight <- function(x, w) sqrt(sum(w * 
            (x - mean.weight(x, w))^2)/sum(w))
        error <- Inf
        iter <- 1
        I <- NROW(x)
        J <- NROW(mu.surrog)
        while (error > tol & iter < max.iter) {
            theta.old <- c(mu.surrog, sd.surrog, 
                pi)
            ff <- vapply(seq(1,J), function(j) pi[j] * 
                dnorm(x, mu.surrog[j], sd.surrog[j]), rep(0, length(x)))
            pi.post <- ff/rowSums(ff)
            pi <- colMeans(pi.post)
            mu.surrog <- vapply(seq(1,J), function(j) mean.weight(x, 
                pi.post[, j]), 0)
            sd.surrog <- vapply(seq(1,J), function(j) sd.weight(x, 
                pi.post[, j]), 0)
            if (var.equal == TRUE) 
                sd.surrog <- rep(sqrt(sum(pi * 
                    sd.surrog^2)), J)
            theta <- c(mu.surrog, sd.surrog, 
                pi)
            error <- sqrt(sum((theta - theta.old)^2)/length(theta))
            if (verbose) 
                cat("-iter.", iter, " error=", 
                    error, "\n")
            iter <- iter + 1
        }
        if (iter == max.iter) 
            stop("Maximum iteration reached")
        out <- list()
        rownames(pi.post) <- rn
        out$z <- pi.post
        out$mu.surrog <- mu.surrog
        out$sd.surrog <- sd.surrog
        out$pi <- pi
        out$loglike <- sum(log(rowSums(ff)))
        out$numpar <- if (var.equal) 
            J * 2 + 1 else J * 3
        out$bic <- -2 * out$loglike + out$numpar * 
            log(I)
        out$aic <- -2 * out$loglike + 2 * out$numpar
        out$x <- x
        out
}