#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。
R の mantelhaen.test()
ファイル名: mantelhaentest.jl 関数名: mantelhaentest
翻訳するときに書いたメモ
R では，exact = TRUE の場合の信頼区間の計算において，uniroot の tol がデフォルトのままの .Machine$double.eps^0.25 になっているので解の精度が不十分である。
Exact conditional test of independence in 2 x 2 x k tables
data: x
S = 16, p-value = 0.03994490358127
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
1.07740143018898 529.83739902335230
sample estimates:
common odds ratio
10.36102408137496
.Machine$double.eps^0.25
[1] 0.0001220703125
なので，小数点以下3桁程度の精度しかない。（信頼区間の推定値としては十分であろうが）
tol = Machine$double.eps にすれば，Julia での解と同じになる。
Exact conditional test of independence in 2 x 2 x k tables
data: x
S = 16, p-value = 0.03994490358127
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
1.077387763518033 531.512782889939672
sample estimates:
common odds ratio
10.36102408137498
なお，mantelhaen.test() 中の d2x2xk() は C によるプログラムを呼ぶことになっているが，以下のような計算をするだけのものである。
d2x2xk = function(K, m, n, t, d) {
y = pmax(0, t - n)
z = pmin(m, t)
tmy = t - y
mpn = m + n
zmy = z - y
tbl = matrix(0, K + 1, sum(zmy) + 1)
pos = 0
tbl[1, 1] = 1
for (i in 1:K) {
for (j in 1:(zmy[i]+1)) {
u = dhyper(tmy[i] - j+1, t[i], mpn[i] - t[i], n[i])
tbl[i + 1, j:(j + pos)] = tbl[i + 1, j:(j + pos)] + tbl[i, 1:(pos + 1)] * u
}
pos = pos + zmy[i]
}
return(tbl[K, ] / sum(tbl[K, ]))
}
==========#
using Match, Rmath, LinearAlgebra, Roots, Printf
function mantelhaentest(x, y = [], z = []; alternative = "twosided", correct = true, exact = false, conflevel = 0.95)
if typeof(x) == Array{Int64, 3}
any(size(x) .< 2) && error("each dimension in table must be >= 2")
else
length(y) == 0 && error("if 'x' is not an array, 'y' must be given")
length(z) == 0 && error("if 'x' is not an array, 'z' must be given")
any(diff([length(x), length(y), length(z)]) .!= 0) && error("'x', 'y', and 'z' must have the same length")
levelsx, tblx = table(x)
levelsy, tbly = table(y)
(length(levelsx) < 2 || length(levelsy) < 2) && error("'x' and 'y' must have at least 2 levels")
x = table(x, y, z)
end
any(sum(x, dims=[1, 2]) .< 2) && error("sample size in each stratum must be > 1")
I, J, K = size(x)
if I == 2 && J == 2
words = (["not equal to", "less than", "greater than"][indexin([alternative], ["twosided", "less", "greater"])])[1]
if !exact
method = "Mantel-Haenszel chi-squared test with continuity correction"
sx = reshape(sum(x, dims=2), I, K)
sy = reshape(sum(x, dims=1), J, K)
n = vec(sum(x, dims=[1, 2]))
delta = sum(x[1, 1, :] .- sx[1, :] .* sy[1, :] ./ n)
yates = correct && abs(delta) >= 0.5 ? 0.5 : 0
denom = sum(vec(prod(vcat(sx, sy), dims=1)) ./ (n .^ 2 .* (n .- 1)))
statistic = (abs(delta) - yates) ^ 2 / denom
df = 1
if alternative == "twosided"
pvalue = pchisq(statistic, df, false)
else
z = sign(delta) * sqrt(statistic)
pvalue = pnorm(z, alternative == "less")
end
sdiag = sum(x[1, 1, :] .* x[2, 2, :] ./ n)
soffd = sum(x[1, 2, :] .* x[2, 1, :] ./ n)
estimate = sdiag / soffd
sd = sqrt(sum((x[1, 1, :] .+ x[2, 2, :]) .* x[1, 1, :] .* x[2, 2, :] ./ n .^ 2) /
(2 * sdiag ^ 2) +
sum(((x[1, 1, :] .+ x[2, 2, :]) .* x[1, 2, :] .* x[2, 1, :] +
(x[1, 2, :] .+ x[2, 1, :]) .* x[1, 1, :] .* x[2, 2, :]) ./ n .^ 2) /
(2 * sdiag * soffd) +
sum((x[1, 2, :] .+ x[2, 1, :]) .* x[1, 2, :] .* x[2, 1, :] ./ n .^ 2) /
(2 * soffd ^ 2))
cint = @match alternative begin
"less" => [0, estimate * exp(qnorm(conflevel) * sd)]
"greater" => [estimate * exp(qnorm(conflevel, false) * sd), Inf]
"twosided" => estimate .* exp.([1, -1] .* qnorm((1 - conflevel) / 2) * sd)
end
println(method)
@printf("Mantel-Haenszel X-squared = %.5f, df = %d, p value = %.5f\n", statistic, df, pvalue)
@printf("alternative hypothesis: true common odds ratio is %s 1\n", words)
@printf("%g %% confidence interval: [%.6f, %.6f]\n", 100conflevel, cint[1], cint[2])
@printf("sample estimates: common odds ratio %g", estimate)
Dict(:statistic => statistic, :df => df, :pvalue => pvalue,
:conflevel => conflevel, :cint => cint, :method => method)
else
function dn2x2xk(ncp)
(length(ncp) == 1 && ncp == 1) && return dc
d = logdc .+ log.(ncp) .* support
d = exp.(d .- maximum(d))
d ./ sum(d)
end
function mn2x2xk(ncp)
ncp == 0 && return lo
ncp == Inf && return hi
sum(support .* dn2x2xk(ncp))
end
function pn2x2xk(q, ncp = 1; uppertail = false)
if length(ncp) == 1 && ncp == 0
uppertail ? float(q <= lo) : float(q >= lo)
elseif length(ncp) == 1 && ncp == Inf
uppertail ? float(q <= hi) : float(q >= hi)
else
d = dn2x2xk(ncp)
if uppertail
select = support .>= q
else
select = support .<= q
end
sum(d[select])
end
end
function mle(x)
x == lo && return 0
x == hi && return Inf
mu = mn2x2xk(1)
if mu > x
find_zero(t -> mn2x2xk(t) - x, (0, 1))
elseif mu < x
1 / find_zero(t -> mn2x2xk(1 / t) - x, (eps(), 1))
else
1
end
end
function ncpU(x, alpha)
x == hi && return Inf
p = pn2x2xk(x, 1)
if p < alpha
find_zero(t -> pn2x2xk(x, t) - alpha, (0, 1))
elseif p > alpha
1 / find_zero(t -> pn2x2xk(x, 1 / t) - alpha, (eps(), 1))
else
1
end
end
function ncpL(x, alpha)
x == lo && return 0
p = pn2x2xk(x, 1; uppertail = true)
if p > alpha
find_zero(t -> pn2x2xk(x, t, uppertail = true) - alpha, (0, 1))
elseif p < alpha
1 / find_zero(t -> pn2x2xk(x, 1 / t, uppertail = true) - alpha, (eps(), 1))
else
1
end
end
method = "Exact conditional test of independence in 2 x 2 x k tables"
mn = reshape(sum(x, dims=1), I, :)
m = mn[1, :]
n = mn[2, :]
t = reshape(sum(x, dims= 2), J, :)[1, :]
s = sum(x[1, 1, :])
lo = sum(pmax(0, t - n))
hi = sum(pmin(m, t))
support = lo:hi
dc = d2x2xk(K, m, n, t, hi - lo + 1)
logdc = log.(dc)
relErr = 1.0000001
alpha2 = (1 - conflevel) / 2
pvalue, cint = @match alternative begin
"less" => (pn2x2xk(s, 1), [0, ncpU(s, 1 - conflevel)])
"greater" => (pn2x2xk(s, 1, uppertail = true), [ncpL(s, 1 - conflevel), Inf])
"twosided" => (sum(dc[dc .<= dc[s - lo + 1] * relErr]), [ncpL(s, alpha2), ncpU(s, alpha2)])
end
statistic = s
estimate = mle(s)
println(method)
@printf("S = %g, p value = %.5f\n", statistic, pvalue)
@printf("alternative hypothesis: true common odds ratio is %s 1\n", words)
@printf("%g %% confidence interval: [%.6f, %.6f]\n", 100conflevel, cint[1], cint[2])
@printf("sample estimates: common odds ratio %g", estimate)
Dict(:statistic => statistic, :pvalue => pvalue,
:conflevel => conflevel, :cint => cint, :method => method)
end
else
df = (I - 1) * (J - 1)
n = zeros(df)
m = zeros(df)
V = zeros(df, df)
for k = 1:K
f = x[:, :, k]
ntot = sum(f)
rowsums = sum(f, dims=2)[1:I-1]
colsums = sum(f, dims=1)[1:J-1]
n .+= vec(f[1:I-1, 1:J-1])
m .+= vec(rowsums .* colsums' ./ ntot)
V .+= kron(
diag(ntot .* colsums, J - 1) .- (colsums .* colsums'),
diag(ntot .* rowsums, I - 1) .- (rowsums .* rowsums')) /
(ntot ^ 2 * (ntot - 1))
end
n = n - m
statistic = n' * (V \ n)
df = df
pvalue = pchisq(statistic, df, false)
method = "Cochran-Mantel-Haenszel test"
println(method)
@printf("Cochran-Mantel-Haenszel M^2 = %.5f, df = %d, p value = %.6f\n", statistic, df, pvalue)
Dict(:statistic => statistic, :df => df, :pvalue => pvalue, :method => method)
end
end
function table(x) # indices が少ないとき
indices = sort(unique(x))
counts = zeros(Int, length(indices))
for i in indexin(x, indices)
counts[i] += 1
end
return indices, counts
end
function table(x, y, z) # 三次元
indicesx = sort(unique(x))
indicesy = sort(unique(y))
indicesz = sort(unique(z))
counts = zeros(Int, length(indicesx), length(indicesy), length(indicesz))
for (i, j, k) in zip(indexin(x, indicesx), indexin(y, indicesy), indexin(z, indicesz))
counts[i, j, k] += 1
end
return indicesx, indicesy, indicesz, counts
end
function diag(x, n)
a = zeros(n, n)
[a[i, i] = x[i] for i = 1:n]
a
end
function adjustlength(x, y)
if length(x) == 1
x = fill(x, length(y))
elseif length(y) == 1
y = fill(y, length(x))
end
length(x) != length(y) && error("length(x) != length(y)")
x, y
end
function pmin(x, y)
x, y = adjustlength(x, y)
[min(x0, y0) for (x0, y0) in zip(x, y)]
end
function pmax(x, y)
x, y = adjustlength(x, y)
[max(x0, y0) for (x0, y0) in zip(x, y)]
end
function d2x2xk(K, m, n, t, d)
y = pmax(0, t - n)
z = pmin(m, t)
tmy = t - y
mpn = m + n
zmy = z - y
tbl = zeros(K + 1, sum(zmy) + 1)
pos = 0
tbl[1, 1] = 1
for i = 1:K
for j = 1:zmy[i]+1
u = dhyper(tmy[i] - j+1, t[i], mpn[i] - t[i], n[i])
tbl[i + 1, j:j + pos] .+= tbl[i, 1:pos + 1] * u
end
pos += zmy[i]
end
tbl[K, :] / sum(tbl[K, :])
end
x = reshape([0, 0, 6, 5, 3, 0, 3, 6, 6, 2, 0, 4, 5, 6, 1, 0, 2,
5, 0, 0], 2, 2, :);
mantelhaentest(x);
#=====
Mantel-Haenszel chi-squared test with continuity correction
Mantel-Haenszel X-squared = 3.92857, df = 1, p value = 0.04747
alternative hypothesis: true common odds ratio is not equal to 1
95 % confidence interval: [1.026713, 47.725133]
sample estimates: common odds ratio 7
=====#
using RCall
R"mantelhaen.test($x)"
#=====
Mantel-Haenszel chi-squared test with continuity correction
data: `#JL`$x
Mantel-Haenszel X-squared = 3.9286, df = 1, p-value = 0.04747
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
1.026713 47.725133
sample estimates:
common odds ratio
7
=====#
mantelhaentest(x, alternative="less");
#=====
Mantel-Haenszel chi-squared test with continuity correction
Mantel-Haenszel X-squared = 3.92857, df = 1, p value = 0.97626
alternative hypothesis: true common odds ratio is less than 1
95 % confidence interval: [0.000000, 35.052454]
sample estimates: common odds ratio 7
=====#
R"mantelhaen.test($x, alternative=\"less\")"
#=====
Mantel-Haenszel chi-squared test with continuity correction
data: `#JL`$x
Mantel-Haenszel X-squared = 3.9286, df = 1, p-value = 0.9763
alternative hypothesis: true common odds ratio is less than 1
95 percent confidence interval:
0.00000 35.05245
sample estimates:
common odds ratio
7
=====#
mantelhaentest(x, alternative="greater");
#=====
Mantel-Haenszel chi-squared test with continuity correction
Mantel-Haenszel X-squared = 3.92857, df = 1, p value = 0.02374
alternative hypothesis: true common odds ratio is greater than 1
95 % confidence interval: [1.397905, Inf]
sample estimates: common odds ratio 7
=====#
R"mantelhaen.test($x, alternative=\"greater\")"
#=====
Mantel-Haenszel chi-squared test with continuity correction
data: `#JL`$x
Mantel-Haenszel X-squared = 3.9286, df = 1, p-value = 0.02374
alternative hypothesis: true common odds ratio is greater than 1
95 percent confidence interval:
1.397905 Inf
sample estimates:
common odds ratio
7
=====#
mantelhaentest(x, exact=true);
#=====
Exact conditional test of independence in 2 x 2 x k tables
S = 16, p value = 0.03994
alternative hypothesis: true common odds ratio is not equal to 1
95 % confidence interval: [1.077388, 531.512783]
sample estimates: common odds ratio 10.361
=====#
R"mantelhaen.test($x, exact=TRUE)"
#=====
Exact conditional test of independence in 2 x 2 x k tables
data: `#JL`$x
S = 16, p-value = 0.03994
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
1.077401 529.837399
sample estimates:
common odds ratio
10.36102
=====#
mantelhaentest(x, exact=true, alternative="less");
#=====
Exact conditional test of independence in 2 x 2 x k tables
S = 16, p value =0.99862
alternative hypothesis: true common odds ratio is less than 1
95 % confidence interval: [0.000000, 261.487879]
sample estimates: common odds ratio 10.361
=====#
R"mantelhaen.test($x, exact=TRUE, alternative=\"less\")"
#=====
Exact conditional test of independence in 2 x 2 x k tables
data: `#JL`$x
S = 16, p-value = 0.9986
alternative hypothesis: true common odds ratio is less than 1
95 percent confidence interval:
0.0000 262.2733
sample estimates:
common odds ratio
10.36102
=====#
mantelhaentest(x, exact=true, alternative="greater");
#=====
Exact conditional test of independence in 2 x 2 x k tables
S = 16, p value = 0.01997
alternative hypothesis: true common odds ratio is greater than 1
95 % confidence interval: [1.384222, Inf]
sample estimates: common odds ratio 10.361
=====#
R"mantelhaen.test($x, exact=TRUE, alternative=\"greater\")"
#=====
Exact conditional test of independence in 2 x 2 x k tables
data: `#JL`$x
S = 16, p-value = 0.01997
alternative hypothesis: true common odds ratio is greater than 1
95 percent confidence interval:
1.384239 Inf
sample estimates:
common odds ratio
10.36102
=====#
#####################################################################################
x = reshape([8, 3, 0, 3, 0, 2, 0, 3, 3, 5, 0, 3,
2, 1, 2, 8, 2, 0, 2, 2, 0, 0, 1, 1,
0, 1, 0, 8, 2, 1, 8, 8, 0, 5, 5, 8,
0, 1, 0, 5, 2, 5, 1, 1, 2, 3, 3, 0,
5, 0, 2, 0, 0, 2, 1, 2, 2, 1, 0, 1,
5, 8, 1, 1, 0, 5, 0, 8, 1, 0, 5, 8], 3, 4, :)
mantelhaentest(x);
#=====
Cochran-Mantel-Haenszel test
Cochran-Mantel-Haenszel M^2 = 33.32493, df = 6, p value = 0.000009
=====#
R"mantelhaen.test($x)"
#=====
Cochran-Mantel-Haenszel test
data: `#JL`$x
Cochran-Mantel-Haenszel M^2 = 33.325, df = 6, p-value = 9.079e-06
=====#
