# R/adtest.R In robCompositions: Compositional Data Analysis

#' Anderson-Darling Normality Tests
#'
#' This function provides three kinds of Anderson-Darling Normality Tests
#' (Anderson and Darling, 1952).
#'
#' Three version of the test are implemented (univariate, angle and radius
#' test) and it depends on the data which test is chosen.
#'
#' If the data is univariate the univariate Anderson-Darling test for normality
#' is applied.
#'
#' If the data is bivariate the angle Anderson-Darling test for normality is
#' performed out.
#'
#' If the data is multivariate the radius Anderson-Darling test for normality
#' is used.
#'
#' If \sQuote{locscatt} is equal to \dQuote{robust} then within the procedure,
#' robust estimates of mean and covariance are provided using \sQuote{covMcd()}
#' from package robustbase.
#'
#' To provide estimates for the corresponding p-values, i.e. to compute the
#' probability of obtaining a result at least as extreme as the one that was
#' actually observed under the null hypothesis, we use Monte Carlo techniques
#' where we check how often the statistic of the underlying data is more
#' extreme than statistics obtained from simulated normal distributed data with
#' the same (column-wise-) mean(s) and (co)variance.
#'
#' @param x either a numeric vector, or a data.frame, or a matrix
#' @param R Number of Monte Carlo simulations to obtain p-values
#' @param locscatt standard for classical estimates of mean and (co)variance.
#' robust for robust estimates using \sQuote{covMcd()} from package robustbase
#' @return \item{statistic }{The result of the corresponding test statistic}
#' \item{method }{The chosen method (univariate, angle or radius)}
#' \item{p.value }{p-value}
#' @author Karel Hron, Matthias Templ
#' @references Anderson, T.W. and Darling, D.A. (1952) Asymptotic theory of
#' certain goodness-of-fit criteria based on stochastic processes. \emph{Annals
#' of Mathematical Statistics}, \bold{23} 193-212.
#' @keywords htest
#' @export
#' @importFrom MASS mvrnorm
#' @importFrom robustbase covMcd
#' @examples
#'
#' data(machineOperators)
#' x <- machineOperators
#'
#DNAME <- deparse(substitute(x))
if( R < 1 ) stop("choose a higher value for R")
if( R < 50 ) warning("maybe, the estimation of the p-value(s) is not accurate; choose a higher value for R")
cv <- function(x, type) {
classical <- function(x){
if( (length(dim(x)) < 1) | is.vector(x) ){
list(mean=mean(as.numeric(x), na.rm=TRUE), varmat=var(x, na.rm=TRUE))
} else {
list(mean=colMeans(x, na.rm=TRUE), varmat=cov(x))
}
}

robust <- function(x){
if( (length(dim(x)) < 1) | is.vector(x) ){
} else {
v <- robustbase::covMcd(x)
list(mean=v$center, varmat=v$cov)
}
}
switch(type,
standard = classical(x),
robust = robust(x))
}

centre <- function(x, type) {
switch(type,
mean = mean(as.numeric(x)),
median = median(as.numeric(x)),
trimmed = mean(as.numeric(x), trim = .1))
}
if(locscatt=="standard") location <- "mean" else location <- "median"

