QSU Research Methods Seminar | 11 October 2022

Mixed-effects models

Mixed-effects models

Mixed-effects models

Mixed-effects models

Mixed-effects models

lme4

MixedModels

jglmm 

~

~

reaction ~ 1 + days

reaction ~ 1 + days + (1 | subj)

reaction ~ 1 + days + (0 + days | subj)

reaction ~ 1 + days + (1 + days | subj)

Usage: setup

using Pkg
Pkg.add("MixedModels")
Pkg.add("DataFrames")
Pkg.add("StatsModels")
options(JULIA_HOME = "/Applications/Julia-1.8.app/Contents/Resources/julia/bin/")
library(jglmm)
jglmm_setup()

Technical note: the first time you fit a model in each session, it will take a long time because Julia is doing just-in-time compilation, subsequent fits will be much faster

Usage: fitting

Linear mixed-effects model

lmer(reaction ~ 1 + days + (1 + days | subj), sleepstudy)

Usage: fitting

Linear mixed-effects model

lmer(reaction ~ 1 + days + (1 + days | subj), sleepstudy)
## Linear mixed model fit by REML ['lmerMod']
## Formula: reaction ~ 1 + days + (1 + days | subj)
##    Data: sleepstudy
## 
## REML criterion at convergence: 1744
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -3.954 -0.463  0.023  0.463  5.179 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  subj     (Intercept) 612.1    24.74        
##           days         35.1     5.92    0.07
##  Residual             654.9    25.59        
## Number of obs: 180, groups:  subj, 18
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   251.41       6.82   36.84
## days           10.47       1.55    6.77
## 
## Correlation of Fixed Effects:
##      (Intr)
## days -0.138

Usage: fitting

Linear mixed-effects model

lmer(reaction ~ 1 + days + (1 + days | subj), sleepstudy)
jglmm(reaction ~ 1 + days + (1 + days | subj), sleepstudy)

Usage: fitting

Linear mixed-effects model

lmer(reaction ~ 1 + days + (1 + days | subj), sleepstudy)
jglmm(reaction ~ 1 + days + (1 + days | subj), sleepstudy)
## Julia Object of type LinearMixedModel{Float64}.
## Linear mixed model fit by maximum likelihood
##  reaction ~ 1 + days + (1 + days | subj)
##    logLik   -2 logLik     AIC       AICc        BIC    
##   -875.9697  1751.9393  1763.9393  1764.4249  1783.0971
## 
## Variance components:
##             Column    Variance Std.Dev.   Corr.
## subj     (Intercept)  565.51067 23.78047
##          days          32.68212  5.71683 +0.08
## Residual              654.94145 25.59182
##  Number of obs: 180; levels of grouping factors: 18
## 
##   Fixed-effects parameters:
## ──────────────────────────────────────────────────
##                 Coef.  Std. Error      z  Pr(>|z|)
## ──────────────────────────────────────────────────
## (Intercept)  251.405      6.63226  37.91    <1e-99
## days          10.4673     1.50224   6.97    <1e-11
## ──────────────────────────────────────────────────

Usage: fitting

Linear mixed-effects model

lmer(reaction ~ 1 + days + (1 + days | subj), sleepstudy)
jglmm(reaction ~ 1 + days + (1 + days | subj), sleepstudy)
## Julia Object of type LinearMixedModel{Float64}.
## Linear mixed model fit by maximum likelihood
##  reaction ~ 1 + days + (1 + days | subj)
##    logLik   -2 logLik     AIC       AICc        BIC    
##   -875.9697  1751.9393  1763.9393  1764.4249  1783.0971
## 
## Variance components:
##             Column    Variance Std.Dev.   Corr.
## subj     (Intercept)  565.51067 23.78047
##          days          32.68212  5.71683 +0.08
## Residual              654.94145 25.59182
##  Number of obs: 180; levels of grouping factors: 18
## 
##   Fixed-effects parameters:
## ──────────────────────────────────────────────────
##                 Coef.  Std. Error      z  Pr(>|z|)
## ──────────────────────────────────────────────────
## (Intercept)  251.405      6.63226  37.91    <1e-99
## days          10.4673     1.50224   6.97    <1e-11
## ──────────────────────────────────────────────────

Technical note:
MixedModels by default uses maximum likelihood estimation with optimizer bobyqa
so the more equivalent call is lmer(REML = FALSE, control = lmerControl(optimizer = "bobyqa"))

Usage: fitting

Generalized linear mixed-effects model (logistic)

Usage: fitting

Generalized linear mixed-effects model (logistic)

