Nothing
# Permutation methods----
#' Permutation p-value
#'
#' Computes a p-value using a permutation procedure.
#'
#' @param pcevObj A pcev object of class \code{PcevClassical} or \code{PcevSingular}
#' \code{PcevBlock}
#' @param shrink Should we use a shrinkage estimate of the residual variance?
#' @param index If \code{pcevObj} is of class \code{PcevBlock}, \code{index} is a
#' vector describing the block to which individual response variables
#' correspond.
#' @param nperm The number of permutations to perform.
#' @param ... Extra parameters.
#' @export
permutePval <- function(pcevObj, ...) UseMethod("permutePval")
#' @rdname permutePval
permutePval.default <- function(pcevObj, ...) {
stop(strwrap("This function should be used with a Pcev object of class
PcevClassical, PcevSingular or PcevBlock"),
call. = FALSE)
}
#' @rdname permutePval
permutePval.PcevClassical <- function(pcevObj, shrink, index, nperm, ...) {
results <- estimatePcev(pcevObj, shrink)
N <- nrow(pcevObj$Y)
PCEV <- pcevObj$Y %*% results$weights
initFit <- lm.fit(pcevObj$X, PCEV)
df1 <- ncol(pcevObj$X) - 1
df2 <- N - ncol(pcevObj$X)
initFstat <- (sum((mean(PCEV) - initFit$fitted.values)^2)/df1)/(sum(initFit$residuals^2)/df2)
initPval <- pf(initFstat, df1, df2, lower.tail = FALSE)
if (nperm) {
permutationPvalues <- replicate(nperm, expr = {
tmp <- pcevObj
tmp$Y <- tmp$Y[sample(N), ]
tmpRes <- try(estimatePcev(tmp, shrink, index), silent = TRUE)
if (inherits(tmpRes, "try-error")) {
return(NA)
} else {
tmpPCEV <- tmp$Y %*% tmpRes$weights
tmpFit <- lm.fit(tmp$X, tmpPCEV)
tmpFstat <- (sum((mean(tmpPCEV) - tmpFit$fitted.values)^2)/df1)/(sum(tmpFit$residuals^2)/df2)
return(pf(tmpFstat, df1, df2, lower.tail = FALSE))
}
})
pvalue <- mean(permutationPvalues < initPval)
} else {
pvalue <- NA
}
results$pvalue <- pvalue
return(results)
}
#' @rdname permutePval
permutePval.PcevBlock <- function(pcevObj, shrink, index, nperm, ...) {
results <- estimatePcev(pcevObj, shrink, index)
N <- nrow(pcevObj$Y)
PCEV <- pcevObj$Y %*% results$weights
initFit <- lm.fit(pcevObj$X, PCEV)
df1 <- ncol(pcevObj$X) - 1
df2 <- N - ncol(pcevObj$X)
initFstat <- (sum((mean(PCEV) - initFit$fitted.values)^2)/df1)/(sum(initFit$residuals^2)/df2)
initPval <- pf(initFstat, df1, df2, lower.tail = FALSE)
if (nperm) {
permutationPvalues <- replicate(nperm, expr = {
tmp <- pcevObj
tmp$Y <- tmp$Y[sample(N), ]
tmpRes <- try(estimatePcev(tmp, shrink, index), silent = TRUE)
if (inherits(tmpRes, "try-error")) {
return(NA)
} else {
tmpPCEV <- tmp$Y %*% tmpRes$weights
tmpFit <- lm.fit(tmp$X, tmpPCEV)
tmpFstat <- (sum((mean(tmpPCEV) - tmpFit$fitted.values)^2)/df1)/(sum(tmpFit$residuals^2)/df2)
return(pf(tmpFstat, df1, df2, lower.tail = FALSE))
}
})
pvalue <- mean(permutationPvalues < initPval)
} else {
pvalue <- NA
}
results$pvalue <- pvalue
return(results)
}
#' @rdname permutePval
permutePval.PcevSingular <- function(pcevObj, shrink, index, nperm, ...) {
results <- estimatePcev(pcevObj, shrink)
N <- nrow(pcevObj$Y)
PCEV <- pcevObj$Y %*% results$weights
initFit <- lm.fit(pcevObj$X, PCEV)
df1 <- ncol(pcevObj$X) - 1
df2 <- N - ncol(pcevObj$X)
initFstat <- (sum((mean(PCEV) - initFit$fitted.values)^2)/df1)/(sum(initFit$residuals^2)/df2)
initPval <- pf(initFstat, df1, df2, lower.tail = FALSE)
if (nperm) {
permutationPvalues <- replicate(nperm, expr = {
tmp <- pcevObj
tmp$Y <- tmp$Y[sample(N), ]
tmpRes <- try(estimatePcev(tmp, shrink, index), silent = TRUE)
if (inherits(tmpRes, "try-error")) {
return(NA)
} else {
tmpPCEV <- tmp$Y %*% tmpRes$weights
tmpFit <- lm.fit(tmp$X, tmpPCEV)
tmpFstat <- (sum((mean(tmpPCEV) - tmpFit$fitted.values)^2)/df1)/(sum(tmpFit$residuals^2)/df2)
return(pf(tmpFstat, df1, df2, lower.tail = FALSE))
}
})
pvalue <- mean(permutationPvalues < initPval)
} else {
pvalue <- NA
}
results$pvalue <- pvalue
return(results)
}
###################
# Wilks methods----
#' Wilks' lambda exact test
#'
#' Computes a p-value using Wilks' Lambda.
#'
#' The null distribution of this test statistic is only known in the case of a
#' single covariate, and therefore this is the only case implemented.
#'
#' @param pcevObj A pcev object of class \code{PcevClassical} or
#' \code{PcevBlock}
#' @param shrink Should we use a shrinkage estimate of the residual variance?
#' @param index If \code{pcevObj} is of class \code{PcevBlock}, \code{index} is
#' a vector describing the block to which individual response variables
#' correspond.
#' @param ... Extra parameters.
#' @export
wilksPval <- function(pcevObj, ...) UseMethod("wilksPval")
#' @rdname wilksPval
wilksPval.default <- function(pcevObj, ...) {
stop(strwrap("This function should be used with a Pcev object of class
PcevClassical, PcevBlock or PcevSingular"),
call. = FALSE)
}
#' @rdname wilksPval
wilksPval.PcevClassical <- function(pcevObj, shrink, index, ...) {
results <- estimatePcev(pcevObj, shrink)
N <- nrow(pcevObj$Y)
p <- ncol(pcevObj$Y)
d <- results$largestRoot
wilksLambda <- (N - p - 1)/p * d
df1 <- p
df2 <- N - p - 1
pvalue <- pf(wilksLambda, df1, df2, lower.tail = FALSE)
results$pvalue <- pvalue
return(results)
}
#' @rdname wilksPval
wilksPval.PcevSingular <- function(pcevObj, shrink, index, ...) {
stop(strwrap("Wilks test is only available for estimation = all."),
call. = FALSE)
}
#' @rdname wilksPval
wilksPval.PcevBlock <- function(pcevObj, shrink, index, ...) {
stop("Only inference = \"permutation\" is available for the block approach.",
call. = FALSE)
}
################################
# Roy's largest root methods----
#' Roy's largest root exact test
#'
#' In the classical domain of PCEV applicability this function uses Johnstone's
#' approximation to the null distribution of ' Roy's Largest Root statistic.
#' It uses a location-scale variant of the Tracy-Widom distribution of order 1.
#'
#' Note that if \code{shrink} is set to \code{TRUE}, the location-scale
#' parameters are estimated using a small number of permutations.
#'
#' @param pcevObj A pcev object of class \code{PcevClassical} or
#' \code{PcevBlock}
#' @param shrink Should we use a shrinkage estimate of the residual variance?
#' @param index If \code{pcevObj} is of class \code{PcevBlock}, \code{index} is
#' a vector describing the block to which individual response variables
#' correspond
#' @param nperm Number of permutations for Tracy-Widom empirical estimate.
#' @param ... Extra parameters.
#' @export
roysPval <- function(pcevObj, ...) UseMethod("roysPval")
#' @rdname roysPval
roysPval.default <- function(pcevObj, ...) {
stop(strwrap("This function should be used with a Pcev object of class
PcevClassical or PcevSingular"),
call. = FALSE)
}
#' @rdname roysPval
roysPval.PcevClassical <- function(pcevObj, shrink, index, ...) {
results <- estimatePcev(pcevObj, shrink)
n <- nrow(pcevObj$Y)
p <- ncol(pcevObj$Y)
q <- ncol(pcevObj$X) - 1
d <- results$largestRoot
params <- JohnstoneParam(p, n - q, q)
if (shrink) {
# Estimate the null distribution using
# permutations and MLE
null_dist <- replicate(25, expr = {
tmp <- pcevObj
tmp$Y <- tmp$Y[sample(n), ]
tmpRes <- try(estimatePcev(tmp, shrink, index), silent = TRUE)
if (inherits(tmpRes, "try-error")) {
return(NA)
} else {
return(tmpRes$largestRoot)
}
})
# Fit a location-scale version of TW distribution
# Note: likelihood may throw warnings at some evaluations,
# which is OK
oldw <- getOption("warn")
options(warn = -1)
res <- optim(c(params[1], params[2]), function(param) logLik(param, log(null_dist)),
control = list(fnscale = -1))
options(warn = oldw)
mu1 <- res$par[1]
sigma1 <- res$par[2]
TW <- (log(d) - mu1)/sigma1
pvalue <- RMTstat::ptw(TW, beta = 1, lower.tail = FALSE, log.p = FALSE)
} else {
# TW <- (log(theta/(1-theta)) - mu)/sigma
TW <- (log(d) - params[1])/params[2]
pvalue <- RMTstat::ptw(TW, beta = 1, lower.tail = FALSE, log.p = FALSE)
}
results$pvalue <- pvalue
return(results)
}
#' @rdname roysPval
roysPval.PcevSingular <- function(pcevObj, shrink, index, nperm, ...) {
results <- estimatePcev(pcevObj, shrink)
n <- nrow(pcevObj$Y)
p <- ncol(pcevObj$Y)
q <- ncol(pcevObj$X)
d <- results$largestRoot
if (missing(nperm)) nperm <- 50
null_dist <- replicate(nperm, expr = {
tmp <- pcevObj
tmp$Y <- tmp$Y[sample(n), ]
tmpRes <- try(estimatePcev(tmp, shrink, index), silent = TRUE)
if (inherits(tmpRes, "try-error")) {
return(NA)
} else {
return(tmpRes$largestRoot)
}
})
# Use method of moments
muTW <- -1.2065335745820
sigmaTW <- sqrt(1.607781034581)
muS <- mean(log(null_dist))
sigmaS <- stats::sd(log(null_dist))
sigma1 <- sigmaS/sigmaTW
mu1 <- muS - sigma1 * muTW
TW <- (log(d) - mu1)/sigma1
pvalue <- RMTstat::ptw(TW, beta = 1, lower.tail = FALSE, log.p = FALSE)
results$pvalue <- pvalue
return(results)
}
#' @rdname roysPval
roysPval.PcevBlock <- function(pcevObj, shrink, index, ...) {
stop("Only inference = \"permutation\" is available for the block approach.",
call. = FALSE)
}
##################
# Print method----
#' @export
print.Pcev <- function(x, ...) {
# Provide a summary of the results
N <- nrow(x$pcevObj$Y)
p <- ncol(x$pcevObj$Y)
q <- ncol(x$pcevObj$X)
if (x$Wilks) {
exact <- "Wilks' lambda test)"
} else {
exact <- "Roy's largest root test)"
}
cat("\nPrincipal component of explained variance\n")
cat("\n", N, "observations,", p, "response variables\n")
cat("\nEstimation method:", x$methods[1])
cat("\nInference method:", x$methods[2])
if (x$methods[2] == "exact") {
cat("\n(performed using", exact)
}
pvalue <- x$pvalue
if (!is.na(pvalue)) {
if (pvalue == 0) {
if (x$methods[2] == "permutation") {
pvalue <- paste0("< ", 1/x$nperm)
}
if (x$methods[1] == "exact") {
pvalue <- "~ 0"
}
}
}
cat("\nP-value obtained:", pvalue, "\n")
if (x$shrink) cat("\nShrinkage parameter rho was estimated at", x$rho, "\n")
cat("\nVariable importance factors")
if (p > 10) cat(" (truncated)\n") else cat("\n")
cat(format(sort(x$VIMP, decreasing = TRUE)[1:10], digits = 3),
"\n\n")
}
########################################################
# Utility functions for null distribution estimation----
dtw_ls <- function(x, mu, sigma, beta=1, log=FALSE) {
x1 <- (x - mu)/sigma
# values for dtw are only available for the
# interval -10 to 6
x1 <- as.numeric(x1 >= -10)*x1 - 10*as.numeric(x1 < -10)
x1 <- as.numeric(x1 <= 6)*x1 + 6*as.numeric(x1 > 6)
return(RMTstat::dtw(x1, beta, log)/sigma)
}
logLik <- function(param, data) {
mu <- param[1]
sigma <- param[2]
data <- data[!is.na(data)]
lL <- sum(log(dtw_ls(data, mu, sigma, log = FALSE)))
return(lL)
}
blockMatrixDiagonal <- function(matrix_list) {
dimensions <- sapply(matrix_list, nrow)
finalDimension <- sum(dimensions)
finalMatrix <- matrix(0, nrow = finalDimension, ncol = finalDimension)
index <- 1
for (k in 1:length(dimensions)) {
finalMatrix[index:(index + dimensions[k] - 1),
index:(index + dimensions[k] - 1)] <- matrix_list[[k]]
index <- index + dimensions[k]
}
return(finalMatrix)
}
JohnstoneParam <- function(p, m, n) {
s_serif <- min(n, p)
m_serif <- 0.5*(abs(n - p) - 1)
n_serif <- 0.5*(abs(m - p) - 1)
N_serif <- 2*(s_serif + m_serif + n_serif) + 1
one_third <- 1/3
gamma <- 2 * asin( sqrt( (s_serif - 0.5)/N_serif ) )
phi <- 2 * asin( sqrt( (s_serif + 2*m_serif + 0.5)/N_serif ) )
mu <- 2 * log(tan( 0.5*(phi + gamma)))
sigma <- (16/N_serif^2)^one_third * (sin(phi + gamma)^2 * sin(phi) * sin(gamma))^(-1*one_third)
return(c(mu, sigma))
}
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.