砂になった人

砂になった人の趣味や勉強についての備忘録

【Julia】ロジスティック回帰をやる

こんにちは。今回は、Juliaでロジスティック回帰を実装してみようという回です。今まではRを使ってきましたが、つい3日前くらいに唐突に「Julia始めてみようかな〜」と思い立って触り出しまして、その練習というか備忘録がてらブログにしてみました。実際にロジスティック回帰をJuliaでやろうとするならばGLMなるパッケージがあるそうですが、よくわからんかった(というよりはドキュメントを読むのが面倒だった)ので自力でやってみました。ただ、ロジットに関する実装系の日本語ブログはほとんどなかったので、後学のためにも記録として残しておきます。

準備

まずは、以下のパッケージを読み込みます。

# Load packages
using Optim # optimization
using Distributions # RNG
using NLSolversBase # calculation of Hessian
using Random 
using LinearAlgebra
# Seed setting
Random.seed!(0)

次に、いくつかの補助的な関数を定義します。頭がRに冒されているので、とりあえずRライクな関数を作りました。

rnorm(N, mu = 0, sd = 1) = rand(Normal(mu, sd), N) # Normal random variates
runif(N, min = 0, max = 1) = rand(Uniform(min, max), N) # Uniform random variates
plogis(x) = exp(x) / (1 + exp(x)) # Logit

それでは、今回使用するデータを擬似的に作成します。今回は二変数による回帰分析用のデータを作成していこうと思います。サイズは1000にして、二つの変数$x_{i1}$と$x_{i2}$はそれぞれ$N(0, 25)$と$U(0, 1)$から生成されるとし、誤差項$\varepsilon_{i}$は$N(0, 1)$から生成します。係数値$\beta$は長さ3の列ベクトルで、$[-3, -1, 2]$と定義します。また、結果変数の$Y_i$は確率$\text{logit}^{-1}(x_i^\top\beta + \varepsilon_i)$のベルヌーイ分布からサンプリングします。

# Define data
## Size
N = 1000
## Coefficient
beta = [-3 -1 2]
## Length of coef
K = length(beta)
## N * 3 matrix of covariates (1st col is intercept)
X = [ones(N) rnorm(N, 0, 5) runif(N)]
## Linear aggregator
y_star = X * beta' + rnorm(N)
## Calculate probability
p = plogis.(y_star)
## Generate outcome variable
Y = rand.(Bernoulli.(p))

対数尤度関数の定義

対数尤度関数loglikを定義します。ロジスティック回帰の場合、尤度関数$\mathcal{L}$はベルヌーイ分布として考えます。よって、対数尤度$l$は

\begin{align} \mathcal{L} &= \prod_{i=1}^N p^{Y_i}(1-p)^{(1 - Y_i)} \\ l &= \log \mathcal{L} = \sum\limits_{i=1}^N Y_i \log(p) + (1 - Y_i) \log(1 - p) \end{align}

となります。しかしながら、JuliaOptim.optimize関数は何だか最小化問題を解くのが通例となってるぽい(?)のでそれに合わせた形で対数尤度関数を定義しなければなりません。とはいえ、単純に符号を反転させればOKです。

# Log-likelihood
function loglik(x, y, b)
    len = length(y)
    prob = plogis.(x * b)
    llike = log.(prob) .* y + log.(ones(len) - prob) .* (ones(len) - y)
    return(-sum(llike)) # Minimization problem
end

最適化

いよいよ、最適化です。その前に、分散共分散行列の計算で用いるへシアンを計算するためにTwiceDifferentiableという関数を実行します。1項目は任意の関数、2項目は初期値、3項目のautodiff = :forwardは自動微分の設定です。その後、optimizeで最適化します。このとき、Nelder-Mead法を用いました。

# For Hessian
fn = TwiceDifferentiable(vars -> loglik(X, Y, vars[1:K]), ones(K); autodiff = :forward);

# Optimize
opt = optimize(fn, ones(K), method = NelderMead())

実行が終わったので、収束しているか確認すると

Optim.converged(opt)
true

無事、収束してました。

結果

早速、得られた計数値を取得し、またへシアンを取得して標準誤差を求めます。へシアンの符号を反転させて逆行列を取ると分散共分散行列が得られるので、その対角成分の平方根を取れば標準誤差が求まります。

# Extract the result
beta_hat = Optim.minimizer(opt)

# Calculate Hessian
H = hessian!(fn, beta_hat)
# Standard error
sd = sqrt.(H  |> inv |> diag)

次にパラメータの95%信頼区間を求めてみましょう。見にくいので、小数点第4位以下を切り捨てます。

