# Content

Using Julia for multilevel modelling

- Demo with stata
- Demo with R
- Demo with Julia

Using Julia for multilevel modelling

- Demo with stata
- Demo with R
- Demo with Julia

- ESS cumulative Dataset
- Non missing countries
- Non missing ESS round

- Linear model
- Predicting global attitudes towards immigration
- Age as a predictor
- Random effect for country
- Random effect for ESS Wave
- Random Intercept + Random Slope

scale : attitudes towards immigration (0-10, 33 values) agea : age (quantitative) pays / cntry : pays essround : vague dbenquC*te

```
mixed scale agea [w=w] || pays: agea || essround: agea
predict yhat , fitted
```

```
(sampling weights assumed)
Obtaining starting values by EM:
Performing gradient-based optimization:
Iteration 0: log pseudolikelihood = -371021.98
Iteration 1: log pseudolikelihood = -371021.98
Computing standard errors:
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
-----------------------------------------------------------
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
------------------------------------------------------------------------------
------------------------------------------------------------------------------
| 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
------------------------------------------------------------------------------
```

`library(lme4)`

```
## Loading required package: Matrix
## Loading required package: Rcpp
```

```
df <- read.csv("/home/me/df.csv")
fm1 <- lmer(scale ~ agea2 + ( 1 + agea2 | cntry ) + ( 1 + agea2 | essround), data=df, weights=w)
```

```
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
```

`summary(fm1)`

```
## 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
##
## 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
##
## 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
```

```
julia -p 4
using DataFrames
using MixedModels
df = readtable("/home/me/df.csv")
fm1 = fit(lmm(scale ~ 1 + agea + (1+ agea |pays) + (1 + agea |essround), df))
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
```

- Stata : 1 min 54 sec
- R : 10 min 35 sec
- Julia : 0 min 3.9 sec

In Stata you get the results by predicting random effects and fitted values (additional computation)

In R and Julia, one only needs to extract fixed and random effects : fixef() and ranef()

nlme (old R library) :

- developed in mid-1990bs by Pinheiro and Bates
- documented in P & B (2000, Springer) and several papers
- implemented in R and C using the .Call interface and Rbs C API.
- fits linear mixed-effects models and nonlinear mixed-effects models
- allows for additional variance and or correlation structures in the conditional distribution of the response, given the random effects.
- multiple formula model specification with PDMat, VarCorr, b& classes.
- allows for multiple nested grouping factors. There is a kludgy method for fitting models with crossed grouping factors but I wouldnbt push it too hard.
- uses S3 methods and class tags (IMO bS3 classesb donbt exist)
- optimization of parameter estimates via EM and Newton-Raphson algorithms, profiled w.r.t. residual variance only
- Internal representation uses relative precision factor and log-Cholesky formulation. Singular covariance matrices correspond to infinite parameter values.
- no longer actively developed. Maintained by R-Core primarily to adjust to changing requirements on R packages
- no explicit methods for fitting GLMMs.
- lme with iterative reweighting is used in the glmmPQL function in the MASS package.

- development began in late 90bs by Bates, Maechler, DebRoy, Sarkar and others. Current development is primarily by Bolker and Walker.
- documented in papers, vignettes. Bates started writing a book but didnbt complete the project - some draft chapters still online.
- implemented in R, C and C++ using facilities provided by Rcpp and Eigen
- fits linear mixed-effects models and GLMMs. Nominally has NLMM capabilities but the implementation is not solid
- single formula model specification. Grouping factors can be nested, partially or fully crossed.
- uses S4 classes and methods, S3 methods, reference classes, Matrix package.
- internal representation uses relative covariance factor and sparse Cholesky factorization. Optimization of profiled log-likelihood uses general nonlinear optimizers that allow for box constraints (in practice only require non-negativity of some of the parameters) for LMMs. GLMMs use PIRLS (penalized iteratively reweighted least squares) to determine the conditional modes of the random effects with Laplace or adaptive Gauss-Hermite approximation to the log-likelihood. Settling on a single robust and reliable optimizer has been difficult.
- all models use a sparse matrix representation to solve the penalized least squares problem
- actively maintained and developed

- development started in 2012 by Bates
- little or no documentation outside of examples
- implemented exclusively in Julia (about 1600 lines of code)
- fits LMMs. Development of GLMM capabilities is planned.
- single formula specification similar to lme4. Because Julia is still a young language not all lme4 formulas will work with MixedModels. (Building the formula language interpretation code from scratch is not easy.) Grouping factors can be nested or crossed (partially or completely).
- uses Julia methods and user-defined types. Julia methods provide for multiple dispatch like S4.
- internal representation (all at the Julia level) provides for several different PLS (penalized least squares) solvers according to the configuration of the grouping factors.
- Models with a nested sequence of grouping factors, including a single grouping factor as a trivial special case, use dense matrix methods and provide an analytic gradient of the profiled log-likelihood. Similarly for models with two crossed or nearly crossed grouping factors (think bsubjectb and bitemb).
- Optimizers from Steven Johnsonbs NLopt Julia package that interfaces to his C library, which is also used in the R packages nloptr and nloptwrap. At present LN_BOBYQA is used for derivative-free optimization and LD_MMA for models that provide a gradient but switching optimizers within the NLopt library is trivial.
- Considerable effort made to reuse storage and cut down on allocation/garbage collection.
- Actively developed. Maintenance hasnbt been too much of an issue because I donbt think there are many users at present.

bBecause I was involved in the development of all three packages I will take the liberty of commenting on the philosophy.b"

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.

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 (my fault entirely) and fragile.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 doesnbt have the infrastructure of R (dataframes, formula language, graphics, b&) but does have a clean implementation of the parts that are available. The most important aspect of Julia is bone languageb. You develop in the same language in which you optimize.

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

Julia remains a work in progress It may not still be production ready But it allows the estimation of complex models a lot easier than alternatives

Solid arguments in favor of Julia

- Itbs free
- Itbs growing
- Important R package developers are working with Julia
- The syntax is inspired by R and rather transparent
- R an Julia can interoperate quite easily