トップ > JuliaLang > 線形混合モデルを（Juliaで）理解したい

線形混合モデルを（Juliaで）理解したい

JuliaLang 研究 R

以前「線形混合モデルを理解したい」という記事を読んで、Rで線形混合モデルを試したことがあるのですが、今回は、Juliaで同じことをやってみたときにハマったところを書き残しておきます。

やりたかったことは、前述の記事からもってきた以下のようなコードです。palmerpenguinsというパッケージに含まれているpenguinsというデータから、NAが含まれるデータを削除して、lme4でモデルを作る、というものです。もちろんRCallを使えば以下のように何も考えずに実行でき、「Julia最高！」となります。

using RCall                            # JuliaからRを呼び出す
R"library(palmerpenguins)"
R"summary(penguins)"
R"library(tidyverse)"                  # Rでパイプ演算子（%>%）を使うために必要
R"usedata <- penguins %>% na.omit()"   # NAデータを削除
R"glmm_model <- lme4::lmer(bill_length_mm ~ flipper_length_mm
    + (1 + flipper_length_mm|species),
    data = usedata)"                   # モデルの作成と分析
R"summary(glmm_model)"
R"lme4::ranef(glmm_model)"             # 回帰係数を表示する

これを、Rを最低限しか使わずにJuliaで実行します。まず、palmerpenguinsパッケージのデータを読み込みます。

using RCall, CSV, DataFrames

R"library(palmerpenguins)"
filepath = R"""palmerpenguins::path_to_file("penguins.csv")"""
penguins = CSV.read(convert(String, filepath), DataFrame; missingstring="NA")
penguins = dropmissing(penguins)

RDatasetsパッケージにはpalmerpenguinsは含まれていません。RDatasetsを使えば同じマシンにインストールされているRのパッケージを参照できるものだと思い込んでいたので、最初にここでハマりました。というわけで、R側にインストールされたパッケージからCSVファイルを読み込む必要があります。Rにインストールされているパッケージを参照して、データファイルへのパスを取ってきてCSV.read()します。path_to_file()の中でダブルクォートを使っているので、ダブルクォート三つで式全体を囲んでいます。式の中のダブルクォートをバックスラッシュでエスケープしてもいいでしょう。

さて、ハマりポイント二つ目が、NAの扱いでした。Juliaでは欠損値をmissingとして扱います。なので、ファイルから読み込むときに欠損値がどのような表記になっているかをmissingstringで指定してあげます。複数の表し方があるときにはmissingstring=["NA", "-999"]といった書き方ができます。

線形混合モデルの分析にはMixedAnova.jlのlme()が使えそうです。このパッケージはGLM.jl、MixedModels.jl、FixedEffectModels.jlなどを統合してANOVAをくっつけたような便利パッケージです。以下のように使いましたが、ANOVAを使わないのであればMixedModels.jlだけでも同じことができるかもしれません。

using MixedAnova
mm_init()

lmm1 = lme(@formula(bill_length_mm ~ flipper_length_mm
    + (1 + flipper_length_mm|species)), penguins)

出力されたのは次のようなモデルです。Rのlme4を使ったものと値がやや違いますが、許容範囲なのかな。Rで実行したときはパラメータが収束していないというwarningが出るので、そのせいかもしれません。とはいえ「Juliaは誤りが多くて使えない（Why I no longer recommend Julia）」なんて意見もあるので気を付けないといけませんね。

Linear mixed model fit by maximum likelihood
 bill_length_mm ~ 1 + flipper_length_mm + (1 + flipper_length_mm | species)
   logLik   -2 logLik     AIC       AICc        BIC    
  -795.0174  1590.0348  1602.0348  1602.2924  1624.8836

Variance components:
               Column        Variance   Std.Dev.    Corr.
species  (Intercept)        143.4020628 11.9750600
         flipper_length_mm    0.0038364  0.0619387 -0.95
Residual                      6.5214560  2.5537141
 Number of obs: 333; levels of grouping factors: 3

  Fixed-effects parameters:
─────────────────────────────────────────────────────────
                       Coef.  Std. Error      z  Pr(>|z|)
─────────────────────────────────────────────────────────
(Intercept)        -0.059232   8.14618    -0.01    0.9942
flipper_length_mm   0.221575   0.0416848   5.32    <1e-06
─────────────────────────────────────────────────────────

係数を得るにはlme4と同様にranef()が使えますが、raneftables()のほうが情報が多くて理解しやすいです。

julia> raneftables(lmm1)
(species = (species = String15["Adelie", "Chinstrap", "Gentoo"],
(Intercept) = [9.780118308163642, 4.4103374755445595, -14.1904557836386],
flipper_length_mm = [-0.06844124963200851, 0.0054061201200947765, 0.06303512951155885]),)

以上、環境はJulia 1.7.3 / Windows 10でした。

julia> versioninfo()
Julia Version 1.7.3
Commit 742b9abb4d (2022-05-06 12:58 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i5-10400 CPU @ 2.90GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, skylake)
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS =