以前「線形混合モデルを理解したい」という記事を読んで、Rで線形混合モデルを試したことがあるのですが、今回は、Juliaで同じことをやってみたときにハマったところを書き残しておきます。
やりたかったことは、前述の記事からもってきた以下のようなコードです。palmerpenguinsというパッケージに含まれているpenguinsというデータから、NAが含まれるデータを削除して、lme4でモデルを作る、というものです。もちろんRCallを使えば以下のように何も考えずに実行でき、「Julia最高！」となります。
using RCall
R"library(palmerpenguins)"
R"summary(penguins)"
R"library(tidyverse)"
R"usedata <- penguins %>% na.omit()"
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 =
