The R Project SVN
R
(root)/branches/tools4pkgs/src/library/stats/tests/ks-test.Rout.save – Rev 71686
Rev 71678 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
R Under development (unstable) (2016-11-22 r71677) -- "Unsuffered Consequences"Copyright (C) 2016 The R Foundation for Statistical ComputingPlatform: x86_64-pc-linux-gnu (64-bit)R is free software and comes with ABSOLUTELY NO WARRANTY.You are welcome to redistribute it under certain conditions.Type 'license()' or 'licence()' for distribution details.R is a collaborative project with many contributors.Type 'contributors()' for more information and'citation()' on how to cite R or R packages in publications.Type 'demo()' for some demos, 'help()' for on-line help, or'help.start()' for an HTML browser interface to help.Type 'q()' to quit R.> #### Some examples of the KS and Wilcoxon tests>> ### ------ Kolmogorov Smirnov (KS) -------------->> ## unrealistic one of PR#14561> ds1 <- c(1.7,2,3,3,4,4,5,5,6,6)> ks.test(ds1, "pnorm", mean = 3.3, sd = 1.55216)One-sample Kolmogorov-Smirnov testdata: ds1D = 0.274, p-value = 0.4407alternative hypothesis: two-sidedWarning message:In ks.test(ds1, "pnorm", mean = 3.3, sd = 1.55216) :ties should not be present for the Kolmogorov-Smirnov test> # how on earth can sigma = 1.55216 be known?>> # R >= 2.14.0 allows the equally invalid> ks.test(ds1, "pnorm", mean = 3.3, sd = 1.55216, exact = TRUE)One-sample Kolmogorov-Smirnov testdata: ds1D = 0.274, p-value = 0.3715alternative hypothesis: two-sidedWarning message:In ks.test(ds1, "pnorm", mean = 3.3, sd = 1.55216, exact = TRUE) :ties should not be present for the Kolmogorov-Smirnov test>> ## Try out the effects of rounding> set.seed(123)> ds2 <- rnorm(1000)> ks.test(ds2, "pnorm") # exact = FALSE is default for n = 1000One-sample Kolmogorov-Smirnov testdata: ds2D = 0.019416, p-value = 0.8452alternative hypothesis: two-sided> ks.test(ds2, "pnorm", exact = TRUE)One-sample Kolmogorov-Smirnov testdata: ds2D = 0.019416, p-value = 0.8379alternative hypothesis: two-sided> ## next two are still close> ks.test(round(ds2, 2), "pnorm")One-sample Kolmogorov-Smirnov testdata: round(ds2, 2)D = 0.019169, p-value = 0.856alternative hypothesis: two-sidedWarning message:In ks.test(round(ds2, 2), "pnorm") :ties should not be present for the Kolmogorov-Smirnov test> ks.test(round(ds2, 2), "pnorm", exact = TRUE)One-sample Kolmogorov-Smirnov testdata: round(ds2, 2)D = 0.019169, p-value = 0.8489alternative hypothesis: two-sidedWarning message:In ks.test(round(ds2, 2), "pnorm", exact = TRUE) :ties should not be present for the Kolmogorov-Smirnov test> # now D has doubled, but p-values remain similar (if very different from ds2)> ks.test(round(ds2, 1), "pnorm")One-sample Kolmogorov-Smirnov testdata: round(ds2, 1)D = 0.03674, p-value = 0.1344alternative hypothesis: two-sidedWarning message:In ks.test(round(ds2, 1), "pnorm") :ties should not be present for the Kolmogorov-Smirnov test> ks.test(round(ds2, 1), "pnorm", exact = TRUE)One-sample Kolmogorov-Smirnov testdata: round(ds2, 1)D = 0.03674, p-value = 0.1311alternative hypothesis: two-sidedWarning message:In ks.test(round(ds2, 1), "pnorm", exact = TRUE) :ties should not be present for the Kolmogorov-Smirnov test>>> ### ------ Wilkoxon (Mann Whitney) -------------->> options(nwarnings = 1000)> (alts <- setNames(, eval(formals(stats:::wilcox.test.default)$alternative)))two.sided less greater"two.sided" "less" "greater"> x0 <- 0:4> (x.set <- list(s0 = lapply(x0, function(m) 0:m),+ s. = lapply(x0, function(m) c(1e-9, seq_len(m)))))$s0$s0[[1]][1] 0$s0[[2]][1] 0 1$s0[[3]][1] 0 1 2$s0[[4]][1] 0 1 2 3$s0[[5]][1] 0 1 2 3 4$s.$s.[[1]][1] 1e-09$s.[[2]][1] 1e-09 1e+00$s.[[3]][1] 1e-09 1e+00 2e+00$s.[[4]][1] 1e-09 1e+00 2e+00 3e+00$s.[[5]][1] 1e-09 1e+00 2e+00 3e+00 4e+00> stats <- setNames(nm = c("statistic", "p.value", "conf.int", "estimate"))>> ## Even with conf.int = TRUE, do not want errors :> RR <-+ lapply(x.set, ## for all data sets+ function(xs)+ lapply(alts, ## for all three alternatives+ function(alt)+ lapply(xs, function(x)+ ## try(+ wilcox.test(x, exact=TRUE, conf.int=TRUE, alternative = alt)+ ## )+ )))There were 52 warnings (use warnings() to see them)> length(ww <- warnings()) # 52 (or 43 for x0 <- 0:3)[1] 52> unique(ww) # 4 different onesWarning messages:1: In wilcox.test.default(x, exact = TRUE, conf.int = TRUE, ... :cannot compute exact p-value with zeroes2: In wilcox.test.default(x, exact = TRUE, conf.int = TRUE, ... :cannot compute exact confidence interval with zeroes3: cannot compute confidence interval when all observations are zero or tied4: In wilcox.test.default(x, exact = TRUE, conf.int = TRUE, ... :requested conf.level not achievable>> cc <- lapply(RR, function(A) lapply(A, function(bb) lapply(bb, class)))> table(unlist(cc))htest30> ## in R <= 3.3.1, with try( .. ) above, we got> ## htest try-error> ## 23 7> uc <- unlist(cc[["s0"]]); noquote(names(uc)[uc != "htest"]) ## these 7 cases :character(0)> ## two.sided1 two.sided2 two.sided3> ## less1 less2> ## greater1 greater2>> ##--- How close are the stats of (0:m) to those of (eps, 1:m) ------------>> ## a version that still works with above try(.) and errors there:> getC <- function(L, C) if(inherits(L,"try-error")) c(L) else L[[C]]> stR <- lapply(stats, function(COMP)+ lapply(RR, function(A)+ lapply(A, function(bb)+ lapply(bb, getC, C=COMP) )))>> ## a) P-value> pv <- stR[["p.value"]]> ## only the first is NaN, all others in [0,1]:> sapply(pv$s0, unlist)two.sided less greater[1,] NaN 1.0000000 1.00000000[2,] 1.0000000 0.9772499 0.50000000[3,] 0.3710934 0.9631809 0.18554668[4,] 0.1814492 0.9693156 0.09072460[5,] 0.1003482 0.9776951 0.05017412> sapply(pv$s., unlist) # not really close, but ..two.sided less greater[1,] 1.0000 1 0.50000[2,] 0.5000 1 0.25000[3,] 0.2500 1 0.12500[4,] 0.1250 1 0.06250[5,] 0.0625 1 0.03125>> pv$s0$two.sided[1] <- 1 ## artificially> stopifnot(all.equal(pv$s0, pv$s., tol = 0.5 + 1e-6), # seen 0.5+ ## "less" are close:+ all.equal(unlist(pv[[c("s0","less")]]),+ unlist(pv[[c("s.","less")]]), tol = 0.03),+ 0 <= unlist(pv), unlist(pv) <= 1) # <- no further NA ..> ## b)> sapply(stR[["statistic"]], unlist)s0 s.two.sided.V 0 1two.sided.V 1 3two.sided.V 3 6two.sided.V 6 10two.sided.V 10 15less.V 0 1less.V 1 3less.V 3 6less.V 6 10less.V 10 15greater.V 0 1greater.V 1 3greater.V 3 6greater.V 6 10greater.V 10 15> ## Conf.int.:> ## c)> sapply(stR[["estimate" ]], unlist)s0 s.two.sided1 NaN 1.0e-09two.sided.midrange 1.0 5.0e-01two.sided.(pseudo)median 1.5 1.0e+00two.sided.(pseudo)median 2.0 1.5e+00two.sided.(pseudo)median 2.5 2.0e+00less1 NaN 1.0e-09less.midrange 1.0 5.0e-01less.(pseudo)median 1.5 1.0e+00less.(pseudo)median 2.0 1.5e+00less.(pseudo)median 2.5 2.0e+00greater1 NaN 1.0e-09greater.midrange 1.0 5.0e-01greater.(pseudo)median 1.5 1.0e+00greater.(pseudo)median 2.0 1.5e+00greater.(pseudo)median 2.5 2.0e+00> ## d) confidence interval> formatCI <- function(ci)+ sprintf("[%g, %g] (%g%%)", ci[[1]], ci[[2]],+ round(100*attr(ci,"conf.level")))> nx <- length(x0)> noquote(vapply(stR[["conf.int"]], function(ss)+ vapply(ss, function(alt) vapply(alt, formatCI, ""), character(nx)),+ matrix("", nx, length(alts)))), , s0two.sided less greater[1,] [NaN, NaN] (0%) [-Inf, NaN] (0%) [NaN, Inf] (0%)[2,] [NaN, NaN] (0%) [-Inf, NaN] (0%) [NaN, Inf] (0%)[3,] [1.5, 1.5] (0%) [-Inf, 1.5] (0%) [1.50001, Inf] (20%)[4,] [1.00007, 2.99993] (60%) [-Inf, 2.00002] (60%) [1.00009, Inf] (80%)[5,] [1.00006, 3.99994] (80%) [-Inf, 3.00002] (80%) [1.00009, Inf] (90%), , s.two.sided less greater[1,] [1e-09, 1e-09] (0%) [-Inf, 1e-09] (50%) [1e-09, Inf] (50%)[2,] [1e-09, 1] (50%) [-Inf, 1] (75%) [1e-09, Inf] (75%)[3,] [1e-09, 2] (75%) [-Inf, 2] (87%) [1e-09, Inf] (87%)[4,] [1e-09, 3] (88%) [-Inf, 3] (95%) [1e-09, Inf] (95%)[5,] [1e-09, 4] (95%) [-Inf, 4] (95%) [1e-09, Inf] (95%)>>> ##-------- 2-sample tests (working unchanged) ------------------>> R2 <- lapply(alts, ## for all three alternatives+ function(alt)+ lapply(seq_along(x0), function(k)+ wilcox.test(x = x.set$s0[[k]], y = x.set$s.[[k]],+ exact=TRUE, conf.int=TRUE, alternative = alt)))There were 27 warnings (use warnings() to see them)> length(w2 <- warnings()) # 27[1] 27> unique(w2) # 3 different onesWarning messages:1: In wilcox.test.default(x = x.set$s0[[k]], y = x.set$s.[[k]], ... :Requested conf.level not achievable2: In wilcox.test.default(x = x.set$s0[[k]], y = x.set$s.[[k]], ... :cannot compute exact p-value with ties3: In wilcox.test.default(x = x.set$s0[[k]], y = x.set$s.[[k]], ... :cannot compute exact confidence intervals with ties>> table(uc2 <- unlist(c2 <- lapply(R2, function(A) lapply(A, class))))htest15> stopifnot(uc2 == "htest")>> stR2 <- lapply(stats,+ function(COMP)+ lapply(R2, function(A) lapply(A, getC, C=COMP)))>> lapply(stats[-3], ## -3: "conf.int" separately+ function(ST) sapply(stR2[[ST]], unlist))$statistictwo.sided less greaterW 0.0 0.0 0.0W 1.5 1.5 1.5W 4.0 4.0 4.0W 7.5 7.5 7.5W 12.0 12.0 12.0$p.valuetwo.sided less greater[1,] 1 0.5 1.0000000[2,] 1 0.5 0.7928919[3,] 1 0.5 0.6734524[4,] 1 0.5 0.6156105[5,] 1 0.5 0.5837406$estimatetwo.sided less greaterdifference in location -1.000000e-09 -1.000000e-09 -1.000000e-09difference in location -4.467848e-05 -4.467848e-05 -4.467848e-05difference in location -4.692131e-05 -4.692131e-05 -4.692131e-05difference in location 1.937902e-05 1.937902e-05 1.937902e-05difference in location 3.741417e-05 3.741417e-05 3.741417e-05>> noquote(sapply(stR2[["conf.int"]], function(.) vapply(., formatCI, "")))two.sided less greater[1,] [-1e-09, -1e-09] (0%) [-Inf, -1e-09] (50%) [-1e-09, Inf] (50%)[2,] [-1, 1] (95%) [-Inf, 1] (95%) [-1, Inf] (95%)[3,] [-2, 2] (95%) [-Inf, 2] (95%) [-2, Inf] (95%)[4,] [-3, 3] (95%) [-Inf, 2.00005] (95%) [-2.00003, Inf] (95%)[5,] [-2.99998, 2.99996] (95%) [-Inf, 2.00005] (95%) [-2.00006, Inf] (95%)>>> proc.time()user system elapsed0.337 0.014 0.340