Score.Test.Interact(SPA3G)
Score.Test.Interact()所属R语言包：SPA3G
Implement the gene-centric gene-gene interaction effect test for H0:
实施基因为中心的基因 - 基因互作效应试验H0：
译者：生物统计家园网 机器人LoveR
描述----------Description----------
Score.Test.Interact returns results of interaction test, including score statistic, p-value, and estimates of variance components.
Score.Test.Interact的互动测试，包括得分统计量，p值，方差分量估计的返回结果。
用法----------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.
K1和K2的矩阵：的elementwise的乘法。
参数：par
numerical vector: initial values of varicance components.
数值向量：初始值varicance组件。
参数：method
the method to be used in maximazing REML. the defalt method is "BFGS". Other options are Average Information "AI" and Fisher Scoreing "FS".
被用来在maximazing REML方法。 defalt方法是BFGS“。其他选项是的平均信息“AI”和Fisher Scoreing“FS”。
参数：test
logical: if TRUE conduct the test.
逻辑：如果为TRUE进行测试。
Details
详细信息----------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.
的初始值（PAR）的长度应该是相同的方差分量估计的数量。下，只能实行零模型（H 有3个方差分量的分数测试。
值----------Value----------
参数：VCs
REML estimats of variance components
方差分量的REML estimats
参数：Fisher.info
fisher information matrix
Fisher信息矩阵
参数：Beta
ML estimate of the overall mean
ML估计的总体平均
参数：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
p-值的测试
实例----------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)
}
}
转载请注明:出自 生物统计家园网(http://www.biostatistic.net)。
注：
注1：为了方便大家学习，本文档为生物统计家园网机器人LoveR翻译而成，仅供个人R语言学习参考使用，生物统计家园保留版权。
注2：由于是机器人自动翻译，难免有不准确之处，使用时仔细对照中、英文内容进行反复理解，可以帮助R语言的学习。
注3：如遇到不准确之处，请在本贴的后面进行回帖，我们会逐渐进行修订。
|
发表于 2012-9-30 12:07:25
