Julia に翻訳--145 重回帰分析

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

重回帰分析
http://aoki2.si.gunma-u.ac.jp/R/mreg.html

ファイル名: mreg.jl　　関数名: mreg

翻訳するときに書いたメモ

出力にこだわらなければ簡単なプログラムなんだけど。

==========#

using Statistics, LinearAlgebra, Rmath, Printf

function mreg(dat)
    n, nc = size(dat)
    nv = nc - 1
    r = cor(dat)
    mall = vec(mean(dat, dims=1))
    m = mall[1:nc-1]
    betas = r[1:nc-1, 1:nc-1] \ r[1:nc-1, nc]
    variance = cov(dat) .* (n - 1)
    prop = diag(variance)
    prop = (prop / prop[nc])[1:nc-1]
    b = betas ./ sqrt.(prop)
    Sr = transpose(variance[:, nc][1:nc-1]) * b
    St = variance[nc, nc]
    Se = St .- Sr
    SS = [Sr, Se, St]
    dfr = nv
    dfe = n - nv - 1
    dft = n - 1
    df = [dfr, dfe, dft]
    Ms = SS ./ df
    f = Ms[1] / Ms[2]
    fvalue = [f, NaN, NaN]
    p = [pf(f, dfr, dfe, false),NaN, NaN]
    b0 = mall[nc] - sum(b .* m)
    b = vcat(b, b0)
    inverse = inv((n - 1) .* cov(dat)[1:nc-1, 1:nc-1])
    SEb = vcat([sqrt(inverse[i, i] * Ms[2]) for i =1:nv],
               sqrt((1 / n + transpose(m) * inverse * m) * Ms[2]))
    tval = b ./ SEb
    pval = pt.(abs.(tval), n - nv - 1, false) * 2
    tolerance = 1 ./ diag(inv(cor(dat)[1:nc-1, 1:nc-1]))
    R2 = 1 - Se / St
    R = sqrt(R2)
    R2s = 1 - Ms[2] / Ms[3]
    loglik = - 0.5 * n * (log(2 * pi) + 1 - log(n) + log(Se))
    AIC = 2 * nc + 2 - 2 * loglik
    printresult(b, SEb, tval, pval, betas, tolerance)
    printanovatable(SS, df, Ms, f, p)
    printstatistics(R2, R, R2s, loglik, AIC)
    Dict(:result => hcat(b, SEb, tval, pval, vcat(betas, NaN), vcat(tolerance, NaN)),
            :anova => hcat(SS, df, Ms, fvalue, p),
            :Rs => [R2, R, R2s, loglik, AIC])
end

diag(a) = [a[i, i] for i in 1:size(a, 1)]

function printresult(b, SEb, tval, pval, betas, tolerance)
    @printf("%6s%12s%12s%12s%12s%18s%12s\n",
            "", "偏回帰係数", "標準誤差", "t 値", "P 値", "標準化偏回帰係数", "トレランス")
    for i = 1:length(betas)
        @printf("%6s%12.5g%12.5g%12.5g%12.5g%18.5g%12.5g\n",
                "Var" * string(i), b[i], SEb[i], tval[i], pval[i], betas[i], tolerance[i])
    end
    @printf("%6s%12.5g%12.5g%12.5g%12.5g\n",
            "定数項", b[end], SEb[end], tval[end], pval[end])
end

function printanovatable(SS, df, Ms, f, p)
    println("回帰の分散分析表")
    @printf("%4s  %10s  %6s  %12s  %12s  %12s\n",
               "", "平方和", "自由度", "平均平方", "F 値", "P 値")
    @printf("回帰  %10.6g  %6d  %12.6g  %12.6g  %12.6g\n",
               SS[1], df[1], Ms[1], f[1], p[1])
    @printf("残差  %10.6g  %6d  %12.6g\n", SS[2], df[2], Ms[2])
    @printf("全体  %10.6g  %6d  %12.6g\n", SS[3], df[3], Ms[3])
end

function printstatistics(R2, R, R2s, loglik, AIC)
    println("重相関係数 = $R")
    println("重相関係数の二乗 = $R2")
    println("自由度調整済重相関係数の二乗 = $R2s")
    println("対数尤度 = $loglik")
    println("AIC = $AIC")
end

dat = [
    1.2 1.9 0.9
    1.6 2.7 1.3
    3.5 3.7 2.0
    4.0 3.1 1.8
    5.6 3.5 2.2
    5.7 7.5 3.5
    6.7 1.2 1.9
    7.5 3.7 2.7
    8.5 0.6 2.1
    9.7 5.1 3.6]

