[R] Anscombe-Glynn, Bonett-Seier, D'Agostino

From: Lukasz Komsta <luke_at_novum.am.lublin.pl>
Date: Sat 30 Apr 2005 - 04:11:38 EST


Dear useRs,

I was searching CRAN for implementation of kurtosis and skewness tests, and found that there is some kind of lack on it.

So, I have written three functions:

  1. Anscombe-Glynn test for kurtosis
  2. Bonett-Seier test based on Geary's kurtosis (which is not widely known, but I was inspired by original paper describing it, found coincidentally in Elsevier database)
  3. D'Agostino test for skewness

These three functions are not enough to make another small package, so I am waiting for ideas about implementing it in some existing package. If there is a need, I will contact maintainer and write manpages with appropriate examples and references.

Regards,

-- 
Lukasz Komsta
Department of Medicinal Chemistry
Medical University of Lublin
6 Chodzki, 20-093 Lublin, Poland
Fax +48 81 7425165


Code:

agostino.test <- function (x, alternative=c("two.sided","less","greater"))
{
     DNAME <- deparse(substitute(x))
     x <- sort(x[complete.cases(x)])
     n <- length(x)

	s <- match.arg(alternative)
	alter <- switch(s, two.sided=0, less=1, greater=2)

     if ((n < 8 || n > 46340))
         stop("sample size must be between 8 and 46340")

	s3 <- (sum((x-mean(x))^3)/n)/(sum((x-mean(x))^2)/n)^(3/2)
	y <- s3*sqrt((n+1)*(n+3)/(6*(n-2)))
	b2 <- 3*(n*n+27*n-70)*(n+1)*(n+3)/((n-2)*(n+5)*(n+7)*(n+9))
	w <- sqrt(-1+sqrt(2*(b2-1)));
	d <- 1/sqrt(log10(w));
	a <- sqrt(2/(w*w-1));
	z <- d*log10(y/a+sqrt((y/a)^2+1));

     pval <- pnorm(z, lower.tail = FALSE)

if (alter == 0) {
	pval <- 2*pval
	if (pval > 1) pval<-2-pval
	alt <- "data have a skewness"

}
else if (alter == 1) { alt <- "data have positive skewness"
}
else { pval <- 1-pval alt <- "data have negative skewness"
}
RVAL <- list(statistic = c(g1 = s3, z = z), p.value = pval, alternative = alt, method = "D'Agostino skewness test", data.name = DNAME) class(RVAL) <- "htest" return(RVAL) } bonett.test <- function (x, alternative=c("two.sided","less","greater")) { DNAME <- deparse(substitute(x)) x <- sort(x[complete.cases(x)]) n <- length(x) s <- match.arg(alternative) alter <- switch(s, two.sided=0, less=1, greater=2) rho <- sqrt(sum((x-mean(x))^2)/n); tau <- sum(abs(x-mean(x)))/n; omega <- 13.29*(log(rho)-log(tau)); z <- sqrt(n+2)*(omega-3)/3.54; pval <- pnorm(z, lower.tail = FALSE) if (alter == 0) { pval <- 2*pval if (pval > 1) pval<-2-pval alt <- "kurtosis is not equal to 3"
}
else if (alter == 1) { alt <- "kurtosis is greater than 3"
}
else { pval <- 1-pval alt <- "kurtosis is lower than 3"
}
RVAL <- list(statistic = c(tau = tau, z = z), alternative = alt, p.value = pval, method = "Bonett-Seier kurtosis test", data.name = DNAME) class(RVAL) <- "htest" return(RVAL) } anscombe.test <- function (x, alternative=c("two.sided","less","greater")) { DNAME <- deparse(substitute(x)) x <- sort(x[complete.cases(x)]) n <- length(x) s <- match.arg(alternative) alter <- switch(s, two.sided=0, less=1, greater=2) b <- n*sum( (x-mean(x))^4 )/(sum( (x-mean(x))^2 )^2); eb2 <- 3*(n-1)/(n+1); vb2 <- 24*n*(n-2)*(n-3)/ ((n+1)^2*(n+3)*(n+5)); m3 <- (6*(n^2-5*n+2)/((n+7)*(n+9)))*sqrt((6*(n+3)*(n+5))/(n*(n-2)*(n-3))); a <- 6+(8/m3)*(2/m3+sqrt(1+4/m3)); xx <- (b-eb2)/sqrt(vb2); z <- ( 1-2/(9*a)-( (1-2/a) / (1+xx*sqrt(2/(a-4))) )^(1/3))/ sqrt(2/(9*a)); pval <- pnorm(z, lower.tail = FALSE) if (alter == 0) { pval <- 2*pval if (pval > 1) pval<-2-pval alt <- "kurtosis is not equal to 3"
}
else if (alter == 1) { alt <- "kurtosis is greater than 3"
}
else { pval <- 1-pval alt <- "kurtosis is lower than 3"
}
RVAL <- list(statistic = c(b2 = b, z = z), p.value = pval, alternative = alt, method = "Anscombe-Glynn kurtosis test", data.name = DNAME) class(RVAL) <- "htest" return(RVAL) } ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Received on Sat Apr 30 04:22:31 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:31:31 EST