--- a +++ b/partyMod/tests/Distributions.R @@ -0,0 +1,117 @@ + +set.seed(290875) +library("party") +if (!require("mvtnorm")) + stop("cannot load package mvtnorm") + + +### get rid of the NAMESPACE +attach(asNamespace("party")) + +### +### +### Regression tests for conditional distributions +### +### functions defined in file `./src/Distributions.c' + +### chisq-distribution of quadratic forms +t <- 2.1 +df <- 2 +storage.mode(t) <- "double" +storage.mode(df) <- "double" +stopifnot(isequal(1 - pchisq(t, df = df), ### P-values!!! + .Call("R_quadformConditionalPvalue", t, df, PACKAGE = "party"))) + +stopifnot(isequal(2*pnorm(-t), + .Call("R_maxabsConditionalPvalue", t, matrix(1), as.integer(1), 0.0, 0.0, 0.0, PACKAGE = "party"))) + + +maxpts <- 25000 +storage.mode(maxpts) <- "integer" +abseps <- 0.0001 +releps <- 0 +tol <- 1e-10 + +a <- 1.96 +b <- diag(2) + +p1 <- .Call("R_maxabsConditionalPvalue", a, b, maxpts, abseps, releps, tol, PACKAGE = "party") +p2 <- pmvnorm(lower = rep(-a,2), upper = rep(a,2), corr = b) +stopifnot(isequal(round(p1, 3), round(1 - p2, 3))) + +b <- diag(4) +p1 <- .Call("R_maxabsConditionalPvalue", a, b, maxpts, abseps, releps, tol, PACKAGE = "party") +p2 <- pmvnorm(lower = rep(-a,4), upper = rep(a,4), corr = b) +stopifnot(isequal(round(p1, 3), round(1 - p2, 3))) + +b <- diag(4) +b[upper.tri(b)] <- c(0.1, 0.2, 0.3) +b[lower.tri(b)] <- t(b)[lower.tri(b)] +p1 <- .Call("R_maxabsConditionalPvalue", a, b, maxpts, abseps, releps, tol, PACKAGE = "party") +p2 <- pmvnorm(lower = rep(-a,4), upper = rep(a,4), corr = b) +stopifnot(isequal(round(p1, 3), round(1 - p2, 3))) + +### Monte-Carlo approximation of P-Values, univariate +mydata = data.frame(y = gl(2, 50), x1 = rnorm(100), + x2 = rnorm(100), x3 = rnorm(100)) +inp <- initVariableFrame(mydata[,"x1",drop = FALSE], trafo = function(data) +ptrafo(data, numeric_trafo = rank)) +resp <- initVariableFrame(mydata[,"y",drop = FALSE], trafo = NULL, response = TRUE) +ls <- new("LearningSample", inputs = inp, responses = resp, + weights = rep(1, inp@nobs), nobs = nrow(mydata), + ninputs = inp@ninputs) +tm <- ctree_memory(ls) +varctrl <- new("VariableControl") +varctrl@teststat <- factor("max", levels = c("max", "quad")) +varctrl@pvalue <- FALSE +gtctrl <- new("GlobalTestControl") +gtctrl@testtype <- factor("MonteCarlo", levels = levels(gtctrl@testtype)) +gtctrl@nresample <- as.integer(19999) + +pvals <- .Call("R_GlobalTest", ls, ls@weights, tm, varctrl, gtctrl, PACKAGE = "party") +wstat <- abs(qnorm(wilcox.test(x1 ~ y, data = mydata, + exact = FALSE, correct = FALSE)$p.value/2)) +wpval <- wilcox.test(x1 ~ y, data = mydata, exact = TRUE)$p.value +stopifnot(isequal(wstat, pvals[[1]])) +stopifnot(abs(wpval - (1 - pvals[[2]])) < 0.01) + +### Monte-Carlo approximations of P-Values, multiple inputs +mydata = data.frame(y = gl(2, 50), x1 = rnorm(100), + x2 = rnorm(100), x3 = rnorm(100)) +inp <- initVariableFrame(mydata[,c("x1", "x2", "x3"), + drop = FALSE], trafo = function(data) +ptrafo(data, numeric_trafo = rank)) +resp <- initVariableFrame(mydata[,"y",drop = FALSE], trafo = NULL, response = TRUE) +ls <- new("LearningSample", inputs = inp, responses = resp, + weights = rep(1, inp@nobs), nobs = nrow(mydata), + ninputs = inp@ninputs) +tm <- ctree_memory(ls) +varctrl <- new("VariableControl") +varctrl@teststat <- factor("max", levels = c("max", "quad")) +varctrl@pvalue <- TRUE +gtctrl <- new("GlobalTestControl") +gtctrl@testtype <- factor("Univariate", levels = levels(gtctrl@testtype)) +gtctrl@nresample <- as.integer(19999) + +pvals <- .Call("R_GlobalTest", ls, ls@weights, tm, varctrl, gtctrl, PACKAGE = "party") +wstat <- c(abs(qnorm(wilcox.test(x1 ~ y, data = mydata, + exact = FALSE, correct = FALSE)$p.value/2)), + abs(qnorm(wilcox.test(x2 ~ y, data = mydata, + exact = FALSE, correct = FALSE)$p.value/2)), + abs(qnorm(wilcox.test(x3 ~ y, data = mydata, + exact = FALSE, correct = FALSE)$p.value/2))) +wpval <- c(wilcox.test(x1 ~ y, data = mydata, + exact = FALSE, correct = FALSE)$p.value, + wilcox.test(x2 ~ y, data = mydata, + exact = FALSE, correct = FALSE)$p.value, + wilcox.test(x3 ~ y, data = mydata, + exact = FALSE, correct = FALSE)$p.value) +stopifnot(isequal(wstat, pvals[[1]])) +stopifnot(isequal(wpval, 1 - pvals[[2]])) + +### Monte-Carlo approximations of P-Values, min-P approach +gtctrl@testtype <- factor("MonteCarlo", levels = levels(gtctrl@testtype)) +gtctrl@nresample <- as.integer(19999) +pvals <- .Call("R_GlobalTest", ls, ls@weights, tm, varctrl, gtctrl, PACKAGE = "party") +stopifnot(isequal(wstat, pvals[[1]])) +stopifnot(all(wpval < (1 - pvals[[2]])))