Programme

Un exemple d’utilisation de Julia pour la modélisation multiniveaux

  1. Demo avec stata
  2. Demo avec R
  3. Demo avec Julia

Données

Modèle

Variables

scale : attitudes vis à vis de l’immigration (0-10, 33 values) agea : age (quantitatif) pays / cntry : pays essround : vague d’enquête

Avec Stata

  1. mixed scale agea [w=w] || pays: agea || essround: agea

  2. predict yhat , fitted

Resultats (1)

Mixed-effects regression
Number of obs = 170359

-----------------------------------------------------------
                |   No. of       Observations per Group
 Group Variable |   Groups  Minimum   Average    Maximum
----------------+------------------------------------------
           pays |   16       7590     10647.4      16201
       essround |   96       1149     1774.6       2858
-----------------------------------------------------------
                                                           
                                                           

Résultats (2)

Wald chi2(1)       =     56.69
Log pseudolikelihood = -371021.98               
Prob > chi2        =    0.0000
(Std. Err. adjusted for 16 clusters in pays)
------------------------------------------------------------------------------
             |               Robust
       scale |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        agea |   -.011606   .0015415    -7.53   0.000    -.0146273   -.0085848
       _cons |   5.823571   .1225367    47.53   0.000     5.583403    6.063738
------------------------------------------------------------------------------

Résultats (3)

  ------------------------------------------------------------------------------
                               |               Robust           
  Random-effects Parameters    |   Estimate   Std. Err.     [95% Conf. Interval]
  -----------------------------+------------------------------------------------
  pays: Independent            |
                     var(agea) |   .0000309   .0000109      .0000154    .0000618
                    var(_cons) |   .2100469   .0697941      .1095167    .4028584
  -----------------------------+------------------------------------------------
  essround: Independent        |
                     var(agea) |   .0000104   2.37e-06      6.70e-06    .0000163
                    var(_cons) |    .045416   .0118672      .0272139    .0757927
  -----------------------------+------------------------------------------------
                 var(Residual) |    3.96908   .1949909      3.604727     4.37026
  ------------------------------------------------------------------------------

Avec R + lme4

  1. library(lme4)

  2. df <- read.csv(“/home/me/df.csv”)

  3. fm1 <- lmer(scale ~ agea2 + ( 1 + agea2 | cntry ) + ( 1 + agea2 | essround), data=df, weights=w)

  4. summary(fm1)

Résultats (1)

Linear mixed model fit by REML ['lmerMod']
Formula: scale ~ agea2 + (1 + agea2 | cntry) + (1 + agea2 | essround)
Data: df
Weights: w
REML criterion at convergence: 823608.8

Scaled residuals: 
   Min      1Q  Median      3Q     Max 
-8.0905 -0.3942  0.0411  0.4416 10.3865 

Résultats (2)

Random effects:
Groups   Name        Variance  Std.Dev.   Corr 
cntry    (Intercept) 2.285e-01 0.4780656      
         agea2       3.572e-05 0.0059769 -0.02  
essround (Intercept) 1.463e-02 0.1209345      
         agea2       6.481e-07 0.0008051  0.16 
Residual             4.131e+00 2.0324689      

Number of obs: 170359, groups:  cntry, 16; essround, 6

Résultats (3)

Fixed effects:
            Estimate Std. Error t value
(Intercept)  5.825191   0.131082   44.44
agea2       -0.011668   0.001591   -7.33
 
Correlation of Fixed Effects:
     (Intr)
     agea2 -0.046

Avec Julia + MixedModels

  1. julia -p 4

  2. using DataFrames
  3. using MixedModels

  4. df = readtable(“/home/me/df.csv”)

  5. fm1 = fit(lmm(scale ~ 1 + agea + (1+ agea |pays) + (1 + agea |essround), df))

Résultats

Linear mixed model fit by maximum likelihood
Formula: scale ~ 1 + agea + ((1 + agea) | pays) + ((1 + agea) | essround)

logLik: -354086.877907, deviance: 708173.755814
Variance components:

               Variance        Std.Dev.     Corr.
pays      0.228742671708289 0.47827050056
          0.000023400329230 0.00483738868  -0.02
essround  0.009791089252406 0.09894993306
          0.000000027877676 0.00016696609   1.00
Residual  3.736088501904096 1.93289640227

Number of obs: 170359; levels of grouping factors: 16, 6
Fixed-effects parameters:
-------------------------------------------
               Estimate  Std.Error  z value
-------------------------------------------
(Intercept)     5.84238   0.126896  46.0406
agea         -0.0116274 0.00123927 -9.38241

Temps de calcul

Différences

Avec Stata, il faut utiliser la commande predict pour obtenir les effets fixes et aléatoires ou encore les valeurs prédites. Cela implique donc d’autres calculs mais les résultats sont immédiatements intégrés au dataframe d’origine.

Avec R et Julia, ces valeurs sont stockées automatiquement dans un objet contenant les résultats. On utilise les extracteurs fixef(), ranef(), fitted() pour les obtenir. Les résultats sont stockés dans un nouvel objet distinct du dataframe d’origine

Revue des principales différences techniques (Douglas Bates)

nlme (old R library) :

LME 4 (recent R library)

MixedModels (Julia)

General comments by Bates :

“Because I was involved in the development of all three packages I will take the liberty of commenting on the philosophy.”

  1. “The purpose of developing nlme was for fitting nonlinear mixed-effects models. Linear mixed-effects models were incorporated mainly as an iterative step in nlme. Numerical methods used are not nearly as sophisticated as those in lme4 and MixedModels.”

  2. “lme4 was developed to provide a use-case for S4 classes and methods. The numerical methods implemented in lme4 are, in my opinion, superior to those in nlme, mainly through the use of the relative covariance factor and the profiled log-likelihood.
    Using C++, Rcpp and RcppEigen was motivated by trying to provide generality and speed. The end result is confusing […] and fragile.”

  3. “MixedModels was developed because I became interested in the Julia language at the same time that I was disillusioned with at least some aspects of R. As a new language Julia doesn’t have the infrastructure of R (dataframes, formula language, graphics, …) but does have a clean implementation of the parts that are available. The most important aspect of Julia is”one language“. You develop in the same language in which you optimize.”

https://stat.ethz.ch/pipermail/r-sig-mixed-models/2014q4/022791.html

Résumé

  1. Julia reste un langage en cours de développement rapide
  2. Il n’est pas pret pour une utilisation quotidienne sur certains aspects
  3. Mais pour les calculs lourds et complexes, il ouvre des possibilités sans précédent

Quelques arguments en faveur de Julia

  1. C’est un langage/logiciel libre (et utilisable pour le secteur privé)
  2. Le langage se développe rapidement
  3. De nombreux chercheurs travaillent avec R et passent à Julia
  4. La syntaxe est voisine de R et de python (relativement transparente)
  5. Julia permet de réutiliser des programmes et fonctions écrits en C, Fortran, R
  6. Julia est compatible avec les outils récents de modélisation avancée (Stan.jl)

Pour aller plus loin

  1. http://julialang.org/
  2. https://github.com/JuliaLang
  3. http://juliacon.org/
  4. http://juliastats.github.io/
  5. google groups [Julia Stats]
  6. Un groupe des utilisateurs Julia à Paris en préparation

Contact

antoine.jardin@sciencespo.fr