R/JGGM.R

In equSA: Learning High-Dimensional Graphical Models

JGGM <- function(data,ALPHA1=0.05,ALPHA2=0.01,structure = "temporal",parallel=FALSE,nCPUs){
## run equSA_sub for each dataset, record the psi score ###
  U <- NULL;
  M=length(data);
if(parallel){
  cores <- nCPUs
  cl<-makeCluster(cores)
  registerDoParallel(cl,cores=cores)
  U <- foreach(i=1:M,.combine="cbind")  %dopar%
  {
    data_i <- data[[i]]
    U1 <- psical(data_i,ALPHA1=ALPHA1)
    if(i==1){
      U1
    }else{
      U1[,3]
    }
  }
  stopImplicitCluster()
  stopCluster(cl)
}else{
  for(i in 1:M)
  {
    data_i <- data[[i]]
    U1 <- psical(data_i,ALPHA1=ALPHA1)
    U<-cbind(U,U1[,3])
  }
  index <- U1[,1:2]
  U<-cbind(index,U)
}
  score_raw <- U


### pre_function ####
Ind <- function(x){
  k1=k2=0;
  for(i in 1:(length(x)-1))
  {
    if(x[i]!=x[i+1]){
      k1=k1+1
    }else{
      k2=k2+1 
    }
  }
  return(c(k1,k2))
}
  
Ind2 <- function(x){
    e_mod <- ifelse(mean(x)>=0.5,1,0)
    cc <- abs(x-e_mod)
    k1 <- sum(cc)
    k2 <- M-k1
    return(c(k1,k2))
}

#### setting ####
p <- dim(data[[1]])[2]
N <- p*(p-1)/2

#### hyperparameter ###

a1 <- 1
b1 <- 10
a2 <- 1
b2 <- 1

##### S #### 
permu <- function(n){
  S <- NULL;
  for(m in 0:n)
  {
    S <- rbind(S,t(apply(combn(1:n,m=m),2,function(cm) replace(rep(0,n),cm,1))))
  }
  return(S)
}
S <- permu(M)


U_initial <- U
psi_final_hat <- NULL;
num <- dim(S)[1]

for(kk in 1:N)
{
  psi <- as.numeric(U[kk,-c(1,2)])
  wi <- NULL;
  psi_com_hat <- NULL;
  for(ite in 1:num)
  {
    n0 <- length(which(S[ite,]==0))
    n1 <- length(which(S[ite,]==1))
    if(structure == "temporal"){
      k1 <- Ind(S[ite,])[1]  ## change
      k2 <- M-1-k1  ## not change
    }else if(structure == "spatial"){
      k1 <- Ind2(S[ite,])[1]  ## change
      k2 <- M-k1  ## not change
    }else{
      stop("please use temporal or spatial structure")
    }
    psi0 <- psi[S[ite,]==0]
    psi1 <- psi[S[ite,]==1]
    ## combined psi score ##
    psi_com <- NULL;
    psi0_com <- sum(psi0)/sqrt(length(psi0))
    psi1_com <- sum(psi1)/sqrt(length(psi1))
    psi_com[S[ite,]==0] <- psi0_com
    psi_com[S[ite,]==1] <- psi1_com
    psi_com_hat <- rbind(psi_com_hat,psi_com)
    #### calculate Prob ####
    p0 <- 1/sqrt(n0)*(1/sqrt(2*pi))^n0*gamma(n0/2-1/2+a2)*(1/2*sum(psi0^2)-1/(2*n0)*(sum(psi0))^2+b2)^(-(n0/2-1/2+a2))
    p1 <- 1/sqrt(n1)*(1/sqrt(2*pi))^n1*gamma(n1/2-1/2+a2)*(1/2*sum(psi1^2)-1/(2*n1)*(sum(psi1))^2+b2)^(-(n1/2-1/2+a2))
    if(n0==0)
    {
      wi[ite] <- beta(k1+a1,k2+b1)*p1
    }else if(n1==0){
      wi[ite] <- beta(k1+a1,k2+b1)*p0
    }else{
      wi[ite] <- beta(k1+a1,k2+b1)*p0*p1
    }
  } 
  Prob <- wi/sum(wi) 
  psi_final <- t(t(psi_com_hat)%*%Prob)
  psi_final_hat <- rbind(psi_final_hat,psi_final)
}
score_com <- cbind(U_initial[,c(1,2)],psi_final_hat)

### combine all psi score ###
psi_tran <- function(z){
  q<-pnorm(-abs(z), log.p=TRUE)
  q<-q+log(2.0)
  s<-qnorm(q,log.p=TRUE)
  s<-(-1)*s
  return(s)
}
tran_psi <- sapply(as.data.frame(score_com[,-(1:2)]),psi_tran)
tran_psi <- cbind(score_com[,1:2],tran_psi)


psi_all <- as.vector(as.matrix(tran_psi[,-c(1,2)]))
num <- 1:length(psi_all)
U <- cbind(num,num,psi_all)
U<-U[order(U[,3]), 1:3]
N<-length(U[,1])
ratio<-ceiling(N/100000)
m<-floor(N/ratio)
m0<-N-m*ratio
s<-sample.int(ratio,m,replace=TRUE)
for(i in 1:length(s)) s[i]<-s[i]+(i-1)*ratio
if(m0>0){
  s0<-sample.int(m0,1)+length(s)*ratio
  s<-c(s,s0)
}
Us<-U[s,]
y <- round(Us,6)
qqqscore <- pcorselR(y,ALPHA2=ALPHA2)


#### cut-off psi score qqqscore ###
A <- array(0,dim = c(M,p,p))
for(k in 1:M)
{
  psik <- tran_psi[,c(1,2,k+2)]
  psik <- psik[psik[,3]>=qqqscore,]
  for(i in 1:nrow(psik)){
    k1<-psik[i,1]
    k2<-psik[i,2]
    A[k,,][k1,k2]<-1
    A[k,,][k2,k1]<-1
  }
}
result <- list()

result$Array = A
result$score.sep = score_raw
result$score.joint = score_com
return(result)
}