[R] D'agostino test

[Gregor.gawron@rmf.ch]

Alexandre,

I found once a code to calculate Shapiro_Francia and D'Agostino tests of normality written long time ago by Peter B. Mandeville. Below his code and references (It is in spanish, I assume):
Gregor


#--------------------------------------------------------------------------------------
I coded the Shapiro-Francia test which handles from 5 to 5000 observations
from 

	Patrick Royston
	1993 A Pocket-Calculator Algorithm for the Shapiro-Francia Test for
	Non-Normality: An Application to Medicine. Statistics in Medicine vol
	12, 181-184
	1991 Estimating Departure from Normality. Statistics in Medicine.
	vol 10, 1283-1293
	1983 A Simple Method for Evaluating the Shapiro-Francia W' Test of 
	Non-Normality. The Statistician 32 (1983) 297-300

and the D'Agostino tests from

	Ralph B. D'Agostino, Albert Belanger, and Ralph B. D'Agostino, Jr.
	1990 A Suggestion for Using Powerful and Informative Tests of Normality.
	The American Statistician, November 1990, Vol. 44, No. 4

for testing normality in R.    

# Lee 1996:33-34
MOMENTS <- function(data,r) sum((data-mean(data))^r)/length(data)

# Peter
SKEW <- function(data) MOMENTS(data,3)/(MOMENTS(data,2)*sqrt(MOMENTS(data,2)))

# Peter
KURTOSIS <- function(data) MOMENTS(data,4)/(MOMENTS(data,2)*MOMENTS(data,2))

# D'Agostino, Belanger and D'Agostino 1990:316-321
DAGOSTINO <- function(data){
    cat("  PRUEBAS DE D'AGOSTINO\n")
    n <- length(data)
    cat("    SIMETRIA\n")
    cat("      COEFICIENTE DE SIMETRIA:",sqrtb1 <- SKEW(data),"\n")
    if(n>8){
        y <- sqrtb1*sqrt((n+1)*(n+3)/(6*(n-2)))
        beta2 <- 3*(n*n+27*n-70)*(n+1)*(n+3)/((n-2)*(n+5)*(n+7)*(n+9))
        w <- sqrt(-1+sqrt(2*(beta2-1)))
        delta <- 1/sqrt(log(w))
        ALPHA <- sqrt(2/(w*w-1))
        cat("      ZETA CALCULADA:",zb1 <-
delta*log(y/ALPHA+sqrt((y/ALPHA)^2+1)),"\n")
        cat("      PROBABILIDAD:",2*(1-pnorm(abs(zb1))),"\n")
    }else
        cat("    LA PRUEBA DE SESGO REQUIERE POR LO MENOS 9 REPETICIONES\n")
    cat("    KURTOSIS\n")
    cat("      COEFICIENTE DE KURTOSIS:",b2 <- KURTOSIS(data),"\n")
    if(n>19){
       meanb2 <- 3*(n-1)/(n+1)
       varb2 <- 24*n*(n-2)*(n-3)/((n+1)*(n+1)*(n+3)*(n+5))
       x <- (b2-meanb2)/sqrt(varb2)
       moment <-
6*(n*n-5*n+2)/((n+7)*(n+9))*sqrt(6*(n+3)*(n+5)/(n*(n-2)*(n-3)))
       a <- 5+8/moment*(2/moment+sqrt(1+4/(moment*moment)))
       cat("      ZETA CALCULADA:",zb2 <-
(1-2/(9*a)-((1-2/a)/(1+x*sqrt(2/(a-4))))^(1/3))/sqrt(2/(9*a)),"\n")
       cat("      PROBABILIDAD:",2*(1-pnorm(abs(zb2))),"\n")
       cat("    OMNIBUS\n")
       cat("      JI-CUADRADA CALCULADA:",k2 <- zb1*zb1+zb2*zb2,"\n")
       cat("      GRADOS DE LIBERTAD: 2\n")
       cat("      PROBABILIDAD:",probji2 <- 1-pchisq(k2,2),"\n")
    }else
       cat("    LAS PRUEBAS DE KURTOSIS Y OMNIBUS REQUIEREN POR LO MENOS 20 REPETICIONES\n")
}

# Royston 1993:183-184
SHAPIRO.FRANCIA <- function(data){
    cat("  PRUEBA DE NORMALIDAD DE SHAPIRO-FRANCIA\n")
    n <- length(data)
    if(n<5 || n>5000)
        cat("    REQUIERE ENTRE 5 Y 5000 REPETICIONES\n")
    else{
        xbar <- mean(data)
        sdata <- sort(data)
        resid <- sdata-xbar
        uniform <- seq(1,n)
        np <- qnorm((uniform-0.375)/(n+0.25))
        cat("    W':",w <-
(sum(np*sdata))^2/(sum(np*np)*sum(resid*resid)),"\n")
        u <- log(n)
        v <- log(u)
        muy <- -1.2725+1.0521*(v-u)
        sigmay <- 1.0308-0.26758*(v+2/u)
        cat("    ZETA CALCULADA:",zeta <- (log(1-w)-muy)/sigmay,"\n")
        cat("    PROBABILIDAD:",probz <- 1-pnorm(zeta),"\n")
    }
}

Peter B.

Peter B. Mandeville                             mandevip at deimos.tc.uaslp.mx
Jefe del Depto. de Informática y Bioestadística rpe1531 at pasteur.fmed.uaslp.mx 
Facultad de Medicine                            Tel: 48 26-23-45 ext. 232
Universidad Autónoma de San Luis Potosí         Fax: 48 28-23-52
Av. V. Carranza 2405
Col. Los Filtros
Apartado Postal 145
San Luis Potosí, S.L.P.
78210 México

#--------------------------------------------------------------------------------------------------


-----Original Message-----
From: r-help-bounces at stat.math.ethz.ch [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Alexandre Bournery
Sent: Montag, 30. August 2004 18:08
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

______________________________________________
R-help at 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


Any information in this communication is confidential and ma...{{dropped}}