R Under development (unstable) (2023-11-26 r85638 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 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.

> 
> set.seed(290875)
> library("party")
Loading required package: grid
Loading required package: mvtnorm
Loading required package: modeltools
Loading required package: stats4
Loading required package: strucchange
Loading required package: zoo

Attaching package: 'zoo'

The following objects are masked from 'package:base':

    as.Date, as.Date.numeric

Loading required package: sandwich
> if (!require("mvtnorm"))
+     stop("cannot load package mvtnorm")
> 
> 
> ### get rid of the NAMESPACE
> attach(list2env(as.list(asNamespace("party"))))
The following objects are masked from package:party:

    cforest, cforest_classical, cforest_control, cforest_unbiased,
    conditionalTree, ctree, ctree_control, edge_simple,
    initVariableFrame, mob, mob_control, node_barplot, node_bivplot,
    node_boxplot, node_density, node_hist, node_inner,
    node_scatterplot, node_surv, node_terminal, nodes, party_intern,
    prettytree, proximity, ptrafo, response, reweight, sctest.mob,
    treeresponse, varimp, varimpAUC, where

> 
> ### 
> ###
> ###    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)))
> 
> stopifnot(isequal(2*pnorm(-t), 
+           .Call(R_maxabsConditionalPvalue, t, matrix(1), as.integer(1), 0.0, 0.0, 0.0)))
> 
> 
> 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)
> 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)
> 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)
> p2 <- pmvnorm(lower = rep(-a,4), upper = rep(a,4), corr = b)
> stopifnot(isequal(round(p1, 3), round(1 - p2, 3)))
> 
> if (FALSE) {
+ ### 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)
+ 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)
+ 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)
+ stopifnot(isequal(wstat, pvals[[1]]))
+ stopifnot(all(wpval < (1 - pvals[[2]])))
+ }
> 
> proc.time()
   user  system elapsed 
   1.76    0.37    2.12