R/GGInterScoreTest2.R

In fhebert/GeneGeneInteractions2: Global tests for the detection of gene-gene interaction effects on binary phenotypes

# GGInterScoreTest2 = function(X1,X2,Y,perm.method="parametric",N=1000){
#   W = cbind(PCsel(X1),PCsel(X2))
#   W = cbind(1,W)
#   n = length(Y)
#   XX1 = matrix(apply(X1,2,function(x){rep(x,ncol(X2))}),nrow=n)
#   XX2 = matrix(rep(X2,ncol(X1)),nrow=n)
#   S1 = as(XX1*XX2,"dgCMatrix")
#   ind = which(colSums(S1)==0)
#   if(length(ind)>0){
#     S1 = S1[,-ind]
#   }
#   X1 = as(X1,"dgCMatrix")
#   X2 = as(X2,"dgCMatrix")
#   X1.1 = X1==1
#   X1.2 = X1==2
#   X2.1 = X2==1
#   X2.2 = X2==2
#   p1 = ncol(X1)
#   p2 = ncol(X2)
#   A1 = as(matrix(apply(cbind(X1.1,X1.2),2,function(x){rep(x,2*p2)}),nrow=n),"dgCMatrix")
#   A2 = as(matrix(rep(cbind(X2.1,X2.2),2*p1),nrow=n),"dgCMatrix")
#   S2 = A1*A2
#   ind = which(colSums(S2)==0)
#   if(length(ind)>0){
#     S2 = S2[,-ind]
#   }
#   res.GLM = GLM(W,Y,rep(0,ncol(W)))
#   hatY0 = res.GLM$probs
#   Z1 = as.vector(crossprod(Y-hatY0,S1))
#   Z2 = as.vector(crossprod(Y-hatY0,S2))
#   s2Y = sum((Y-hatY0)^2)/length(Y)
#   #######
#   # W = as.matrix(W)
#   WS1 = crossprod(W,S1)
#   A1 = as.matrix(crossprod(S1)-t(WS1)%*%solve(crossprod(W))%*%WS1)
#   WS2 = crossprod(W,S2)
#   A2 = as.matrix(crossprod(S2)-t(WS2)%*%solve(crossprod(W))%*%WS2)
#   # A = (t(W)%*%W-t(W)%*%X%*%solve(t(X)%*%X)%*%t(X)%*%W)
#   #######
#   Gamma1 = s2Y*A1
#   Gamma2 = s2Y*A2
#   d1 = sqrt(diag(Gamma1))
#   d2 = sqrt(diag(Gamma2))
#   Z1 = Z1/d1
#   Z2 = Z2/d2
#   D1 = matrix(rep(d1,ncol(S1)),ncol=ncol(S1))
#   D2 = matrix(rep(d2,ncol(S2)),ncol=ncol(S2))
#   Sigma1 = Gamma1/(D1*t(D1))
#   Sigma2 = Gamma2/(D2*t(D2))
#   if(perm.method=="parametric"){
#     Y0 = ParametricResampling(hatY0,N)
#     hatY0 = PermGLMProbs(W,Y0,rep(0,ncol(W)))
#   }else{
#     Y0 = sapply(1:N,function(i){sample(Y)})
#     hatY0 = PermGLMProbs(W,Y0,rep(0,ncol(W)))
#   }
#   d0 = Y0-hatY0
#   Z01 = as.matrix(crossprod(d0,S1))
#   Z02 = as.matrix(crossprod(d0,S2))
#   # Z0 = t(d0)%*%W
#   sY0 = sqrt(crossprod(rep(1,length(Y)),d0^2)[1,]/length(Y))
#   # sY0 = sqrt((rep(1,length(Y))%*%(d0^2))[1,]/length(Y))
#   dA1 = sqrt(diag(A1))
#   dA2 = sqrt(diag(A2))
#   Z01 = Z01/(matrix(sY0,ncol=1)%*%matrix(dA1,nrow=1))
#   Z01 = as.matrix(Z01)
#   Z02 = Z02/(matrix(sY0,ncol=1)%*%matrix(dA2,nrow=1))
#   Z02 = as.matrix(Z02)
#   diag(Sigma1) = diag(Sigma2) = 1
#   ind = which(colSums(is.na(Sigma1))==(ncol(Sigma1)-1))
#   if(length(ind)>0){
#     Sigma1 = Sigma1[-ind,-ind]
#     Z1 = Z1[-ind]
#     Z01 = Z01[,-ind]
#   }
#   ind = which(colSums(is.na(Sigma2))==(ncol(Sigma2)-1))
#   if(length(ind)>0){
#     Sigma2 = Sigma2[-ind,-ind]
#     Z2 = Z2[-ind]
#     Z02 = Z02[,-ind]
#   }
#   return(list(Z1=Z1,Z01=Z01,Sigma1=Sigma1,Z2=Z2,Z02=Z02,Sigma2=Sigma2))
# }