http://districtdatalabs.silvrback.com/intro-to-r-for-microsoft-excel-users
であるが。ggplot2 と dplyr をお勧めなのはいいとして，
diamonds <- data.frame(ggplot2::diamonds)
diamonds$size[diamonds$carat < 0.5] <- "Small"
diamonds$size[diamonds$carat >=0.5 & diamonds$carat < 1] <- "Medium"
diamonds$size[diamonds$carat >= 1] <- "Large"
というのを例示しているが，すじが悪い。
diamonds$size2 <- cut(diamonds$carat, c(0, 0.5, 1, 99999), right=FALSE, labels=c("Small", "Medium", "Large"))
のほうを勧めるほうがよいと思われる。
ifelse もお勧めではない。
Brunner-Munzel 検定 ッテ，もはや，参照する価値もない記事です。パチもん検定です。
後半に追記あり（2021/10/13)
http://oku.edu.mie-u.ac.jp/~okumura/stat/brunner-munzel.html
で紹介されているが，結構時間がかかる。
> system.time(print(okumura(x, y)))
[1] 0.008037645
ユーザ システム 経過
1898.105 10.102 1923.545
brunner.munzel.test の必要な部分だけ取り出して，2倍強の速度になった。
> system.time(print(brunner.munzel.permutation.test(x, y)$p.value))
[1] 0.008037645
ユーザ システム 経過
788.666 5.037 794.799
brunner.munzel.permutation.test = function(x, y) {
BM = function(X) {
x = xandy[X]
y = xandy[-X]
r1 = rank(x)
r2 = rank(y)
r = rank(c(x, y))
m1 = mean(r[1:n1])
m2 = mean(r[n1 + 1:n2])
v1 = sum((r[1:n1] - r1 - m1 + (n1 + 1)/2)^2)/(n1 - 1)
v2 = sum((r[n1 + 1:n2] - r2 - m2 + (n2 + 1)/2)^2)/(n2 - 1)
(m2 - m1)/sqrt(n1 * v1 + n2 * v2)
}
x = na.omit(x)
y = na.omit(y)
n1 = length(x)
n2 = length(y)
xandy = c(x, y)
structure(list(method = "Brunner-Munzel Permutation Test",
p.value = mean(abs(combn(n1+n2, n1, BM)) >= abs(BM(1:n1))),
data.name = paste(deparse(substitute(x)), "and", deparse(substitute(y)))),
class = "htest")
}
okumura = function(x, y) {
bm = brunner.munzel.test(x, y)$statistic
n1 = length(x)
n2 = length(y)
N = n1 + n2
xandy = c(x, y)
foo = function(X) {
brunner.munzel.test(xandy[X], xandy[-X])$statistic
}
z = combn(1:N, n1, foo)
mean(abs(z) >= abs(bm))
}
########################################################
Julia で書いたら，マシンも違うが，約10秒（少なくとも35倍ほど速い）。
using Combinatorics # combinations
using StatsBase # tiedrank
using Statistics # mean
function brunner_munzel_permutation_test(x, y)
function BM(n1, n2, n, x, y)
r1 = tiedrank(x)
r2 = tiedrank(y)
r = tiedrank(vcat(x, y))
m1 = mean(r[1:n1])
m2 = mean(r[n1 + 1:n])
v1 = sum((r[1:n1] .- r1 .- m1 .+ (n1 + 1)/2).^2)/(n1 - 1)
v2 = sum((r[n1 + 1:n] .- r2 .- m2 .+ (n2 + 1)/2).^2)/(n2 - 1)
abs(m2 - m1)/sqrt(n1 * v1 + n2 * v2)
end
all_perm(x, n) = Iterators.product((x for i = 1:n)...)
x = [z for z in x if !ismissing(z)]; # 欠損値 missing を除くデータ
y = [z for z in y if !ismissing(z)]; # 欠損値 missing を除くデータ
n1 = length(x)
n2 = length(y)
xandy = vcat(x, y);
n = n1 + n2
s0 = BM(n1, n2, n, x, y)
results = Float64[]
for i in combinations(1:n, n1)
flag = falses(n)
flag[i] .= true
x = xandy[flag]
y = xandy[.!flag]
append!(results, BM(n1, n2, n, x, y))
end
mean(results .>= s0 - eps())
end
