Fitting Mixed Effects Models – Python, Julia or R?

By: Dean Markwick's Blog -- Julia

Re-posted from: https://dm13450.github.io/2022/01/06/Mixed-Models-Benchmarking.html

I’m benchmarking how long it takes to fit a mixed
effects model using lme4 in R, statsmodels in Python, plus
showing how MixedModels.jl in Julia is also a viable option.


Enjoy these types of posts? Then you should sign up for my newsletter. It’s a short monthly recap of anything and everything I’ve found interesting recently plus
any posts I’ve written. So sign up and stay informed!






Data science is always up for debating whether R or Python is the better
language when it comes to analysing some data. Julia has been
making its case as a viable alternative as well. In most cases you can
perform a task in all three with a little bit of syntax
adjustment, so no need for any real commitment. Yet, I’ve
recently been having performance issues in R with a mixed model.

I have a dataset with 1604 groups for a random effect that has been
grinding to a halt when fitting in R using lme4. The team at lme4
have a vignette titled
performance tips
which at the bottom suggests using Julia to speed things up. So I’ve
taken it upon myself to benchmark the basic model-fitting performances
to see if there is a measurable difference. You can use this post as
an example of fitting a mixed effects model in Python, R and Julia.

The Setup

In our first experiment, I am using the palmerspenguins dataset to
fit a basic linear model. I’ve followed the most basic method in all
three languages, using what the first thing in Google
displays.

The dataset has 333 observations with 3 groups for the random effects
parameter.

I’ll make sure all the parameters are close across the three
languages before benchmarking the performance, again, using what
Google says is the best approach to time some code.

I’m running everything on a Late 2013 Macbook. 2.6GHz i5 with 16GB of
RAM. I’m more than happy to repeat this on an M1 Macbook if someone is
willing to sponsor the purchase!

Now onto the code.

Mixed Effects Models in R with lme4

We will start with R as that is where the dataset comes from. Loading
up the palmerspenguins package and filtering out the NaN in the
relevant columns provide a consistent dataset for the other
languages.

  • R – 4.1.0
  • lme4 – 1.1-27.1
require(palmerpenguins)
require(lme4)
require(microbenchmark)

testData <- palmerpenguins::penguins %>%
                    drop_na(sex, species, island, bill_length_mm)
lmer(bill_length_mm ~ (1 | species) + sex + 1,  testData)
  • Intercept – 43.211
  • sex:male – 3.694
  • species variance – 29.55

This testData gets saved as a csv for the rest of the languages
to read and use in the benchmarking.

To benchmark the function I use the
microbenchmark
package. It is an excellent and lightweight way of quickly working out
how long a function takes to run.

microbenchmark(lmer(bill_length_mm ~ (1 | species) + sex + 1,  testData), times = 1000)

This outputs the relevant quantities of the benchmarking times.

Mixed Effects Models Python with statsmodels

  • Python 3.9.4
  • statmodels 0.13.1
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
import timeit as tt

modelData = pd.read_csv("~/Downloads/penguins.csv")

md = smf.mixedlm("bill_length_mm  ~ sex + 1", modelData, groups=modelData["species"])
mdf = md.fit(method=["lbfgs"])

Which gives us parameter values:

  • Intercept – 43.211
  • sex:male – 3.694
  • Species variance: 29.496

When we benchmark the code, we define a specific function and
repeatedly run it 10000 times. This is all contained in the timeit
module, part of the Python standard library.

def run():
    md = smf.mixedlm("bill_length_mm  ~ sex + 1", modelData, groups=modelData["species"])
    mdf = md.fit(method=["lbfgs"])

times = tt.repeat('run()', repeat = 10000, setup = "from __main__ import run", number = 1)

We’ll be taking the mean, median and, range of the times array.

Mixed Effects Models in Julia with MixedModels.jl

  • Julia 1.6.0
  • MixedModels 4.5.0

Julia follows the R syntax very closely, so this needs little
explanation.

using DataFrames, DataFramesMeta, CSV, MixedModels
using BenchmarkTools

modelData = CSV.read("/Users/deanmarkwick/Downloads/penguins.csv",
                                     DataFrame)

m1 = fit(MixedModel, @formula(bill_length_mm ~ 1 + (1|species) + sex), modelData)
  • Intercept – 43.2
  • sex:male coefficient – 3.694
  • group variance of – 19.68751

We use the
BenchmarkTools.jl
package to run the function 10,000 times.

@benchmark fit(MixedModel, @formula(bill_length_mm ~ 1 + (1|species) + sex), $modelData)

