Kui mängu tulevad juhuslikud faktorid on R-is põhitööriistaks paketi lme4 poolt pakutav (kuivõrd see pakett on praegu alles arenemisjärgus siis nn. mittestabiilne, ent rohkemate võimalustega pakett lme4a võib sageli võimaldada teha ka keerukamaid asju). Põhifunktsioon selles paketis kannab nime lmer (või selle üldistus glmer, millest tuleb juttu järgmises punktis). Siiski tasub ära mainida, et paketist nlme pärinev funktsioon lme (samuti segamudelite sobitamieks) omab mõningaid võimalusi, mida lmer-is veel pole (näiteks eelmises punktis kirjeldatud korrelatsioonistruktuuride määramine).

> library(lme4)
> m6=lmer(weight~Diet+(I(Time-10.5)|Chick),data=ChickWeight)

> m7=lmer(weight~Diet+(1|Chick)+(I(Time-10.5)-1|Chick),data=ChickWeight)

Funktsiooni lmer süntaks on muidu traditsioniline, ent juhuslikud faktorid paigutatakse sulgude sisse, püstkriipsust vasakul näidatakse tunnust, mille parameeterist juttu on ja paremal tunnust, mille alusel toimub jagamine. Seega siis mudel m6 on selline, kus iga kaal sõltub lisaks saadud toidutüübile veel ka ajast, kusjuures igal objektil on nii erinev vabaliige kui ka erinev sõltuvus ajast. Mudel m7 erineb eelnimetatust vaid seeläbi, et eeldatakse juhuslike efektide omavahelist sõltumatust (ehk nt konkreetse objekti suur vabaliikmele vastav juhuslik efekt ei anna alust arvata nagu ka ajast sõluv efekt peaks olema suur).

Järgnevad mudelite väljatrükid:

> summary(m6)
Linear mixed model fit by REML
Formula: weight ~ Diet + (I(Time - 10.5) | Chick)
Data: ChickWeight
AIC BIC logLik deviance REMLdev
4909 4944 -2447 4905 4893
Random effects:
Groups Name Variance Std.Dev. Corr
Chick (Intercept) 4326.675 65.7775
I(Time - 10.5) 84.727 9.2047 0.999
Residual 163.460 12.7851
Number of obs: 578, groups: Chick, 50

Fixed effects:
Estimate Std. Error t value
(Intercept) 54.923 1.390 39.52
Diet2 2.854 2.365 1.21
Diet3 1.996 2.365 0.84
Diet4 9.274 2.368 3.92

Correlation of Fixed Effects:
(Intr) Diet2 Diet3
Diet2 -0.588
Diet3 -0.588 0.345
Diet4 -0.587 0.345 0.345
> summary(m7)
Linear mixed model fit by REML
Formula: weight ~ Diet + (1 | Chick) + (I(Time - 10.5) - 1 | Chick)
Data: ChickWeight
AIC BIC logLik deviance REMLdev
5036 5067 -2511 5045 5022
Random effects:
Groups Name Variance Std.Dev.
Chick (Intercept) 592.130 24.3337
Chick I(Time - 10.5) 84.737 9.2053
Residual 163.486 12.7862
Number of obs: 578, groups: Chick, 50

Fixed effects:
Estimate Std. Error t value
(Intercept) 99.231 5.640 17.594
Diet2 19.811 9.612 2.061
Diet3 38.976 9.612 4.055
Diet4 32.748 9.614 3.406

Correlation of Fixed Effects:
(Intr) Diet2 Diet3
Diet2 -0.587
Diet3 -0.587 0.344
Diet4 -0.587 0.344 0.344

Juhuslike efektide prognoosid väljastab funktsioon ranef. Mudelisse oluliste argumentide valikul saab kasutada funktsiooni anova. Kuivõrd mittetasakaaluliste andmestike korral ei ole klassikalised lähendid testidena kasutatavad, siis leitakse paketis lme4 (juhuslikud) usaldusvahemikud parameetrite järeljaotusest simuleerides. Selleks kasutatakse kahe funktsiooni abi, ent esialgu on selline simuleerimine võimalik vaid mudeli korral, kus juhuslikud efektid ei ole korreleeritud ning ka siis tuleb juhuslike efektide  usaldusvahemikesse suhtuda mõningase ettevaatlikusega (paketis lme4a võimaldavad usaldusvahemikke leida funktsioonid profile ja confint).

> HPDinterval(mcmcsamp(m7,n=2000))
$fixef
                lower     upper
(Intercept) 94.050477 105.99396
Diet2        8.713968  28.66470
Diet3       28.276681  48.28054
Diet4       22.194475  42.24714
attr(,"Probability")
[1] 0.95

$ST
         lower     upper
[1,] 0.7168018 0.9168519
[2,] 0.4592102 0.6633417
attr(,"Probability")
[1] 0.95

$sigma
        lower    upper
[1,] 13.93005 16.02898
attr(,"Probability")
[1] 0.95