mreg(dat)
# 偏回帰係数    標準誤差        t 値        P 値  標準化偏回帰係数  トレランス
# Var1     0.20462   0.0075643      27.051  2.4192e-08           0.67067     0.98358
# Var2     0.28663    0.010802      26.536  2.7638e-08           0.65793     0.98358
# 定数項     0.14918    0.054506      2.7368     0.02905
# 回帰の分散分析表
#   平方和  自由度      平均平方          F 値          P 値
# 回帰     6.67164       2       3.33582       823.477   4.93187e-09
# 残差   0.0283563       7     0.0040509
# 全体         6.7       9      0.744444
# 重相関係数 = 0.9978816139144453
# 重相関係数の二乗 = 0.9957677153884982
# 自由度調整済重相関係数の二乗 = 0.9945584912137834
# 対数尤度 = 15.13806917315012
# AIC = -22.27613834630024

Julia に翻訳--144 リッカート尺度

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

リッカート尺度
http://aoki2.si.gunma-u.ac.jp/R/likert.html

ファイル名: likert.jl　　関数名: likert

翻訳するときに書いたメモ

==========#

using Rmath

function likert(dat)
    resp = dat / sum(dat)
    cum = cumsum(resp)
    result =(dnorm.(qnorm.(vcat(0, cum[1:end-1]))) .- dnorm.(qnorm.(cum))) ./ resp
    println(result)
    Dict(:result => result)
end

dat = [7, 18, 34, 26, 15]
likert(dat)
# [-1.9181130652891667, -1.0194925450770291, -0.20873614063905502, 0.5984157124911674, 1.5543918350245474]

Julia に翻訳--143 クロンバックのα信頼性係数

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

クロンバックのα信頼性係数
http://aoki2.si.gunma-u.ac.jp/R/alpha.html

ファイル名: alpha.jl　　関数名: alpha

翻訳するときに書いたメモ

出力にはこだわらない。
必要なら，戻り値を使って。

==========#

using Statistics

function alpha(x; detail = false)
    nr, k = size(x)
    (k > 1 && nr > 1) || error("parameter error")
    α = alpha0(x)
    println("α = $α")
    if detail == true && k >= 3
        R2 = 1 .- 1 ./ diag(inv(cor(x)))
        α2 = similar(R2)
        cor2 = similar(R2)
        for i in 1:k
            x2 = x[:, [i != j for j in 1:k]]
            α2[i] = alpha0(x2)
            cor2[i] = cor(x[:, i], vec(sum(x2, dims=2)))
        end
        println("それぞれの変数を除いたときの α\n$α2")
        println("それぞれの変数とその変数を除いたときの合計値との相関係数\n$cor2")
        println("それぞれの変数をその変数以外の変数で予測したときの決定係数\n$R2")
        Dict(:alpha => α, :alpha2 => α2, :R2 => R2, :cor2 => cor2)
    end
end

function alpha0(x)
    k = size(x, 2)
    VarCovMat = cov(x)
    Sy2 = sum(VarCovMat)
    Sj2 = sum(diag(VarCovMat))
    k / (k - 1) * (1 - Sj2 / Sy2)
end

diag(a) = [a[i, i] for i in 1:size(a, 1)]

x = [
    49 44 37 54
    36 36 36 42
    42 51 45 35
    47 67 54 40
    54 59 68 54
    51 55 59 67
    45 36 48 46
    72 49 50 58];

alpha(x)
alpha(x, detail = true)
# α = 0.7403365794138637
# それぞれの変数を除いたときの α
# [0.6657399977072106, 0.7307558528428093, 0.592778807780536, 0.72184]
# それぞれの変数とその変数を除いたときの合計値との相関係数
# [0.5599374682949105, 0.4444039960775981, 0.6801937600489253, 0.4589177292647371]
# それぞれの変数をその変数以外の変数で予測したときの決定係数
# [0.4298971517659005, 0.551857280000196, 0.6071655584759322, 0.4821633678144336]

Julia に翻訳--142 スミルノフ・グラブス検定

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

スミルノフ・グラブス検定
http://aoki2.si.gunma-u.ac.jp/R/SG.html

ファイル名: sg.jl　　関数名: sg

翻訳するときに書いたメモ

==========#

using Statistics, Rmath