# 95% CI
q_ci = quantile(Normal(0, 1), [0.975]) # 1.96
upr = beta_hat + sd .* q_ci # Upper
upr = round.(upr, digits = 3)
lwr = beta_hat - sd .* q_ci # Lower
lwr = round.(lwr, digits = 3)

以上をまとめた表を出力してみます(DataFrame形式ですが)。

using DataFrames
re(x, y) = repeat([x], y)
ci = map(string, re("[", 3), lwr, re(", ", 3), upr, re("]", 3))
var_name = ["Intercept", "x1", "x2"]
rename(DataFrame(hcat(var_name, beta_hat, ci, sd), :auto), 
        :x1 => "Variable", :x2 => "Est.", :x3 => "95% CI", :x4 => "Std.Error")
3×4 DataFrame
 Row │ Variable   Est.       95% CI            Std.Error 
     │ Any        Any        Any               Any       
─────┼───────────────────────────────────────────────────
   1 │ Intercept  -2.90563   [-3.475, -2.336]  0.29058
   2 │ x1         -0.839131  [-0.953, -0.726]  0.0578705
   3 │ x2         2.25445    [1.451, 3.058]    0.409948

なかなか良い推定値なのではないでしょうか。実際の値である$[-3, -1, 2]$に近く、しっかりと実装できてそうです。しかしながら、少々不安なのでRでも確かめることにしました。まず、Juliaで使用したデータを書き出しちゃいましょう。

using CSV
d = rename(DataFrame(hcat(Y, X), :auto),
            :x1 => :Y, :x2 => :Intercept, 
            :x3 => :x1, :x4 => :x2) 
d |> CSV.write("d.csv", delim = ',', writeheader = true)

それでは、実際にglmでロジットを回してみます。

d <- read.csv('d.csv')
summary(glm(Y ~ x1 + x2, data = d, family = 'binomial'))
Call:
glm(formula = Y ~ x1 + x2, family = "binomial", data = d)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.4804  -0.3232  -0.0565   0.2520   2.6856  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.90562    0.29058  -9.999  < 2e-16 ***
x1          -0.83913    0.05787 -14.500  < 2e-16 ***
x2           2.25441    0.40995   5.499 3.81e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1312.48  on 999  degrees of freedom
Residual deviance:  515.87  on 997  degrees of freedom
AIC: 521.87

Number of Fisher Scoring iterations: 7

グッドですね。ほとんど推定値がドンピシャでした。きちんと動かせましたね。

予測確率の計算

最後に、せっかくなので変数$x_{2i}$を用いた予測確率を計算していきたいと思います。ここでのアプローチは、King, Tomz & Wittenberg (2000, AJPS)*1が提案したシュミレーションを用いたprocedureであるCLARIFYを利用して、予測確率を計算します。手順は

  1. モデルの推定値・分散を用いて多変量正規分布からランダムに変数をサンプリング
  2. 興味がある説明変数以外は中央値(または平均値)に固定した擬似的な説明変数の行列を作成
  3. 1と2から予測確率の行列を計算し、点推定値と95%信頼区間を求める

です。そこまで難しくはないのでちゃっちゃっとやります。

# MVnormal rng
rmvnorm(N, mu, sigma) = rand(MvNormal(mu, sigma), N)
b_rng = rmvnorm(N, beta_hat, sd.^2)
# Generate hypothetical `x2`
x2_hyp = collect(0:0.01:1)
# Generate hypothetical covariate
X_hyp = [ones(length(x2_hyp)) re(median(X[:, 2]), length(x2_hyp))  x2_hyp]
# Calculate probability and CI
EY = plogis.(X_hyp * b_rng)
lower = zeros(size(EY)[1])
med = zeros(size(EY)[1])
upper = zeros(size(EY)[1])
for l in 1:size(EY)[1]
    lower[l] = quantile(EY[l, :], 0.025)
    med[l] = median(EY[l, :])
    upper[l] = quantile(EY[l, :], 0.975)
end
# Plot
using Plots
plot(x2_hyp, med,
        ribbon = (med .- lower, upper .- med),
        legend=:bottomright, fa = .3,
        label = false, xlabel = "x2", ylabel = "Predicted probabilities",
        dpi = 300)

こちらも良い感じですね。すっきりと実装から可視化までスムーズにできたと思います。

おわりに

今のところ、正直言ってJuliaのご利益がいまいち分かってませんw まだRでいいかなぁとか思ったり。

*1:King, G., Tomz, M., & Wittenberg, J. (2000). Making the most of statistical analyses: Improving interpretation and presentation. American journal of political science, 347-361.