### Tritiated Water Diffusion Across Human Chorioamnion
### Hollander & Wolfe (1999), Table 4.1, page 110
water_transfer <- data.frame(
pd = c(0.80, 0.83, 1.89, 1.04, 1.45, 1.38, 1.91, 1.64, 0.73, 1.46,
1.15, 0.88, 0.90, 0.74, 1.21),
age = factor(c(rep("At term", 10), rep("12-26 Weeks", 5))))
### Wilcoxon-Mann-Whitney test, cf. Hollander & Wolfe (1999), page 111
### exact p-value and confidence interval for the difference in location
### (At term - 12-26 Weeks)
wt <- wilcox_test(pd ~ age, data = water_transfer,
distribution = "exact", conf.int = TRUE)
print(wt)
### extract observed Wilcoxon statistic, i.e, the sum of the
### ranks for age = "12-26 Weeks"
statistic(wt, "linear")
### its expectation
expectation(wt)
### and variance
covariance(wt)
### and, finally, the exact two-sided p-value
pvalue(wt)
### Confidence interval for difference (12-26 Weeks - At term)
wilcox_test(pd ~ age, data = water_transfer,
xtrafo = function(data)
trafo(data, factor_trafo = function(x) as.numeric(x == levels(x)[2])),
distribution = "exact", conf.int = TRUE)
### Permutation test, asymptotic p-value
oneway_test(pd ~ age, data = water_transfer)
### approximate p-value (with 99\% confidence interval)
pvalue(oneway_test(pd ~ age, data = water_transfer,
distribution = approximate(B = 9999)))
### exact p-value
pt <- oneway_test(pd ~ age, data = water_transfer, distribution = "exact")
pvalue(pt)
### plot density and distribution of the standardised
### test statistic
layout(matrix(1:2, nrow = 2))
s <- support(pt)
d <- sapply(s, function(x) dperm(pt, x))
p <- sapply(s, function(x) pperm(pt, x))
plot(s, d, type = "S", xlab = "Teststatistic", ylab = "Density")
plot(s, p, type = "S", xlab = "Teststatistic", ylab = "Cumm. Probability")
### Length of YOY Gizzard Shad from Kokosing Lake, Ohio,
### sampled in Summer 1984, Hollander & Wolfe (1999), Table 6.3, page 200
YOY <- data.frame(length = c(46, 28, 46, 37, 32, 41, 42, 45, 38, 44,
42, 60, 32, 42, 45, 58, 27, 51, 42, 52,
38, 33, 26, 25, 28, 28, 26, 27, 27, 27,
31, 30, 27, 29, 30, 25, 25, 24, 27, 30),
site = factor(c(rep("I", 10), rep("II", 10),
rep("III", 10), rep("IV", 10))))
### Kruskal-Wallis test, approximate exact p-value
kw <- kruskal_test(length ~ site, data = YOY,
distribution = approximate(B = 9999))
kw
pvalue(kw)
### Nemenyi-Damico-Wolfe-Dunn test (joint ranking)
### Hollander & Wolfe (1999), page 244
### (where Steel-Dwass results are given)
if (require("multcomp")) {
NDWD <- oneway_test(length ~ site, data = YOY,
ytrafo = function(data) trafo(data, numeric_trafo = rank),
xtrafo = function(data) trafo(data, factor_trafo = function(x)
model.matrix(~x - 1) %*% t(contrMat(table(x), "Tukey"))),
teststat = "maxtype", distribution = approximate(B = 90000))
### global p-value
print(pvalue(NDWD))
### sites (I = II) != (III = IV) at alpha = 0.01 (page 244)
print(pvalue(NDWD, method = "single-step"))
}
Run the code above in your browser using DataLab