Nothing
#'@title Do a permutation test for functional regression
#' @description Performs a permutation F test (Ramsay, Hooker, and Graves, 2009, p. 145)
#' for the significance of a functional covariate, and a permutation
#' likelihood ratio test. The permutation test function currently doesn't allow models
#' with multiple functional covariates, but
#' subject-level covariates are allowed.
#' @param object An object of class funreg
#' @param num.permute The number of permutations to use. Ramsay, Hooker and Graves (2009)
#' recommended ``several hundred'' (p. 145), but for a quicker initial look
#' it might suffice to use 100.
#' @param seed An optional random number seed.
#' @references
#' Onghena, P., & May, R. B. (1995). Pitfalls in computing and interpreting
#' randomization test p values: A commentary on Chen and Dunlap.
#' Behavior Research Methods, Instruments, & Computers, 27(3), 408-411. DOI: 10.3758/BF03200438.
#'
#' Phipson, Belinda and Smyth, Gordon K. (2010) Permutation P-values
#' Should Never Be Zero: Calculating Exact P-values When Permutations
#' Are Randomly Drawn. Statistical Applications
#' in Genetics and Molecular Biology: Vol. 9: Iss. 1, Article 39. DOI: 10.2202/1544-6115.1585.
#'
#' Ramsay, J. O., Hooker, G., & Graves, S. (2009). Functional
#' data analysis with R and MATLAB. NY: Springer.
#'
#' Sen, S. (2014) Permutation Tests.
#'
#' @return Returns a list with several components. First,
#' \code{pvalue.F} is the p-value for the F test. Second, \code{conf.int.for.pvalue.F} is
#' the confidence interval for estimating the p-value that would
#' be obtained from the dataset as \code{num.permute} approached infinity.
#' The idea of a confidence interval for a p-value is explained further by
#' Sen (2013), with a STATA example.
#' Third, \code{orig.F} is the F statistic calculated on the original
#' dataset. Last, \code{permuted.F} is the vector of F statistics calculated
#' on each of the random permuted datasets. Also included are \code{pvalue.LR},
#' \code{conf.int.for.pvalue.LR}, \code{orig.LR}, \code{permuted.LR} for
#' the permutation test with a likelihood ratio statistic.
#' A more conservative alternative formula for the p-value is used in
#' \code{pvalue.F.better} and \code{pvalue.LR.better}.
#' It is not obvious whether to define the p-value as the
#' proportion of permuted datasets with statistics less than or equal to
#' the original, or simply less than the original. This should usually not
#' matter, as a tie is not likely. We made the arbitrary decision to use
#' the former here because it was presented in this way in the
#' Wikipedia article for permutation tests.
#' The conservative alternative formula is
#' the number of less extreme permuted datasets plus one,
#' over the total number of datasets plus one. Adding one to the numerator
#' and denominator is suggested by some authors, partly in order to prevent
#' a nonsensical zero p-value (Onghena & May, 1995; Phipson, Belinda & Smyth, 2010).
#' @importFrom stats fitted var logLik binom.test
#'@export
funreg.permutation <- function(object,num.permute=500,seed=NULL) {
stopifnot(class(object)=="funreg");
p <- num.functional.covs.in.model(object);
if (p>1) {
stop(paste("funreg.permutation currently does not work for",
"multiple functional covariates."));
}
if (!is.null(seed)) {set.seed(seed);}
result <- object;
id <- result$data$id;
y <- result$data$response;
s <- result$data$other.covariates;
x <- result$data$x;
time <- result$data$time;
short.id <- result$subject.info[,"id"];
short.y <- result$subject.info[,"response"];
short.s <- NULL;
if (ncol(result$subject.info)>3) {
short.s <- as.matrix(result$subject.info[,4:ncol(result$subject.info),drop=FALSE]);
stopifnot(length(short.y)==nrow(short.s));
}
stopifnot(length(short.y)==length(short.id));
stopifnot(identical(as.numeric(fitted(object)$id),as.numeric(short.id)));
short.yhat <- fitted(object)$fitted;
Fstat.orig <- var(short.yhat)/mean((short.y-short.yhat)^2);
Fstat.permuted <- rep(NA,num.permute);
loglikH1 <- logLik(result$intermediate.work.results$fit.model);
loglikH1.permuted <- rep(NA,num.permute);
for (this.permutation in 1:num.permute) {
# Permute the covariates and responses once;
shuffled.indices <- sample(1:length(short.y));
shuffled.short.y <- short.y[shuffled.indices]; # do random permutation;
shuffled.y <- NA*result$data$response;
if (!is.null(short.s)) {
shuffled.short.s <- short.s[shuffled.indices,,drop=FALSE];
shuffled.s <- NA*s;
} else {
shuffled.short.s <- NULL;
shuffled.s <- NULL;
}
for (i in short.id) {
shuffled.y[which(id==i)] <- shuffled.short.y[which(short.id==i)];
if (!is.null(short.s)) {
for (j in 1:ncol(shuffled.short.s)) {
shuffled.s[which(id==i),j] <-
shuffled.short.s[which(short.id==i),j];
}
}
}
result.permuted <- redo.funreg(object=result,
id=id,
response=shuffled.y,
time=time,
other.covariates=shuffled.s,
x=x);
Fstat.permuted[this.permutation] <- var(fitted(result.permuted)$fitted)/
(mean((fitted(result.permuted)$response-fitted(result.permuted)$fitted)^2));
loglikH1.permuted[this.permutation] <-
logLik(result.permuted$intermediate.work.results$fit.model);
}
permutation.pvalue.F <- mean(Fstat.orig<=Fstat.permuted);
permutation.pvalue.F.better <- (sum(Fstat.orig<=Fstat.permuted)+1)/num.permute;
conf.int.for.pvalue.F <- binom.test(x=sum(Fstat.orig<=Fstat.permuted),
n=length(Fstat.permuted))$conf.int;
permutation.pvalue.LR <- mean(loglikH1<=loglikH1.permuted);
permutation.pvalue.LR.better <- (sum(loglikH1<=loglikH1.permuted)+1)/num.permute;
conf.int.for.pvalue.LR <- binom.test(x=sum(loglikH1<=loglikH1.permuted),
n=length(loglikH1.permuted))$conf.int;
cat(paste("Permutation F test two-sided p-value: ",permutation.pvalue.F.better),"\n");
cat(paste("Permutation LR test two-sided p-value: ",permutation.pvalue.LR.better),"\n");
return(list(pvalue.F=permutation.pvalue.F,
conf.int.for.pvalue.F=conf.int.for.pvalue.F,
pvalue.F.better=permutation.pvalue.F.better,
orig.F=Fstat.orig,
permuted.F=Fstat.permuted,
pvalue.LR=permutation.pvalue.LR,
pvalue.LR.better=permutation.pvalue.LR.better,
conf.int.for.pvalue.LR=conf.int.for.pvalue.LR,
orig.loglikH1=loglikH1,
permuted.loglikH1 = loglikH1.permuted));
}
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.