### 1-dim:
if((length(dim(x)) < 1) | is.vector(x) ){
x <- as.vector(x)
n <- length(x)
if (n < 8 )
stop("sample size must be greater than 7")
stat=function(x, N=n, location="mean"){
x <- sort(x[complete.cases(x)])
estCv <- cv(x, locscatt)
p <- pnorm((x - estCv$mean)/sqrt(estCv$varmat))
h <- (2 * seq(1:N) - 1) * (log(p) + log(1 - rev(p)))
A <- (25/N^2-4/N-1)*(centre(h, location)+N)
}
A=stat(x, location=location)
estCv <- cv(x, locscatt)
mv <- estCv$mean varmat <- estCv$varmat
n=length(c(x))
p <- sapply(X=1:R, FUN=function(X,...){ stat(rnorm(n, mv, sqrt(varmat)), location=location) })
pvalue=mean(p>A)
RVAL <- list(statistic = c(A = A), method = "A-D univariate normality test",
p.value=pvalue)
} else if (dim(x)[2] == 2){
### 2-dim:
n <- nrow(x)
stat=function(x, N=n, location="mean"){
estCv <- cv(x, locscatt)
varmat <- estCv$varmat mu <- estCv$mean
u <- (1/sqrt(det(varmat)))*((x[,1]-mu[1])*sqrt(varmat[2,2])-(x[,2]-mu[2])*(varmat[1,2]/sqrt(varmat[2,2])))
v <- (x[,2]-mu[2])/sqrt(varmat[2,2])
teta <- atan(v/u)+(1-sign(u))*pi/2+(1+sign(u))*(1-sign(v))*pi/2
z=teta/(2*pi)
p=sort(z)
h <- (2 * seq(1:N) - 1) * (log(p) + log(1 - rev(p)))
A <- -N-centre(h, location)
}
A=stat(x, location=location)
estCv <- cv(x, locscatt)
varmat <- estCv$varmat mv <- estCv$mean
p <- sapply(X=1:R, FUN=function(X,...){ stat(mvrnorm(n, mv, sqrt(varmat)), location=location) })
pvalue=mean(p>A)
RVAL <- list(statistic = c(A = A), method = "A-D bivariate normality test",
p.value=pvalue)
} else {
### >= 3:
n <- nrow(x)
stat=function(x, N=n, location="mean"){
estCv <- cv(x, locscatt)
#varmat <- var(x)
#mu <- apply(x,2,mean)
u <- mahalanobis(x, center=estCv$mean, cov=estCv$varmat)
z <- pchisq(u,ncol(x))
p=sort(z)
h <- (2 * seq(1:N) - 1) * (log(p) + log(1 - rev(p)))
A <- -N-centre(h, location)
#par(mfrow=c(1,2)); plot(p) ; plot(h)
}
A=stat(x, location=location)
#print(paste("A =", A))
estCv <- cv(x, locscatt)
# estCv$mean=colMeans(x) #weg #estCv$varmat=var(x) #weg
p=numeric(R)
n=nrow(x)
p <- sapply(X=1:R, FUN=function(X,...){ stat(mvrnorm(n, estCv$mean, estCv$varmat), location=location) })
pvalue=mean(p>A)
RVAL <- list(statistic = A, method = "A-D radius test",
p.value=pvalue)
}

class(RVAL) <- "htest"
return(RVAL)

}

######################################################################
#Cramer-vom Mises
#
#cvmtest=function (x)
#{
#    DNAME <- deparse(substitute(x))
#    x <- sort(x[complete.cases(x)])
#    n <- length(x)
#    if (n < 8)
#        stop("sample size must be greater than 7")
#    p <- pnorm((x - mean(x))/sd(x))
#    W <- (1/(12 * n) + sum((p - (2 * seq(1:n) - 1)/(2 * n))^2))*((2*n+1)/(2*n))
#    RVAL <- list(statistic = c(W = W), method = "Cramer-von Mises normality test",
#        data.name = DNAME)
#    class(RVAL) <- "htest"
#    return(RVAL)
#}
#
#######################################################################
##Watson test
#
#wattest=function (x)
#{
#    DNAME <- deparse(substitute(x))
#    x <- sort(x[complete.cases(x)])
#    n <- length(x)
#    if (n < 8)
#        stop("sample size must be greater than 7")
#    p <- pnorm((x - mean(x))/sd(x))
#    W <- (1/(12 * n) + sum((p - (2 * seq(1:n) - 1)/(2 * n))^2))*((2*n+1)/(2*n))
#    WW <- W-((2*n+1)/2)*(mean(p)-1/2)^2
#    RVAL <- list(statistic = c(WW = WW), method = "Watson normality test",
#        data.name = DNAME)
#    class(RVAL) <- "htest"
#    return(RVAL)
#}


## Try the robCompositions package in your browser

Any scripts or data that you put into this service are public.

robCompositions documentation built on Aug. 25, 2023, 5:13 p.m.