[R] D'agostino test

[Douglas G. Scofield]

> -----Original Message-----
> From: Alexandre Bournery [mailto:bournery at mnhn.fr] 
> Sent: Monday, August 30, 2004 11:08 AM
> To: R-help at stat.math.ethz.ch
> Subject: [R] D'agostino test
> 
> Hi, Does anyone know if the D'agostino test is available with R ?
> Alex
> 
> 
> 

Hi, 

See below, taken pretty much verbatim from Zar (1999)

Douglas G. Scofield, PhD        Department of Biology 
scofield at bio.indiana.edu        Indiana University
off: (812) 856-0115             1001 E. 3rd St.
fax: (812) 855-6705             Bloomington, IN  47405-3700
cell: (786) 514-9141
 
---------------------------------------------------------------

dagostino.pearson.test <- function(x) {
    # from Zar (1999), implemented by Doug Scofield,
scofield at bio.indiana.edu
    DNAME <- deparse(substitute(x))
    n <- length(x)
    x2 <- x * x
    x3 <- x * x2
    x4 <- x * x3
    # compute Z_g1
    k3 <- ((n*sum(x3)) - (3*sum(x)*sum(x2)) + (2*(sum(x)^3)/n)) /
((n-1)*(n-2))
    g1 <- k3 / sqrt(var(x)^3)
    sqrtb1 <- ((n - 2)*g1) / sqrt(n*(n - 1))
    A <- sqrtb1 * sqrt(((n + 1)*(n + 3)) / (6*(n - 2)))
    B <- (3*(n*n + 27*n - 70)*(n+1)*(n+3)) / ((n-2)*(n+5)*(n+7)*(n+9))
    C <- sqrt(2*(B - 1)) - 1
    D <- sqrt(C)
    E <- 1 / sqrt(log(D))
    F <- A / sqrt(2/(C - 1))
    Zg1 <- E * log(F + sqrt(F*F + 1))
    # compute Z_g2
    G <- (24*n*(n-2)*(n-3)) / (((n+1)^2)*(n+3)*(n+5))
    k4 <- (((n*n*n + n*n)*sum(x4)) - (4*(n*n + n)*sum(x3)*sum(x)) -
(3*(n*n - n)*sum(x2)^2) + (12*n*sum(x2)*sum(x)^2) - (6*sum(x)^4)) /
(n*(n-1)*(n-2)*(n-3))
    g2 <- k4 / var(x)^2
    H <- ((n-2)*(n-3)*abs(g2)) / ((n+1)*(n-1)*sqrt(G))
    J <- ((6*(n*n - 5*n + 2)) / ((n+7)*(n+9))) * sqrt((6*(n+3)*(n+5)) /
(n*(n-2)*(n-3)))
    K <- 6 + (8/J)*(2/J + sqrt(1 + 4/(J*J)))
    L <- (1 - 2/K) / (1 + H*sqrt(2/(K-4)))
    Zg2 <- (1 - 2/(9*K) - (L^(1/3))) / (sqrt(2/(9*K)))
    K2 <- Zg1*Zg1 + Zg2*Zg2
    pk2 <- pchisq(K2, 2, lower.tail=FALSE)
    RVAL <- list(statistic = c(K2 = K2), p.value = pk2, method =
"D'Agostino-Pearson normality test\n\nK2 is distributed as Chi-squared
with df=2", alternative = "distribution is not normal", data.name =
DNAME)
    class(RVAL) <- "htest"
    return(RVAL)
}