glmer vs julia vs INLA

R / modelos mixtos / julia

Hablábamos el otro día mi amigo Carlos y yo sobre los modelos mixtos y el uso de lme4, Stan o INLA. Total, que el problema es que queríamos un atajo que permitiera tener una estimación de los efectos aleatorios en un tiempo menor que lo que queda hasta el fin del universo.

Pues nada, investigando vi que existía una librería en Julia llamada MixedModels y que es del autor de lme4 así que me puse a probar a ver si es verdad el lema de Julia, “tan rápido como C, tan fácil como Python”.

Vamos a la prueba.

En primer lugar cargamos nuestros conocidos datos de la Encuesta de Población Activa, que suelo utilizar por ser muy sencilla de entender y por tener un volumen de datos considerable.

## Libreria
library(tidyverse)
## ── Attaching packages ──────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.1.0       ✔ purrr   0.3.1  
## ✔ tibble  2.0.1       ✔ dplyr   0.8.0.1
## ✔ tidyr   0.8.3       ✔ stringr 1.4.0  
## ✔ readr   1.3.1       ✔ forcats 0.4.0
## ── Conflicts ─────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(MicroDatosEs)
# install.packages("JuliaCall")
library(JuliaCall)

##  cargar datos

# fpath_dropbox

# enlace_dropbox <- "https://www.dropbox.com/s/h8am8g2yk3dq1y2/md_EPA_2018T4.txt?dl=0"


fpath <- "~/Dropbox/Public/datos_4t18/md_EPA_2018T4.txt"

epa <- epa2005(fpath)
## Warning: 607 parsing failures.
##  row     col           expected actual                                            file
## 1641 REPAIRE 1/0/T/F/TRUE/FALSE    115 '~/Dropbox/Public/datos_4t18/md_EPA_2018T4.txt'
## 1963 REGEST  1/0/T/F/TRUE/FALSE    115 '~/Dropbox/Public/datos_4t18/md_EPA_2018T4.txt'
## 1963 REPAIRE 1/0/T/F/TRUE/FALSE    115 '~/Dropbox/Public/datos_4t18/md_EPA_2018T4.txt'
## 2093 REPAIRE 1/0/T/F/TRUE/FALSE    115 '~/Dropbox/Public/datos_4t18/md_EPA_2018T4.txt'
## 2168 REGEST  1/0/T/F/TRUE/FALSE    310 '~/Dropbox/Public/datos_4t18/md_EPA_2018T4.txt'
## .... ....... .................. ...... ...............................................
## See problems(...) for more details.
names(epa) <- tolower(names(epa))
epa <- subset(epa, select = c(prov, edad, nforma, aoi))

recodificacion <- function(dat) {
  dat$aoi[grepl("^Inactivos", dat$aoi)] <- "i"
  dat$aoi[grepl("[O-o]cupados", dat$aoi)] <- "o"
  dat$aoi[grepl("^Parados", dat$aoi)] <- "p"
  
  dat$aoi <- factor(dat$aoi)
  dat$nforma3 <- dat$nforma
  dat$nforma3[dat$nforma == "Analfabetos" |
                dat$nforma == "Educación primaria" |
                dat$nforma == "Educación primaria incompleta"] <-
    "Est primarios o menos"
  dat$nforma3[dat$nforma == "Educación superior"] <-
    "Est. Universitarios"
  dat$nforma3[dat$nforma == "Primera etapa de educación secundaria" |
                dat$nforma == "Segunda etapa de educación secundaria, orientación general" |
                dat$nforma == "Segunda etapa de educación secundaria, orientación profesional"] <-
    "Est. Secundarios"
  
  dat$nforma3 <- factor(dat$nforma3)
  
  dat$gedad <- dat$edad
  dat$gedad[dat$edad == "de 0 A 4 años" |
              dat$edad == "de 5 A 9 años" |
              dat$edad == "de 10 A 15 años"] <- "15 años o menos "
  dat$gedad[dat$edad == "de 16 A 19 años" |
              dat$edad == "de 20 A 24 años" |
              dat$edad == "de 25 A 29 años" |
              dat$edad == "de 30 A 34 años"] <- "De 16 a 34"
  
  dat$gedad[dat$edad == "de 35 A 39 años" |
              dat$edad == "de 40 A 44 años" |
              dat$edad == "de 45 A 49 años" |
              dat$edad == "de 50 A 54 años"] <- "De 35 a 54"
  
  dat$gedad[dat$edad == "de 55 A 59 años" |
              dat$edad == "de 60 A 64 años" |
              dat$edad == "65 o más años"] <- "55 o más"
  
  
  dat$gedad <-
    factor(dat$gedad,
           levels = c("15 años o menos ", "De 16 a 34", "De 35 a 54", "55 o más")
    )
  
  dat
}

epa <- recodificacion(epa)

# eliminar menores de 16 años  
epa <- epa[epa$gedad != "15 años o menos ", ]

# eliminar inactivos
epa <- epa[epa$aoi != "i", ]

epa$paro <- ifelse(epa$aoi=="p", 1,0)

No voy a contar lo que es un modelo mixto, ni coger datos de train y test, tan sólo ejecutar modelos simples y ver qué tarda más, y menos.

Modelo con glmer

Uso nAGQ=0 y optimizer = "nloptwrap" porque he leído por ahí que es lo más rápido.

Pego el cacho dónde lo he leído

" the option “nAGQ=0” tells glmer to ignore estimating those integrals. For some models this has minimal impact on parameter estimates, and this NCAA hockey model is one of those. The second option tells glmer to fit using the “nloptwrap” optimizer (there are several other optimizers available, too), which tends to be faster than the default optimization method."