Usage: fitting

Generalized linear mixed-effects model (logistic)

glmer(incid / hsz ~ period + (1 | herd), data = cbpp,
      weights = cbpp$hsz, family = "binomial")

Usage: fitting

Generalized linear mixed-effects model (logistic)

glmer(incid / hsz ~ period + (1 | herd), data = cbpp,
      weights = cbpp$hsz, family = "binomial")
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: incid/hsz ~ period + (1 | herd)
##    Data: cbpp
## Weights: cbpp$hsz
## 
##      AIC      BIC   logLik deviance df.resid 
##      194      204      -92      184       51 
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -2.382 -0.789 -0.203  0.514  2.879 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  herd   (Intercept) 0.412    0.642   
## Number of obs: 56, groups:  herd, 15
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -1.398      0.231   -6.05  1.5e-09 ***
## period2       -0.992      0.303   -3.27  0.00107 ** 
## period3       -1.128      0.323   -3.49  0.00047 ***
## period4       -1.580      0.422   -3.74  0.00018 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr) perid2 perid3
## period2 -0.363              
## period3 -0.340  0.280       
## period4 -0.260  0.213  0.198

Usage: fitting

Generalized linear mixed-effects model (logistic)

glmer(incid / hsz ~ period + (1 | herd), data = cbpp,
      weights = cbpp$hsz, family = "binomial")
jglmm(incid / hsz ~ period + (1 | herd), data = cbpp,
      weights = cbpp$hsz, family = "binomial")

Usage: fitting

Generalized linear mixed-effects model (logistic)

glmer(incid / hsz ~ period + (1 | herd), data = cbpp,
      weights = cbpp$hsz, family = "binomial")
jglmm(incid / hsz ~ period + (1 | herd), data = cbpp,
      weights = cbpp$hsz, family = "binomial")
## Julia Object of type GeneralizedLinearMixedModel{Float64, Binomial{Float64}}.
## Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 1)
##   :(incid / hsz) ~ 1 + period + (1 | herd)
##   Distribution: Binomial{Float64}
##   Link: LogitLink()
## 
##    logLik    deviance     AIC       AICc        BIC    
##    -92.0263   100.0959   194.0526   195.2526   204.1793
## 
## Variance components:
##         Column   Variance Std.Dev. 
## herd (Intercept)  0.412501 0.642262
## 
##  Number of obs: 56; levels of grouping factors: 15
## 
## Fixed-effects parameters:
## ───────────────────────────────────────────────────
##                  Coef.  Std. Error      z  Pr(>|z|)
## ───────────────────────────────────────────────────
## (Intercept)  -1.39852     0.227891  -6.14    <1e-09
## period: 2    -0.992353    0.305386  -3.25    0.0012
## period: 3    -1.12868     0.326049  -3.46    0.0005
## period: 4    -1.58032     0.428794  -3.69    0.0002
## ───────────────────────────────────────────────────

Usage: fitting

Generalized linear mixed-effects model (logistic)

jglmm(incid / hsz ~ period + (1 | herd), data = cbpp,
      weights = cbpp$hsz, family = "binomial", link = "logit",
      contrasts = list(period = "effects"))

Usage: fitting

Generalized linear mixed-effects model (logistic)

jglmm(incid / hsz ~ period + (1 | herd), data = cbpp,
      weights = cbpp$hsz, family = "binomial", link = "logit",
      contrasts = list(period = "effects"))
## Julia Object of type GeneralizedLinearMixedModel{Float64, Binomial{Float64}}.
## Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 1)
##   :(incid / hsz) ~ 1 + period + (1 | herd)
##   Distribution: Binomial{Float64}
##   Link: LogitLink()
## 
##    logLik    deviance     AIC       AICc        BIC    
##    -92.0263   100.0959   194.0526   195.2526   204.1793
## 
## Variance components:
##         Column   Variance Std.Dev. 
## herd (Intercept)  0.412525 0.642281
## 
##  Number of obs: 56; levels of grouping factors: 15
## 
## Fixed-effects parameters:
## ─────────────────────────────────────────────────────
##                   Coef.  Std. Error       z  Pr(>|z|)
## ─────────────────────────────────────────────────────
## (Intercept)  -2.32386      0.221802  -10.48    <1e-24
## period: 2    -0.0670095    0.232325   -0.29    0.7730
## period: 3    -0.203318     0.246292   -0.83    0.4091
## period: 4    -0.655016     0.315527   -2.08    0.0379
## ─────────────────────────────────────────────────────

Usage: fitting

Generalized linear mixed-effects model (logistic)

