Score.Test.Interact: Implement the gene-centric gene-gene interaction effect test...

In SPA3G: SPA3G: R package for the method of Li and Cui (2012)

Description

Score.Test.Interact returns results of interaction test, including score statistic, p-value, and estimates of variance components.

Usage

Score.Test.Interact(Y, K1, K2, K3, par, method = "BFGS", test = TRUE)

Arguments

Y	numerical vector: quantitative phenotypes.
K1	matrix: kernel matrix of the first gene.
K2	matrix: kernel matrix of the second gene.
K3	matrix: elementwise multiplication of K1 and K2.
par	numerical vector: initial values of varicance components.
method	the method to be used in maximazing REML. the defalt method is "BFGS". Other options are Average Information "AI" and Fisher Scoreing "FS".
test	logical: if TRUE conduct the test.

Details

The length of the initial values (par) should be the same as the number of variance components you intend to estimate. And the score test can only be implemented under the null model (H0:) which has 3 variance components.

Value

VCs	REML estimats of variance components
Fisher.info	fisher information matrix
Beta	ML estimate of the overall mean
restricted.logLik	restricted log-likelihood
Score	score statistic
df	estimated degree of freedom for the scaled chi-square
scale	estimated scale parameter for the scaled chi-square
p.value	p-value of the test

Examples

## The function is currently defined as
function (Y, K1, K2, K3, par, method = "BFGS", test = TRUE) 
{
    p <- length(par)
    if (p != 3 & test == TRUE) 
        cat("Error: Not matched initial values!")
    theta.new <- par
    theta.old <- rep(0, p)
    X <- matrix(1, n, 1)
    Vs <- array(0, c(n, n, 4))
    Vs[, , 1] <- diag(1, n)
    Vs[, , 2] <- K1
    Vs[, , 3] <- K2
    Vs[, , 4] <- K3
    Sigma <- 0
    for (i in 1:p) {
        Sigma <- Sigma + theta.new[i] * Vs[, , i]
    }
    W <- solve(Sigma)
    R <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% 
        W
    kk <- g.old <- 0
    tt <- c()
    while (sum(abs(theta.new - theta.old)) > 1e-05 & kk < 100) {
        if (method == "BFGS") {
            s <- theta.new - theta.old
            theta.old <- theta.new
            g <- c()
            for (i in 1:p) {
                g[i] <- -t(Y) %*% R %*% Vs[, , i] %*% R %*% Y + 
                  TT(R, Vs[, , i])
            }
            delta <- g - g.old
            g.old <- g
            if (kk == 0 | t(s) %*% delta <= 0) {
                AI <- matrix(0, p, p)
                for (i in 1:p) {
                  for (j in i:p) {
                    AI[i, j] <- AI[j, i] <- t(Y) %*% R %*% Vs[, 
                      , i] %*% R %*% Vs[, , j] %*% R %*% Y
                  }
                }
                H_inv <- solve(AI)
            }
            else {
                rho <- c(1/(t(delta) %*% s))
                H_inv <- (diag(1, p) - (s %*% t(delta)) * rho) %*% 
                  H_inv %*% (diag(1, p) - rho * delta %*% t(s)) + 
                  rho * s %*% t(s)
            }
        }
        if (method == "AI") {
            theta.old <- theta.new
            g <- c()
            for (i in 1:p) {
                g[i] <- t(Y) %*% R %*% Vs[, , i] %*% R %*% Y - 
                  TT(R, Vs[, , i])
            }
            H <- matrix(0, p, p)
            for (i in 1:p) {
                for (j in i:p) {
                  H[i, j] <- H[j, i] <- -t(Y) %*% R %*% Vs[, 
                    , i] %*% R %*% Vs[, , j] %*% R %*% Y
                }
            }
            H_inv <- solve(H)
        }
        if (method == "FS") {
            theta.old <- theta.new
            g <- c()
            for (i in 1:p) {
                g[i] <- t(Y) %*% R %*% Vs[, , i] %*% R %*% Y - 
                  TT(R, Vs[, , i])
            }
            H <- matrix(0, p, p)
            for (i in 1:p) {
                AA <- R %*% Vs[, , i]
                for (j in i:p) {
                  BB <- R %*% Vs[, , j]
                  H[i, j] <- H[j, i] <- -TRACE(AA %*% BB)
                }
            }
            H_inv <- solve(H)
        }
        theta.new <- theta.old - H_inv %*% (g)
        alpha <- 0.5
        while (length(which(theta.new < 0)) > 0 & alpha > 1e-08) {
            theta.new <- theta.old - alpha * H_inv %*% (g)
            alpha <- alpha/2
        }
        theta.new[which(theta.new < 0)] <- 0
        Sigma.new <- 0
        for (i in 1:p) {
            Sigma.new <- Sigma.new + theta.new[i] * Vs[, , i]
        }
        W.new <- solve(Sigma.new)
        R <- W.new - W.new %*% X %*% solve(t(X) %*% W.new %*% 
            X) %*% t(X) %*% W.new
        kk <- kk + 1
    }
    a1 <- R %*% Vs[, , 1]
    a2 <- R %*% Vs[, , 2]
    a3 <- R %*% Vs[, , 3]
    a4 <- R %*% Vs[, , 4]
    b11 <- TT(a1, a1)
    b12 <- TT(a1, a2)
    b13 <- TT(a1, a3)
    b14 <- TT(a1, a4)
    b22 <- TT(a2, a2)
    b23 <- TT(a2, a3)
    b24 <- TT(a2, a4)
    b33 <- TT(a3, a3)
    b34 <- TT(a3, a4)
    b44 <- TT(a4, a4)
    if (test == FALSE) {
        eigen.sigma <- eigen(Sigma.new)
        lR <- -(sum(log(eigen.sigma$values)) + log(det(t(X) %*% 
            W.new %*% X)) + t(Y) %*% R %*% Y)/2
        H <- matrix(c(b11, b12, b13, b14, b12, b22, b23, b24, 
            b13, b23, b33, b34, b14, b24, b34, b44), 4, 4)/2
        beta <- solve(t(X) %*% W.new %*% X) %*% t(X) %*% W.new %*% 
            Y
        object <- list(VCs = theta.new, fisher.info = H, Beta = beta, 
            restricted.logLik = lR)
        return(object)
    }
    if (test == TRUE) {
        eigen.sigma <- eigen(Sigma.new)
        lR <- -(sum(log(eigen.sigma$values)) + log(det(t(X) %*% 
            W.new %*% X)) + t(Y) %*% R %*% Y)/2
        W0 <- W.new
        beta <- solve(t(X) %*% W0 %*% X) %*% t(X) %*% W0 %*% 
            Y
        Q <- t(Y - X %*% beta) %*% W0 %*% K3 %*% W0 %*% (Y - 
            X %*% beta)/2
        e <- TT(R, K3)/2
        Its <- c(b14, b24, b34)
        Iss <- matrix(c(b11, b12, b13, b12, b22, b23, b13, b23, 
            b33), 3, 3)
        Itt <- (b44 - Its %*% solve(Iss) %*% Its)/2
        k <- Itt/e/2
        v = 2 * e^2/Itt
        pvalue <- pchisq(Q/k, df = v, lower.tail = F)
        object <- list(VCs = theta.new, fisher.info = Iss/2, 
            Beta = beta, restricted.logLik = lR, Score = Q, df = v, 
            scale = k, p.value = pvalue)
        class(object) <- "Score Test: tau3=0"
        return(object)
    }
  }

SPA3G documentation built on May 2, 2019, 11:12 a.m.