[R] Thank you!

Kenneth Cabrera krcabrer at epm.net.co
Wed Feb 23 01:21:49 CET 2000


Thank you very much for your help, your fuction runs on R without
any change, and the results are the same that your obtain with Splus,
but there still a small difference.

In SAS
W=0.960439

With your shapiro.wilk.test
W=0.9606107

But the p-values are almost the same.

Thank you very much for your help.

> Here is a function for the Shapiro-wilk test that I obtained from StatLib.
> Using this function in Splus, with  your sample, it gives the same values
> as the SAS routine. I'll copy it into R later today and see what happens.
>
> Cheers,
>
> Anne
>
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> Anne E. York
> National Marine Mammal Laboratory
> Seattle WA 98115-0070  USA
> e-mail: anne.york at noaa.gov
> Voice: +1 206-526-4039
> Fax: +1 206-526-6615
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>
> Shapiro.Wilk.test(test1)
> $W:
> [1] 0.9606107
>
> $n:
> [1] 76
>
> $p:
> [1] 0.06176898
>
> +++++++++++++++++++++++++++++++++++++++++++
> "Shapiro.Wilk.test" <- function(x)
> {
> #
> #This function is an S version of the procedure described by
> # J. P. Royston (1982) in "An Extension of Shapiro and Wilk's W Test
> # for Normality to Large Samples" from Applied Statistics, 31 no.2
> #  pp. 115-124.
> #
>         n <- length(x)
>         index <- 1:n
>         m <- qnorm((index - 0.375)/(n + 0.25))
>         y <- sort(x)
>         mu <- mean(y)
>         SSq <- sum((y - mu)^2)
>         astar <- 2 * m
>         ends <- c(1, n)
>         astar.p <- astar[ - ends]
>         if(n <= 20)
>                 m <- n - 1
>         else m <- n
>         if(m < 20)
>                 aa <- gamma(0.5 * (m + 1))/(sqrt(2) * gamma(0.5 * m + 1))
>         else {
>                 f1 <- (6 * m + 7)/(6 * m + 13)
>                 f2 <- exp(1)/(m + 2)
>                 f3 <- (m + 1)/(m + 2)
>                 f3 <- f3^(m - 2)
>                 aa <- f1 * sqrt(f2 * f3)
>         }
>         astar.1 <- (aa * sum(astar.p^2))/(1 - 2 * aa)
>         astar.1 <- sqrt(astar.1)
>         astar[1] <-  - astar.1
>         astar[n] <- astar.1
>         A <- astar/sqrt(sum(astar^2))
>         W <- (sum(A * y)^2)/SSq
>         if(n <= 20) {
>                 u <- log(n) - 3
>                 lambda <- 0.118898 + 0.133414 * u + 0.327907 * u^2
>                 logmu <- -0.37542 - 0.492145 * u - 1.124332 * u^2 -
> 0.199422 *
>                         u^3
>                 logsigma <- -3.15805 + 0.729399 * u + 3.01855 * u^2 +
> 1.558776
> *
>                         u^3
>         }
>         if(n > 20) {
>                 u <- log(n) - 5
>                 lambda <- 0.480385 + 0.318828 * u - 0.0241665 * u^3 +
>                         0.00879701 * u^4 + 0.002989646 * u^5
>                 logmu <- -1.91487 - 1.37888 * u - 0.04183209 * u^2 +
> 0.1066339
> *
>                         u^3 - 0.03513666 * u^4 - 0.01504614 * u^5
>                 logsigma <- -3.73538 - 1.015807 * u - 0.331885 * u^2 +
>                         0.1773538 * u^3 - 0.01638782 * u^4 - 0.03215018 *
> u^5 +
>
>                         0.003852646 * u^6
>         }
>         mu <- exp(logmu)
>         sigma <- exp(logsigma)
>         y <- (1 - W)^lambda
>         z <- (y - mu)/sigma
>         p <- 1 - pnorm(z)
>         if(n < 7) {
>                 warning("n is too small for this program to correctly
> estimate
> p\n"
>                         )
>                 p <- NA
>         }
>         if(n > 2000) {
>                 warning("n is too large for this program to correctly
> estimate
> p\n"
>                         )
>                 p <- NA
>         }
>         out <- list(W = W, n = n, p = p)
>         out
> }
>
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>
> On Mon, 21 Feb 2000, Kenneth Cabrera wrote:
>
> >Hello:
> >
> >I want to ask about the accuracy of the Shapiro-Wilk's test.
> >
> >I use this short program in SAS
> >____________________________________________________________________________
> >data test1;
> > input x @@;
> > cards;
> > 1.00 1.70 2.13 2.13 2.03 2.50 2.00 2.87 2.40 2.20 1.47 1.70 1.70 1.50 1.80
> > 1.03 0.63 0.60 1.87 1.13 1.10 1.27 0.83 0.67 1.73 2.23 2.50 1.60 1.97 2.17
> > 2.10 0.90 0.80 2.23 0.10 0.43 0.83 0.10 0.40 0.60 1.67 1.13 1.53 1.47 0.67
> > 0.50 1.03 1.33 1.73 1.27 0.90 1.70 2.17 0.70 0.90 0.70 1.07 0.23 0.57 0.90
> > 0.67 1.30 1.03 0.33 0.70 1.47 1.53 1.07 0.60 0.40 0.27 1.53 1.43 2.13 0.87
> > 1.13
> > ;
> >run;
> >proc univariate data=test1 normal;
> >run;
> >_____________________________________________________________________________
> >
> >And I obtain the following result for the Shapiro-Wilk's test
> >
> > W:Normal   0.960439  Pr<W  0.0602
> >
> >When I use the same data in R, with the following statement:
> >
> >___________________________________________________________________________________
> >
> >test1<-c(1.00,1.70,2.13,2.13,2.03,2.50,2.00,2.87,2.40,2.20,1.47,1.70,1.70,1.50,1.80,
> >
> > 1.03,0.63,0.60,1.87,1.13,1.10,1.27,0.83,0.67,1.73,2.23,2.50,1.60,1.97,2.17,
> > 2.10,0.90,0.80,2.23,0.10,0.43,0.83,0.10,0.40,0.60,1.67,1.13,1.53,1.47,0.67,
> > 0.50,1.03,1.33,1.73,1.27,0.90,1.70,2.17,0.70,0.90,0.70,1.07,0.23,0.57,0.90,
> > 0.67,1.30,1.03,0.33,0.70,1.47,1.53,1.07,0.60,0.40,0.27,1.53,1.43,2.13,0.87,
> > 1.13)
> >shapiro.test(test1)
> >____________________________________________________________________________________
> >
> >I obtain the following answer:
> >
> >       Shapiro-Wilk normality test
> >
> >data:  test1
> >W = 0.9733, p-value = 0.1083
> >
> >The rest of the statistics are the same, mean, median, sd, etc.
> >But I don't understand the difference in the Shapiro-Wilk's test of
> >normality.
> >
> >Thank you very much for your help.
> >
> >Kenneth Cabrera
> >Universidad Nacional de Colombia, Sede Medellin
> >Facultad de Ciencias
> >ICNE
> >Instituto de Ciencias Naturales y Ecologia
> >krcabrer at epm.net.co
> >
> >
> >-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
> >r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
> >Send "info", "help", or "[un]subscribe"
> >(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
> >_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
> >

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._



More information about the R-help mailing list