Examples

Oats random effect model

This example uses data of an agricultural trial where different varieties of oats were grown with different amounts of nitrogen fertilizer, and the yield was measured for each plot. The pots were grouped into blocks that are treated as random effects because they are geographically contiguous and are separated from each other. So we expect plots within blocks to be more alike than those in different blocks.

First, import the data:

library(MASS)
data(oats)
oats$nitrogen = substr(oats$N, 1, 3) |> as.numeric()

lme4 model:

library(lme4)

Loading required package: Matrix

lme_oats = lmer(Y ~ nitrogen + V + (1|B), data=oats)
summary(lme_oats)

Linear mixed model fit by REML ['lmerMod']
Formula: Y ~ nitrogen + V + (1 | B)
   Data: oats

REML criterion at convergence: 588

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.84069 -0.80849  0.04022  0.70484  2.22148 

Random effects:
 Groups   Name        Variance Std.Dev.
 B        (Intercept) 245.0    15.65   
 Residual             234.7    15.32   
Number of obs: 72, groups:  B, 6

Fixed effects:
            Estimate Std. Error t value
(Intercept)   82.400      7.516  10.964
nitrogen      73.667      8.075   9.123
VMarvellous    5.292      4.423   1.196
VVictory      -6.875      4.423  -1.554

Correlation of Fixed Effects:
            (Intr) nitrgn VMrvll
nitrogen    -0.322              
VMarvellous -0.294  0.000       
VVictory    -0.294  0.000  0.500

INLA model:

library(INLA)

Loading required package: sp

This is INLA_23.09.09 built 2023-10-16 17:35:11 UTC.
 - See www.r-inla.org/contact-us for how to get help.

inla_oats = inla(Y ~ nitrogen + V + f(B, model='iid'), data=oats)
summary(inla_oats)

Call:
   c("inla.core(formula = formula, family = family, contrasts = contrasts, 
   ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " 
   scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, 
   ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = 
   verbose, ", " lincomb = lincomb, selection = selection, control.compute 
   = control.compute, ", " control.predictor = control.predictor, 
   control.family = control.family, ", " control.inla = control.inla, 
   control.fixed = control.fixed, ", " control.mode = control.mode, 
   control.expert = control.expert, ", " control.hazard = control.hazard, 
   control.lincomb = control.lincomb, ", " control.update = 
   control.update, control.lp.scale = control.lp.scale, ", " 
   control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, 
   ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = 
   num.threads, ", " keep = keep, working.directory = working.directory, 
   silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug = 
   debug, .parent.frame = .parent.frame)" ) 
Time used:
    Pre = 1.68, Running = 0.984, Post = 0.0312, Total = 2.7 
Fixed effects:
              mean     sd 0.025quant 0.5quant 0.975quant   mode kld
(Intercept) 84.805  5.288     74.485   84.779     95.269 84.780   0
nitrogen    65.559 10.556     44.584   65.640     86.068 65.638   0
VMarvellous  5.222  5.961     -6.495    5.223     16.935  5.223   0
VVictory    -6.722  5.961    -18.432   -6.724      4.998 -6.724   0

Random effects:
  Name    Model
    B IID model

Model hyperparameters:
                                            mean      sd 0.025quant 0.5quant
Precision for the Gaussian observations    0.002    0.00      0.002    0.002
Precision for B                         7174.354 1978.87   3633.184 7053.491
                                        0.975quant     mode
Precision for the Gaussian observations   3.00e-03    0.002
Precision for B                           1.12e+04 6931.326

Marginal log-Likelihood:  -343.99 
 is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')