function sg(x)
    n = length(x)
    t = abs.(extrema(x) .- mean(x)) ./ std(x)
    df = n - 2
    p = @. n * pt(sqrt(df / ((n - 1) ^ 2 / t ^ 2 / n - 1)), df, false)
    p = [min(p0, 1) for p0 in p]
    println("min: $(minimum(x))")
    println("t = $(t[1]),  df = $df,  pvalue = $(p[1])")
    println("max: $(maximum(x))")
    println("t = $(t[2]),  df = $df,  pvalue = $(p[2])")
    Dict(:t => t, :df => df, :p => p)
end

x = [133, 134, 134, 134, 135, 135, 139, 140, 140, 140,
    141, 142, 142, 144, 144, 147, 147, 149, 150, 164]

sg(x)
# min: 133
# t = 1.1724347855313042,  df = 18,  pvalue = 1.0
# max: 164
# t = 3.0052064042928888,  df = 18,  pvalue = 0.004864034251605727

 

Julia に翻訳--141 Breslow-Day 検定

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

Breslow-Day 検定
http://aoki2.si.gunma-u.ac.jp/R/BD-test.html

ファイル名: bdtest.jl　　関数名: bdtest

翻訳するときに書いたメモ

@. なんか書かなくてもよいようにしてほしいなあ

==========#

using Rmath

function bdtest(m)
    d1, d2, k = size(m)
    Nk = vec(sum(m, dims=[1,2])) # vex() がないと，1x1x8 の3次元行列になる
    psiMH = sum(m[1, 1, :] .* m[2, 2, :] ./ Nk) / sum(m[2, 1, :] .* m[1, 2, :] ./ Nk)
    nk1 = m[1, 1, :] .+ m[1, 2, :]
    nk2 = m[2, 1, :] .+ m[2, 2, :]
    xkt = m[1, 1, :] .+ m[2, 1, :]
    a1 = psiMH - 1
    @. b1 = -psiMH * (nk1 + xkt) - nk2 + xkt
    @. c1 = psiMH * nk1 * xkt
    @. e = (-b1 - sqrt(b1 ^ 2 - 4 * a1 * c1)) / (2 * a1)
    @. v = 1 / (1 / e + 1 / (nk1 - e) + 1 / (xkt - e) + 1 / (nk2 - xkt + e))
    chisqBD = sum((m[1, 1, :] .- e) .^ 2 ./ v) # この行は @. を使うとエラーになる
    df = k - 1
    p = pchisq(chisqBD, df, false)
    println("chisq. = $chisqBD,  df = $df,  p value = $p")
    Dict(:chisq => chisqBD, :df => df, :pvalue => p)
end

m = [1, 8, 0, 10,
     2, 47, 0, 60,
     3, 29, 1, 13,
     3, 17, 0, 7,
     13, 45, 1, 45,
     13, 26, 0, 18,
     8, 15, 0, 10,
     11, 4, 0, 4 ];
m = reshape(m, 2, 2, :)

bdtest(m)
# chisq. = 8.943202421541203,  df = 7,  p value = 0.25675965693687036

 

Julia の小ネタ--014 三次元配列の作り方

以下のような 2×3×4 行列を作る方法は 2 通りある。

2×3×4 Array{Int64,3}:
[:, :, 1] =
 1  3  5
 2  4  6

[:, :, 2] =
 7   9  11
 8  10  12

[:, :, 3] =
 13  15  17
 14  16  18

[:, :, 4] =
 19  21  23
 20  22  24

第1の方法

Julia は R や FORTRAN と同じように，列優先（次元表示の左にある添え字が最も速く変化する）なので，列優先で 1 次元配列（ベクトル）を用意して，reshape() で次元を定義し直す。

reshape([1,2,3,4,5,6,7,8,9,10,
         11,12,13,14,15,16,17,18,19,20,
         21,22,23,24], (2, 3, 4))

この場合は，以下のように書くこともできる。

reshape(1:24, 2, 3, :)

第2の方法

行優先で第1，第2次元の配列を定義し，cat で連結するときに dims=3 を指定する。

cat([1 3 5; 2 4 6], [7 9 11; 8 10 12],
    [13 15 17; 14 16 18], [19 21 23; 20 22 24], dims=3)

 

Julia に翻訳--140 Chow 検定

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

Chow 検定
http://aoki2.si.gunma-u.ac.jp/R/Chow.html

ファイル名: chow.jl　　関数名: chow

翻訳するときに書いたメモ

==========#

using Statistics, Rmath