As a side note, if you run this benchmarking code in a Jupyter
notebook, you get this beautiful output from the BenchmarkTools
package. Gives you a lovely overview of all the different metrics and
the distribution on the running times.

Julia benchmarking screenshot

Timing Results

All the parameters are close enough, how about the running times?

In milliseconds:

Language Mean Median Min Max
Julia 0.482 0.374 0.320 34
Python 340 260 19 1400
R 29.5 24.5 20.45 467

Julia blows both Python and R out of the water. About 60 times
faster.

I don’t think Python is that slow in practice, I think it is more of
an artefact of the benchmarking code that doesn’t behave in the
same way as Julia and R.

What About Bigger Data and More Groups?

What if we increased the scale of the problem and also the number of
groups in the random effects parameters?

I’ll now fit a Poisson mixed model to some football data. I’ll
be modeling the goals scored by each team as a Poisson variable, with
a fixed effect of whether the team played at home or not and random
effects for the team and another random effect of the opponent.

This new data set is from
football-data.co.uk and has 98,242
rows with 151 groups in the random effects. Much bigger than the
Palmer Penguins dataset.

Now, poking around the statsmodels documentation, there doesn’t
appear to be a way to fit this model in a frequentist
way. The closest is the PoissonBayesMixedGLM, which isn’t comparable
to the R/Julia methods. So in this case we will be dropping Python
from the analysis. If I’m wrong, please let me know in the comments
below and I’ll add it to the benchmarking.

With generalised linear models in both R and Julia, there are
additional parameters to help speed up the fitting but at the expense of
the parameter accuracy. I’ll be testing these parameters to judge how
much of a tradeoff there is between speed and accuracy.

R

The basic mixed-effects generalised linear model doesn’t change much from the
above in R.

glmer(Goals ~ Home + (1 | Team) + (1 | Opponent), data=modelData, family="poisson")

The documentation states that you can pass nAGQ=0 to speed up the
fitting process but might lose some accuracy. So our fast version of
this model is simply:

glmer(Goals ~ Home + (1 | Team) + (1 | Opponent), data=modelData, family="poisson", nAGQ = 0)

Julia

Likewise for Julia hardly any difference in fitting this type of Poisson model.

fit(MixedModel, @formula(Goals ~ Home + (1 | Team) + (1 | Opponent)), footballData, Poisson())

And even mode simply, there is a fast parameter to use which speeds
up the fitting.

fit(MixedModel, @formula(Goals ~ Home + (1 | Team) + (1 | Opponent)), footballData, Poisson(), fast = true)

Big Data Results

Let’s check the fitted coefficients.

Method Intercept Home \(\sigma _\text{Team}\) \(\sigma_\text{Opponent}\)
R slow 0.1345 0.2426 0.2110 0.2304
R fast 0.1369 0.2426 0.2110 0.2304
Julia slow 0.13455 0.242624 0.211030 0.230415
Julia fast 0.136924 0.242625 0.211030 0.230422

The parameters are all very similar, showing that for this parameter
specification the different speed flags do not change the coefficient
results, which is good. But for any specific model, you should verify
on a subsample at least to make sure the flags don’t change anything.

Now, what about speed.

Language Additional Parameter Mean Median Min Max
Julia 11.151 10.966 9.963 16.150
Julia fast=true 5.94 5.924 4.98 8.15
R 35.4 33.12 24.33 66.48
R nAGQ=0 8.06 7.99 7.37 9.56

So setting fast=true gives a 2x speed boost in Julia which is
nice. Likewise, setting nAGQ=0 in R improves the speed by almost 3x over
the default. Julia set to fast = true is the quickest, but I’m
surprised that R can get close with its speed-up parameter.

Conclusion

If you are fitting a large mixed-effects model with lots of groups
hopefully, this convinces you that Julia is the way forward. The syntax
for fitting the model is equivalent, so you can do all your
preparation in R before importing the data into Julia to do the model
fitting.

Introduction to Julia

By: Jaafar Ballout

Re-posted from: https://www.supplychaindataanalytics.com/introduction-to-julia/

After having introduced in one of my previous posts optimization and linear programming, I explain the basics of the Julia programming language in this article. The article represents a tutorial and is based on the official Julia documentation. My tutorial covers the key aspects that I will later use, in upcoming blog posts, for solving supply chain and operation research problems (e.g. network problems). My tutorial also highlights some major syntactic and functional differences when compared to other popular programming languages (namely Python and Matlab).

Defining variables in Julia

A variable is a name associated with a value and saved in the computer memory. Assigning a value for the variable is done using the = operator. Unicode can be used as a variable name. Most Julia editors support LaTeX syntax which can be used to create Unicode characters. The double quotes ” are used for strings while the single quote ‘ is used for a character.

