using RCall, DataFrames
@rimport base as R
R.Sys_setenv(LANG = "en") # こうしておかないとR側でエラーが出たときにハングする
RObject{LglSxp} [1] TRUE
# 以下はJulia側でggplot2を便利に使うための設定
rcalljl_options(; kwargs...) = rcopy(RCall.rcall_p(:options, rcalljl_options=Dict(kwargs)))
function rplotsize(w, h)
if RCall.ijulia_mime == MIME("image/svg+xml")
rcalljl_options(; width=w/100, height=h/100)
else
rcalljl_options(; width=w, height=h)
end
end
function rplotpng(; kwargs...)
RCall.ijulia_setdevice(MIME("image/png"); kwargs...)
RCall.ijulia_mime, rcalljl_options()
end
function rplotsvg(; kwargs...)
RCall.ijulia_setdevice(MIME("image/svg+xml"); kwargs...)
RCall.ijulia_mime, rcalljl_options()
end
# Avoid svg display bug
using Base64
function RCall.ijulia_displayfile(m::MIME"image/svg+xml", f)
open(f) do io
svg = read(io, String)
base64 = base64encode(svg)
html = """<img src="data:$m;base64,$base64" />"""
display(MIME("text/html"), html)
end
end
@rlibrary ggplot2
rplotsvg()
(MIME type image/svg+xml, OrderedCollections.OrderedDict{Symbol, Any}(:rcalljl_options => OrderedCollections.OrderedDict{Symbol, Any}(:height => 5, :width => 6)))
rplotsize(640, 400) # グラフの表示サイズを設定
t = range(-10, 10, length=1000)
data = DataFrame(t = t, x = cos.(t), y = sin.(t))
ggplot(data=data, aes(x=:t, y=:x)) +
geom_line(color="red") +
geom_line(aes(y=:y), color="blue")
RObject{VecSxp}
@rlibrary stats # これでRのfisher_testなどをJuliaで使用可能
A = [
10 10
7 27
]
fisher_result_R = fisher_test(A) # P値が5%未満なのに, 95%信頼区間が1を含む
RObject{VecSxp} Fisher's Exact Test for Count Data data: structure(c(10L, 7L, 10L, 27L), .Dim = c(2L, 2L)) p-value = 0.03516 alternative hypothesis: true odds ratio is not equal to 1 95 percent confidence interval: 0.9836015 15.3827687 sample estimates: odds ratio 3.751532
p-value = 0.03516 が5%未満なのに、95%信頼区間 [0.9836015, 15.3827687] が1を含んでいる。
fisher_result = rcopy(fisher_result_R) # RのデータをJuliaのデータに変換9
OrderedCollections.OrderedDict{Symbol, Any} with 7 entries: :p_value => 0.0351564 :conf_int => [0.983602, 15.3828] :estimate => 3.75153 :null_value => 1.0 :alternative => "two.sided" :method => "Fisher's Exact Test for Count Data" :data_name => "structure(c(10L, 7L, 10L, 27L), .Dim = c(2L, 2L))"
R言語では以下のように結果が自動的に表示されて便利である.
Fisher's Exact Test for Count Data
data: structure(c(10L, 7L, 10L, 27L), .Dim = c(2L, 2L))
p-value = 0.03516
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.9836015 15.3827687
sample estimates:
odds ratio
3.751532
このような表示をJuliaで実現するには以下のようにすればよい.
# fisher_testの自前での実装 (エラーチェックの類に欠けた手抜きであることに注意!)
module My # 使い捨てモジュール
using Distributions, Roots
⪅(x ,y) = x < y || x ≈ y
or(a, b, c, d) = a*d/(b*c)
or(x::AbstractVecOrMat) = or(x...)
pval(d::FisherNoncentralHypergeometric, k) =
sum(pdf(d, j) for j in support(d) if pdf(d,j) ⪅ pdf(d, k))
pval(a, b, c, d, ω=1.0) =
pval(FisherNoncentralHypergeometric(a+b, c+d, a+c, ω), a)
pval(A::AbstractVecOrMat, ω=1.0) = pval(A..., ω)
ci(a, b, c, d, α=0.05) =
exp.(find_zeros(t -> pval(a, b, c, d, exp(t)) - α, -1e2, 1e2))
ci(A::AbstractVecOrMat, α=0.05) = ci(A..., α)
struct FisherTest{D, T}
data::D
ω::T
α::T
odds_ratio::T
p_value::T
conf_int::Vector{T}
end
Base.show(io::IO, x::FisherTest) = print(io, """
Fisher's Exact Test for Count Data
data: $(x.data)
p-value: $(x.p_value)
alternative hypothesis: true odds ratio is not equal to $(x.ω)
$(100(1 - x.α)) percent confidence interval:
$(x.conf_int)
odds ratio: $(x.odds_ratio)
""")
function fisher_test(data; ω=1.0, α = 0.05)
odds_ratio = or(data)
p_value = pval(data, ω)
conf_int = ci(data, α)
FisherTest(data, ω, α, odds_ratio, p_value, conf_int)
end
end
Main.My
My_fisher_result = My.fisher_test(A)
Fisher's Exact Test for Count Data data: [10 10; 7 27] p-value: 0.03515636840648691 alternative hypothesis: true odds ratio is not equal to 1.0 95.0 percent confidence interval: [1.0691205185478725, 13.492616894184634] odds ratio: 3.857142857142857
P値と95%信頼区間が整合的になっている!
dump(My_fisher_result)
Main.My.FisherTest{Matrix{Int64}, Float64} data: Array{Int64}((2, 2)) [10 10; 7 27] ω: Float64 1.0 α: Float64 0.05 odds_ratio: Float64 3.857142857142857 p_value: Float64 0.03515636840648691 conf_int: Array{Float64}((2,)) [1.0691205185478725, 13.492616894184634]