thr3ads.net - R help - [R] D'agostino test [Aug 2004]

If this information is useful, please help other people find it:
Share via:

Alexandre Bournery

2004-Aug-30 16:08 UTC

[R] D'agostino test

Hi, Does anyone know if the D'agostino test is available with R ?
Alex

Gregor.gawron@rmf.ch

2004-Aug-31 07:50 UTC

head link

[R] D'agostino test

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}}

Douglas G. Scofield

2004-Aug-31 13:29 UTC

head link

[R] D'agostino test

> -----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)
}

Diethelm Wuertz

2004-Sep-01 07:00 UTC

head link

[R] D'agostino test

Several versions of the D'Agostino Test are implemented in
"fBasics" from Rmetrics beside many other tests for normality.

Diethelm Wuertz
www.Rmetrics.org


Alexandre Bournery wrote:
> 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
>

Reasonably Related Threads

Search for more possibly parallel threads

R help - Aug 2004 - D'agostino test

[R] D'agostino test

[R] D'agostino test

[R] D'agostino test

[R] D'agostino test

Reasonably Related Threads