Nothing
# R Header
### Copyright (C) 2017-2023 Torsten Hothorn
###
### This file is part of the 'libcoin' R add-on package.
###
### 'libcoin' is free software: you can redistribute it and/or modify
### it under the terms of the GNU General Public License as published by
### the Free Software Foundation, version 2.
###
### 'libcoin' is distributed in the hope that it will be useful,
### but WITHOUT ANY WARRANTY; without even the implied warranty of
### MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
### GNU General Public License for more details.
###
### You should have received a copy of the GNU General Public License
### along with 'libcoin'. If not, see <http://www.gnu.org/licenses/>.
###
###
### DO NOT EDIT THIS FILE
###
### Edit 'libcoin.w' and run 'nuweb -r libcoin.w'
# LinStatExpCov
LinStatExpCov <-
function(X, Y, ix = NULL, iy = NULL, weights = integer(0),
subset = integer(0), block = integer(0), checkNAs = TRUE,
varonly = FALSE, nresample = 0, standardise = FALSE,
tol = sqrt(.Machine$double.eps))
{
if (missing(X) && !is.null(ix) && is.null(iy)) {
X <- ix
ix <- NULL
}
if (missing(X)) X <- integer(0)
## <FIXME> for the time being only!!! </FIXME>
## if (length(subset) > 0) subset <- sort(subset)
if (is.null(ix) && is.null(iy))
.LinStatExpCov1d(X = X, Y = Y,
weights = weights, subset = subset,
block = block, checkNAs = checkNAs,
varonly = varonly, nresample = nresample,
standardise = standardise, tol = tol)
else if (!is.null(ix) && !is.null(iy))
.LinStatExpCov2d(X = X, Y = Y, ix = ix, iy = iy,
weights = weights, subset = subset,
block = block, checkNAs = checkNAs,
varonly = varonly, nresample = nresample,
standardise = standardise, tol = tol)
else
stop("incorrect call to ", sQuote("LinStatExpCov()"))
}
# LinStatExpCov1d
.LinStatExpCov1d <-
function(X, Y, weights = integer(0), subset = integer(0), block = integer(0),
checkNAs = TRUE, varonly = FALSE, nresample = 0, standardise = FALSE,
tol = sqrt(.Machine$double.eps))
{
if (NROW(X) != NROW(Y))
stop("dimensions of X and Y don't match")
N <- NROW(X)
if (is.integer(X)) {
if (is.null(attr(X, "levels")) || checkNAs) {
rg <- range(X)
if (anyNA(rg))
stop("no missing values allowed in X")
stopifnot(rg[1] > 0) # no missing values allowed here!
if (is.null(attr(X, "levels")))
attr(X, "levels") <- seq_len(rg[2])
}
}
if (is.factor(X) && checkNAs)
stopifnot(!anyNA(X))
# Check weights, subset, block
if (is.null(weights)) weights <- integer(0)
if (length(weights) > 0) {
if (!((N == length(weights)) && all(weights >= 0)))
stop("incorrect weights")
if (checkNAs) stopifnot(!anyNA(weights))
}
if (is.null(subset)) subset <- integer(0)
if (length(subset) > 0 && checkNAs) {
rs <- range(subset)
if (anyNA(rs)) stop("no missing values allowed in subset")
if (!((rs[2] <= N) && (rs[1] >= 1L)))
stop("incorrect subset")
}
if (is.null(block)) block <- integer(0)
if (length(block) > 0) {
if (!((N == length(block)) && is.factor(block)))
stop("incorrect block")
if (checkNAs) stopifnot(!anyNA(block))
}
if (checkNAs) {
# Handle Missing Values
ms <- !complete.cases(X, Y)
if (all(ms))
stop("all observations are missing")
if (any(ms)) {
if (length(subset) > 0) {
if (all(subset %in% which(ms)))
stop("all observations are missing")
subset <- subset[!(subset %in% which(ms))]
} else {
subset <- seq_len(N)[-which(ms)]
}
}
}
ret <- .Call(R_ExpectationCovarianceStatistic, X, Y, weights, subset,
block, as.integer(varonly), as.double(tol))
ret$varonly <- as.logical(ret$varonly)
ret$Xfactor <- as.logical(ret$Xfactor)
if (nresample > 0) {
ret$PermutedLinearStatistic <-
.Call(R_PermutedLinearStatistic, X, Y, weights, subset,
block, as.double(nresample))
if (standardise)
ret$StandardisedPermutedLinearStatistic <-
.Call(R_StandardisePermutedLinearStatistic, ret)
}
class(ret) <- c("LinStatExpCov1d", "LinStatExpCov")
ret
}
# LinStatExpCov2d
.LinStatExpCov2d <-
function(X = numeric(0), Y, ix, iy, weights = integer(0), subset = integer(0),
block = integer(0), checkNAs = TRUE, varonly = FALSE, nresample = 0,
standardise = FALSE, tol = sqrt(.Machine$double.eps))
{
IF <- function(x) is.integer(x) || is.factor(x)
if (!((length(ix) == length(iy)) && IF(ix) && IF(iy)))
stop("incorrect ix and/or iy")
N <- length(ix)
# Check ix
if (is.null(attr(ix, "levels"))) {
rg <- range(ix)
if (anyNA(rg))
stop("no missing values allowed in ix")
stopifnot(rg[1] >= 0)
attr(ix, "levels") <- seq_len(rg[2])
} else {
## lev can be data.frame (see inum::inum)
lev <- attr(ix, "levels")
if (!is.vector(lev)) lev <- seq_len(NROW(lev))
attr(ix, "levels") <- lev
if (checkNAs) stopifnot(!anyNA(ix))
}
# Check iy
if (is.null(attr(iy, "levels"))) {
rg <- range(iy)
if (anyNA(rg))
stop("no missing values allowed in iy")
stopifnot(rg[1] >= 0)
attr(iy, "levels") <- seq_len(rg[2])
} else {
## lev can be data.frame (see inum::inum)
lev <- attr(iy, "levels")
if (!is.vector(lev)) lev <- seq_len(NROW(lev))
attr(iy, "levels") <- lev
if (checkNAs) stopifnot(!anyNA(iy))
}
if (length(X) > 0) {
if (!(NROW(X) == (length(attr(ix, "levels")) + 1) &&
all(complete.cases(X))))
stop("incorrect X")
}
if (!(NROW(Y) == (length(attr(iy, "levels")) + 1) &&
all(complete.cases(Y))))
stop("incorrect Y")
# Check weights, subset, block
if (is.null(weights)) weights <- integer(0)
if (length(weights) > 0) {
if (!((N == length(weights)) && all(weights >= 0)))
stop("incorrect weights")
if (checkNAs) stopifnot(!anyNA(weights))
}
if (is.null(subset)) subset <- integer(0)
if (length(subset) > 0 && checkNAs) {
rs <- range(subset)
if (anyNA(rs)) stop("no missing values allowed in subset")
if (!((rs[2] <= N) && (rs[1] >= 1L)))
stop("incorrect subset")
}
if (is.null(block)) block <- integer(0)
if (length(block) > 0) {
if (!((N == length(block)) && is.factor(block)))
stop("incorrect block")
if (checkNAs) stopifnot(!anyNA(block))
}
ret <- .Call(R_ExpectationCovarianceStatistic_2d, X, ix, Y, iy,
weights, subset, block, as.integer(varonly), as.double(tol))
ret$varonly <- as.logical(ret$varonly)
ret$Xfactor <- as.logical(ret$Xfactor)
if (nresample > 0) {
ret$PermutedLinearStatistic <-
.Call(R_PermutedLinearStatistic_2d, X, ix, Y, iy, block, nresample, ret$Table)
if (standardise)
ret$StandardisedPermutedLinearStatistic <-
.Call(R_StandardisePermutedLinearStatistic, ret)
}
class(ret) <- c("LinStatExpCov2d", "LinStatExpCov")
ret
}
# vcov LinStatExpCov
vcov.LinStatExpCov <-
function(object, ...)
{
if (object$varonly)
stop("cannot extract covariance matrix")
drop(.Call(R_unpack_sym, object$Covariance, NULL, 0L))
}
# doTest
### note: lower = FALSE => p-value; lower = TRUE => 1 - p-value
doTest <-
function(object, teststat = c("maximum", "quadratic", "scalar"),
alternative = c("two.sided", "less", "greater"), pvalue = TRUE,
lower = FALSE, log = FALSE, PermutedStatistics = FALSE,
minbucket = 10L, ordered = TRUE, maxselect = object$Xfactor,
pargs = GenzBretz())
{
teststat <- match.arg(teststat, choices = c("maximum", "quadratic", "scalar"))
if (!any(teststat == c("maximum", "quadratic", "scalar")))
stop("incorrect teststat")
alternative <- alternative[1]
if (!any(alternative == c("two.sided", "less", "greater")))
stop("incorrect alternative")
if (maxselect)
stopifnot(object$Xfactor)
if (teststat == "quadratic" || maxselect) {
if (alternative != "two.sided")
stop("incorrect alternative")
}
test <- which(c("maximum", "quadratic", "scalar") == teststat)
if (test == 3) {
if (length(object$LinearStatistic) != 1)
stop("scalar test statistic not applicable")
test <- 1L # scalar is maximum internally
}
alt <- which(c("two.sided", "less", "greater") == alternative)
if (!pvalue && (NCOL(object$PermutedLinearStatistic) > 0))
object$PermutedLinearStatistic <- matrix(NA_real_, nrow = 0, ncol = 0)
if (!maxselect) {
if (teststat == "quadratic") {
ret <- .Call(R_QuadraticTest, object, as.integer(pvalue), as.integer(lower),
as.integer(log), as.integer(PermutedStatistics))
} else {
ret <- .Call(R_MaximumTest, object, as.integer(alt), as.integer(pvalue),
as.integer(lower), as.integer(log), as.integer(PermutedStatistics),
as.integer(pargs$maxpts), as.double(pargs$releps),
as.double(pargs$abseps))
if (teststat == "scalar") {
var <- if (object$varonly) object$Variance else object$Covariance
ret$TestStatistic <- object$LinearStatistic - object$Expectation
ret$TestStatistic <-
if (var > object$tol) ret$TestStatistic / sqrt(var) else NaN
}
}
} else {
ret <- .Call(R_MaximallySelectedTest, object, as.integer(ordered), as.integer(test),
as.integer(minbucket), as.integer(lower), as.integer(log))
}
if (!PermutedStatistics) ret$PermutedStatistics <- NULL
ret
}
# Contrasts
lmult <-
function(x, object)
{
stopifnot(!object$varonly)
stopifnot(is.numeric(x))
if (is.vector(x)) x <- matrix(x, nrow = 1)
P <- object$dimension[1]
stopifnot(ncol(x) == P)
Q <- object$dimension[2]
ret <- object
xLS <- x %*% matrix(object$LinearStatistic, nrow = P)
xExp <- x %*% matrix(object$Expectation, nrow = P)
xExpX <- x %*% matrix(object$ExpectationX, nrow = P)
if (Q == 1) {
xCov <- tcrossprod(x %*% vcov(object), x)
} else {
zmat <- matrix(0, nrow = P * Q, ncol = nrow(x))
mat <- rbind(t(x), zmat)
mat <- mat[rep.int(seq_len(nrow(mat)), Q - 1),, drop = FALSE]
mat <- rbind(mat, t(x))
mat <- matrix(mat, ncol = Q * nrow(x))
mat <- t(mat)
xCov <- tcrossprod(mat %*% vcov(object), mat)
}
if (!is.matrix(xCov)) xCov <- matrix(xCov)
if (length(object$PermutedLinearStatistic) > 0) {
xPS <- apply(object$PermutedLinearStatistic, 2, function(y)
as.vector(x %*% matrix(y, nrow = P)))
if (!is.matrix(xPS)) xPS <- matrix(xPS, nrow = 1)
ret$PermutedLinearStatistic <- xPS
}
ret$LinearStatistic <- as.vector(xLS)
ret$Expectation <- as.vector(xExp)
ret$ExpectationX <- as.vector(xExpX)
ret$Covariance <- as.vector(xCov[lower.tri(xCov, diag = TRUE)])
ret$Variance <- diag(xCov)
ret$dimension <- c(NROW(x), Q)
ret$Xfactor <- FALSE
if (length(object$StandardisedPermutedLinearStatistic) > 0)
ret$StandardisedPermutedLinearStatistic <-
.Call(R_StandardisePermutedLinearStatistic, ret)
ret
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.