Skip to content
Incomplete sheet

This sheet is incomplete and could use some attention. Please submit code snippet suggestions as an issue or PR here.

nlme

Useful links:

  • https://quantdev.ssri.psu.edu/sites/qdev/files/ILD_Ch06_2017_MLMwithHeterogeneousVariance.html

Code 

library(nlme)

Create

Action Code Details
Random intercept model with subject as grouping factor 
lme(fixed = y ~ 1, random = ~ 1 | subject)

Specify

Further specification options for lme()

Action Code Details
Random intercept per group 
random = ~ 1 | subject
Random intercept per group 
random = list(subject = pdSymm(form = ~ 1))
Random intercept, with different between-subject SD per secondary grouping 
random = list(subject = pdDiag(form = ~ diagnosis))
Within-subject error per secondary grouping 
weights = varIdent(form = ~ 1 | diagnosis)
With AR-1 correlation structure 
correlation = corAR1(form = ~ 1 | subject)

Model fit

Action Code Details
Log-likelihood 
logLik(model)
AIC 
AIC(model)
BIC 
BIC(model)
ICC 
performance::icc(model)
Proportion of the total variance explained by the between-cluster variance. Higher is better. 0.1 is considered a meaningful threshold.
Marginal pseudo R-squared 
r.squaredGLMM(model)[1, 'R2m']
Conditional pseudo R-squared 
r.squaredGLMM(model)[1, 'R2c']
Plot residuals against predicted 
plot(model)
Plot random effects 
plot(ranef(model))
Q-Q plot 
car::qqPlot(residuals(model))

Extract info from a fitted model.

Action Code Details
Fixed effects coefficients 
fixef(model)
Fixed effects coefficient intervals 
intervals(model)
Variance and SD summary 
VarCorr(model)
Variance-covariance matrix 
vcov(model)
Within-group variance 
sigma(model)^2
Within-group SD 
sigma(model)
Within-subject SD for group $GROUP (when using weights=varIdent) 
sigma(model) * coef(
    model$modelStruct$varStruct,
    unconstrained = FALSE,
    allCoef = TRUE
)['$GROUP']
Between-group variance 
VarCorr(model)['Residual', 'Variance']
Between-group SD 
VarCorr(model)['Residual', 'StdDev']

Extract effect size

Action Code Details
Cohen's d for main effects 
EMAtools::lme.dscore(model, type = 'nlme')
Cohen's d for group difference 
emg = emmeans::emmeans(model, ~ diagnosis)
emmeans::eff_size(
    emg,
    sigma = sigma(model),
    edf = data.frame(pairs(emg))$df[1]
)
Use with care
Cohen's f for main effects 
effectsize::cohens_f(model)
d = 2 * f

Hypothesis testing

Action Code Details
Coefficient significance 
car::Anova(model)
Test for nested model improvement (likelihood ratio test) 
performance::test_lrt(model0, model1)
Test for nested model improvement (likelihood ratio test) 
anova(model0, model1)