【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$は
となります。しかしながら、Julia
のOptim.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から予測確率の行列を計算し、点推定値と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
でいいかなぁとか思ったり。