chow = function(dat1, dat2)
    ess12 = ess(dat1) + ess(dat2)
    essc = ess(vcat(dat1, dat2))
    nr, df1 = size(dat1)
    df2 = nr + size(dat2, 1) - 2 * df1
    f = (essc - ess12) * df2 / (df1 * ess12)
    p = pf(f, df1, df2, false)
    println("F = $f,  df1 = $df1,  df2 = $df2,  p value = $p")
    Dict(:F => f, :df1 => df1, :df2 => df2, :pvalue => p)
end

function ess(dat) # 回帰分析を行い，誤差変動を返す関数
    r = cor(dat)
    betas = r[1:end-1, 1:end-1] \ r[1:end-1, end]
    SS = cov(dat, corrected=false) .* size(dat, 1)
    prop = [SS[i, i] for i in 1:size(SS, 1)]
    Sr = sum(SS[1:end-1, end] .* (betas .* sqrt.(prop[end] ./ prop[1:end-1])))
    St = SS[end, end]
    St - Sr
end

dat1 = [
    1.2 1.9 0.9
    1.6 2.7 1.3
    3.5 3.7 2.0
    4.0 3.1 1.8
    5.6 3.5 2.2
    5.7 7.5 3.5
    6.7 1.2 1.9
    7.5 3.7 2.7
    8.5 0.6 2.1
    9.7 5.1 3.6
    ];

dat2 = [
    1.4 1.3 0.5
    1.5 2.3 1.3
    3.1 3.2 2.5
    4.4 3.6 1.1
    5.1 3.1 2.8
    5.2 7.3 3.3
    6.5 1.5 1.3
    7.8 3.2 2.2
    8.1 0.1 2.8
    9.5 5.6 3.9
    ];

chow(dat1, dat2)
# F = 0.07154752544573709,  df1 = 3,  df2 = 14,  p value = 0.9742318542171906

Julia に翻訳--139 分散・共分散行列の同等性の検定

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

分散・共分散行列の同等性の検定
http://aoki2.si.gunma-u.ac.jp/R/eq-cov.html

ファイル名: boxm.jl　　関数名: boxm

翻訳するときに書いたメモ

関数を入れ子にして一行で書くのもよいが，for ループで書くとわかりやすいというメリットもある。

==========#

using Statistics, LinearAlgebra, Rmath

function boxm(x, gvar)
    nv = size(x, 2)
    nv >= 2 || error("分散共分散行列を計算する変数は2個以上必要です")
    index, ni = table(gvar)
    n = length(gvar)
    g = length(ni)
    sumlogdetSi = 0
    S = zeros(nv, nv)
    for (i, gvari) in enumerate(index)
        y = x[gvar .== gvari, :]
        covy = cov(y)
        sumlogdetSi += (ni[i] - 1) * log(det(covy))
        S .+= (ni[i] - 1) .* covy
    end
    S ./= n - g
    M = (n - g) * log(det(S)) - sumlogdetSi
    df1 = (g - 1) * nv * (nv + 1) / 2
    rho = 1 - (2 * nv ^ 2 + 3 * nv - 1) / (6 * (nv + 1) * (g - 1)) * (sum(1 ./ (ni .- 1)) - 1 / (n - g))
    tau = (nv - 1) * (nv + 2) / (6 * (g - 1)) * (sum(1 ./ (ni .- 1) .^ 2) - 1 / (n - g) ^ 2)
    df2 = (df1 + 2) / abs(tau - (1 - rho) ^ 2)
    gamma = (rho - df1 / df2) / df1
    F = M * gamma
    p = pf(F, df1, df2, false)
    println("M = $M\nF = $F,  df1 = $df1,  df2 = $df2\np value = $p")
    Dict(:M => M, :F => F, :df1 => df1, :df2 => df2, :pvalue => p)
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 rep(x, n::Array{Int64,1})
    length(x) == length(n) || error("length(x) wasn't length(n)")
    vcat([repeat([x[i]], n[i]) for i = 1:length(x)]...)
end

x = [
   2.9 161.7 120.8
   2.3 114.8  85.2
   2.0 128.4  92.0
   3.2 149.2  97.3
   2.7 126.0  81.1
   4.4 133.8 107.6
   4.1 161.3 114.0
   2.1 111.5  77.3
   4.8 198.7 172.9
   3.6 199.3 157.9
   2.0 188.4 152.7
   4.9 183.6 164.2
   3.9 173.5 172.2
   4.4 184.9 163.2 ];
gvar = rep(1:2, [8, 6]);

boxm(x, gvar)
# M = 9.730421574785467
# F = 1.1535477771673381,  df1 = 6.0,  df2 = 795.1946969489453
# p value = 0.32937843690375435