jlme

Julia (mixed-effects) regression modelling from R

https://github.com/yjunechoe/jlme

Repository

Julia (mixed-effects) regression modelling from R

README.Rmd

---
output: github_document
---



```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%"
)
```

# jlme


[![CRAN status](https://www.r-pkg.org/badges/version/jlme)](https://CRAN.R-project.org/package=jlme)
[![Lifecycle: experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://lifecycle.r-lib.org/articles/stages.html#experimental)
[![R-CMD-check](https://github.com/yjunechoe/jlme/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/yjunechoe/jlme/actions/workflows/R-CMD-check.yaml)
[![test-coverage](https://github.com/yjunechoe/jlme/actions/workflows/test-coverage.yaml/badge.svg)](https://github.com/yjunechoe/jlme/actions/workflows/test-coverage.yaml)
[![codecov](https://codecov.io/gh/yjunechoe/jlme/graph/badge.svg?token=d2bLJunNd5)](https://app.codecov.io/gh/yjunechoe/jlme/)


Julia (mixed-effects) regression modelling from R. Powered by the [`{JuliaConnectoR}`](https://github.com/stefan-m-lenz/JuliaConnectoR) R package and Julia libraries [GLM](https://github.com/JuliaStats/GLM.jl), [StatsModels](https://github.com/JuliaStats/StatsModels.jl), and [MixedModels](https://github.com/JuliaStats/MixedModels.jl).

## Zero-setup test drive

As of March 2025, **Google Colab** supports Julia. This means `{jlme}` *just works* out of the box. Try it out in a [demo notebook](https://colab.research.google.com/drive/1ZKF7SrsMYugzY4g5Bd3tLuAk9VPOcP0P?usp=sharing).

## Installation

You can install the development version of `{jlme}` from [GitHub](https://github.com/yjunechoe/jlme) with:

``` r
# install.packages("remotes")
remotes::install_github("yjunechoe/jlme")
```

`{jlme}` is experimental and under active development: see [NEWS.md](https://github.com/yjunechoe/jlme/blob/main/NEWS.md#jlme-development-version) for the latest updates.

## Setup

```{r example}
library(jlme)
jlme_setup()
```

Using `{jlme}` requires a prior installation of the **Julia programming language**, which can be downloaded from either the [official website](https://julialang.org/) or using the command line utility [juliaup](https://github.com/JuliaLang/juliaup).

If you are encountering issues with setting up Julia, please make sure that you're have the latest version (>=1.1.4) of the [`{JuliaConnectoR}`](https://github.com/stefan-m-lenz/JuliaConnectoR) package installed and see `` ?JuliaConnectoR::`Julia-Setup` `` for troubleshooting.

## Usage (table of contents)

-  [Fit models](#fit-models)
-  [Diagnose models](#diagnose-models)
-  [Assess uncertainty](#assess-uncertainty)
-  [Julia interoperability](#julia-interoperability)
-  [Tips and tricks](#tips-and-tricks)
-  [Acknowledgments](#acknowledgments)


## Fit models

[↑Back to table of contents](#usage-table-of-contents)

Once set up, `(g)lm()` and `(g)lmer` complements in Julia are available via `jlm()` and `jlmer()`, respectively.

### Fixed effects models

`jlm()` with `lm()`/`glm()` syntax:

```{r}
# lm(mpg ~ hp, mtcars)
jlm(mpg ~ hp, mtcars)
```

Contrasts in factor columns are preserved: 

```{r}
x <- mtcars

# Sum code `am`
x$am_sum <- factor(x$am)
contrasts(x$am_sum) <- contr.sum(2)
# Helmert code `cyl`
x$cyl_helm <- factor(x$cyl)
contrasts(x$cyl_helm) <- contr.helmert(3)
colnames(contrasts(x$cyl_helm)) <- c("4vs6", "4&6vs8")

jlm(mpg ~ am_sum + cyl_helm, x)
```

### Mixed effects models

`jlmer()` with `lmer()`/`glmer()` syntax:

```{r}
# lme4::lmer(Reaction ~ Days + (Days | Subject), lme4::sleepstudy)
jlmer(Reaction ~ Days + (Days | Subject), lme4::sleepstudy, REML = TRUE)
```

```{r}
# lme4::glmer(r2 ~ Anger + Gender + (1 | id), lme4::VerbAgg, family = "binomial")
jlmer(r2 ~ Anger + Gender + (1 | id), lme4::VerbAgg, family = "binomial")
```


## Diagnose models

[↑Back to table of contents](#usage-table-of-contents)

Supports `{broom}`-style `tidy()` and `glance()` methods for Julia regression models.

### Summarize model fit

Get information about model components with `tidy()`

```{r}
# Note that MixedModels/`jlmer()` defaults to ML (REML=false)
jmod <- jlmer(Reaction ~ Days + (Days | Subject), lme4::sleepstudy)
tidy(jmod)
```

Get goodness-of-fit measures of a model with `glance()`

```{r}
glance(jmod)
```

### Inspect model objects

Check singular fit

```{r}
issingular(jmod)
```

List all properties of a MixedModel object (properties are accessible via `$`)

```{r}
propertynames(jmod)
```

Check optimization summary

```{r}
jmod$optsum
```


## Assess uncertainty

[↑Back to table of contents](#usage-table-of-contents)

Functions `parametricbootstrap()` and `profilelikelihood()` can be used to assess the variability of parameter estimates.

### Parametric bootstrap

Experimental support for [`MixedModels.parametricbootstrap`](https://juliastats.org/MixedModels.jl/stable/bootstrap/) via `parametricbootstrap()`:

```{r}
samp <- parametricbootstrap(jmod, nsim = 100L, seed = 42L)
samp
```

```{r}
tidy(samp)
```

### Profiling

Experimental support for [`MixedModels.profile`](https://juliastats.org/MixedModels.jl/stable/api/#MixedModels.profile-Tuple%7BLinearMixedModel%7D) via `profilelikelihood()`:

```{r}
prof <- profilelikelihood(jmod)
prof
```

```{r}
tidy(prof)
```


## Julia interoperability

[↑Back to table of contents](#usage-table-of-contents)

Functions `jl_get()` and `jl_put()` transfers data between R and Julia.

### Bring Julia objects into R

Example 1: extract PCA of random effects and return as an R list:

```{r}
jmod$rePCA
```

```{r}
jl_get(jmod$rePCA)
```

Example 2: extract fitlog and plot

```{r fitlog}
fitlog <- jl_get(jl("refit!(deepcopy(x); thin=1)", x = jmod)$optsum$fitlog)
thetas <- t(sapply(fitlog, `[[`, 1))
matplot(thetas, type = "o", xlab = "iterations")
```

### Julia session

See information about the running Julia environment (e.g., the list of loaded Julia libraries) with `jlme_status()`:

```{r}
jlme_status()
```

On setup, `{jlme}` loads [GLM](https://github.com/JuliaStats/GLM.jl), [StatsModels](https://github.com/JuliaStats/StatsModels.jl), and [MixedModels](https://github.com/JuliaStats/MixedModels.jl), as well as those specified in `jlme_setup(add)`. Other libraries such as `Random` (required for `parametricbootstrap()`) are loaded on an as-needed basis.

### More with `{JuliaConnectoR}`

While users will typically not need to interact with `{JuliaConnectoR}` directly, it may be useful for extending `{jlme}` features with other packages in the Julia modelling ecosystem. A simple way to do that is to use `juliaImport()`, which creates makeshift bindings to any Julia library.

Here's an example replicating a workflow using [`Effects.empairs`](https://beacon-biosignals.github.io/Effects.jl/dev/emmeans/) for post-hoc pairwise comparisons:

```{r}
# New model: 2 (M/F) by 3 (curse/scold/shout) factorial
jmod2 <- jlmer(
  r2 ~ Gender * btype + (1 | id),
  data = lme4::VerbAgg,
  family = "binomial"
)
jmod2
```

```{r}
library(JuliaConnectoR)
# Install `Effects.jl` (takes a minute)
juliaEval('using Pkg; Pkg.add("Effects")')
Effects <- juliaImport("Effects")
# Call `Effects.empairs` using R syntax `Effects$empairs()`
pairwise <- Effects$empairs(jmod2, dof = glance(jmod2)$df.residual)
pairwise
```

Note that Julia `DataFrame` objects such as the one above can be collected into an R data frame using `as.data.frame()`. This lets you, for example, apply p-value corrections using the familiar `p.adjust()` function in R, though the [option to do that](https://beacon-biosignals.github.io/Effects.jl/dev/emmeans/#Multiple-Comparisons-Correction) exists in Julia as well.

```{r}
pairwise_df <- as.data.frame(pairwise)
cbind(
  pairwise_df[, 1:2],
  round(pairwise_df[, 3:4], 2),
  pvalue = format.pval(p.adjust(pairwise_df[, 7], "bonferroni"), 1)
)
```


## Tips and tricks

[↑Back to table of contents](#usage-table-of-contents)

### Displaying MixedModels

MixedModels.jl supports various [display formats](https://juliastats.org/MixedModels.jl/stable/mime/) for mixed-effects models which are available in `{jlme}` via the `format` argument of `print()`:

```{r, results="asis"}
# Rendered via {knitr} with chunk option `results="asis"`
print(jmod, format = "markdown")
```

### Data type conversion

Be sure to pass integers (vs. doubles) to Julia functions that expect Integer type:

```{r}
jl_put(1)
jl_put(1L)
```

The `Dict` (dictionary) data type is common and relevant for modelling workflows in Julia. `{jlme}` offers `jl_dict()` as an opinionated constructor for that usecase:

```{r}
# Note use of `I()` to protect length-1 values as scalar
jl_dict(a = 1:2, b = "three", c = I(4.5))
```

See `` ?JuliaConnectoR::`JuliaConnectoR-package` `` for a comprehensive list of data type conversion rules.

### Performance (linear algebra backend)

Using [`MKL.jl`](https://github.com/JuliaLinearAlgebra/MKL.jl) or [`AppleAccelerate.jl`](https://github.com/JuliaLinearAlgebra/AppleAccelerate.jl) may improve model fitting performance (but see the system requirements first). This should be supplied to the `add` argument to `jlme_setup()`, to ensure that they are loaded first, prior to other packages.

```{r, eval = FALSE}
# Not run
jlme_setup(add = "MKL", restart = TRUE)
jlme_status() # Should now see MKL loaded here
```

### Performance (data transfer)

If the data is large, a sizable overhead will come from transferring the data from R to Julia. If you are also looking to fit many models to the same data, you should first filter to keep only the columns you need and then use `jl_data()` to send the data to Julia. The Julia object can then be used in place of an R data frame.

```{r}
data_r <- mtcars

# Keep only columns you need + convert with `jl_data()`
data_julia <- jl_data(data_r[, c("mpg", "am")])

jlm(mpg ~ am, data_julia)
```

If your data has custom contrasts, you can use `jl_contrasts()` to also convert that to Julia first before passing it to the model.

```{r}
data_r$am <- as.factor(data_r$am)
contrasts(data_r$am) <- contr.sum(2)

data_julia <- jl_data(data_r[, c("mpg", "am")])
contrasts_julia <- jl_contrasts(data_r)

jlm(mpg ~ am, data_julia, contrasts = contrasts_julia)
```

### Just learn Julia

If you spend non-negligible time fitting regression models for your work, please just [learn Julia](https://julialang.org/learning/)! It's a great high-level language that feels close to R in syntax and its REPL-based workflow.


## Acknowledgments

[↑Back to table of contents](#usage-table-of-contents)

- The [JuliaConnectoR](https://github.com/stefan-m-lenz/JuliaConnectoR) R package for powering the R interface to Julia.

- The [Julia](https://julialang.org/) packages [GLM.jl](https://github.com/JuliaStats/GLM.jl), [StatsModels](https://github.com/JuliaStats/StatsModels.jl), and [MixedModels.jl](https://github.com/JuliaStats/MixedModels.jl) for interfaces to and implementations of (mixed effects) regression models.

- Other R interfaces to MixedModels. It has come to my attention after working on this for a while and publishing it to CRAN that there were prior efforts by [Mika Braginsky](https://github.com/mikabr) in [`{jglmm}`](https://github.com/mikabr/jglmm) and by [Philip Alday](https://github.com/palday) in [`{jlme}`](https://github.com/palday/jlme). These packages/scripts predate `{JuliaConnectoR}` (which I find to be leaner and more robust) and are instead built on top of [`{JuliaCall}`](https://github.com/JuliaInterop/JuliaCall).

```{r, echo = FALSE}
stop_julia()
```

Owner

GitHub Events

Committers

Last synced: about 1 year ago

Issues and Pull Requests

Last synced: about 1 year ago

Packages