Here are some examples demonstrating the above explanation.

In[]:
x = 7 # variable name: x; variable value = 7

Out[]:
7
In[]:
y = 10 # variable name: y; variable value = 10

Out[]:
10

Variable names are case-sensitive and have no semantic meaning.

In[]:
Y = 100 

Out[]:
100
In[]:
n = "Jaafar" # variable name: n; variable value = "Jaafar" which is  a string

Out[]:
"Jaafar"
In[]:
Letter = 'a' # character

Out[]:
'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)
In[]:
α = 1 # unicode is easier in Julia: write \alpha and press tab 

Out[]:
1
In[]:
😠=0 # \:angry: and then press <tab>

Out[]: 
0

Integer and floating point numbers

Integers and floating points are the building blocks of mathematical operations. Using the typeof command in Julia I can find the type of any pre-defined variable.

In[]:
x = 1; #semicolon is to avoid printing the variables
y = 10000;
z = 1.1;
d = 5e-3

println("x is ", typeof(x))
println("y is ", typeof(y))
println("z is ", typeof(z))
println("d is ", typeof(d))

Out[]:
x is Int64
y is Int64
z is Float64
d is Float64

Using the Sys.WORD_SIZE, Julia’s internal variable, I can indicate whether the targetted system is 32-bit or 64-bit.

In[]:
# 64-bit system
Sys.WORD_SIZE

Out[]:
64

Main mathematical operations in Julia

The arithmetic operations are similar to Python except for the power. For power operations, Julia uses ^ instead of ** (as known from Python).

Boolean operations are as follows:

  • a && b: a and b
  • a || b: a or b
  • !a: negation
In[]:
a = 1;
b = 3;

# Addition
a + b;

# Subtraction
a - b;

# times
a*b;

# divison
a/b;

# Power 
a^b; # this is different from python where ** is used to raise a to the bth power

#Updating operators 
a+=1; # a = a + 1
b*=2; # b = b * 2

Vectorized operators are very important in linear algebra. Dot operators are used for arrays where elementary operations are performed.

In[]:
[1,2,3].^1 # [1^1,2^1,3^1]

Out[]:
3-element Vector{Int64}:
 1
 2
 3

Basic collections

Tuples, named tuples, and dictionaries

  • Tuples: Ordered immutable collections of elements
  • NamedTuples: Exactly like tuples but also assign a name for each variable
  • Dictionaries: Unordered mutable collections of pairs: key-value
In[]:
# Tuple 
favoritesongs = ("outnumbered", "Power Over Me", "Bad Habits") # elements have the same type

# Tuple
favoritethings= ("yellow", 'j', pi) # elements with different types

# NamedTuple
favoritesongs_named = (a = "outnumbered", b = "Power Over Me", c = "Bad Habits") # it is between a Tuple and Dictionary

# Dictionary
myDict = Dict("name" => "Jaafar", "age" => "twenty", "hobby"=> "biking")

Out[]:
Dict{String, String} with 3 entries:
  "name"  => "Jaafar"
  "hobby" => "biking"
  "age"   => "twenty"
In[]:
# Tuple access
favoritesongs[1] # indexing starts by 1 not 0 (unlike Python)

# NamedTuple access
favoritesongs_named[1] # accessed by index
favoritesongs_named.a # accessed by key

# Dictionary access
myDict["name"] # call the kay to output the value in the dictionary 

Out[]:
"Jaafar"

Vectors, arrays, and matrices

Like any numerical computation language, Julia provides an easy way to handle matrices and their corresponding operations. Unlike Matlab, Julia arrays are indexed with square brackets, A[i,j]. However, similarly to Matlab, indexing starts using one, not zero, which makes it more convenient especially in loops later by using f(i) instead of f(i-1).

  • array: ordered and mutable collection of items of the same type
  • vector: array of dimension one
  • matrix: array of dimension two
  • tensor: array of n-dimension (usually 3 and above)
In[]:
#vector
array_V = [1, 2, 3, 4, 5, 6, 7, 8, 9] # acts as a vector of one-dimension
typeof(array_V)

Out[]:
Vector{Int64} (alias for Array{Int64, 1})
In[]:
#Matrix
array_M = [1 2 3; 4 5 6; 7 8 9] # acts as a matrix of two-dimension

Out[]:
3×3 Matrix{Int64}:
 1  2  3
 4  5  6
 7  8  9
In[]:
# Random vector
vec = rand(9)