jglmm(incid / hsz ~ period + (1 | herd), data = cbpp,
      weights = cbpp$hsz, family = "binomial",
      contrasts = list(period = "effects"))
   family %in% c("bernoulli", "binomial", "gamma", "normal", "poisson")
     link %in% c("cauchit", "cloglog", "identity", "inverse", "logit",
                 "log", "probit", "sqrt")
contrasts %in% c("dummy", "effects", "helmert")

Usage: extraction

sm <- jglmm(reaction ~ 1 + days + (1 + days | subj), sleepstudy)
glance(sm)

Usage: extraction

Tidy form of estimated effects:

tidy(sm)

Usage: extraction

Observations augmented with fitted values and residuals:

augment(sm)

Usage: extraction

Conditional modes of random effects:

ranef(sm)$subj

Usage: extraction

Other methods:

nobs()         # number of observations
sigma()        # sigma (residual standard deviation)
logLik()       # log likelihood and degrees of freedom
extractAIC()   # AIC and BIC
deviance()     # deviance

fixef()        # fixed effects estimates
fitted()       # fitted values
vcov()         # variance-covariance matrix
hatvalues()    # diagonal elements of the hat matrix

predict()      # predicted values
simulate()     # simulated values

Implementation

Communicate between R and Julia using JuliaCall package

  • julia_setup() to start Julia process
  • julia_library() to load library in Julia
  • julia_assign("variable", value) to pass value from R to Julia
  • julia_eval("command") to run command in Julia and return value to R

Implementation: setup

jglmm_setup <- function() {
  julia_setup()
  julia_library("MixedModels")
  julia_library("DataFrames")
  julia_library("StatsModels")
}

Implementation: fitting

jglmm <- function(formula, data, family = "normal", link = NULL,
                  weights = NULL, contrasts = NULL, REML = FALSE) {}
  • pass formula, data, weights to Julia using julia_assign()
  • choose between MixedModels.LinearMixedModel and MixedModels.GeneralizedLinearMixedModel:
    • if family == "normal" & link == "identity" then LinearMixedModel
    • otherwise GeneralizedLinearMixedModel
  • convert arguments from R to Julia:
    • family: "binomial" to Binomial()
    • link: "logit" to LogitLink()
    • contrasts: list(period = "effects") to Dict(:period => EffectsCoding())
  • fit model:
    • model <- julia_eval("fit(MixedModels.GeneralizedLinearMixedModel, formula, data, wts = weights, Binomial(), LogitLink(), contrasts = Dict(:period => EffectsCoding()))")
  • return list(formula, data, model) with class "jglmm"

Implementation: extraction

logLik(object, ...)
#' @importFrom stats logLik
#' @param object An object of class `jglmm`, as returned by `jglmm()`.
#' @param ... Optional additional arguments, currently none are used.
#' @export
logLik.jglmm <- function(object, ...) {
  julia_assign("model", object$model)
  loglik <- julia_eval("loglikelihood(model);")
  attr(loglik, "df") <- julia_eval("dof(model);")
  class(loglik) <- "logLik"
  loglik
}
model_lmer <- lmer(reaction ~ 1 + days + (1 + days | subj), sleepstudy)
logLik(model_lmer)
## 'log Lik.' -872 (df=6)
model_jglmm <- jglmm(reaction ~ 1 + days + (1 + days | subj), sleepstudy)
logLik(model_jglmm)
## 'log Lik.' -876 (df=6)

Implementation: extraction

#' @importFrom generics augment
#' @param x An object of class `jglmm`, as returned by `jglmm()`.
#' @param ... Optional additional arguments, currently none are used.
#' @export
augment.jglmm <- function(x, ...) {
  julia_assign("model", x$model)
  fits <- julia_eval("fitted(model)")
  resids <- julia_eval("residuals(model)")
  x$data |> dplyr::mutate(.fitted = fits, .resid = resids)
}

Implementation

NAMESPACE:

export(jglmm)
export(jglmm_setup)
S3method(augment,jglmm)
S3method(deviance,jglmm)
S3method(extractAIC,jglmm)
S3method(fitted,jglmm)
S3method(fixef,jglmm)
S3method(glance,jglmm)
S3method(hatvalues,jglmm)
S3method(logLik,jglmm)
S3method(nobs,jglmm)
S3method(predict,jglmm)
S3method(print,jglmm)
S3method(ranef,jglmm)
S3method(sigma,jglmm)
S3method(simulate,jglmm)
S3method(summary,jglmm)
S3method(tidy,jglmm)
S3method(vcov,jglmm)