# Discrepancy between lme4 and GLM.jl

**URL:** <https://discourse.julialang.org/t/discrepancy-between-lme4-and-glm-jl/85601>\
**Category:** Machine Learning\
**Tags:** statistics, linear-regression\
**Created:** [August 11, 2022, 2:55am UTC](https://discourse.julialang.org/t/discrepancy-between-lme4-and-glm-jl/85601 "2022-08-11T02:55:42Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![kevbonham](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kevbonham/32/216165_2.png) [@kevbonham](https://discourse.julialang.org/u/kevbonham)\
**Post date:** [August 11, 2022, 2:55am UTC](https://discourse.julialang.org/t/discrepancy-between-lme4-and-glm-jl/85601/1 "2022-08-11T02:55:42Z")

</div>

There I was, happily running thousands of linear models with GLM.jl, but the results didn’t seem quite right. For one thing, all of my p-values were exceptionally low (\>90% of features significant, even after FDR correction). For another, all of the coeficients were negative (which is unexpected). I ran my data through a standard tool in my field, that’s built on top of `lme4` in R, and got **very** different results - frankly more like I’d expect.

I did some digging, and indeed, the results with the exact same data and exact same model are giving me different results.

First a verification, here I’m running a very simple model with the data, `y ~ x`, and both GLM.jl and lme4 give the exact same result:

> **Code in julia (using RCall for the R bits)**
>
> ```julia
> using CSV
> using DataFrames
> using DataFrames.PrettyTables
> using GLM
> using RCall
> 
> dat = CSV.read("/home/kevin/Downloads/dat.csv", DataFrame)
> 
> modjl = lm(@formula(y ~ x), dat)
> modjldf = DataFrame(coeftable(modjl))
> pretty_table(select(modjldf, 1:5))
> 
> R"library('lme4')"
> @rput dat
> R"modR <- lm(y ~ x, dat)"
> R"summary(modR)"
> 
> ```

GLM results:

```julia
┌─────────────┬──────────┬────────────┬──────────┬─────────────┐
│ Name │ Coef. │ Std. Error │ t │ Pr(>|t|) │
│ String │ Float64 │ Float64 │ Float64 │ Float64 │
├─────────────┼──────────┼────────────┼──────────┼─────────────┤
│ (Intercept) │ 104.43 │ 1.64392 │ 63.5247 │ 4.0734e-169 │
│ x │ 0.248086 │ 0.254027 │ 0.976613 │ 0.329598 │
└─────────────┴──────────┴────────────┴──────────┴─────────────┘

```

`lme4` results:

```julia
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 104.4295 1.6439 63.525 <2e-16 ***
x 0.2481 0.2540 0.977 0.33

```

But when I throw in a covariate

> **Code**
>
> ```julia
> modjl2 = lm(@formula(y ~ x + c1), dat)
> modjl2df = DataFrame(coeftable(modjl2))
> pretty_table(select(modjl2df, 1:5))
> 
> R"modR2 <- lm(y ~ x + c1, dat)"
> R"summary(modR2)"
> 
> ```

GLM.jl:

```julia
┌─────────────┬─────────────┬────────────┬───────────┬──────────────┐
│ Name │ Coef. │ Std. Error │ t │ Pr(>|t|) │
│ String │ Float64 │ Float64 │ Float64 │ Float64 │
├─────────────┼─────────────┼────────────┼───────────┼──────────────┤
│ (Intercept) │ 0.0 │ NaN │ NaN │ NaN │
│ x │ -9.48368 │ 0.272352 │ -34.8213 │ 3.99587e-104 │
│ c1 │ -4.08581e-8 │ 2.57778e-7 │ -0.158501 │ 0.874176 │
└─────────────┴─────────────┴────────────┴───────────┴──────────────┘

```

lme4:

```julia
Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.147e+02 7.949e+00 14.425 < 2e-16 ***
x 4.362e-01 7.181e-01 0.607 0.544088
c1 -6.721e-07 2.006e-07 -3.351 0.000915 ***

```

Notice in the GLM model, the inflated (and sign-reversed) coefficient, and the wildly low p-value. Also note the intercept going to 0 and all the other stats are NaN… FWIW - this only happens with certain variables. The data table also contains `c2`, which is a categorical variable (encoded with integers and a few missing values and if I only do `y ~ x + c2`, there’s no discrepancy.

Data table [is here](https://gist.github.com/kescobo/6853789f94cf83b0b80737c4b4700de1) if anyone wants to play. I’m stumped…

---

<div class="post-metadata">

**Author:** ![kevbonham](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kevbonham/32/216165_2.png) [@kevbonham](https://discourse.julialang.org/u/kevbonham)\
**Post date:** [August 11, 2022, 3:13am UTC](https://discourse.julialang.org/t/discrepancy-between-lme4-and-glm-jl/85601/2 "2022-08-11T03:13:04Z")

</div>

Oh FFS, it’s the size of the numbers in c1… (on the order of 1e7). If I just divide them by 1e6, the discrepancy goes away. This is a bug, right?

here’s a MWE:

> **Code**
>
> ```julia
> using DataFrames
> using DataFrames.PrettyTables
> using GLM
> using RCall
> using Distributions
> 
> df = DataFrame(y = rand(Normal(), 200), x = rand(Normal(), 200), c1 = rand(Normal(1e7, 1e6), 200))
> df.c2 = df.c1 ./ 1e6
> 
> modjl = DataFrame(coeftable(lm(@formula(y ~ x + c1), df)))
> modjl2 = DataFrame(coeftable(lm(@formula(y ~ x + c2), df)))
> 
> @rput df
> R"modR <- lm(y ~ x + c1, df)"
> R"modR2 <- lm(y ~ x + c2, df)"
> 
> pretty_table(select(modjl, 1:5))
> pretty_table(select(modjl2, 1:5))
> R"summary(modR)"
> R"summary(modR2)"
> 
> ```

GLM:

```julia
┌─────────────┬─────────────┬────────────┬───────────┬──────────┐
│ Name │ Coef. │ Std. Error │ t │ Pr(>|t|) │
│ String │ Float64 │ Float64 │ Float64 │ Float64 │
├─────────────┼─────────────┼────────────┼───────────┼──────────┤
│ (Intercept) │ 0.0 │ NaN │ NaN │ NaN │
│ x │ 0.0885528 │ 0.0666665 │ 1.32829 │ 0.18561 │
│ c1 │ -4.56781e-9 │ 6.54316e-9 │ -0.698105 │ 0.485931 │
└─────────────┴─────────────┴────────────┴───────────┴──────────┘
┌─────────────┬───────────┬────────────┬───────────┬──────────┐
│ Name │ Coef. │ Std. Error │ t │ Pr(>|t|) │
│ String │ Float64 │ Float64 │ Float64 │ Float64 │
├─────────────┼───────────┼────────────┼───────────┼──────────┤
│ (Intercept) │ -0.180284 │ 0.663792 │ -0.271597 │ 0.786216 │
│ x │ 0.088447 │ 0.0668241 │ 1.32358 │ 0.187177 │
│ c2 │ 0.0131475 │ 0.0655554 │ 0.200556 │ 0.841253 │
└─────────────┴───────────┴────────────┴───────────┴──────────┘

```

lme4:

```julia
Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.803e-01 6.638e-01 -0.272 0.786
x 8.845e-02 6.682e-02 1.324 0.187
c1 1.315e-08 6.556e-08 0.201 0.841

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.18028 0.66379 -0.272 0.786
x 0.08845 0.06682 1.324 0.187
c2 0.01315 0.06556 0.201 0.841

```

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [August 11, 2022, 4:37am UTC](https://discourse.julialang.org/t/discrepancy-between-lme4-and-glm-jl/85601/3 "2022-08-11T04:37:10Z")

</div>

> [@kevbonham](#):
>
> This is a bug, right?

I’d say yes for sure. File a bug!

---

<div class="post-metadata">

**Author:** ![kevbonham](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kevbonham/32/216165_2.png) [@kevbonham](https://discourse.julialang.org/u/kevbonham)\
**Post date:** [August 11, 2022, 2:03pm UTC](https://discourse.julialang.org/t/discrepancy-between-lme4-and-glm-jl/85601/4 "2022-08-11T14:03:21Z")

</div>

It looks like it’s this issue: [https://github.com/JuliaStats/GLM.jl/issues/426](https://github.com/JuliaStats/GLM.jl/issues/426)

Using `dropcollinear=false` also fixes it.

---

<div class="post-metadata">

**Author:** ![palday](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/palday/32/12640_2.png) [@palday](https://discourse.julialang.org/u/palday)\
**Post date:** [August 11, 2022, 3:26pm UTC](https://discourse.julialang.org/t/discrepancy-between-lme4-and-glm-jl/85601/5 "2022-08-11T15:26:23Z")

</div>

The core issue here is the underlying numerics, which is a challenging topic. Numerical stability is unfortunately not scale invariant. I wish more introductory stats courses taught practical aspects of statistical computation like rescaling instead of spending time focusing on arcana from pre-computer-era statistical practice. 😕

Just FYI: you’ve loaded but are not using `lme4` in the R section. `lm` is part of base R (technically part of the `stats` package, but that package is autoloaded at R startup).

`lmer` and `glmer` are the key functions from `lme4` and the comparable Julia package is `MixedModels.jl`.

---

<div class="post-metadata">

**Author:** ![kevbonham](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kevbonham/32/216165_2.png) [@kevbonham](https://discourse.julialang.org/u/kevbonham)\
**Post date:** [August 11, 2022, 3:44pm UTC](https://discourse.julialang.org/t/discrepancy-between-lme4-and-glm-jl/85601/6 "2022-08-11T15:44:08Z")

</div>

> [@palday](#):
>
> Just FYI: you’ve loaded but are not using `lme4` in the R section. `lm` is part of base R (technically part of the `stats` package, but that package is autoloaded at R startup).

Can you tell I don’t use R much 😳, thanks for the heads up!

---

<div class="post-metadata">

**Author:** ![Eric](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eric/32/25565_2.png) [@Eric](https://discourse.julialang.org/u/Eric)\
**Post date:** [August 15, 2022, 7:45pm UTC](https://discourse.julialang.org/t/discrepancy-between-lme4-and-glm-jl/85601/7 "2022-08-15T19:45:17Z")

</div>

Yes, definitely with the collinear option.  
It is a problem mainly as it makes one think that numerical issue is the underlying cause here. And while indeed the data could be rescaled more meaningfully, this is not the main issue. In my view, the main issue in your model is with the variable _x_, which does not bring anything to your model (a model with only an intercept and no variable will work just as well as one with only the _x_ variable; the R² value will give you hint too).  
Now other solutions like R, SAS, and Julia (with for instance, the \ operator) will give you accurate coefficients, so I don’t think a numerical issue is at play in this specific case.

---

<div class="post-metadata">

**Author:** ![palday](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/palday/32/12640_2.png) [@palday](https://discourse.julialang.org/u/palday)\
**Post date:** [November 1, 2022, 4:18pm UTC](https://discourse.julialang.org/t/discrepancy-between-lme4-and-glm-jl/85601/8 "2022-11-01T16:18:30Z")

</div>

@Eric numerical stability is a property of the method and how it interacts with the data. Other solutions in R, SAS and Julia uses different methods, so saying that they work is a statement about the stability of those methods and not about the current default Cholesky method used in GLM.jl. Collinearity is well defined in theory, but less well defined in practice, where it is sensitive to rounding error and the numerical tolerances used in e.g. the pivoting methods used to handle it. That’s what’s happening with the intercept in the one model: it’s being pivoted out as part of the collinearity handling.

Relatedly, this is why GLM.jl will switch to using a QR decomposition method as its default solve method soon: we’ve had a number of issues reported by people running into numerical issues on the Cholesky decomposition.