Fit a (mixed-effects) regression model in Julia

Fit a (mixed-effects) regression model in Julia

Usage

jlm(formula, data, family = "gaussian", contrasts = jl_contrasts(data), ...)

jlmer(
  formula,
  data,
  family = NULL,
  contrasts = jl_contrasts(data),
  ...,
  fit = TRUE,
  optsum = list(),
  progress = interactive()
)

Arguments

formula

A formula written in Julia syntax. Can be a string or a language object.

data

A data frame

family

A distribution family

contrasts

A Julia dictionary of contrasts Inferred from data by default.

...

Additional arguments to the fit() function called in Julia

fit

Whether to fit the model. If FALSE, returns the unfit model object.

optsum

A list of values to set for the optimizer. See $optsum of unfit model for possible options.

progress

Whether to print model fitting progress. Defaults to interactive()

Value

A julia model object of class jlme

Examples

# \donttest{
jlme_setup(restart = TRUE)
#> Starting Julia (v1.10.5) ...
#> Successfully set up Julia connection. (13s)

# Fixed effects models
lm(mpg ~ hp, mtcars)
#> 
#> Call:
#> lm(formula = mpg ~ hp, data = mtcars)
#> 
#> Coefficients:
#> (Intercept)           hp  
#>    30.09886     -0.06823  
#> 
jlm(mpg ~ hp, mtcars)
#> <Julia object of type StatsModels.TableRegressionModel>
#> 
#> mpg ~ 1 + hp
#> 
#> ────────────────────────────────────────────────────────────────────────────
#>                   Coef.  Std. Error      z  Pr(>|z|)   Lower 95%   Upper 95%
#> ────────────────────────────────────────────────────────────────────────────
#> (Intercept)  30.0989      1.63392    18.42    <1e-75  26.8964     33.3013
#> hp           -0.0682283   0.0101193  -6.74    <1e-10  -0.0880617  -0.0483948
#> ──────────────────────────────────────────────────────────────────────────── 

# Auto-handling of contrasts
x <- mtcars
x$cyl_helm <- factor(x$cyl)
contrasts(x$cyl_helm) <- contr.helmert(3)
colnames(contrasts(x$cyl_helm)) <- c("4vs6", "4&6vs8")
lm(mpg ~ cyl_helm, x)
#> 
#> Call:
#> lm(formula = mpg ~ cyl_helm, data = x)
#> 
#> Coefficients:
#>    (Intercept)    cyl_helm4vs6  cyl_helm4&6vs8  
#>         20.502          -3.460          -2.701  
#> 
jlm(mpg ~ cyl_helm, x)
#> <Julia object of type StatsModels.TableRegressionModel>
#> 
#> mpg ~ 1 + cyl_helm
#> 
#> ─────────────────────────────────────────────────────────────────────────────
#>                      Coef.  Std. Error      z  Pr(>|z|)  Lower 95%  Upper 95%
#> ─────────────────────────────────────────────────────────────────────────────
#> (Intercept)       20.5022     0.593528  34.54    <1e-99   19.3389    21.6655
#> cyl_helm: 4vs6    -3.46039    0.779174  -4.44    <1e-05   -4.98754   -1.93324
#> cyl_helm: 4&6vs8  -2.70108    0.387175  -6.98    <1e-11   -3.45993   -1.94223
#> ───────────────────────────────────────────────────────────────────────────── 

# Mixed effects models
library(lme4)
#> Loading required package: Matrix

glmer(r2 ~ Anger + Gender + (1 | id), VerbAgg, family = "binomial")
#> Generalized linear mixed model fit by maximum likelihood (Laplace
#>   Approximation) [glmerMod]
#>  Family: binomial  ( logit )
#> Formula: r2 ~ Anger + Gender + (1 | id)
#>    Data: VerbAgg
#>       AIC       BIC    logLik  deviance  df.resid 
#>  9504.505  9532.240 -4748.253  9496.505      7580 
#> Random effects:
#>  Groups Name        Std.Dev.
#>  id     (Intercept) 1.059   
#> Number of obs: 7584, groups:  id, 316
#> Fixed Effects:
#> (Intercept)        Anger      GenderM  
#>    -1.10090      0.04626      0.26002  
jlmer(r2 ~ Anger + Gender + (1 | id), VerbAgg, family = "binomial")
#> <Julia object of type GeneralizedLinearMixedModel>
#> 
#> Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 1)
#>   r2 ~ 1 + Anger + Gender + (1 | id)
#>   Distribution: Bernoulli{Float64}
#>   Link: LogitLink()
#> 
#>    logLik    deviance     AIC       AICc        BIC    
#>  -4748.2525  9496.5050  9504.5050  9504.5102  9532.2401
#> 
#> Variance components:
#>       Column   VarianceStd.Dev.
#> id (Intercept)  1.12074 1.05865
#> 
#>  Number of obs: 7584; levels of grouping factors: 316
#> 
#> Fixed-effects parameters:
#> ────────────────────────────────────────────────────
#>                   Coef.  Std. Error      z  Pr(>|z|)
#> ────────────────────────────────────────────────────
#> (Intercept)  -1.10115     0.280681   -3.92    <1e-04
#> Anger         0.0462741   0.0134906   3.43    0.0006
#> Gender: M     0.260057    0.153847    1.69    0.0910
#> ──────────────────────────────────────────────────── 

# Set optimizer options via `optsum`
jlmer(
  r2 ~ Anger + Gender + (1 | id), VerbAgg, family = "binomial",
  optsum = list(
    optimizer = jl(":LN_NELDERMEAD"),
    maxfeval = 10L
  )
)
#> <Julia object of type GeneralizedLinearMixedModel>
#> 
#> Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 1)
#>   r2 ~ 1 + Anger + Gender + (1 | id)
#>   Distribution: Bernoulli{Float64}
#>   Link: LogitLink()
#> 
#>    logLik    deviance     AIC       AICc        BIC    
#>  -4749.2651  9498.5303  9506.5303  9506.5356  9534.2655
#> 
#> Variance components:
#>       Column   VarianceStd.Dev.
#> id (Intercept)  1.0 1.0
#> 
#>  Number of obs: 7584; levels of grouping factors: 316
#> 
#> Fixed-effects parameters:
#> ────────────────────────────────────────────────────
#>                   Coef.  Std. Error      z  Pr(>|z|)
#> ────────────────────────────────────────────────────
#> (Intercept)  -0.851285     0.267259  -3.19    0.0014
#> Anger         0.0353915    0.012845   2.76    0.0059
#> Gender: M     0.204939     0.14652    1.40    0.1619
#> ──────────────────────────────────────────────────── 

stop_julia()
# }