将混合效应模型公式从R（lme4）改写为Julia

Question

我希望能够在Julia中获得与R中lme4库的lmer函数相同的结果。 请在下面找到一个使用R内置的mtcars数据集的示例

library(lme4)
data<-mtcars    
data$vs<-as.factor(data$vs)
data$am<-as.factor(data$am)
data$gear<-as.factor(data$gear)
str(data)
model <- lmer(mpg ~ cyl:gear + hp:am + (1|gear:am), data = data)

我在Julia的MixedModels包中找到了lmm()函数，它应该能够运行相同的结果，但我不知道如何使用lmm()重写lmer()函数的第一个参数中的公式，特别是交互（:）运算符。

如果您能提供一个简短的示例，我将不胜感激。

- JoannaW

---

这可能会有所帮助 https://github.com/dmbates/MixedModels.jl - MLavoie

交互操作符只是 *（这会给您主效应和交互作用）或 &（如果您只想要交互作用）。 - mmagnuski

Answer 1

这里R和Julia模型之间的对应关系似乎并不完全相同，还存在不同数值算法的问题。但是我尝试使用MixedModels重新创建相同的模型：

using RDatasets
using MixedModels

mtcars = dataset("datasets","mtcars")
mtcars[:AM] = PooledDataArray(mtcars[:AM])
mtcars[:Gear] = PooledDataArray(mtcars[:Gear])
mtcars[:GearAM] = PooledDataArray(collect(zip(mtcars[:Gear],mtcars[:AM])))
m = fit!(lmm(MPG ~ 1 + Gear + AM + Cyl + HP + (1|GearAM),mtcars))

手动创建混合效应列很麻烦，也许有更好的方法。请注意R和Julia之间系数命名的差异。两者都有6个固定效应系数。

在我的电脑上，Julia的解决方案似乎不同，但实现了更好的对数似然值。随机效应预计会比较弱，因为其变量已经存在于固定效应中（它仅解释Gear和AM之间的依赖关系），并且只有32个数据点。

希望这能帮到您，如果您有更好的理解，可以在另一个答案或评论中添加。

- Dan Getz

---

我需要将它转换为PooledDataArray，还是只需转换为分类变量？ - skan

Answer 2

Tl;DR：如果使用REML=true和@formula(MPG ~ 1 + Cyl & Gear + HP & AM + (1 | Gear & AM)，则可以完全复现。

更长的答案：支持(1 | Gear & AM)需要MixedModels.jl的3.0版本或更高版本，该版本截至目前（2021年9月1日）未发布，因此您需要使用pkg"add MixedModels#master"进行安装。这是在Julia 1.5.0上运行的。

using RDatasets, MixedModels
mtcars = RDatasets.dataset("datasets","mtcars")
categorical!(mtcars, [:AM, :Gear])
m = fit(MixedModel, 
        @formula(MPG ~ Cyl & Gear + HP & AM + (1 | Gear & AM)), 
        mtcars, 
        REML=true)

产生以下输出：

Linear mixed model fit by REML
 MPG ~ 1 + Cyl & Gear + HP & AM + (1 | Gear & AM)
 REML criterion at convergence: 168.11435769585046

Variance components:
             Column    Variance  Std.Dev. 
Gear & AM (Intercept)  16.913543 4.1126078
Residual                7.988539 2.8264004
 Number of obs: 32; levels of grouping factors: 4

  Fixed-effects parameters:
──────────────────────────────────────────────────────
                    Coef.  Std. Error      z  Pr(>|z|)
──────────────────────────────────────────────────────
(Intercept)    34.3507      3.50385     9.80    <1e-21
Cyl & Gear: 3  -1.08837     0.781649   -1.39    0.1638
Cyl & Gear: 4  -1.56833     0.875625   -1.79    0.0733
Cyl & Gear: 5  -0.167413    1.42839    -0.12    0.9067
HP & AM: 0     -0.0412826   0.0215381  -1.92    0.0553
HP & AM: 1     -0.0613247   0.0298067  -2.06    0.0396
──────────────────────────────────────────────────────

与您发布的 R 代码相比：

Linear mixed model fit by REML ['lmerMod']
Formula: mpg ~ cyl:gear + hp:am + (1 | gear:am)
   Data: data

REML criterion at convergence: 168.1

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.42280 -0.54922 -0.07898  0.54892  2.09217 

Random effects:
 Groups   Name        Variance Std.Dev.
 gear:am  (Intercept) 16.914   4.113   
 Residual              7.989   2.826   
Number of obs: 32, groups:  gear:am, 4

Fixed effects:
            Estimate Std. Error t value
(Intercept) 34.35066    3.50385   9.804
cyl:gear3   -1.08837    0.78165  -1.392
cyl:gear4   -1.56833    0.87562  -1.791
cyl:gear5   -0.16741    1.42839  -0.117
hp:am0      -0.04128    0.02154  -1.917
hp:am1      -0.06132    0.02981  -2.057

不修改数据

您还可以通过为变量提供对比来控制哪些变量被视为分类变量：

m2 = fit(MixedModel,
         @formula(MPG ~ Cyl & Gear + HP & AM + (1 | Gear & AM)),
         mtcars,
         REML=true,
         contrasts = Dict(:AM => DummyCoding(), :Gear => DummyCoding()))

说明

在R和Julia（使用StatsModels.jl）中指定回归公式有两个主要区别。

1. 在Julia中，您使用&来创建术语之间的交互作用，但在R中您使用:。
2. 在Julia中，您必须明确说明您正在使用@formula宏定义公式。在R中，通过“非标准评估”的魔力，任何事物都可能是宏。

因此，您模型的等效Julia公式为

@formula(mpg ~ cyl&gear + hp&am + (1|gear&am))

使用 categorical! 将 DataFrame 中的列转换为分类变量，请参阅 DataFrames.jl 文档。

- Dave Kleinschmidt