R/hessianPoisson.R

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

hessianPoisson <- function(beta, y, w, X, 
        variant) {
        n <- NROW(w)
        J <- NCOL(w)
        K <- NCOL(X)
        y <- vapply(seq(1,J), function(j) y, rep(0, n))
        eta <- vapply(seq(1,J), function(j) {
            Xj <- cbind(X[, , j])
            betaj <- beta[, j, drop = FALSE]
            Xj %*% betaj
        }, rep(0, n))
        lambda <- exp(eta)  ##
        h <- w * exp(-exp(eta)) * exp(y * eta)
        
        g <- rowSums(h)
        dlogP <- function(k) {
            if (variant[k]) 
                colSums(dh(k)/g) else sum(rowSums(dh(k))/g)
        }
        dh <- function(k) h * (y - lambda) * 
            X[, k, ]  #######
        S <- NULL
        for (k in seq_len(K)) {
            Sk <- dlogP(k)
            names(Sk) <- rep(paste("var", k, 
                sep = ""), length(Sk))
            S <- c(S, Sk)
        }
        
        d2h <- function(k, kprime, j, jprime) X[, 
            k, j] * dh(kprime)[, jprime] * (y - 
            lambda)[, j] - X[, k, j] * X[, kprime, 
            jprime] * h[, j] * lambda[, jprime]  #######
        d2logP <- function(k, kprime, variantk, 
            variantkprime) {
            if (variantk & variantkprime) {
                res <- matrix(0, J, J)
                for (j in seq_len(J)) {
                    for (jprime in seq_len(J)) {
                    dhk <- dh(k)[, j]
                    dhkprime <- dh(kprime)[, 
                        jprime]
                    d2hkkprime <- d2h(k, kprime, 
                        j, jprime)
                    if (j != jprime) 
                        res[j, jprime] <- sum(-dhk * 
                        dhkprime/g^2)
                    if (j == jprime) 
                        res[j, jprime] <- sum((d2hkkprime * 
                        g - dhk * dhkprime)/g^2)
                    }
                }
            }
            if (variantk & !variantkprime) {
                res <- matrix(NA, J, 1)
                for (j in seq_len(J)) {
                    dhk <- dh(k)[, j]
                    dhkprime <- rowSums(dh(kprime))
                    d2hkkprime <- d2h(k, kprime, 
                    j, j)
                    res[j, 1] <- sum((d2hkkprime * 
                    g - dhk * dhkprime)/g^2)
                }
                res <- cbind(res)
            }
            if (!variantk & variantkprime) {
                res <- matrix(NA, 1, J)
                for (j in seq_len(J)) {
                    dhk <- rowSums(dh(k))
                    dhkprime <- dh(kprime)[, 
                    j]
                    d2hkkprime <- d2h(k, kprime, 
                    j, j)
                    res[1, j] <- sum((d2hkkprime * 
                    g - dhk * dhkprime)/g^2)
                }
                res <- rbind(res)
            }
            if (!variantk & !variantkprime) {
                dhk <- rowSums(dh(k))
                dhkprime <- rowSums(dh(kprime))
                d2hkkprime <- rowSums(d2h(k, 
                    kprime, seq_len(J), seq_len(J)))
                res <- cbind(sum((d2hkkprime * 
                    g - dhk * dhkprime)/g^2))
            }
            res
        }
        H <- NULL
        for (k in seq_len(K)) {
            Hfilak <- NULL
            for (kprime in seq_len(K)) {
                res <- d2logP(k, kprime, variant[k], 
                    variant[kprime])
                colnames(res) <- rep(paste("var", 
                    kprime, sep = ""), NCOL(res))
                rownames(res) <- rep(paste("var", 
                    k, sep = ""), NROW(res))
                Hfilak <- cbind(Hfilak, res)
            }
            H <- rbind(H, Hfilak)
        }
        return(list(S = S, H = H))
}