library(tictoc)
library(lme4)
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following object is masked from 'package:tidyr':
## 
##     expand
tic()
  fit_mixto <- glmer(paro ~ (1 | prov) + (1 | gedad) + (1 | nforma3),
    family = binomial,
    data = epa,
    nAGQ = 0L,
    control = glmerControl(optimizer = "nloptwrap")
  )
toc()
## 10.903 sec elapsed

Summary del modelo

summary(fit_mixto)
## Generalized linear mixed model fit by maximum likelihood (Adaptive
##   Gauss-Hermite Quadrature, nAGQ = 0) [glmerMod]
##  Family: binomial  ( logit )
## Formula: paro ~ (1 | prov) + (1 | gedad) + (1 | nforma3)
##    Data: epa
## Control: glmerControl(optimizer = "nloptwrap")
## 
##      AIC      BIC   logLik deviance df.resid 
##  57620.1  57657.0 -28806.1  57612.1    74564 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.2062 -0.4286 -0.3446 -0.2478  5.7528 
## 
## Random effects:
##  Groups  Name        Variance Std.Dev.
##  prov    (Intercept) 0.1324   0.3638  
##  gedad   (Intercept) 0.1879   0.4334  
##  nforma3 (Intercept) 0.3759   0.6131  
## Number of obs: 74568, groups:  prov, 52; gedad, 3; nforma3, 3
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.6461     0.4367  -3.769 0.000164 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Modelo con MixedModels en julia

Fijándome en este post del creador de la librería y en la docu de la librería, me instalé la librería JuliaCall que permite usar julia dentro de R

library(JuliaCall)
julia_setup()
## Julia version 1.0.3 at location /media/hd2/Descargas/julia/julia-1.0.3/bin will be used.
## Loading setup script for JuliaCall...
## Finish loading setup script for JuliaCall.
julia_library("MixedModels")
julia_assign("form", formula(fit_mixto))
julia_assign("epa", epa)

Ponemos opción fast=true en GeneralizedLinearMixedModels porque la documentación dice que es lo más rápido. Veamos el tiempo que tarda

julia_eval("@elapsed gm1 = fit!(GeneralizedLinearMixedModel(form, epa, Bernoulli()), fast=true)")
## [1] 38.94062

¡Ojo!, utilizando modelos más complejos, glmer pudo hacerlos pero MixedModels en julia, no. La primera vez que julia ejecuta el modelo tarda como unos 40 segundos, porque tiene que compilar no sé qué historia, las sucesivas veces tarda entorno a 9.5 segundos. No está mal, unos 11 segundos de glmer vs 9.5 de MixedModels

Summary del modelo

julia_eval("gm1")
## Julia Object of type GeneralizedLinearMixedModel{Float64}.
## Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 1)
##   Formula: paro ~ (1 | prov) + (1 | gedad) + (1 | nforma3)
##   Distribution: Bernoulli{Float64}
##   Link: LogitLink()
## 
##   Deviance: 57612.1234
## 
## Variance components:
##              Column    Variance   Std.Dev.  
##  prov    (Intercept)  0.13269447 0.36427252
##  gedad   (Intercept)  0.19182336 0.43797643
##  nforma3 (Intercept)  0.38060583 0.61693260
## 
##  Number of obs: 74568; levels of grouping factors: 52, 4, 3
## 
## Fixed-effects parameters:
##              Estimate Std.Error  z value P(>|z|)
## (Intercept)  -1.64609  0.440015 -3.74097  0.0002

Modelo con INLA

Con INLA se pueden especificar muchas cosas, voy a poner un modelo sencillo, seguro que se puede optimizar un montón, pero no me apetece ahora mismo.

library(INLA)
## Loading required package: sp
## Loading required package: splines
family1 <- "binomial"

formula1 <- paro ~ f(prov, model = "iid") +
  f(gedad, model = "iid") + f(nforma3, model = "iid")

tic()

m_inla <- inla(formula1,
  family = family1,
  data = epa,
  control.compute = list(dic = FALSE, cpo = FALSE),
  control.inla = list(int.strategy = "eb")
)

toc()
## 41.776 sec elapsed

Ojo, que a veces el modelo de INLA da un crash.

Resumen del modelo

summary(m_inla)
## 
## Call:
## c("inla(formula = formula1, family = family1, data = epa, control.compute = list(dic = FALSE, ",  "    cpo = FALSE), control.inla = list(int.strategy = \"eb\"))" )
## 
## Time used:
##  Pre-processing    Running inla Post-processing           Total 
##          1.2621         40.1902          0.3095         41.7618 
## 
## Fixed effects:
##                mean     sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept) -1.6466 0.3516    -2.3369  -1.6466    -0.9569 -1.6466   0
## 
## Random effects:
## Name   Model
##  prov   IID model 
## gedad   IID model 
## nforma3   IID model 
## 
## Model hyperparameters:
##                       mean    sd      0.025quant 0.5quant 0.975quant
## Precision for prov     8.0160  1.7225  5.1678     7.8371  11.8852   
## Precision for gedad   11.2787  9.1258  2.1820     8.7681  35.2338   
## Precision for nforma3  5.2391  4.3753  0.9665     4.0213  16.7312   
##                       mode   
## Precision for prov     7.4915
## Precision for gedad    5.3011
## Precision for nforma3  2.3706
## 
## Expected number of effective parameters(std dev): 52.89(0.00)
## Number of equivalent replicates : 1409.83 
## 
## Marginal Likelihood:  -28830.08

Pues al menos con este conjunto de datos y este modelo me quedaría con glmer en vez de la implementación con julia o INLA

 
comments powered by Disqus