R Under development (unstable) (2024-09-06 r87103 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > if(requireNamespace("dagitty")) { + library(pcalg) + (doExtras <- pcalg:::doExtras()) + + ## Minimalistic CRAN checks + + ## Test 1 ############################ + ## Test that "no adjustment set" and "empty adjustment set" are distinguished properly + x <- 1; y <- 2 + cpdag <- matrix(c(0,1,1,0),2,2) ## 1 --- 2 => no adj set + dag <- matrix(c(0,1,0,0),2,2) ## 1 --> 2 => empty adj set + + adjC <- adjustment(amat = cpdag, amat.type = "cpdag", x = 1, y = 2, set.type = "canonical") + adjD <- adjustment(amat = dag, amat.type = "dag", x = 1, y = 2, set.type = "canonical") + adjP <- adjustment(amat = dag, amat.type = "pdag", x = 1, y = 2, set.type = "canonical") + + stopifnot(!identical(adjC, adjD), identical(adjD, adjP) ) + + ## Test 2 ############################### + gacVSadj <- function(amat, x, y ,z, V, type) { + ## gac(z) is TRUE IFF z is returned by adjustment() + ## x,y,z: col positions as used in GAC + ## Result: TRUE is result is equal + typeDG <- switch(type, + dag = "dag", + cpdag = "cpdag", + mag = "mag", + pag = "pag") + gacRes <- gac(amat,x,y, z, type)$gac + adjRes <- adjustment(amat = amat, amat.type = typeDG, x = x, y = y, set.type = "all") + if (gacRes) { ## z is valid adj set + res <- any(sapply(adjRes, function(xx) setequal(z, xx))) + } else { ## z is not valid adj set + res <- all(!sapply(adjRes, function(xx) setequal(z, xx))) + } + res + } + + xx <- TRUE + + ## CPDAG 1: Paper Fig 1 + mFig1 <- matrix(c(0,1,1,0,0,0, 1,0,1,1,1,0, 0,0,0,0,0,1, + 0,1,1,0,1,1, 0,1,0,1,0,1, 0,0,0,0,0,0), 6,6) + type <- "cpdag" + x <- 3; y <- 6 + + V <- as.character(1:ncol(mFig1)) + rownames(mFig1) <- colnames(mFig1) <- V + + xx <- xx & gacVSadj(mFig1,x,y, z=c(2,4), V=V, type) + xx <- xx & gacVSadj(mFig1,x,y, z=c(4,5), V=V, type) + + type <- "pag" + mFig3a <- matrix(c(0,1,0,0, 1,0,1,1, 0,1,0,1, 0,1,1,0), 4,4) + V <- as.character(1:ncol(mFig3a)) + rownames(mFig3a) <- colnames(mFig3a) <- V + xx <- xx & gacVSadj(mFig3a, x=2, y=4, z=NULL, V=V, type) + + ## DAG 1 from Marloes' Talk + mMMd1 <- matrix(c(0,1,0,1,0,0, 0,0,1,0,1,0, 0,0,0,0,0,1, + 0,0,0,0,0,0, 0,0,0,0,0,0, 0,0,0,0,0,0),6,6) + V <- as.character(1:ncol(mMMd1)) + rownames(mMMd1) <- colnames(mMMd1) <- V + + type <- "dag" + x <- 1; y <- 3 + xx <- xx & gacVSadj(mMMd1, x,y, z=NULL, V=V, type) + xx <- xx & gacVSadj(mMMd1, x,y, z= 2, V=V, type) + + if (!xx) { + stop("Problem when testing function gacVSadj.") + } else { + message("OK, no issues were found.") + } + + ############################################################ + ## Extensive checks + ############################################################ + if (doExtras) { + cat("doExtras is ON\n") + ## Test that "no adjustment set" and "empty adjustment set" are distinguished properly + x <- 1; y <- 2 + cpdag <- matrix(c(0,1,1,0),2,2) ## 1 --- 2 => no adj set + dag <- matrix(c(0,1,0,0),2,2) ## 1 --> 2 => empty adj set + + adjC <- adjustment(amat = cpdag, amat.type = "cpdag", x = 1, y = 2, set.type = "canonical") + adjD <- adjustment(amat = dag, amat.type = "dag", x = 1, y = 2, set.type = "canonical") + adjP <- adjustment(amat = dag, amat.type = "pdag", x = 1, y = 2, set.type = "canonical") + + stopifnot(!identical(adjC, adjD), identical(adjD, adjP) ) + + adjCAll <- adjustment(amat = cpdag, amat.type = "cpdag", x = 1, y = 2, set.type = "all") + adjDAll <- adjustment(amat = dag, amat.type = "dag", x = 1, y = 2, set.type = "all") + adjPAll <- adjustment(amat = dag, amat.type = "pdag", x = 1, y = 2, set.type = "all") + + stopifnot( !identical(adjCAll, adjDAll), identical(adjDAll, adjPAll) ) + + adjCMin <- adjustment(amat = cpdag, amat.type = "cpdag", x = 1, y = 2, set.type = "minimal") + adjDMin <- adjustment(amat = dag, amat.type = "dag", x = 1, y = 2, set.type = "minimal") + adjPMin <- adjustment(amat = dag, amat.type = "pdag", x = 1, y = 2, set.type = "minimal") + + stopifnot( !identical(adjCMin, adjDMin), identical(adjDMin, adjPMin) ) + + + ##################################################################################### + ## Test 1: Compare CPDAG and PDAG implementation and validate all sets using gac() + ##################################################################################### + nreps <- 100 + simRes <- data.frame(setType = rep(NA, nreps), id = rep(NA,nreps), + rtCPDAG = rep(NA,nreps), rtPDAG = rep(NA, nreps), + nSet = rep(NA, nreps), gacCheck = rep(NA, nreps)) + proc.time() + for (i in 1:nreps) { + cat("i = ",i,"\n") + ## generate a graph + seed <- i + set.seed(seed) + p <- sample(x=5:10, size = 1) + prob <- sample(x=3:7/10, size = 1) + g <- pcalg:::randomDAG(p, prob) ## true DAG + cpdag <- dag2cpdag(g) + cpdag.mat <- t(as(cpdag,"matrix")) ## has correct encoding + + ## define input + amat <- cpdag.mat + x <- sample(x = 1:p, size = 1) + y <- sample(x = setdiff(1:p,x), size = 1) + set.type <- sample(x = c("all", "minimal"), size = 1) + simRes$setType[i] <- set.type + + ## run both methods + simRes$rtCPDAG[i] <- system.time(res1 <- adjustment(amat = amat, amat.type = "cpdag", x = x, y = y, set.type = set.type))[3] + simRes$rtPDAG[i] <- system.time(res2 <- adjustment(amat = amat, amat.type = "pdag", x = x, y = y, set.type = set.type))[3] + simRes$nSet[i] <- length(res1) + + if (length(res1) == 0) { + res1 <- vector("list", 0) + } + if (length(res2) == 0) { + res2 <- vector("list", 0) + } + ## compare results + simRes$id[i] <- identical(res1,res2) + + ## compare results with gac() based on "pdag" + if (length(res2) > 0) { + gc <- TRUE + for (j in 1:length(res2)) { + gc <- gc & gac(amat = amat, x = x, y = y, z = res2[[j]], type = "cpdag")$gac + } + simRes$gacCheck[i] <- gc + } + + } + proc.time() + + summary(simRes) + table(is.na(simRes$gacCheck), simRes$nSet == 0) + + ################################################ + ## Test 2: Check using predefined graphs + ################################################ + gacVSadj <- function(amat, x, y ,z, V, type) { + ## gac(z) is TRUE IFF z is returned by adjustment() + ## x,y,z: col positions as used in GAC + ## Result: TRUE is result is equal + typeDG <- switch(type, + dag = "dag", + cpdag = "cpdag", + mag = "mag", + pag = "pag") + gacRes <- gac(amat,x,y, z, type)$gac + adjRes <- adjustment(amat = amat, amat.type = typeDG, x = x, y = y, set.type = "all") + if (gacRes) { ## z is valid adj set + res <- any(sapply(adjRes, function(xx) setequal(z, xx))) + } else { ## z is not valid adj set + res <- all(!sapply(adjRes, function(xx) setequal(z, xx))) + } + res + } + + xx <- TRUE + ################################################## + ## DAG / CPDAG + ################################################## + ## CPDAG 1: Paper Fig 1 + mFig1 <- matrix(c(0,1,1,0,0,0, 1,0,1,1,1,0, 0,0,0,0,0,1, + 0,1,1,0,1,1, 0,1,0,1,0,1, 0,0,0,0,0,0), 6,6) + type <- "cpdag" + x <- 3; y <- 6 + + V <- as.character(1:ncol(mFig1)) + rownames(mFig1) <- colnames(mFig1) <- V + + xx <- xx & gacVSadj(mFig1,x,y, z=c(2,4), V=V, type) + xx <- xx & gacVSadj(mFig1,x,y, z=c(4,5), V=V, type) + xx <- xx & gacVSadj(mFig1,x,y, z=c(4,2,1), V=V, type) + xx <- xx & gacVSadj(mFig1,x,y, z=c(4,5,1), V=V, type) + xx <- xx & gacVSadj(mFig1,x,y, z=c(4,2,5), V=V, type) + xx <- xx & gacVSadj(mFig1,x,y, z=c(4,2,5,1), V=V, type) + xx <- xx & gacVSadj(mFig1,x,y, z= 2, V=V, type) + xx <- xx & gacVSadj(mFig1,x,y, z= NULL, V=V, type) + + ## CPDAG 2: Paper Fig 5a + mFig5a <- matrix(c(0,1,0,0,0, 1,0,1,0,0, 0,0,0,0,1, 0,0,1,0,0, 0,0,0,0,0), 5,5) + V <- as.character(1:ncol(mFig5a)) + rownames(mFig5a) <- colnames(mFig5a) <- V + + type <- "cpdag" + x <- c(1,5); y <- 4 + xx <- xx & gacVSadj(mFig5a, x,y, z=c(2,3), V=V, type) + xx <- xx & gacVSadj(mFig5a, x,y, z= 2, V=V, type) + + ## DAG 1 from Marloes' Talk + mMMd1 <- matrix(c(0,1,0,1,0,0, 0,0,1,0,1,0, 0,0,0,0,0,1, + 0,0,0,0,0,0, 0,0,0,0,0,0, 0,0,0,0,0,0),6,6) + V <- as.character(1:ncol(mMMd1)) + rownames(mMMd1) <- colnames(mMMd1) <- V + + type <- "dag" + x <- 1; y <- 3 + xx <- xx & gacVSadj(mMMd1, x,y, z=NULL, V=V, type) + xx <- xx & gacVSadj(mMMd1, x,y, z= 2, V=V, type) + xx <- xx & gacVSadj(mMMd1, x,y, z= 4, V=V, type) + xx <- xx & gacVSadj(mMMd1, x,y, z= 5, V=V, type) + xx <- xx & gacVSadj(mMMd1, x,y, z= 6, V=V, type) + xx <- xx & gacVSadj(mMMd1, x,y, z=c(4,5), V=V, type) + + ## DAG 2 from Marloes' Talk + mMMd2 <- matrix(c(0,1,0,1,0,0, 0,0,0,0,0,0, 0,1,0,0,1,0, + 0,0,0,0,1,0, 0,0,0,0,0,1, 0,0,0,0,0,0), 6,6) + V <- as.character(1:ncol(mMMd2)) + rownames(mMMd2) <- colnames(mMMd2) <- V + + type <- "dag" + x <- 4; y <- 6 + xx <- xx & gacVSadj(mMMd2, x,y, z= 1, V=V, type) + xx <- xx & gacVSadj(mMMd2, x,y, z= 3, V=V, type) + xx <- xx & gacVSadj(mMMd2, x,y, z= 5, V=V, type) + xx <- xx & gacVSadj(mMMd2, x,y, z=c(1,5), V=V, type) + xx <- xx & gacVSadj(mMMd2, x,y, z=c(1,2), V=V, type) + xx <- xx & gacVSadj(mMMd2, x,y, z=c(1,3), V=V, type) + xx <- xx & gacVSadj(mMMd2, x,y, z= 2, V=V, type) + + ################################################## + ## PAG + ################################################## + type <- "pag" + mFig3a <- matrix(c(0,1,0,0, 1,0,1,1, 0,1,0,1, 0,1,1,0), 4,4) + V <- as.character(1:ncol(mFig3a)) + rownames(mFig3a) <- colnames(mFig3a) <- V + xx <- xx & gacVSadj(mFig3a, x=2, y=4, z=NULL, V=V, type) + + mFig3b <- matrix(c(0,2,0,0, 3,0,3,3, 0,2,0,3, 0,2,2,0), 4,4) + V <- as.character(1:ncol(mFig3b)) + rownames(mFig3b) <- colnames(mFig3b) <- V + xx <- xx & gacVSadj(mFig3b, x=2, y=4, z=NULL, V=V, type) + + mFig3c <- matrix(c(0,3,0,0, 2,0,3,3, 0,2,0,3, 0,2,2,0), 4,4) + V <- as.character(1:ncol(mFig3c)) + rownames(mFig3c) <- colnames(mFig3c) <- V + xx <- xx & gacVSadj(mFig3c, x=2, y=4, z=NULL, V=V, type) + + mFig4a <- matrix(c(0,0,1,0,0,0, 0,0,1,0,0,0, 2,2,0,3,3,2, + 0,0,2,0,2,2, 0,0,2,1,0,2, 0,0,1,3,3,0), 6,6) + V <- as.character(1:ncol(mFig4a)) + rownames(mFig4a) <- colnames(mFig4a) <- V + xx <- xx & gacVSadj(mFig4a, x=3, y=4, z=NULL, V=V, type) + xx <- xx & gacVSadj(mFig4a, x=3, y=4, z= 6, V=V, type) + xx <- xx & gacVSadj(mFig4a, x=3, y=4, z=c(1,6), V=V, type) + xx <- xx & gacVSadj(mFig4a, x=3, y=4, z=c(2,6), V=V, type) + xx <- xx & gacVSadj(mFig4a, x=3, y=4, z=c(1,2,6), V=V, type) + + mFig4b <- matrix(c(0,0,1,0,0,0, 0,0,1,0,0,0, 2,2,0,0,3,2, + 0,0,0,0,2,2, 0,0,2,3,0,2, 0,0,2,3,2,0), 6,6) + V <- as.character(1:ncol(mFig4b)) + rownames(mFig4b) <- colnames(mFig4b) <- V + xx <- xx & gacVSadj(mFig4b, x=3, y=4, z=NULL, V=V, type) + xx <- xx & gacVSadj(mFig4b, x=3, y=4, z= 6, V=V, type) + xx <- xx & gacVSadj(mFig4b, x=3, y=4, z=c(5,6), V=V, type) + + mFig5b <- matrix(c(0,1,0,0,0,0,0, 2,0,2,3,0,3,0, 0,1,0,0,0,0,0, 0,2,0,0,3,0,0, + 0,0,0,2,0,2,3, 0,2,0,0,2,0,0, 0,0,0,0,2,0,0), 7,7) + V <- as.character(1:ncol(mFig5b)) + rownames(mFig5b) <- colnames(mFig5b) <- V + xx <- xx & gacVSadj(mFig5b, x=c(2,7), y=6, z=NULL, V=V, type) + xx <- xx & gacVSadj(mFig5b, x=c(2,7), y=6, z=c(4,5), V=V, type) + xx <- xx & gacVSadj(mFig5b, x=c(2,7), y=6, z=c(4,5,1), V=V, type) + xx <- xx & gacVSadj(mFig5b, x=c(2,7), y=6, z=c(4,5,3), V=V, type) + xx <- xx & gacVSadj(mFig5b, x=c(2,7), y=6, z=c(1,3,4,5), V=V, type) + + ## PAG in Marloes' talk + mMMp <- matrix(c(0,0,0,3,2,0,0, 0,0,0,0,1,0,0, 0,0,0,0,1,0,0, 2,0,0,0,0,3,2, + 3,2,2,0,0,0,3, 0,0,0,2,0,0,0, 0,0,0,2,2,0,0), 7,7) + V <- as.character(1:ncol(mMMp)) + rownames(mMMp) <- colnames(mMMp) <- V + + x <- c(5,6); y <- 7 + xx <- xx & gacVSadj(mMMp, x,y, z=NULL, V=V, type) + xx <- xx & gacVSadj(mMMp, x,y, z= 1, V=V, type) + xx <- xx & gacVSadj(mMMp, x,y, z= 4, V=V, type) + xx <- xx & gacVSadj(mMMp, x,y, z= 2, V=V, type) + xx <- xx & gacVSadj(mMMp, x,y, z= 3, V=V, type) + xx <- xx & gacVSadj(mMMp, x,y, z=c(2,3), V=V, type) + xx <- xx & gacVSadj(mMMp, x,y, z=c(1,4), V=V, type) + xx <- xx & gacVSadj(mMMp, x,y, z=c(1,4,2), V=V, type) + xx <- xx & gacVSadj(mMMp, x,y, z=c(1,4,3), V=V, type) + xx <- xx & gacVSadj(mMMp, x,y, z=c(1,4,2,3), V=V, type) + + ################################################## + ## V=V, type = "pag" -- Tests from Ema + ################################################## + type <- "pag" + pag.m <- readRDS(system.file(package="pcalg", "external", "gac-pags.rds")) + m1 <- pag.m[["m1"]] + V <- colnames(m1) + x <- 6; y <- 9 + xx <- xx & gacVSadj(m1,x,y, z=NULL, V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=1, V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=2, V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=3, V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=4, V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=c(2,3), V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=c(2,3,8), V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=c(2,3,7,8), V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=c(2,3,5,8), V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=c(2,3,5,7,8), V=V, type) + + x <- c(6,8); y <- 9 + xx <- xx & gacVSadj(m1,x,y, z=NULL, V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=1, V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=2, V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=3, V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=4, V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=c(2,3), V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=c(2,3,4), V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=c(2,3,7), V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=c(2,3,5), V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=c(2,3,5,7), V=V, type) + + x <- 3; y <- 1 + xx <- xx & gacVSadj(m1,x,y, z=NULL, V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=2, V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=4, V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=5, V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=6, V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=c(2,6), V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=c(2,8), V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=c(2,7,8), V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=c(2,5,8), V=V, type) + xx <- xx & gacVSadj(m1,x,y, z=c(2,5,7,8), V=V, type) + + m2 <- pag.m[["m2"]] + V <- colnames(m2) + x <- 3; y <-1 + xx <- xx & gacVSadj(m2,x,y, z=NULL, V=V, type) + xx <- xx & gacVSadj(m2,x,y, z=2, V=V, type) + xx <- xx & gacVSadj(m2,x,y, z=4, V=V, type) + xx <- xx & gacVSadj(m2,x,y, z=c(2,8), V=V, type) + xx <- xx & gacVSadj(m2,x,y, z=8, V=V, type) + xx <- xx & gacVSadj(m2,x,y, z=9, V=V, type) + xx <- xx & gacVSadj(m2,x,y, z=c(2,8,9), V=V, type) + xx <- xx & gacVSadj(m2,x,y, z=c(2,5), V=V, type) + + x <- c(3,9); y <- 1 + xx <- xx & gacVSadj(m2,x,y, z=NULL, V=V, type) + xx <- xx & gacVSadj(m2,x,y, z=2, V=V, type) + xx <- xx & gacVSadj(m2,x,y, z=4, V=V, type) + xx <- xx & gacVSadj(m2,x,y, z=c(2,8), V=V, type) + xx <- xx & gacVSadj(m2,x,y, z=8, V=V, type) + xx <- xx & gacVSadj(m2,x,y, z=9, V=V, type) + xx <- xx & gacVSadj(m2,x,y, z=c(2,8,9), V=V, type) + xx <- xx & gacVSadj(m2,x,y, z=c(2,5), V=V, type) + + m3 <- pag.m[["m3"]] + V <- colnames(m3) + x <- 1; y <- 9 + xx <- xx & gacVSadj(m3,x,y, z=NULL, V=V, type) + xx <- xx & gacVSadj(m3,x,y, z=2, V=V, type) + xx <- xx & gacVSadj(m3,x,y, z=3, V=V, type) + xx <- xx & gacVSadj(m3,x,y, z=5, V=V, type) + xx <- xx & gacVSadj(m3,x,y, z=7, V=V, type) + xx <- xx & gacVSadj(m3,x,y, z=8, V=V, type) + xx <- xx & gacVSadj(m3,x,y, z=c(2,3), V=V, type) + xx <- xx & gacVSadj(m3,x,y, z=c(5,7), V=V, type) + + x <- 1; y <- 8 + xx <- xx & gacVSadj(m3,x,y, z=NULL, V=V, type) + xx <- xx & gacVSadj(m3,x,y, z=2, V=V, type) + xx <- xx & gacVSadj(m3,x,y, z=3, V=V, type) + xx <- xx & gacVSadj(m3,x,y, z=5, V=V, type) + xx <- xx & gacVSadj(m3,x,y, z=7, V=V, type) + xx <- xx & gacVSadj(m3,x,y, z=9, V=V, type) + xx <- xx & gacVSadj(m3,x,y, z=c(2,3), V=V, type) + xx <- xx & gacVSadj(m3,x,y, z=c(5,9), V=V, type) + + if (!xx) { + stop("Problem when testing function gacVSadj.") + } else { + message("OK, no issues were found.") + } + + ################################################## + ## given same graph, type=cpdag and type=pdag + ## should give same canonical set + ################################################## + m <- rbind(c(0,1,0,0,0,0), + c(1,0,1,0,0,0), + c(0,1,0,0,0,0), + c(0,0,0,0,0,0), + c(0,1,1,1,0,0), + c(1,0,1,1,1,0)) + colnames(m) <- rownames(m) <- as.character(1:6) + + ## You can see that the current adjustment function outputs different sets + ## if type = "cpdag" or type = "pdag" which shouldn't happen + ## because it is the same graph: + res1 <- adjustment(m,amat.type="cpdag",2,4,set.type="canonical") + res2 <- adjustment(m,amat.type="pdag",2,4,set.type="canonical") + + if (!all.equal(res1, res2)) { + stop("Canonical set is not the same for type=cpdag and type=pdag\n") + } + + } ## doExtras + else { + cat("doExtras is OFF\n") + } + + } # only if CRAN package {dagitty} is available Loading required namespace: dagitty OK, no issues were found. doExtras is OFF > > proc.time() user system elapsed 2.12 0.18 2.10