Out[]:
9-element Vector{Float64}:
 0.7130265942088201
 0.9545688377050932
 0.7878361868436774
 0.4973658015754845
 0.44265779030703434
 0.01870528656705095
 0.010563833645745424
 0.8906392694739755
 0.5416448302194592
In[]:
# Random matrix 
mat = rand(3,3)

Out[]:
3×3 Matrix{Float64}:
 0.412231  0.0180507  0.862113
 0.534452  0.711949   0.541887
 0.52126   0.894952   0.443401
In[]:
# Random tensor
ten = rand(3,3,3) # three-dimenisonal 

Out[]:
3×3×3 Array{Float64, 3}:
[:, :, 1] =
 0.517095    0.976259  0.114393
 0.00295048  0.759259  0.302369
 0.988611    0.688391  0.438473

[:, :, 2] =
 0.163933  0.138108  0.770564
 0.899507  0.109004  0.577751
 0.63999   0.280642  0.751499

[:, :, 3] =
 0.361409  0.575224  0.525733
 0.858351  0.586987  0.638436
 0.101579  0.447222  0.364909

It is important to note that ranges act like vector. However, specifying a range is easier to code. The syntax is as follows: range_name = start:step:end

In[]:
# Range
r = 1:1:9 

Out[]:
1:1:9
In[]:
collect(r) # transofrm the range output to a vector output (better output)

Out[]:
9-element Vector{Int64}:
 1
 2
 3
 4
 5
 6
 7
 8
 9

More on indices and ranges in Julia

Working with matrices and arrays, in general, requires good command of indexing and slicing operations in Julia. This can be done more easily using the ranges. This is due to their compact code.

In[]:
# Define an array
name_letters = ['j','a','a','f','a','r']

# Index the array
name_letters[1] # returns j

# slice the array using a range: start:end
name_letters[1:3] # returns jaa

# slice the array using a range with step of 2: start:step:end
name_letters[1:3:6] # returns jf

Out[]:
2-element Vector{Char}:
 'j': ASCII/Unicode U+006A (category Ll: Letter, lowercase)
 'f': ASCII/Unicode U+0066 (category Ll: Letter, lowercase)

Printing output in Julia

Although printing the output of the code is simple it is important for e.g. debugging the code. Print commands are also helpful to readers that are trying to understand the function and purpose of the code.

In[]:
# Print on new line
println("Jaafar") ;
println("Ballout");

# Print one same line
print("Jaafar");
print(" Ballout");

Out[]:
Jaafar
Ballout
Jaafar Ballout

Specifying loops and conditions in Julia

Both loops and conditions in Julia require an end command, unlike Python where indentation is enough to end the if-condition, loop, or function-definition.

In[]:
#If condition
if length("Jaafar") > length("Ballout")
    print("Your first name is bigger than your last name")
elseif length("Jaafar") == length("Ballout")
    println("First name and last name have same number of characters")
else
    println("Your first name is smaller than your last name")
end

Out[]:
Your first name is smaller than your last name

The code block below prints each character in my name into a new line.

In[]:
name_letters = ['j','a','a','f','a','r']
# For loop
for i in 1:length(name_letters)
    println(name_letters[i])
end

Out[]:
j
a
a
f
a
r

The code below finds the location or index of a character in my name.

In[]:
#If condition in For Loop
for i in 1:length(name_letters)
    if name_letters[i] == 'a'
        println(i)
    end
end

Out[]:
2
3
5

Defining functions in Julia

Functions, e.g. routines or methods, can be defined in Julia. The linear optimization model, below, is from my previous blog post. In the coding example that follows, I define the objective function as a function in Julia.

Here is how I can define the objective function in Julia:

In[]:
function max(x,y)
    return 5x + 4y # no need to write the multiplication * sign; Julia will understand
end

Out[]:
max (generic function with 1 method)
In[]:
# Optimal soultion of the constrained problem above from previous post 
z = max(3.75,1.25) # optimal value
print(z)

Out[]:
23.75

Import files into Julia

Importing files is very import especially in supply chain management and logistics problems. Because I might use .csv files in future posts, explaining how to import these files in Julia is necessary at this stage. I will be using Pandas.jl which is a Julia interface to the excellent Pandas package in Python.

In[]:
using Pandas
In[]:
df_list = Pandas.read_csv("https://gist.githubusercontent.com/brooksandrew/e570c38bcc72a8d102422f2af836513b/raw/89c76b2563dbc0e88384719a35cba0dfc04cd522/edgelist_sleeping_giant.csv");

I will introduce useful packages like DataFrames.jl and PyPlots in Julia in future work. These packages are very useful for solving network problems as known from supply chain management and operations research.

The post Introduction to Julia appeared first on Supply Chain Data Analytics.