| Line 19... |
Line 19... |
| 19 |
ks.test(round(ds2, 2), "pnorm")
|
19 |
ks.test(round(ds2, 2), "pnorm")
|
| 20 |
ks.test(round(ds2, 2), "pnorm", exact = TRUE)
|
20 |
ks.test(round(ds2, 2), "pnorm", exact = TRUE)
|
| 21 |
# now D has doubled, but p-values remain similar (if very different from ds2)
|
21 |
# now D has doubled, but p-values remain similar (if very different from ds2)
|
| 22 |
ks.test(round(ds2, 1), "pnorm")
|
22 |
ks.test(round(ds2, 1), "pnorm")
|
| 23 |
ks.test(round(ds2, 1), "pnorm", exact = TRUE)
|
23 |
ks.test(round(ds2, 1), "pnorm", exact = TRUE)
|
| - |
|
24 |
summary(w1 <- warnings())
|
| - |
|
25 |
|
| 24 |
|
26 |
|
| 25 |
|
- |
|
| 26 |
### ------ Wilkoxon (Mann Whitney) --------------
|
27 |
### ------ Wilcoxon (Mann Whitney) --------------
|
| 27 |
|
28 |
|
| 28 |
options(nwarnings = 1000)
|
29 |
options(nwarnings = 1000)
|
| 29 |
(alts <- setNames(, eval(formals(stats:::wilcox.test.default)$alternative)))
|
30 |
(alts <- setNames(, eval(formals(stats:::wilcox.test.default)$alternative)))
|
| 30 |
x0 <- 0:4
|
31 |
x0 <- setNames(0:4, paste0("n_", 1L + 0:4))
|
| 31 |
(x.set <- list(s0 = lapply(x0, function(m) 0:m),
|
32 |
str(x.set <- list(s0 = lapply(x0, function(m) 0:m),
|
| 32 |
s. = lapply(x0, function(m) c(1e-9, seq_len(m)))))
|
33 |
s. = lapply(x0, function(m) c(1e-9, seq_len(m)))))
|
| 33 |
stats <- setNames(nm = c("statistic", "p.value", "conf.int", "estimate"))
|
34 |
stats <- setNames(nm = c("statistic", "p.value", "conf.int", "estimate"))
|
| 34 |
|
- |
|
| 35 |
## Even with conf.int = TRUE, do not want errors :
|
35 |
dig.r <- c(Inf, 20, 15, 10, 7, 5); names(dig.r) <- paste0("d.rank=", dig.r)
|
| 36 |
RR <-
|
36 |
RL <- lapply(dig.r, function(dig.rank)
|
| 37 |
lapply(x.set, ## for all data sets
|
37 |
lapply(x.set, ## for all data sets
|
| 38 |
function(xs)
|
38 |
function(xs)
|
| 39 |
lapply(alts, ## for all three alternatives
|
39 |
lapply(alts, ## for all three alternatives
|
| 40 |
function(alt)
|
40 |
function(alt)
|
| 41 |
lapply(xs, function(x)
|
41 |
lapply(xs, function(x)
|
| 42 |
## try(
|
42 |
## try( nolonger: Even with conf.int = TRUE, do not want errors :
|
| 43 |
wilcox.test(x, exact=TRUE, conf.int=TRUE, alternative = alt)
|
43 |
wilcox.test(x, exact=TRUE, conf.int=TRUE, alternative = alt,
|
| - |
|
44 |
digits.rank = dig.rank)
|
| 44 |
## )
|
45 |
## )
|
| 45 |
)))
|
46 |
))))
|
| 46 |
length(ww <- warnings()) # 52 (or 43 for x0 <- 0:3)
|
47 |
length(ww <- warnings()) # 1 ; was 52 (or 43 for x0 <- 0:3)
|
| - |
|
48 |
if(!identical(w1, ww)) # typically they *are* identical ==> no new warnings
|
| - |
|
49 |
summary(ww)
|
| 47 |
unique(ww) # 4 different ones
|
50 |
## unique(ww) # (were 4 different ones)
|
| - |
|
51 |
|
| - |
|
52 |
## RL is very nested tree of "htest" results
|
| - |
|
53 |
str(RL, max.level = 2)
|
| - |
|
54 |
str(RL[["d.rank=7"]]$s.$two.sided, max.level = 1) # 5 x {list of 9, class "htest"}
|
| - |
|
55 |
str(RL[["d.rank=7"]]$s.$two.sided [1:2], give.attr = FALSE)
|
| 48 |
|
56 |
|
| 49 |
cc <- lapply(RR, function(A) lapply(A, function(bb) lapply(bb, class)))
|
57 |
cc <- lapply(RL, function(A) lapply(A, function(B) lapply(B, function(bb) lapply(bb, class))))
|
| 50 |
table(unlist(cc))
|
58 |
##N table(unlist(cc))
|
| 51 |
## in R <= 3.3.1, with try( .. ) above, we got
|
59 |
## in R <= 3.3.1, with try( .. ) above and only digits.rank=Inf, we got
|
| 52 |
## htest try-error
|
60 |
## htest try-error
|
| 53 |
## 23 7
|
61 |
## 23 7
|
| - |
|
62 |
stopifnot(identical("htest", unique(unlist(cc))))
|
| 54 |
uc <- unlist(cc[["s0"]]); noquote(names(uc)[uc != "htest"]) ## these 7 cases :
|
63 |
##N uc <- unlist(cc[["s0"]]); noquote(names(uc)[uc != "htest"]) ## these 7 cases :
|
| 55 |
## two.sided1 two.sided2 two.sided3
|
64 |
## two.sided1 two.sided2 two.sided3
|
| 56 |
## less1 less2
|
65 |
## less1 less2
|
| 57 |
## greater1 greater2
|
66 |
## greater1 greater2
|
| 58 |
|
67 |
|
| 59 |
##--- How close are the stats of (0:m) to those of (eps, 1:m) ------------
|
68 |
##--- How close are the stats of (0:m) to those of (eps, 1:m) ------------
|
| 60 |
|
69 |
|
| 61 |
## a version that still works with above try(.) and errors there:
|
70 |
## Contained errors in older versions:
|
| 62 |
getC <- function(L, C) if(inherits(L,"try-error")) c(L) else L[[C]]
|
- |
|
| 63 |
stR <- lapply(stats, function(COMP)
|
71 |
stR <- lapply(stats, function(COMP)
|
| 64 |
lapply(RR, function(A)
|
72 |
lapply(RL, function(A)
|
| 65 |
lapply(A, function(bb)
|
73 |
lapply(A, function(B)
|
| - |
|
74 |
lapply(B, function(bb)
|
| 66 |
lapply(bb, getC, C=COMP) )))
|
75 |
lapply(bb, `[[`, COMP) ))))
|
| - |
|
76 |
if(interactive()) {
|
| - |
|
77 |
str(stR, max.level = 2)
|
| - |
|
78 |
str(stR$p.value$`d.rank=7`, give.attr = FALSE)
|
| - |
|
79 |
}
|
| 67 |
|
80 |
|
| - |
|
81 |
## short cut & speedup:
|
| - |
|
82 |
unList <- function(x) .Internal(unlist(x, FALSE, TRUE))
|
| - |
|
83 |
## recurs. use.nms
|
| 68 |
## a) P-value
|
84 |
## a) P-value s
|
| 69 |
pv <- stR[["p.value"]]
|
85 |
pv <- stR$p.value
|
| - |
|
86 |
|
| - |
|
87 |
str(apv <- simplify2array(
|
| - |
|
88 |
lapply(pv, \(L0) simplify2array(
|
| - |
|
89 |
lapply(L0, \(L) simplify2array(lapply(L, unList))))))
|
| - |
|
90 |
)
|
| - |
|
91 |
## num [1:5, 1:3, 1:2, 1:6] 1 1 0.5 0.25 0.125 ...
|
| - |
|
92 |
## - attr(*, "dimnames")=List of 4
|
| - |
|
93 |
## ..$ : chr [1:5] "n_1" "n_2" "n_3" "n_4" ...
|
| 70 |
## only the first is NaN, all others in [0,1]:
|
94 |
## ..$ : chr [1:3] "two.sided" "less" "greater"
|
| - |
|
95 |
## ..$ : chr [1:2] "s0" "s."
|
| - |
|
96 |
## ..$ : chr [1:6] "d.rank=Inf" "d.rank=20" "d.rank=15" "d.rank=10" ...
|
| - |
|
97 |
object.size(apv) # 3144 bytes
|
| - |
|
98 |
object.size(dpv <- as.data.frame.table(apv)) # 8176 bytes
|
| - |
|
99 |
str(dpv)
|
| - |
|
100 |
str(ps. <- apv[,,"s.",])
|
| 71 |
sapply(pv$s0, unlist)
|
101 |
str(ps0 <- apv[,,"s0",])
|
| - |
|
102 |
(epsC <- .Machine$double.eps) # 2.2e-16
|
| - |
|
103 |
allEQ <- function(x,y) all.equal(x, y, tolerance = 8*epsC)
|
| - |
|
104 |
stopifnot(exprs = {
|
| 72 |
sapply(pv$s., unlist) # not really close, but ..
|
105 |
0 <= unlist(apv) ; unlist(apv) <= 1 # (no NA, ..)
|
| - |
|
106 |
## ps0 does *not* vary with digits.rank :
|
| - |
|
107 |
identical(dim(print(ups0 <- drop(unique(ps0, MARGIN = 3)))), c(5L, 3L))
|
| - |
|
108 |
is.function(uNam <- function(x) `names<-`(x, NULL))
|
| - |
|
109 |
allEQ(uNam(ups0[,"two.sided"]), c(1, 2^(0:-3)))
|
| - |
|
110 |
allEQ(uNam(ups0[, "greater" ]), 2^(0:-4))
|
| - |
|
111 |
allEQ(uNam(ups0[, "less" ]), c(1,1,1,1,1))
|
| - |
|
112 |
## NB: This has changed from R 4.6.0 (now using digits.zap = digits.rank):
|
| - |
|
113 |
## ps., eps=1e-9: two possible results, boundary for 1e-9 between digits.rank = 7 and 10
|
| - |
|
114 |
identical(dim(print(ups. <- unique(ps., MARGIN = 3))), c(5L, 3L, 2L))
|
| - |
|
115 |
allEQ(ps.[,,"d.rank=Inf"], ps.[,,"d.rank=10"])
|
| - |
|
116 |
is.function(uNam <- function(x) `dimnames<-`(x, NULL))
|
| - |
|
117 |
allEQ(uNam(ups.[,"two.sided", ]), 2^cbind(0:-4, c(0L,0:-3)))
|
| - |
|
118 |
allEQ(uNam(ups.[, "greater" , ]), cbind(2^-(1:5), 2^-c(1,1:4)))
|
| - |
|
119 |
allEQ(uNam(ups.[, "less" , ]), cbind(1, rep(1, 5)))
|
| - |
|
120 |
})
|
| - |
|
121 |
|
| 73 |
|
122 |
|
| 74 |
pv$s0$two.sided[1] <- 1 ## artificially
|
123 |
## b) Statistic - Estimate
|
| - |
|
124 |
|
| 75 |
stopifnot(all.equal(pv$s0, pv$s., tolerance = 0.5 + 1e-6), # seen 0.5
|
125 |
(astat <- simplify2array(
|
| 76 |
## "less" are close:
|
126 |
lapply(stR[["statistic"]],
|
| 77 |
all.equal(unlist(pv[[c("s0","less")]]),
|
127 |
\(L0) simplify2array(
|
| - |
|
128 |
lapply(L0, \(L) simplify2array(sapply(L, unList)))))))
|
| - |
|
129 |
|
| - |
|
130 |
(aEst <- simplify2array(
|
| 78 |
unlist(pv[[c("s.","less")]]),
|
131 |
lapply(stR[["estimate"]],
|
| 79 |
tolerance = 0.03),
|
132 |
\(L0) simplify2array(
|
| 80 |
0 <= unlist(pv), unlist(pv) <= 1) # <- no further NA ..
|
133 |
lapply(L0, \(L) simplify2array(lapply(L, unList)))))))
|
| 81 |
## b)
|
134 |
str(aEst)
|
| - |
|
135 |
## num [1:5, 1:3, 1:2, 1:6] 0 0.5 1 1.5 2 0 0.5 1 1.5 2 ...
|
| 82 |
sapply(stR[["statistic"]], unlist)
|
136 |
## - attr(*, "dimnames")=List of 4
|
| - |
|
137 |
## ..$ : chr [1:5] "n_1.(pseudo)median" "n_2.(pseudo)median" "n_3.(pseudo)median" "n_4.(pseudo)median" ...
|
| - |
|
138 |
## ..$ : chr [1:3] "two.sided" "less" "greater"
|
| - |
|
139 |
## ..$ : chr [1:2] "s0" "s."
|
| - |
|
140 |
## ..$ : chr [1:6] "d.rank=Inf" "d.rank=20" "d.rank=15" "d.rank=10" ...
|
| - |
|
141 |
|
| 83 |
## Conf.int.:
|
142 |
## Confidence Interval :
|
| 84 |
## c)
|
143 |
## c) C.I. extraction
|
| 85 |
sapply(stR[["estimate" ]], unlist)
|
144 |
CI <- stR[["conf.int"]]
|
| - |
|
145 |
str(CI, max.level = 2)
|
| - |
|
146 |
|
| - |
|
147 |
if(FALSE) { # printing chCI below is much nicer
|
| - |
|
148 |
str(aCI <- simplify2array(lapply(CI,
|
| - |
|
149 |
function(A) simplify2array(lapply(A,
|
| - |
|
150 |
function(B) simplify2array(lapply(B,
|
| - |
|
151 |
function(C) t(simplify2array(C)))))))))
|
| - |
|
152 |
## num [1:5, 1:2, 1:3, 1:2, 1:6] 0 0 0 0 0 ...
|
| - |
|
153 |
## - attr(*, "dimnames")=List of 5
|
| - |
|
154 |
## ..$ : chr [1:5] "n_1" "n_2" "n_3" "n_4" ...
|
| - |
|
155 |
## ..$ : NULL
|
| - |
|
156 |
## ..$ : chr [1:3] "two.sided" "less" "greater"
|
| - |
|
157 |
## ..$ : chr [1:2] "s0" "s."
|
| - |
|
158 |
## ..$ : chr [1:6] "d.rank=Inf" "d.rank=20" "d.rank=15" "d.rank=10" ...
|
| - |
|
159 |
aCI # has *lost* the conf.level which is varying !!
|
| - |
|
160 |
}
|
| - |
|
161 |
|
| 86 |
## d) confidence interval
|
162 |
## d) confidence interval __printing__
|
| 87 |
formatCI <- function(ci)
|
163 |
formatCI <- function(ci)
|
| 88 |
sprintf("[%g, %g] (%g%%)", ci[[1]], ci[[2]],
|
164 |
sprintf("[%g, %g] (%g%%)", ci[[1]], ci[[2]],
|
| 89 |
round(100*attr(ci,"conf.level")))
|
165 |
round(100*attr(ci,"conf.level")))
|
| 90 |
nx <- length(x0)
|
166 |
nx <- length(x0)
|
| - |
|
167 |
Fval <- array("", dim = c(nx, length(alts), length(x.set)))
|
| 91 |
noquote(vapply(stR[["conf.int"]], function(ss)
|
168 |
(chCI <- noquote(vapply(CI, function(D)
|
| - |
|
169 |
vapply(D, function(ss)
|
| 92 |
vapply(ss, function(alt) vapply(alt, formatCI, ""), character(nx)),
|
170 |
vapply(ss, function(alt) vapply(alt, formatCI, ""), character(nx)),
|
| 93 |
matrix("", nx, length(alts))))
|
171 |
Fval[,,1]),
|
| - |
|
172 |
Fval)
|
| - |
|
173 |
))
|
| 94 |
|
174 |
|
| 95 |
|
175 |
|
| 96 |
##-------- 2-sample tests (working unchanged) ------------------
|
176 |
##-------- 2-sample tests (working unchanged) ------------------
|
| 97 |
|
177 |
|
| 98 |
R2 <- lapply(alts, ## for all three alternatives
|
178 |
R2 <- lapply(alts, ## for all three alternatives
|
| 99 |
function(alt)
|
179 |
function(alt)
|
| 100 |
lapply(seq_along(x0), function(k)
|
180 |
lapply(seq_along(x0), function(k)
|
| 101 |
wilcox.test(x = x.set$s0[[k]], y = x.set$s.[[k]],
|
181 |
wilcox.test(x = x.set$s0[[k]], y = x.set$s.[[k]],
|
| 102 |
exact=TRUE, conf.int=TRUE, alternative = alt)))
|
182 |
exact=TRUE, conf.int=TRUE, alternative = alt)))
|
| 103 |
length(w2 <- warnings()) # 27
|
183 |
length(w2 <- warnings()) # 1 (was '27')
|
| 104 |
unique(w2) # 3 different ones
|
184 |
if(!identical(w2, w1)) # typically they *are* identical ==> no new warnings
|
| - |
|
185 |
summary(w2)
|
| - |
|
186 |
|
| 105 |
|
187 |
|
| 106 |
table(uc2 <- unlist(c2 <- lapply(R2, function(A) lapply(A, class))))
|
188 |
table(uc2 <- unlist(c2 <- lapply(R2, function(A) lapply(A, class))))
|
| 107 |
stopifnot(uc2 == "htest")
|
189 |
stopifnot(uc2 == "htest")
|
| 108 |
|
190 |
|
| 109 |
stR2 <- lapply(stats,
|
191 |
stR2 <- lapply(stats,
|
| 110 |
function(COMP)
|
192 |
function(COMP)
|
| 111 |
lapply(R2, function(A) lapply(A, getC, C=COMP)))
|
193 |
lapply(R2, function(A) lapply(A, `[[`, COMP)))
|
| 112 |
|
194 |
|
| 113 |
lapply(stats[-3], ## -3: "conf.int" separately
|
195 |
lapply(stats[-3], ## -3: "conf.int" separately
|
| 114 |
function(ST) sapply(stR2[[ST]], unlist))
|
196 |
function(ST) sapply(stR2[[ST]], unlist))
|
| 115 |
|
197 |
|
| 116 |
noquote(sapply(stR2[["conf.int"]], function(.) vapply(., formatCI, "")))
|
198 |
noquote(sapply(stR2[["conf.int"]], function(.) vapply(., formatCI, "")))
|
The R Project SVN
R