R/EdgeOrientation.R
In audreyqyfu/mrpc: PC Algorithm with the Principle of Mendelian Randomization
#This is the step 2 of MRPC to direction determination of the undirected edges.
EdgeOrientation <- function (gInput, GV, suffStat, FDR, alpha, indepTest,
FDRcontrol = c("LOND", "ADDIS", "NONE"), tau = 0.5, lambda = 0.25,
verbose = FALSE) {
g <- as(gInput@graph, "matrix") # g ia an adjacency from undirected graph (skleton)
g1 <- g
p <- nrow(g)
# tarmat (adjacency matrix for directed graph) is updated in every step,
# and contains the output of final topology
tarmat <- matrix(0, nrow(g),ncol(g)) #same row and column from g
rownames(tarmat) <- rownames(g) #same row names from g
colnames(tarmat) <- colnames(g) #same column names from g
# extract all edges
g[lower.tri(g)] <- 0 # Use only upper triangular because g is symmetric matrix
edges <- which (g==1, arr.ind = TRUE)
#Step-1 start
if (GV>0) {
# identify edges involving Vs
edgesWithBothVs <- edges[which (edges[,1] <= GV & edges[,2] <= GV), ]
edgesWithVFirst <- edges[which (edges[,1] <= GV & edges[,2] > GV), ]
edgesWithVSecond <- edges[which (edges[,1] > GV & edges[,2] <= GV), ]
# assign 1s to corresponding edges in tarmat
# edges between two Vs are undirected (or bidirected)
if (length (edgesWithBothVs) > 2) {
tarmat[edgesWithBothVs] <- 1
tarmat[edgesWithBothVs[,2:1]] <- 1
} else {
tarmat[edgesWithBothVs[1], edgesWithBothVs[2]] <- 1
tarmat[edgesWithBothVs[2], edgesWithBothVs[1]] <- 1
}
# Edges involving one V go from V to the other node
# Need to distinguish whether edgesWithVFirst or edgesWithVSecond is
# a matrix or a vector.
if (length (edgesWithVFirst) > 2) {
tarmat[edgesWithVFirst] <- 1
} else {
tarmat[edgesWithVFirst[1], edgesWithVFirst[2]] <- 1
}
if (length(edgesWithVSecond) > 2) {
tarmat[edgesWithVSecond[,2:1]] <- 1
} else {
tarmat[edgesWithVSecond[2], edgesWithVSecond[1]] <- 1
}
}
#Step-1 end
#Step-2 start
#Start to orient v-structures
if (verbose)
cat("\n V-structures are as follows :\n")
# extract edges involving at least one gene node
#edgesWithGs <- edges[which (edges[,1]>GV | edges[,2]>GV), ]
#ind <- edgesWithGs
m <- gInput@test #current test
Alpha <- gInput@alpha #alpha
pval <- gInput@pval
R <- gInput@R #decision of test
ind <- which(g1 == 1, arr.ind = TRUE) #Pullout the all relation in adjacency matrix from undirected graph
V <- colnames(g)
# Update the addis vectors/values output from the ModiSkeleton function.
w0 <- tau * lambda * FDR/2
alphai <- gInput@alphai
kappai <- gInput@kappai
kappai_star <- gInput@kappai_star
K <- gInput@K
Ci <- gInput@Ci
Ci_sum <- sum(Ci)
Si <- gInput@Si
Si_sum <- sum(Si)
gammai_sum <- gInput@gammai_sum
Ci_plus <- gInput@Ci_plus
gammai <- gInput@gammai
normalizer <- gInput@normalizer
exponent <- gInput@exponent
for (i in seq_len(nrow(ind))) {
x <- ind[i, 1]
y <- ind[i, 2]
allZ <- setdiff(which(g1[y, ] == 1), x)
for (z in allZ) {
# Triplet x-y-z is directed x-->y<--z if x and z conditionally dependent given y
if ((g1[x, z] == 0 & g1[x, y] == 1 & g1[y, z] == 1) & !(tarmat[y, x] ==1) & !(tarmat[z, y] ==1) & !(tarmat[y, z] ==1) &
!(y %in% gInput@sepset[[x]][[z]] || y %in% gInput@sepset[[z]][[x]]))
{
m <- m+1
# Increase the length of the R, pval, and other addis vectors if
# their length is greater than or equal to the current iteration.
if (m >= length(R)) {
R <- c(R, numeric(m / 2))
pval <- c(pval, numeric(m / 2))
kappai <- c(kappai, numeric(m / 2))
kappai_star <- c(kappai_star, numeric(m / 2))
Ci <- c(Ci, numeric(m / 2))
Si <- c(Si, numeric(m / 2))
Ci_plus <- c(Ci_plus, numeric(m / 2))
gammai <- c(gammai, numeric(m / 2))
alphai <- c(alphai, numeric(m / 2))
}
if(indepTest=="gaussCItest") #if indepTest=gaussCItest
{
pval[m] <- gaussCItest(x, z, y, suffStat)
}
if(indepTest=="disCItest") #if indepTest=disCItest
{
pval[m] <- disCItest(x, z, y, suffStat) #additional
}
#pval=disCItest(x, z, y, suffStat) #additional conditional test
if (FDRcontrol == 'LOND') {
# Calculate alpha using the LOND method.
alphai[m] <- SeqFDR(m,FDR,a=2,R)
Alpha <- alphai[m]
# Calculate the alpha value using the ADDIS algorithm
} else if (FDRcontrol == 'ADDIS') {
# Run ADDIS on the current iteration and update all the vectors
# and other values.
run_addis <- addis(alpha = FDR,
tau = tau,
lambda = lambda,
iter = m,
w0 = w0,
pval = pval,
alphai = alphai,
gammai = gammai,
kappai = kappai,
kappai_star = kappai_star,
K = K,
Ci = Ci,
Si = Si,
Ri = R,
Ci_plus = Ci_plus,
Ci_sum = Ci_sum,
Si_sum = Si_sum,
gammai_sum = gammai_sum,
normalizer = normalizer,
exponent = exponent)
# Update all values and vectors output from the addis function.
alphai[[m]] <- run_addis[[1]]
gammai[[m]] <- run_addis[[2]]
K <- run_addis[[3]]
R[[m]] <- run_addis[[5]]
Si[[m - 1]] <- run_addis[[6]]
Ci[[m - 1]] <- run_addis[[7]]
Ci_sum <- run_addis[[9]]
Si_sum <- run_addis[[10]]
gammai_sum <- run_addis[[12]]
# Only update the kappai and Ci_plus vectors if K is greater
# than one. If K is zero then the first element in kappai will
# remain zero until the first rejection.
if (K != 0) {
kappai[[K]] <- run_addis[[4]]
kappai_star[[K]] <- run_addis[[11]]
Ci_plus[1:K] <- run_addis[[8]]
}
# Update the Alpha value with the current alphai value.
Alpha <- alphai[[m]]
} else if (FDRcontrol == 'NONE') {
Alpha <- alpha
}
if (verbose){
cat("x=", x, " y=", z, " S=", y,"\n")
cat("Test number =", m, "\n")
cat("Additional pval value =", pval[m], "\n")
cat("Alpha value =", Alpha, "\n")
}
if (pval[m] <= Alpha) { #Reject H0 (H0:nodes are independent)
R[m] <- 1
if (verbose) {
cat(V[x], "->", V[y], "<-", V[z], "\n") #Printout the v-structures
cat("Since pval<Alpha,additional test is rejected;", "Nodes", V[x] ,"and" ,V[z] ,"are dependent given", V[y], "\n")
}
tarmat[x, y] <- tarmat[z, y] <- 1 #directed x-->y<--z
}
else {
R[m] <- 0 #Accept H0
#tarmat[x, y] <- tarmat[y, z] <- 1 #directed x-->y-->z
if (verbose)
cat("Since pval>Alpha,additional test is accepted;", "Nodes", V[x] ,"and" ,V[z] ,"are independent given", V[y], "\n")
}
}
}
}
#Step-2 end to orient v-structures
#Step-3
#Orient remaining edges, weather gMR is applicable or not based on the following process:
#Pullout all edges from g1 matrix and define 1st and 2nd node (goal to make a triplet for the remaining edges)
#extract all edges with involving those 1st and 2nd nodes (from row and column of g)
#ignore the nodes that already have direction
#form a triplet
#determine the direction using test results from part 1 based on gMR
#Repeat until all undirected edges to directed
# m <- m
# Alpha <- Alpha
# R <- R
#start
#when data contain genetic variants
if (any(tarmat == 1) & GV>0) #if at least one edge directed already made so far
{
WW1 <- unique(which(tarmat==1,arr.ind = T)[,2]) #pullout the all canidate genes for v-structure
#if(WW1<GV){
#WW1=WW1[-c(1:GV)] #ignor canidate genes for GV
WW1 <- setdiff(WW1,1:GV)
#}
if(length(WW1)!=0)
{
#WW1=setdiff(WW1,1:GV)
for (v1 in 1:length(WW1))
{
WW2 <- unique(which(tarmat[,WW1[v1]]==1,arr.ind = T)) #edge between v-structure
WW3 <- unique(which(g1[,WW1[v1]]==1,arr.ind = T)) #edges between canidate genes and others
WW3 <- setdiff(WW3,WW2) #ignore if already direction
if(length(WW3)!=0)
{
for (v2 in 1:length(WW3))
{
x <- WW2[1]
y <- WW1[v1]
z <- WW3[v2]
if ((!is.na(x) & !is.na(y) & !is.na(z)) & (x!=0 & y!=0 & z!=0) & (x!=y & y!=z & z!=x))
{
#Case-2: If,y and z are adjacent, x and z conditionally independent given y,
#then the edge direction will be y-->z
if (g1[y, z] == 1 & tarmat[y, z]!=1 & tarmat[z, y]!=1 & ((y %in% gInput@sepset[[x]][[z]]) || (y %in% gInput@sepset[[z]][[x]])))
{
tarmat[y, z] <- 1
tarmat[z, y] <- 0
}
#Case-3: If,y and z are adjacent, x and z conditionally dependent given y,
#then the edge direction will be z-->y.
if (g1[y, z] == 1 & g1[x, z] != 1 & tarmat[x,y]==1 & tarmat[y, x]!=1 & tarmat[y, z]!=1 & tarmat[z, y]!=1 &!(y %in% gInput@sepset[[x]][[z]]) & !(y %in% gInput@sepset[[z]][[x]]))
{
m <- m+1
# Increase the length of the R, pval, and other addis vectors if
# their length is greater than or equal to the current iteration.
if (m >= length(R)) {
R <- c(R, numeric(m / 2))
pval <- c(pval, numeric(m / 2))
kappai <- c(kappai, numeric(m / 2))
kappai_star <- c(kappai_star, numeric(m / 2))
Ci <- c(Ci, numeric(m / 2))
Si <- c(Si, numeric(m / 2))
Ci_plus <- c(Ci_plus, numeric(m / 2))
gammai <- c(gammai, numeric(m / 2))
alphai <- c(alphai, numeric(m / 2))
}
if(indepTest=="gaussCItest") #if indepTest=gaussCItest for continuous data
{
pval[m] <- gaussCItest(x, z, y, suffStat) #additional pval
}
if(indepTest=="disCItest") #if indepTest=disCItest for discrete data
{
pval[m] <- disCItest(x, z, y, suffStat) #additional pval
}
if (FDRcontrol == 'LOND') {
# Calculate alpha using the LOND method.
alphai[m] <- SeqFDR(m,FDR,a=2,R)
Alpha <- alphai[m]
# Calculate the alpha value using the ADDIS algorithm
} else if (FDRcontrol == 'ADDIS') {
# Run ADDIS on the current iteration and update all the vectors
# and other values.
run_addis <- addis(alpha = FDR,
tau = tau,
lambda = lambda,
iter = m,
w0 = w0,
pval = pval,
alphai = alphai,
gammai = gammai,
kappai = kappai,
kappai_star = kappai_star,
K = K,
Ci = Ci,
Si = Si,
Ri = R,
Ci_plus = Ci_plus,
Ci_sum = Ci_sum,
Si_sum = Si_sum,
gammai_sum = gammai_sum,
normalizer = normalizer,
exponent = exponent)
# Update all values and vectors output from the addis function.
alphai[[m]] <- run_addis[[1]]
gammai[[m]] <- run_addis[[2]]
K <- run_addis[[3]]
R[[m]] <- run_addis[[5]]
Si[[m - 1]] <- run_addis[[6]]
Ci[[m - 1]] <- run_addis[[7]]
Ci_sum <- run_addis[[9]]
Si_sum <- run_addis[[10]]
gammai_sum <- run_addis[[12]]
# Only update the kappai and Ci_plus vectors if K is greater
# than one. If K is zero then the first element in kappai will
# remain zero until the first rejection.
if (K != 0) {
kappai[[K]] <- run_addis[[4]]
kappai_star[[K]] <- run_addis[[11]]
Ci_plus[1:K] <- run_addis[[8]]
}
# Update the Alpha value with the current alphai value.
Alpha <- alphai[[m]]
} else if (FDRcontrol == 'NONE') {
Alpha <- alpha
}
if (verbose){
cat("Additional pval value =", pval[m], "\n")
cat("Alpha value =", Alpha, "\n")
}
if (pval[m] <= Alpha) { #Reject H0 (H0:nodes are independent)
R[m] <- 1
if (verbose) {
V <- colnames(g)
cat("\n", V[x], "->", V[y], "<-", V[z], "\n") #Printout the v-structures
cat("Since pval<Alpha,additional test is rejected", "Nodes", V[x] ,"and" ,V[z] ,"are dependent given", V[y], "\n")
}
tarmat[z, y] <- 1 #directed z-->y
}
else {
tarmat[y, z] <- 1
R[m] <- 0 #Accept H0
}
}
#Case-4:If, y & z have relation and x & y conditionally dependent given z and x & z conditionally independent given y,
#then edge direction will be z<-->y
if (g1[y, z] == 1 & g1[x, z] == 1 & g1[x, y] == 1 & !(z %in% gInput@sepset[[x]][[y]]) &!(y %in% gInput@sepset[[x]][[z]]))
{
tarmat[y, z] <- 1
tarmat[z, y] <- 1
}
}
}
}
}
#if(length(WW3)!=0)
#{
#Orient remaining edges,
T1 <- which(tarmat==1,arr.ind = T)
G1 <- which(g==1,arr.ind = T)
D1 <- setdiff(G1,T1)
#
if(length(D1)!=0)
{
for (d1 in 1:length(D1))
{
Rem1_row <- which(g[,D1[d1]]==1,arr.ind = T)
Rem1_col <- which(g[D1[d1],]==1,arr.ind = T)
Rem1 <- c(Rem1_row,Rem1_col)
Rem2 <- which(tarmat[,Rem1]==1,arr.ind = T)
if(length(Rem1)!=0 & length(Rem2)!=0)
{
for (d11 in 1:length(Rem1))
{
x <- Rem2[d11]
y <- Rem1[d11]
z <- D1[d1]
}
if ((!is.na(x) & !is.na(y) & !is.na(z)) & (x!=0 & y!=0 & z!=0) & (x!=y & y!=z & z!=x))
{
#Case-2: If,y and z are adjacent, x and z conditionally independent given y,
#then the edge direction will be y-->z
if (g1[y, z] == 1 & tarmat[y, z]!=1 & tarmat[z, y]!=1 & ((y %in% gInput@sepset[[x]][[z]]) || (y %in% gInput@sepset[[z]][[x]])))
{
tarmat[y, z] <- 1
tarmat[z, y] <- 0
}
#Case-3: If,y and z are adjacent, x and z conditionally dependent given y,
#then the edge direction will be z-->y.
if (g1[y, z] == 1 & g1[x, z] != 1 & tarmat[x,y]==1 & tarmat[y, x]!=1 & tarmat[y, z]!=1 & tarmat[z, y]!=1 &!(y %in% gInput@sepset[[x]][[z]]) & !(y %in% gInput@sepset[[z]][[x]]))
{
m <- m+1
# Increase the length of the R, pval, and other addis vectors if
# their length is greater than or equal to the current iteration.
if (m >= length(R)) {
R <- c(R, numeric(m / 2))
pval <- c(pval, numeric(m / 2))
kappai <- c(kappai, numeric(m / 2))
kappai_star <- c(kappai_star, numeric(m / 2))
Ci <- c(Ci, numeric(m / 2))
Si <- c(Si, numeric(m / 2))
Ci_plus <- c(Ci_plus, numeric(m / 2))
gammai <- c(gammai, numeric(m / 2))
alphai <- c(alphai, numeric(m / 2))
}
if(indepTest=="gaussCItest") #if indepTest=gaussCItest
{
pval[m] <- gaussCItest(x, z, y, suffStat)
}
if(indepTest=="disCItest") #if indepTest=gaussCItest
{
pval[m] <- disCItest(x, z, y, suffStat) #additional
}
if (FDRcontrol == 'LOND') {
# Calculate alpha using the LOND method.
alphai[m] <- SeqFDR(m,FDR,a=2,R)
Alpha <- alphai[m]
# Calculate the alpha value using the ADDIS algorithm
} else if (FDRcontrol == 'ADDIS') {
# Run ADDIS on the current iteration and update all the vectors
# and other values.
run_addis <- addis(alpha = FDR,
tau = tau,
lambda = lambda,
iter = m,
w0 = w0,
pval = pval,
alphai = alphai,
gammai = gammai,
kappai = kappai,
kappai_star = kappai_star,
K = K,
Ci = Ci,
Si = Si,
Ri = R,
Ci_plus = Ci_plus,
Ci_sum = Ci_sum,
Si_sum = Si_sum,
gammai_sum = gammai_sum,
normalizer = normalizer,
exponent = exponent)
# Update all values and vectors output from the addis function.
alphai[[m]] <- run_addis[[1]]
gammai[[m]] <- run_addis[[2]]
K <- run_addis[[3]]
R[[m]] <- run_addis[[5]]
Si[[m - 1]] <- run_addis[[6]]
Ci[[m - 1]] <- run_addis[[7]]
Ci_sum <- run_addis[[9]]
Si_sum <- run_addis[[10]]
gammai_sum <- run_addis[[12]]
# Only update the kappai and Ci_plus vectors if K is greater
# than one. If K is zero then the first element in kappai will
# remain zero until the first rejection.
if (K != 0) {
kappai[[K]] <- run_addis[[4]]
kappai_star[[K]] <- run_addis[[11]]
Ci_plus[1:K] <- run_addis[[8]]
}
# Update the Alpha value with the current alphai value.
Alpha <- alphai[[m]]
} else if (FDRcontrol == 'NONE') {
Alpha <- alpha
}
if (pval[m] <= Alpha) { #Reject H0 (H0:nodes are independent)
R[m] <- 1
if (verbose) {
V <- colnames(g)
cat("\n", V[x], "->", V[y], "<-", V[z], "\n") #Printout the v-structures
cat("Since pval<Alpha,additional test is rejected", "Nodes", V[x] ,"and" ,V[z] ,"are dependent given", V[y], "\n")
}
tarmat[z, y] <- 1 #directed z-->y
}
else {
tarmat[y, z] <- 1
R[m] <- 0 #Accept H0
}
}
#Case-4:If, y & z have relation and x & y conditionally dependent given z and x & z conditionally independent given y,
#then edge direction will be z<-->y
if (g1[y, z] == 1 & g1[x, z] == 1 & g1[x, y] == 1 & !(z %in% gInput@sepset[[x]][[y]]) &!(y %in% gInput@sepset[[x]][[z]]))
{
tarmat[y, z] <- 1
tarmat[z, y] <- 1
}
}
}
}
}
#Check the reamning edge orientation
ind <- which(g1 == 1, arr.ind = TRUE) #Pullout the all relation in adjacency matrix from undirected graph
for (i in seq_len(nrow(ind))) {
x <- ind[i, 1]
y <- ind[i, 2]
if(tarmat[x, y]==0 & tarmat[y, x]==0) #bidirected if still no edge in tarmat
{
tarmat[x, y] <- 1
tarmat[y, x] <- 1
}
}
}
}#end when data contain genetic variants
#Start to orient remaining edges, when data not necessary to contain genetic variants.
#First identify all the canidate genes for v-structure
#Then identify the edges (undirected) with canidate genes
#ignore the nodes that already have direction
#Make a triplet
#determine the direction using test results from part 1 based on gMR
#Repeat until all undirected edges to directed
#Start when data not necessary to contain genetic variants.
if (any(tarmat == 1) & GV==0) #if at least one edge directed already made so far
{
WW1 <- unique(which(tarmat == 1,arr.ind = T)[,2]) #pullout the all canidate genes for v-structure
#WW1=WW1[-c(1:GV)]
if(length(WW1)!=0)
{
#WW1=setdiff(WW1,1:GV)
for (v1 in 1:length(WW1))
{
WW2 <- unique(which(tarmat[,WW1[v1]]==1,arr.ind = T)) #edge between v-structure
WW3 <- unique(which(g1[,WW1[v1]]==1,arr.ind = T)) #edges between canidate genes and others
WW3 <- setdiff(WW3,WW2) #ignore if already direction
if(length(WW3)!=0)
{
for (v2 in 1:length(WW3))
{
x <- WW2[1]
y <- WW1[v1]
z <- WW3[v2]
if ((!is.na(x) & !is.na(y) & !is.na(z)) & (x!=0 & y!=0 & z!=0) & (x!=y & y!=z & z!=x))
{
#Case-2: If,y and z are adjacent, x and z conditionally independent given y,
#then the edge direction will be y-->z
if (g1[y, z] == 1 & tarmat[y, z]!=1 & tarmat[z, y]!=1 & ((y %in% gInput@sepset[[x]][[z]]) || (y %in% gInput@sepset[[z]][[x]])))
{
tarmat[y, z] <- 1
tarmat[z, y] <- 0
}
#Case-3: If,y and z are adjacent, x and z conditionally dependent given y,
#then the edge direction will be z-->y.
if (g1[y, z] == 1 & g1[x, z] != 1 & tarmat[x,y]==1 & tarmat[y, x]!=1 & tarmat[y, z]!=1 & tarmat[z, y]!=1 &!(y %in% gInput@sepset[[x]][[z]]) & !(y %in% gInput@sepset[[z]][[x]]))
{
m <- m+1
# Increase the length of the R, pval, and other addis vectors if
# their length is greater than or equal to the current iteration.
if (m >= length(R)) {
R <- c(R, numeric(m / 2))
pval <- c(pval, numeric(m / 2))
kappai <- c(kappai, numeric(m / 2))
kappai_star <- c(kappai_star, numeric(m / 2))
Ci <- c(Ci, numeric(m / 2))
Si <- c(Si, numeric(m / 2))
Ci_plus <- c(Ci_plus, numeric(m / 2))
gammai <- c(gammai, numeric(m / 2))
alphai <- c(alphai, numeric(m / 2))
}
if(indepTest=="gaussCItest") #if indepTest=gaussCItest
{
pval[m] <- gaussCItest(x, z, y, suffStat)
}
if(indepTest=="disCItest") #if indepTest=disCItest
{
pval[m] <- disCItest(x, z, y, suffStat) #additional
}
if (FDRcontrol == 'LOND') {
# Calculate alpha using the LOND method.
alphai[m] <- SeqFDR(m,FDR,a=2,R)
Alpha <- alphai[m]
# Calculate the alpha value using the ADDIS algorithm
} else if (FDRcontrol == 'ADDIS') {
# Run ADDIS on the current iteration and update all the vectors
# and other values.
run_addis <- addis(alpha = FDR,
tau = tau,
lambda = lambda,
iter = m,
w0 = w0,
pval = pval,
alphai = alphai,
gammai = gammai,
kappai = kappai,
kappai_star = kappai_star,
K = K,
Ci = Ci,
Si = Si,
Ri = R,
Ci_plus = Ci_plus,
Ci_sum = Ci_sum,
Si_sum = Si_sum,
gammai_sum = gammai_sum,
normalizer = normalizer,
exponent = exponent)
# Update all values and vectors output from the addis function.
alphai[[m]] <- run_addis[[1]]
gammai[[m]] <- run_addis[[2]]
K <- run_addis[[3]]
R[[m]] <- run_addis[[5]]
Si[[m - 1]] <- run_addis[[6]]
Ci[[m - 1]] <- run_addis[[7]]
Ci_sum <- run_addis[[9]]
Si_sum <- run_addis[[10]]
gammai_sum <- run_addis[[12]]
# Only update the kappai and Ci_plus vectors if K is greater
# than one. If K is zero then the first element in kappai will
# remain zero until the first rejection.
if (K != 0) {
kappai[[K]] <- run_addis[[4]]
kappai_star[[K]] <- run_addis[[11]]
Ci_plus[1:K] <- run_addis[[8]]
}
# Update the Alpha value with the current alphai value.
Alpha <- alphai[[m]]
} else if (FDRcontrol == 'NONE') {
Alpha <- alpha
}
if (verbose) {
cat("Additional pval value =", pval[m], "\n")
cat("Alpha value =", Alpha, "\n")
}
#Alpha=SeqFDR(m,FDR,a=2,R) #Alpha valued from sequential FDR test
if (pval[m] <= Alpha) { #Reject H0 (H0:nodes are independent)
R[m] <- 1
if (verbose) {
V <- colnames(g)
cat("\n", V[x], "->", V[y], "<-", V[z], "\n") #Printout the v-structures
cat("Since pval<Alpha,additional test is rejected", "Nodes", V[x] ,"and" ,V[z] ,"are dependent given", V[y], "\n")
}
tarmat[z, y] <- 1 #directed z-->y
}
else {
tarmat[y, z] <- 1
R[m] <- 0 #Accept H0
}
}
#Case-4:If, y & z have relation and x & y conditionally dependent given z and x & z conditionally independent given y,
#then edge direction will be z<-->y
if (g1[y, z] == 1 & g1[x, z] == 1 & g1[x, y] == 1 & !(z %in% gInput@sepset[[x]][[y]]) &!(y %in% gInput@sepset[[x]][[z]]))
{
tarmat[y, z] <- 1
tarmat[z, y] <- 1
}
}
}
}
}
#if(length(WW3)!=0)
#{
#Orient remaining edges,
T1 <- which(tarmat==1,arr.ind = T)
G1 <- which(g==1,arr.ind = T)
D1 <- setdiff(G1,T1)
#
if(length(D1)!=0)
{
for (d1 in 1:length(D1))
{
Rem1 <- which(g[,D1[d1]]==1,arr.ind = T)
Rem2 <- which(tarmat[,Rem1]==1,arr.ind = T)
if(length(Rem1)!=0 & length(Rem2)!=0)
{
for (d11 in 1:length(Rem1))
{
x <- Rem2[d11]
y <- Rem1[d11]
z <- D1[d1]
}
if ((!is.na(x) & !is.na(y) & !is.na(z)) & (x!=0 & y!=0 & z!=0) & (x!=y & y!=z & z!=x))
{
#Case-2: If,y and z are adjacent, x and z conditionally independent given y,
#then the edge direction will be y-->z
if (g1[y, z] == 1 & tarmat[y, z]!=1 & tarmat[z, y]!=1 & ((y %in% gInput@sepset[[x]][[z]]) || (y %in% gInput@sepset[[z]][[x]])))
{
tarmat[y, z] <- 1
tarmat[z, y] <- 0
}
#Case-3: If,y and z are adjacent, x and z conditionally dependent given y,
#then the edge direction will be z-->y.
if (g1[y, z] == 1 & g1[x, z] != 1 & tarmat[x,y]==1 & tarmat[y, x]!=1 & tarmat[y, z]!=1 & tarmat[z, y]!=1 &!(y %in% gInput@sepset[[x]][[z]]) & !(y %in% gInput@sepset[[z]][[x]]))
{
m <- m+1
# Increase the length of the R, pval, and other addis vectors if
# their length is greater than or equal to the current iteration.
if (m >= length(R)) {
R <- c(R, numeric(m / 2))
pval <- c(pval, numeric(m / 2))
kappai <- c(kappai, numeric(m / 2))
kappai_star <- c(kappai_star, numeric(m / 2))
Ci <- c(Ci, numeric(m / 2))
Si <- c(Si, numeric(m / 2))
Ci_plus <- c(Ci_plus, numeric(m / 2))
gammai <- c(gammai, numeric(m / 2))
alphai <- c(alphai, numeric(m / 2))
}
if(indepTest=="gaussCItest") #if indepTest=gaussCItest
{
pval[m] <- gaussCItest(x, z, y, suffStat)
}
if(indepTest=="disCItest") #if indepTest=gaussCItest
{
pval[m] <- disCItest(x, z, y, suffStat) #additional
}
if (FDRcontrol == 'LOND') {
# Calculate alpha using the LOND method.
alphai[m] <- SeqFDR(m,FDR,a=2,R)
Alpha <- alphai[m]
# Calculate the alpha value using the ADDIS algorithm
} else if (FDRcontrol == 'ADDIS') {
# Run ADDIS on the current iteration and update all the vectors
# and other values.
run_addis <- addis(alpha = FDR,
tau = tau,
lambda = lambda,
iter = m,
w0 = w0,
pval = pval,
alphai = alphai,
gammai = gammai,
kappai = kappai,
kappai_star = kappai_star,
K = K,
Ci = Ci,
Si = Si,
Ri = R,
Ci_plus = Ci_plus,
Ci_sum = Ci_sum,
Si_sum = Si_sum,
gammai_sum = gammai_sum,
normalizer = normalizer,
exponent = exponent)
# Update all values and vectors output from the addis function.
alphai[[m]] <- run_addis[[1]]
gammai[[m]] <- run_addis[[2]]
K <- run_addis[[3]]
R[[m]] <- run_addis[[5]]
Si[[m - 1]] <- run_addis[[6]]
Ci[[m - 1]] <- run_addis[[7]]
Ci_sum <- run_addis[[9]]
Si_sum <- run_addis[[10]]
gammai_sum <- run_addis[[12]]
# Only update the kappai and Ci_plus vectors if K is greater
# than one. If K is zero then the first element in kappai will
# remain zero until the first rejection.
if (K != 0) {
kappai[[K]] <- run_addis[[4]]
kappai_star[[K]] <- run_addis[[11]]
Ci_plus[1:K] <- run_addis[[8]]
}
# Update the Alpha value with the current alphai value.
Alpha <- alphai[[m]]
} else if (FDRcontrol == 'NONE') {
Alpha <- alpha
}
if (pval[m] <= Alpha) { #Reject H0 (H0:nodes are independent)
R[m] <- 1
if (verbose) {
V <- colnames(g)
cat("\n", V[x], "->", V[y], "<-", V[z], "\n") #Printout the v-structures
cat("Since pval<Alpha,additional test is rejected", "Nodes", V[x] ,"and" ,V[z] ,"are dependent given", V[y], "\n")
}
tarmat[z, y] <- 1 #directed z-->y
}
else {
tarmat[y, z] <- 1
R[m] <- 0 #Accept H0
}
}
#Case-4:If, y & z have relation and x & y conditionally dependent given z and x & z conditionally independent given y,
#then edge direction will be z<-->y
if (g1[y, z] == 1 & g1[x, z] == 1 & g1[x, y] == 1 & !(z %in% gInput@sepset[[x]][[y]]) &!(y %in% gInput@sepset[[x]][[z]]))
{
tarmat[y, z] <- 1
tarmat[z, y] <- 1
}
}
}
}
}
#Check the reaming edges
ind <- which(g1 == 1, arr.ind = TRUE) #Pullout the all relation in adjacency matrix from undirected graph
for (i in seq_len(nrow(ind))) {
x <- ind[i, 1]
y <- ind[i, 2]
if(tarmat[x, y]==0 & tarmat[y, x]==0) #bidirected if still no edge in tarmat
{
tarmat[x, y] <- 1
tarmat[y, x] <- 1
}
}
}
}#end when data not necessary to contain genetic variants.
#Step 4
#Produce as a same skeleton if no edges involves with genetic variants and no v-structures
if(all(tarmat==0))
{
tarmat <- g1
}
#Step 5
#If found the direction with only genetic variants, no v-structures and
#all the other nodes still undirected, then remaining edges will be bidirected
if(GV>0 & (any(tarmat[1:GV,]==1) || any(tarmat[,1:GV]==1)) & all(tarmat[-c(1:GV),-c(1:GV)]==0))
{
tarmat1 <- g1
tarmat1[1:GV,] <- tarmat[1:GV,]
tarmat1[,1:GV] <- tarmat[,1:GV]
tarmat <- tarmat1
}
gInput@graph <- as(tarmat, "graphNEL")
gInput@R <- R[1:m]
gInput@K <- K
gInput@pval <- pval[1:m]
gInput@kappai <- kappai[1:K]
gInput@kappai_star <- kappai_star[1:K]
gInput@Ci <- Ci[1:m]
gInput@Si <- Si[1:m]
gInput@Ci_plus <- Ci_plus[1:K]
gInput@gammai <- gammai[1:m]
gInput@gammai_sum <- gammai_sum
gInput
}
audreyqyfu/mrpc documentation built on April 17, 2022, 7:35 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#This is the step 2 of MRPC to direction determination of the undirected edges.
EdgeOrientation <- function (gInput, GV, suffStat, FDR, alpha, indepTest,
FDRcontrol = c("LOND", "ADDIS", "NONE"), tau = 0.5, lambda = 0.25,
verbose = FALSE) {
g <- as(gInput@graph, "matrix") # g ia an adjacency from undirected graph (skleton)
g1 <- g
p <- nrow(g)
# tarmat (adjacency matrix for directed graph) is updated in every step,
# and contains the output of final topology
tarmat <- matrix(0, nrow(g),ncol(g)) #same row and column from g
rownames(tarmat) <- rownames(g) #same row names from g
colnames(tarmat) <- colnames(g) #same column names from g
# extract all edges
g[lower.tri(g)] <- 0 # Use only upper triangular because g is symmetric matrix
edges <- which (g==1, arr.ind = TRUE)
#Step-1 start
if (GV>0) {
# identify edges involving Vs
edgesWithBothVs <- edges[which (edges[,1] <= GV & edges[,2] <= GV), ]
edgesWithVFirst <- edges[which (edges[,1] <= GV & edges[,2] > GV), ]
edgesWithVSecond <- edges[which (edges[,1] > GV & edges[,2] <= GV), ]
# assign 1s to corresponding edges in tarmat
# edges between two Vs are undirected (or bidirected)
if (length (edgesWithBothVs) > 2) {
tarmat[edgesWithBothVs] <- 1
tarmat[edgesWithBothVs[,2:1]] <- 1
} else {
tarmat[edgesWithBothVs[1], edgesWithBothVs[2]] <- 1
tarmat[edgesWithBothVs[2], edgesWithBothVs[1]] <- 1
}
# Edges involving one V go from V to the other node
# Need to distinguish whether edgesWithVFirst or edgesWithVSecond is
# a matrix or a vector.
if (length (edgesWithVFirst) > 2) {
tarmat[edgesWithVFirst] <- 1
} else {
tarmat[edgesWithVFirst[1], edgesWithVFirst[2]] <- 1
}
if (length(edgesWithVSecond) > 2) {
tarmat[edgesWithVSecond[,2:1]] <- 1
} else {
tarmat[edgesWithVSecond[2], edgesWithVSecond[1]] <- 1
}
}
#Step-1 end
#Step-2 start
#Start to orient v-structures
if (verbose)
cat("\n V-structures are as follows :\n")
# extract edges involving at least one gene node
#edgesWithGs <- edges[which (edges[,1]>GV | edges[,2]>GV), ]
#ind <- edgesWithGs
m <- gInput@test #current test
Alpha <- gInput@alpha #alpha
pval <- gInput@pval
R <- gInput@R #decision of test
ind <- which(g1 == 1, arr.ind = TRUE) #Pullout the all relation in adjacency matrix from undirected graph
V <- colnames(g)
# Update the addis vectors/values output from the ModiSkeleton function.
w0 <- tau * lambda * FDR/2
alphai <- gInput@alphai
kappai <- gInput@kappai
kappai_star <- gInput@kappai_star
K <- gInput@K
Ci <- gInput@Ci
Ci_sum <- sum(Ci)
Si <- gInput@Si
Si_sum <- sum(Si)
gammai_sum <- gInput@gammai_sum
Ci_plus <- gInput@Ci_plus
gammai <- gInput@gammai
normalizer <- gInput@normalizer
exponent <- gInput@exponent
for (i in seq_len(nrow(ind))) {
x <- ind[i, 1]
y <- ind[i, 2]
allZ <- setdiff(which(g1[y, ] == 1), x)
for (z in allZ) {
# Triplet x-y-z is directed x-->y<--z if x and z conditionally dependent given y
if ((g1[x, z] == 0 & g1[x, y] == 1 & g1[y, z] == 1) & !(tarmat[y, x] ==1) & !(tarmat[z, y] ==1) & !(tarmat[y, z] ==1) &
!(y %in% gInput@sepset[[x]][[z]] || y %in% gInput@sepset[[z]][[x]]))
{
m <- m+1
# Increase the length of the R, pval, and other addis vectors if
# their length is greater than or equal to the current iteration.
if (m >= length(R)) {
R <- c(R, numeric(m / 2))
pval <- c(pval, numeric(m / 2))
kappai <- c(kappai, numeric(m / 2))
kappai_star <- c(kappai_star, numeric(m / 2))
Ci <- c(Ci, numeric(m / 2))
Si <- c(Si, numeric(m / 2))
Ci_plus <- c(Ci_plus, numeric(m / 2))
gammai <- c(gammai, numeric(m / 2))
alphai <- c(alphai, numeric(m / 2))
}
if(indepTest=="gaussCItest") #if indepTest=gaussCItest
{
pval[m] <- gaussCItest(x, z, y, suffStat)
}
if(indepTest=="disCItest") #if indepTest=disCItest
{
pval[m] <- disCItest(x, z, y, suffStat) #additional
}
#pval=disCItest(x, z, y, suffStat) #additional conditional test
if (FDRcontrol == 'LOND') {
# Calculate alpha using the LOND method.
alphai[m] <- SeqFDR(m,FDR,a=2,R)
Alpha <- alphai[m]
# Calculate the alpha value using the ADDIS algorithm
} else if (FDRcontrol == 'ADDIS') {
# Run ADDIS on the current iteration and update all the vectors
# and other values.
run_addis <- addis(alpha = FDR,
tau = tau,
lambda = lambda,
iter = m,
w0 = w0,
pval = pval,
alphai = alphai,
gammai = gammai,
kappai = kappai,
kappai_star = kappai_star,
K = K,
Ci = Ci,
Si = Si,
Ri = R,
Ci_plus = Ci_plus,
Ci_sum = Ci_sum,
Si_sum = Si_sum,
gammai_sum = gammai_sum,
normalizer = normalizer,
exponent = exponent)
# Update all values and vectors output from the addis function.
alphai[[m]] <- run_addis[[1]]
gammai[[m]] <- run_addis[[2]]
K <- run_addis[[3]]
R[[m]] <- run_addis[[5]]
Si[[m - 1]] <- run_addis[[6]]
Ci[[m - 1]] <- run_addis[[7]]
Ci_sum <- run_addis[[9]]
Si_sum <- run_addis[[10]]
gammai_sum <- run_addis[[12]]
# Only update the kappai and Ci_plus vectors if K is greater
# than one. If K is zero then the first element in kappai will
# remain zero until the first rejection.
if (K != 0) {
kappai[[K]] <- run_addis[[4]]
kappai_star[[K]] <- run_addis[[11]]
Ci_plus[1:K] <- run_addis[[8]]
}
# Update the Alpha value with the current alphai value.
Alpha <- alphai[[m]]
} else if (FDRcontrol == 'NONE') {
Alpha <- alpha
}
if (verbose){
cat("x=", x, " y=", z, " S=", y,"\n")
cat("Test number =", m, "\n")
cat("Additional pval value =", pval[m], "\n")
cat("Alpha value =", Alpha, "\n")
}
if (pval[m] <= Alpha) { #Reject H0 (H0:nodes are independent)
R[m] <- 1
if (verbose) {
cat(V[x], "->", V[y], "<-", V[z], "\n") #Printout the v-structures
cat("Since pval<Alpha,additional test is rejected;", "Nodes", V[x] ,"and" ,V[z] ,"are dependent given", V[y], "\n")
}
tarmat[x, y] <- tarmat[z, y] <- 1 #directed x-->y<--z
}
else {
R[m] <- 0 #Accept H0
#tarmat[x, y] <- tarmat[y, z] <- 1 #directed x-->y-->z
if (verbose)
cat("Since pval>Alpha,additional test is accepted;", "Nodes", V[x] ,"and" ,V[z] ,"are independent given", V[y], "\n")
}
}
}
}
#Step-2 end to orient v-structures
#Step-3
#Orient remaining edges, weather gMR is applicable or not based on the following process:
#Pullout all edges from g1 matrix and define 1st and 2nd node (goal to make a triplet for the remaining edges)
#extract all edges with involving those 1st and 2nd nodes (from row and column of g)
#ignore the nodes that already have direction
#form a triplet
#determine the direction using test results from part 1 based on gMR
#Repeat until all undirected edges to directed
# m <- m
# Alpha <- Alpha
# R <- R
#start
#when data contain genetic variants
if (any(tarmat == 1) & GV>0) #if at least one edge directed already made so far
{
WW1 <- unique(which(tarmat==1,arr.ind = T)[,2]) #pullout the all canidate genes for v-structure
#if(WW1<GV){
#WW1=WW1[-c(1:GV)] #ignor canidate genes for GV
WW1 <- setdiff(WW1,1:GV)
#}
if(length(WW1)!=0)
{
#WW1=setdiff(WW1,1:GV)
for (v1 in 1:length(WW1))
{
WW2 <- unique(which(tarmat[,WW1[v1]]==1,arr.ind = T)) #edge between v-structure
WW3 <- unique(which(g1[,WW1[v1]]==1,arr.ind = T)) #edges between canidate genes and others
WW3 <- setdiff(WW3,WW2) #ignore if already direction
if(length(WW3)!=0)
{
for (v2 in 1:length(WW3))
{
x <- WW2[1]
y <- WW1[v1]
z <- WW3[v2]
if ((!is.na(x) & !is.na(y) & !is.na(z)) & (x!=0 & y!=0 & z!=0) & (x!=y & y!=z & z!=x))
{
#Case-2: If,y and z are adjacent, x and z conditionally independent given y,
#then the edge direction will be y-->z
if (g1[y, z] == 1 & tarmat[y, z]!=1 & tarmat[z, y]!=1 & ((y %in% gInput@sepset[[x]][[z]]) || (y %in% gInput@sepset[[z]][[x]])))
{
tarmat[y, z] <- 1
tarmat[z, y] <- 0
}
#Case-3: If,y and z are adjacent, x and z conditionally dependent given y,
#then the edge direction will be z-->y.
if (g1[y, z] == 1 & g1[x, z] != 1 & tarmat[x,y]==1 & tarmat[y, x]!=1 & tarmat[y, z]!=1 & tarmat[z, y]!=1 &!(y %in% gInput@sepset[[x]][[z]]) & !(y %in% gInput@sepset[[z]][[x]]))
{
m <- m+1
# Increase the length of the R, pval, and other addis vectors if
# their length is greater than or equal to the current iteration.
if (m >= length(R)) {
R <- c(R, numeric(m / 2))
pval <- c(pval, numeric(m / 2))
kappai <- c(kappai, numeric(m / 2))
kappai_star <- c(kappai_star, numeric(m / 2))
Ci <- c(Ci, numeric(m / 2))
Si <- c(Si, numeric(m / 2))
Ci_plus <- c(Ci_plus, numeric(m / 2))
gammai <- c(gammai, numeric(m / 2))
alphai <- c(alphai, numeric(m / 2))
}
if(indepTest=="gaussCItest") #if indepTest=gaussCItest for continuous data
{
pval[m] <- gaussCItest(x, z, y, suffStat) #additional pval
}
if(indepTest=="disCItest") #if indepTest=disCItest for discrete data
{
pval[m] <- disCItest(x, z, y, suffStat) #additional pval
}
if (FDRcontrol == 'LOND') {
# Calculate alpha using the LOND method.
alphai[m] <- SeqFDR(m,FDR,a=2,R)
Alpha <- alphai[m]
# Calculate the alpha value using the ADDIS algorithm
} else if (FDRcontrol == 'ADDIS') {
# Run ADDIS on the current iteration and update all the vectors
# and other values.
run_addis <- addis(alpha = FDR,
tau = tau,
lambda = lambda,
iter = m,
w0 = w0,
pval = pval,
alphai = alphai,
gammai = gammai,
kappai = kappai,
kappai_star = kappai_star,
K = K,
Ci = Ci,
Si = Si,
Ri = R,
Ci_plus = Ci_plus,
Ci_sum = Ci_sum,
Si_sum = Si_sum,
gammai_sum = gammai_sum,
normalizer = normalizer,
exponent = exponent)
# Update all values and vectors output from the addis function.
alphai[[m]] <- run_addis[[1]]
gammai[[m]] <- run_addis[[2]]
K <- run_addis[[3]]
R[[m]] <- run_addis[[5]]
Si[[m - 1]] <- run_addis[[6]]
Ci[[m - 1]] <- run_addis[[7]]
Ci_sum <- run_addis[[9]]
Si_sum <- run_addis[[10]]
gammai_sum <- run_addis[[12]]
# Only update the kappai and Ci_plus vectors if K is greater
# than one. If K is zero then the first element in kappai will
# remain zero until the first rejection.
if (K != 0) {
kappai[[K]] <- run_addis[[4]]
kappai_star[[K]] <- run_addis[[11]]
Ci_plus[1:K] <- run_addis[[8]]
}
# Update the Alpha value with the current alphai value.
Alpha <- alphai[[m]]
} else if (FDRcontrol == 'NONE') {
Alpha <- alpha
}
if (verbose){
cat("Additional pval value =", pval[m], "\n")
cat("Alpha value =", Alpha, "\n")
}
if (pval[m] <= Alpha) { #Reject H0 (H0:nodes are independent)
R[m] <- 1
if (verbose) {
V <- colnames(g)
cat("\n", V[x], "->", V[y], "<-", V[z], "\n") #Printout the v-structures
cat("Since pval<Alpha,additional test is rejected", "Nodes", V[x] ,"and" ,V[z] ,"are dependent given", V[y], "\n")
}
tarmat[z, y] <- 1 #directed z-->y
}
else {
tarmat[y, z] <- 1
R[m] <- 0 #Accept H0
}
}
#Case-4:If, y & z have relation and x & y conditionally dependent given z and x & z conditionally independent given y,
#then edge direction will be z<-->y
if (g1[y, z] == 1 & g1[x, z] == 1 & g1[x, y] == 1 & !(z %in% gInput@sepset[[x]][[y]]) &!(y %in% gInput@sepset[[x]][[z]]))
{
tarmat[y, z] <- 1
tarmat[z, y] <- 1
}
}
}
}
}
#if(length(WW3)!=0)
#{
#Orient remaining edges,
T1 <- which(tarmat==1,arr.ind = T)
G1 <- which(g==1,arr.ind = T)
D1 <- setdiff(G1,T1)
#
if(length(D1)!=0)
{
for (d1 in 1:length(D1))
{
Rem1_row <- which(g[,D1[d1]]==1,arr.ind = T)
Rem1_col <- which(g[D1[d1],]==1,arr.ind = T)
Rem1 <- c(Rem1_row,Rem1_col)
Rem2 <- which(tarmat[,Rem1]==1,arr.ind = T)
if(length(Rem1)!=0 & length(Rem2)!=0)
{
for (d11 in 1:length(Rem1))
{
x <- Rem2[d11]
y <- Rem1[d11]
z <- D1[d1]
}
if ((!is.na(x) & !is.na(y) & !is.na(z)) & (x!=0 & y!=0 & z!=0) & (x!=y & y!=z & z!=x))
{
#Case-2: If,y and z are adjacent, x and z conditionally independent given y,
#then the edge direction will be y-->z
if (g1[y, z] == 1 & tarmat[y, z]!=1 & tarmat[z, y]!=1 & ((y %in% gInput@sepset[[x]][[z]]) || (y %in% gInput@sepset[[z]][[x]])))
{
tarmat[y, z] <- 1
tarmat[z, y] <- 0
}
#Case-3: If,y and z are adjacent, x and z conditionally dependent given y,
#then the edge direction will be z-->y.
if (g1[y, z] == 1 & g1[x, z] != 1 & tarmat[x,y]==1 & tarmat[y, x]!=1 & tarmat[y, z]!=1 & tarmat[z, y]!=1 &!(y %in% gInput@sepset[[x]][[z]]) & !(y %in% gInput@sepset[[z]][[x]]))
{
m <- m+1
# Increase the length of the R, pval, and other addis vectors if
# their length is greater than or equal to the current iteration.
if (m >= length(R)) {
R <- c(R, numeric(m / 2))
pval <- c(pval, numeric(m / 2))
kappai <- c(kappai, numeric(m / 2))
kappai_star <- c(kappai_star, numeric(m / 2))
Ci <- c(Ci, numeric(m / 2))
Si <- c(Si, numeric(m / 2))
Ci_plus <- c(Ci_plus, numeric(m / 2))
gammai <- c(gammai, numeric(m / 2))
alphai <- c(alphai, numeric(m / 2))
}
if(indepTest=="gaussCItest") #if indepTest=gaussCItest
{
pval[m] <- gaussCItest(x, z, y, suffStat)
}
if(indepTest=="disCItest") #if indepTest=gaussCItest
{
pval[m] <- disCItest(x, z, y, suffStat) #additional
}
if (FDRcontrol == 'LOND') {
# Calculate alpha using the LOND method.
alphai[m] <- SeqFDR(m,FDR,a=2,R)
Alpha <- alphai[m]
# Calculate the alpha value using the ADDIS algorithm
} else if (FDRcontrol == 'ADDIS') {
# Run ADDIS on the current iteration and update all the vectors
# and other values.
run_addis <- addis(alpha = FDR,
tau = tau,
lambda = lambda,
iter = m,
w0 = w0,
pval = pval,
alphai = alphai,
gammai = gammai,
kappai = kappai,
kappai_star = kappai_star,
K = K,
Ci = Ci,
Si = Si,
Ri = R,
Ci_plus = Ci_plus,
Ci_sum = Ci_sum,
Si_sum = Si_sum,
gammai_sum = gammai_sum,
normalizer = normalizer,
exponent = exponent)
# Update all values and vectors output from the addis function.
alphai[[m]] <- run_addis[[1]]
gammai[[m]] <- run_addis[[2]]
K <- run_addis[[3]]
R[[m]] <- run_addis[[5]]
Si[[m - 1]] <- run_addis[[6]]
Ci[[m - 1]] <- run_addis[[7]]
Ci_sum <- run_addis[[9]]
Si_sum <- run_addis[[10]]
gammai_sum <- run_addis[[12]]
# Only update the kappai and Ci_plus vectors if K is greater
# than one. If K is zero then the first element in kappai will
# remain zero until the first rejection.
if (K != 0) {
kappai[[K]] <- run_addis[[4]]
kappai_star[[K]] <- run_addis[[11]]
Ci_plus[1:K] <- run_addis[[8]]
}
# Update the Alpha value with the current alphai value.
Alpha <- alphai[[m]]
} else if (FDRcontrol == 'NONE') {
Alpha <- alpha
}
if (pval[m] <= Alpha) { #Reject H0 (H0:nodes are independent)
R[m] <- 1
if (verbose) {
V <- colnames(g)
cat("\n", V[x], "->", V[y], "<-", V[z], "\n") #Printout the v-structures
cat("Since pval<Alpha,additional test is rejected", "Nodes", V[x] ,"and" ,V[z] ,"are dependent given", V[y], "\n")
}
tarmat[z, y] <- 1 #directed z-->y
}
else {
tarmat[y, z] <- 1
R[m] <- 0 #Accept H0
}
}
#Case-4:If, y & z have relation and x & y conditionally dependent given z and x & z conditionally independent given y,
#then edge direction will be z<-->y
if (g1[y, z] == 1 & g1[x, z] == 1 & g1[x, y] == 1 & !(z %in% gInput@sepset[[x]][[y]]) &!(y %in% gInput@sepset[[x]][[z]]))
{
tarmat[y, z] <- 1
tarmat[z, y] <- 1
}
}
}
}
}
#Check the reamning edge orientation
ind <- which(g1 == 1, arr.ind = TRUE) #Pullout the all relation in adjacency matrix from undirected graph
for (i in seq_len(nrow(ind))) {
x <- ind[i, 1]
y <- ind[i, 2]
if(tarmat[x, y]==0 & tarmat[y, x]==0) #bidirected if still no edge in tarmat
{
tarmat[x, y] <- 1
tarmat[y, x] <- 1
}
}
}
}#end when data contain genetic variants
#Start to orient remaining edges, when data not necessary to contain genetic variants.
#First identify all the canidate genes for v-structure
#Then identify the edges (undirected) with canidate genes
#ignore the nodes that already have direction
#Make a triplet
#determine the direction using test results from part 1 based on gMR
#Repeat until all undirected edges to directed
#Start when data not necessary to contain genetic variants.
if (any(tarmat == 1) & GV==0) #if at least one edge directed already made so far
{
WW1 <- unique(which(tarmat == 1,arr.ind = T)[,2]) #pullout the all canidate genes for v-structure
#WW1=WW1[-c(1:GV)]
if(length(WW1)!=0)
{
#WW1=setdiff(WW1,1:GV)
for (v1 in 1:length(WW1))
{
WW2 <- unique(which(tarmat[,WW1[v1]]==1,arr.ind = T)) #edge between v-structure
WW3 <- unique(which(g1[,WW1[v1]]==1,arr.ind = T)) #edges between canidate genes and others
WW3 <- setdiff(WW3,WW2) #ignore if already direction
if(length(WW3)!=0)
{
for (v2 in 1:length(WW3))
{
x <- WW2[1]
y <- WW1[v1]
z <- WW3[v2]
if ((!is.na(x) & !is.na(y) & !is.na(z)) & (x!=0 & y!=0 & z!=0) & (x!=y & y!=z & z!=x))
{
#Case-2: If,y and z are adjacent, x and z conditionally independent given y,
#then the edge direction will be y-->z
if (g1[y, z] == 1 & tarmat[y, z]!=1 & tarmat[z, y]!=1 & ((y %in% gInput@sepset[[x]][[z]]) || (y %in% gInput@sepset[[z]][[x]])))
{
tarmat[y, z] <- 1
tarmat[z, y] <- 0
}
#Case-3: If,y and z are adjacent, x and z conditionally dependent given y,
#then the edge direction will be z-->y.
if (g1[y, z] == 1 & g1[x, z] != 1 & tarmat[x,y]==1 & tarmat[y, x]!=1 & tarmat[y, z]!=1 & tarmat[z, y]!=1 &!(y %in% gInput@sepset[[x]][[z]]) & !(y %in% gInput@sepset[[z]][[x]]))
{
m <- m+1
# Increase the length of the R, pval, and other addis vectors if
# their length is greater than or equal to the current iteration.
if (m >= length(R)) {
R <- c(R, numeric(m / 2))
pval <- c(pval, numeric(m / 2))
kappai <- c(kappai, numeric(m / 2))
kappai_star <- c(kappai_star, numeric(m / 2))
Ci <- c(Ci, numeric(m / 2))
Si <- c(Si, numeric(m / 2))
Ci_plus <- c(Ci_plus, numeric(m / 2))
gammai <- c(gammai, numeric(m / 2))
alphai <- c(alphai, numeric(m / 2))
}
if(indepTest=="gaussCItest") #if indepTest=gaussCItest
{
pval[m] <- gaussCItest(x, z, y, suffStat)
}
if(indepTest=="disCItest") #if indepTest=disCItest
{
pval[m] <- disCItest(x, z, y, suffStat) #additional
}
if (FDRcontrol == 'LOND') {
# Calculate alpha using the LOND method.
alphai[m] <- SeqFDR(m,FDR,a=2,R)
Alpha <- alphai[m]
# Calculate the alpha value using the ADDIS algorithm
} else if (FDRcontrol == 'ADDIS') {
# Run ADDIS on the current iteration and update all the vectors
# and other values.
run_addis <- addis(alpha = FDR,
tau = tau,
lambda = lambda,
iter = m,
w0 = w0,
pval = pval,
alphai = alphai,
gammai = gammai,
kappai = kappai,
kappai_star = kappai_star,
K = K,
Ci = Ci,
Si = Si,
Ri = R,
Ci_plus = Ci_plus,
Ci_sum = Ci_sum,
Si_sum = Si_sum,
gammai_sum = gammai_sum,
normalizer = normalizer,
exponent = exponent)
# Update all values and vectors output from the addis function.
alphai[[m]] <- run_addis[[1]]
gammai[[m]] <- run_addis[[2]]
K <- run_addis[[3]]
R[[m]] <- run_addis[[5]]
Si[[m - 1]] <- run_addis[[6]]
Ci[[m - 1]] <- run_addis[[7]]
Ci_sum <- run_addis[[9]]
Si_sum <- run_addis[[10]]
gammai_sum <- run_addis[[12]]
# Only update the kappai and Ci_plus vectors if K is greater
# than one. If K is zero then the first element in kappai will
# remain zero until the first rejection.
if (K != 0) {
kappai[[K]] <- run_addis[[4]]
kappai_star[[K]] <- run_addis[[11]]
Ci_plus[1:K] <- run_addis[[8]]
}
# Update the Alpha value with the current alphai value.
Alpha <- alphai[[m]]
} else if (FDRcontrol == 'NONE') {
Alpha <- alpha
}
if (verbose) {
cat("Additional pval value =", pval[m], "\n")
cat("Alpha value =", Alpha, "\n")
}
#Alpha=SeqFDR(m,FDR,a=2,R) #Alpha valued from sequential FDR test
if (pval[m] <= Alpha) { #Reject H0 (H0:nodes are independent)
R[m] <- 1
if (verbose) {
V <- colnames(g)
cat("\n", V[x], "->", V[y], "<-", V[z], "\n") #Printout the v-structures
cat("Since pval<Alpha,additional test is rejected", "Nodes", V[x] ,"and" ,V[z] ,"are dependent given", V[y], "\n")
}
tarmat[z, y] <- 1 #directed z-->y
}
else {
tarmat[y, z] <- 1
R[m] <- 0 #Accept H0
}
}
#Case-4:If, y & z have relation and x & y conditionally dependent given z and x & z conditionally independent given y,
#then edge direction will be z<-->y
if (g1[y, z] == 1 & g1[x, z] == 1 & g1[x, y] == 1 & !(z %in% gInput@sepset[[x]][[y]]) &!(y %in% gInput@sepset[[x]][[z]]))
{
tarmat[y, z] <- 1
tarmat[z, y] <- 1
}
}
}
}
}
#if(length(WW3)!=0)
#{
#Orient remaining edges,
T1 <- which(tarmat==1,arr.ind = T)
G1 <- which(g==1,arr.ind = T)
D1 <- setdiff(G1,T1)
#
if(length(D1)!=0)
{
for (d1 in 1:length(D1))
{
Rem1 <- which(g[,D1[d1]]==1,arr.ind = T)
Rem2 <- which(tarmat[,Rem1]==1,arr.ind = T)
if(length(Rem1)!=0 & length(Rem2)!=0)
{
for (d11 in 1:length(Rem1))
{
x <- Rem2[d11]
y <- Rem1[d11]
z <- D1[d1]
}
if ((!is.na(x) & !is.na(y) & !is.na(z)) & (x!=0 & y!=0 & z!=0) & (x!=y & y!=z & z!=x))
{
#Case-2: If,y and z are adjacent, x and z conditionally independent given y,
#then the edge direction will be y-->z
if (g1[y, z] == 1 & tarmat[y, z]!=1 & tarmat[z, y]!=1 & ((y %in% gInput@sepset[[x]][[z]]) || (y %in% gInput@sepset[[z]][[x]])))
{
tarmat[y, z] <- 1
tarmat[z, y] <- 0
}
#Case-3: If,y and z are adjacent, x and z conditionally dependent given y,
#then the edge direction will be z-->y.
if (g1[y, z] == 1 & g1[x, z] != 1 & tarmat[x,y]==1 & tarmat[y, x]!=1 & tarmat[y, z]!=1 & tarmat[z, y]!=1 &!(y %in% gInput@sepset[[x]][[z]]) & !(y %in% gInput@sepset[[z]][[x]]))
{
m <- m+1
# Increase the length of the R, pval, and other addis vectors if
# their length is greater than or equal to the current iteration.
if (m >= length(R)) {
R <- c(R, numeric(m / 2))
pval <- c(pval, numeric(m / 2))
kappai <- c(kappai, numeric(m / 2))
kappai_star <- c(kappai_star, numeric(m / 2))
Ci <- c(Ci, numeric(m / 2))
Si <- c(Si, numeric(m / 2))
Ci_plus <- c(Ci_plus, numeric(m / 2))
gammai <- c(gammai, numeric(m / 2))
alphai <- c(alphai, numeric(m / 2))
}
if(indepTest=="gaussCItest") #if indepTest=gaussCItest
{
pval[m] <- gaussCItest(x, z, y, suffStat)
}
if(indepTest=="disCItest") #if indepTest=gaussCItest
{
pval[m] <- disCItest(x, z, y, suffStat) #additional
}
if (FDRcontrol == 'LOND') {
# Calculate alpha using the LOND method.
alphai[m] <- SeqFDR(m,FDR,a=2,R)
Alpha <- alphai[m]
# Calculate the alpha value using the ADDIS algorithm
} else if (FDRcontrol == 'ADDIS') {
# Run ADDIS on the current iteration and update all the vectors
# and other values.
run_addis <- addis(alpha = FDR,
tau = tau,
lambda = lambda,
iter = m,
w0 = w0,
pval = pval,
alphai = alphai,
gammai = gammai,
kappai = kappai,
kappai_star = kappai_star,
K = K,
Ci = Ci,
Si = Si,
Ri = R,
Ci_plus = Ci_plus,
Ci_sum = Ci_sum,
Si_sum = Si_sum,
gammai_sum = gammai_sum,
normalizer = normalizer,
exponent = exponent)
# Update all values and vectors output from the addis function.
alphai[[m]] <- run_addis[[1]]
gammai[[m]] <- run_addis[[2]]
K <- run_addis[[3]]
R[[m]] <- run_addis[[5]]
Si[[m - 1]] <- run_addis[[6]]
Ci[[m - 1]] <- run_addis[[7]]
Ci_sum <- run_addis[[9]]
Si_sum <- run_addis[[10]]
gammai_sum <- run_addis[[12]]
# Only update the kappai and Ci_plus vectors if K is greater
# than one. If K is zero then the first element in kappai will
# remain zero until the first rejection.
if (K != 0) {
kappai[[K]] <- run_addis[[4]]
kappai_star[[K]] <- run_addis[[11]]
Ci_plus[1:K] <- run_addis[[8]]
}
# Update the Alpha value with the current alphai value.
Alpha <- alphai[[m]]
} else if (FDRcontrol == 'NONE') {
Alpha <- alpha
}
if (pval[m] <= Alpha) { #Reject H0 (H0:nodes are independent)
R[m] <- 1
if (verbose) {
V <- colnames(g)
cat("\n", V[x], "->", V[y], "<-", V[z], "\n") #Printout the v-structures
cat("Since pval<Alpha,additional test is rejected", "Nodes", V[x] ,"and" ,V[z] ,"are dependent given", V[y], "\n")
}
tarmat[z, y] <- 1 #directed z-->y
}
else {
tarmat[y, z] <- 1
R[m] <- 0 #Accept H0
}
}
#Case-4:If, y & z have relation and x & y conditionally dependent given z and x & z conditionally independent given y,
#then edge direction will be z<-->y
if (g1[y, z] == 1 & g1[x, z] == 1 & g1[x, y] == 1 & !(z %in% gInput@sepset[[x]][[y]]) &!(y %in% gInput@sepset[[x]][[z]]))
{
tarmat[y, z] <- 1
tarmat[z, y] <- 1
}
}
}
}
}
#Check the reaming edges
ind <- which(g1 == 1, arr.ind = TRUE) #Pullout the all relation in adjacency matrix from undirected graph
for (i in seq_len(nrow(ind))) {
x <- ind[i, 1]
y <- ind[i, 2]
if(tarmat[x, y]==0 & tarmat[y, x]==0) #bidirected if still no edge in tarmat
{
tarmat[x, y] <- 1
tarmat[y, x] <- 1
}
}
}
}#end when data not necessary to contain genetic variants.
#Step 4
#Produce as a same skeleton if no edges involves with genetic variants and no v-structures
if(all(tarmat==0))
{
tarmat <- g1
}
#Step 5
#If found the direction with only genetic variants, no v-structures and
#all the other nodes still undirected, then remaining edges will be bidirected
if(GV>0 & (any(tarmat[1:GV,]==1) || any(tarmat[,1:GV]==1)) & all(tarmat[-c(1:GV),-c(1:GV)]==0))
{
tarmat1 <- g1
tarmat1[1:GV,] <- tarmat[1:GV,]
tarmat1[,1:GV] <- tarmat[,1:GV]
tarmat <- tarmat1
}
gInput@graph <- as(tarmat, "graphNEL")
gInput@R <- R[1:m]
gInput@K <- K
gInput@pval <- pval[1:m]
gInput@kappai <- kappai[1:K]
gInput@kappai_star <- kappai_star[1:K]
gInput@Ci <- Ci[1:m]
gInput@Si <- Si[1:m]
gInput@Ci_plus <- Ci_plus[1:K]
gInput@gammai <- gammai[1:m]
gInput@gammai_sum <- gammai_sum
gInput
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.