• Introduction to linear mixed models

• The R package skewlmm

• Challenge 1: within subject serial correlation

• Challenge 2: outliers

• Challenge 3: skewness

• Challenge 4: censored/missing data

Goal: Evaluate the impact of different high-fat diets on mice, two kinds of supplementation were considered, in addition to a standard high-fat diet (control).

Available in the R package skewlmm.

The weight of 52 mice was evaluated at baseline (week 0), and at weeks 1,2,…,10 of treatment.

data(miceweight)
head(miceweight, 4)

##   treat mouseID week weight
## 1     C     101    0   33.8
## 2     C     101    1   37.2
## 3     C     101    2   34.8
## 4     C     101    3   38.3

Linear mixed-effects models (LMMs) are an important class of statistical models that can be used to analyze correlated data. They are particularly useful for handling:

• repeated measures data;

• clustered data;

• longitudinal studies.

Such data are encountered in a variety of fields including biostatistics, public health, psychometrics, educational measurement, sociology and geophysics.

LMMs were introduced by Laird & Ware (1982)Laird, N. M., & Ware, J. H. (1982). Random-effects models for longitudinal data. Biometrics, 963-974. and allow to incorporate both “fixed” effects and “random” effects, as well as (possibly correlated) errors.

In the model, \[ \mathbf{Y}_{i}= \mathbf{X}_i\boldsymbol{\beta} + \mathbf{Z}_i \mathbf{b}_i + \boldsymbol{\epsilon}_{i}, \]

we typically assume that the random effects are independent of the error terms and that

\[ \mathbf{b}_i \sim N(\mathbf{0}, \mathbf{D}), \] \[ \boldsymbol\epsilon_i \sim N(\mathbf{0}, \boldsymbol{\Omega}_i). \]

This assumption yields a covariance structure for the outcomes

\[ \text{Var}(\mathbf{Y}_i) = \boldsymbol{\Sigma}_i= \mathbf{Z}_i \mathbf{DZ}_i^\top + \boldsymbol{\Omega}_i. \]

The expected model for subject \(i\) is

\[ E(\mathbf{Y}_i|\mathbf{X}_i, \mathbf{b}_i) = \mathbf{X}_i\boldsymbol\beta + \mathbf{Z}_i\mathbf{b}_i. \]

This interpretation of the mean given the subject-specific effect is called “conditional” interpretation.

The average over the random effects leads to a population-averaged “marginal” interpretation:

\[ E(E(\mathbf{Y}_i|\mathbf{X}_i,\mathbf{b}_i)) = \mathbf{X}_i\boldsymbol\beta. \]

So the fixed effects (\(\boldsymbol\beta\) coefficients) are also the estimated values in the population.

Our model is

\[ \mathbf{Y}_{i}= \mathbf{X}_i\boldsymbol{\beta} + \mathbf{Z}_i \mathbf{b}_i + \boldsymbol{\epsilon}_{i}. \]

If we fit a model with only random intercepts, we have

\[ \mathbf{Y}_{i}= \mathbf{X}_i\boldsymbol\beta + b_{0i} + \boldsymbol\epsilon_{i}. \]

We assume that \(b_{0i} \sim N(0, \sigma^2_b)\) and \(\boldsymbol\epsilon_i \sim N(\mathbf{0}, \sigma^2_e\mathbf{I})\), where

• \(\sigma^2_e\) is the amount of random variability (of the residuals),

• \(\sigma^2_b\) is the amount of variability in the population around the population average line.

Two main popular functions to fit LMMs in R are:

• lme() in the nlme package: supports several random effects and error level dependence structures;

• lmer() in the lme4 package: has more efficient linear algebra tools and is more efficient for fitting models with crossed random effects, but does not support special dependence structures.

Fitting the model \[ \begin{aligned} \mathbf{y}_{i} = & (b_{0i} +\beta_0+ \beta_1\,\text{treat}_{1i}+ \beta_2\,\text{treat}_{2i})\mathbf{1}_{n_i} \\ & + (b_{1i}+\beta_3+\beta_4\,\text{treat}_{1i}+\beta_5\,\text{treat}_{2i}) {\bf t}_i +\mathbb{\epsilon}_{i} \end{aligned} \] where \(y_{i}\) denotes the weight in grams (g) and \(t_{ij}\) denotes the numbers of weeks of diet, for \(t_{i1}=5, t_{i2}=6, \ldots, t_{i6}=10\), \(i = 1, \ldots, 52\). Also, \(t_{ij}\) is (week - 7.5).

Obs.: Time variable was centered in 0 and the first four weeks were excluded as no treatment was received during that time.

randomEffects = ranef(fnlme)

The skewlmm packageSchumacher, F. L., Matos, L. A., Valeriano, K. A. L., & Lachos, V. H. (2023). skewlmm: Scale Mixture of Skew-Normal Linear Mixed Models. R package version 1.1.0. contains flexible and robust models for longitudinal data. It can be installed from GitHub or its released version from CRAN as follows:

remotes::install_github("fernandalschumacher/skewlmm")
#or
install.packages("skewlmm")

Its three main functions are smn.lmm, smsn.lmm, and smn.clmm:

library(skewlmm)
smn.lmm(data, formFixed, groupVar, ...)
smsn.lmm(data, formFixed, groupVar, ...)
smn.clmm(data, formFixed, groupVar, ci, lcl, ucl, ...)

fskewlmm <- smn.lmm(formFixed = weight ~ treat*weekt, formRandom = ~ weekt, 
                   groupVar = "mouseID", data = miceweight)
summary(fskewlmm)

## Linear mixed models with distribution norm and dependency structure UNC 
## Call:
## smn.lmm(data = miceweight, formFixed = weight ~ treat * weekt, 
##     groupVar = "mouseID", formRandom = ~weekt)
## 
## Distribution norm
## Random effects: 
##   Formula: ~weekt
##   Structure:  
##   Estimated variance (D):
##             (Intercept)     weekt
## (Intercept)   22.952931 1.3137366
## weekt          1.313737 0.8853784
## 
## Fixed effects: weight ~ treat * weekt
## with approximate confidence intervals
##                    Value Std.error CI 95% lower CI 95% upper
## (Intercept)   43.9939394 1.6480925    40.763737   47.2241414
## treatT1        0.9155844 1.9111436    -2.830188    4.6613571
## treatT2       -0.9389394 1.9806544    -4.820951    2.9430719
## weekt          1.7589610 0.3746267     1.024706    2.4932159
## treatT1:weekt -0.1703896 0.4492493    -1.050902    0.7101229
## treatT2:weekt -0.5881039 0.4344989    -1.439706    0.2634984
## 
## Dependency structure: UNC
##   Estimate(s):
##   sigma2 
## 7.579361 
## 
## Model selection criteria:
##    logLik      AIC      BIC
##  -862.978 1745.957 1783.387
## 
## Number of observations: 312 
## Number of groups: 52

fskewlmm

## Linear mixed models with distribution norm and dependency structure UNC 
## Call:
## smn.lmm(data = miceweight, formFixed = weight ~ treat * weekt, 
##     groupVar = "mouseID", formRandom = ~weekt)
## 
## Fixed: weight ~ treat * weekt
## Random:
##   Formula: ~weekt
##   Structure: General positive-definite 
##   Estimated variance (D):
##             (Intercept)     weekt
## (Intercept)   22.952931 1.3137366
## weekt          1.313737 0.8853784
## 
## Estimated parameters:
##      (Intercept) treatT1 treatT2  weekt treatT1:weekt treatT2:weekt sigma2
##          43.9939  0.9156 -0.9389 1.7590       -0.1704       -0.5881 7.5794
## s.e.      1.6481  1.9111  1.9807 0.3746        0.4492        0.4345 0.5441
##      Dsqrt11 Dsqrt12 Dsqrt22
##       4.7854  0.2306  0.9123
## s.e.  0.7351  0.1937  0.2028
## 
## Model selection criteria:
##    logLik      AIC      BIC
##  -862.978 1745.957 1783.387
## 
## Number of observations: 312 
## Number of groups: 52

acfresid(fskewlmm, maxLag = 4) %>% kable(digits = 4)

acfresid(fskewlmm, maxLag = 4, calcCI = TRUE) %>% plot()

The following methods/functions are available:

• The misspecification of the distributional assumption may result in invalid statistical inferences, especially when dealing with heavy tails and skewness

• For instance, substantial bias in the ML estimates of regression parameters can result when the random-effects distribution is misspecified (Drikvandi et al., 2017)Drikvandi, R., Verbeke, G., & Molenberghs, G. (2017). Diagnosing misspecification of the random-effects distribution in mixed models. Biometrics, 73(1), 63-71.

• In longitudinal studies, repeated measures are collected over time and hence the error term can be serially correlated (not necessarily handled by the random effects specification)

The following dependence structures can be helpful to model serial correlation:

1. Uncorrelated (UNC): \(\boldsymbol{\Omega}_i = \sigma^2_e \mathbf{I}_{n_i}\).

2. Autoregressive dependence of order \(p\) (AR(\(p\))): \[\boldsymbol{\Omega}_i = \sigma^2_e \mathbf{R}_i = \frac{\sigma^2_e}{1 - \phi_1 \rho_1 - \ldots - \phi_p \rho_p} [\rho_{|r-s|}],\] where \(\rho_1, \ldots, \rho_p\) are the theoretical autocorrelations of the process and functions of \(\boldsymbol{\phi} = (\phi_1, \ldots, \phi_p)^\top\), and they satisfy the Yule-Walker equations.

3. Damped exponential correlation (DEC): \[\boldsymbol{\Omega}_i = \sigma_e^2 \mathbf{R}_i = \sigma_e^2 \left[ \phi_1^{|t_{ij} - t_{ik}|^{\phi_2}} \right], \quad 0 < \phi_1 < 1, \phi_2 > 0.\]

We can fit an LMM with various correlation structures with the nlme package by specifying the correlation argument:

# Autoregressive dependence of order 1
fnlmeAR1 = update(fnlme, correlation = corAR1(), method = "ML",
                  control = lmeControl(maxIter = 200, msMaxIter = 200, 
                                       msMaxEval = 100, opt='optim'))

And with the skewlmm package using the depStruct argument:

# Autoregressive dependence of order 1
fskewlmmAR1 = update(fskewlmm, depStruct = "ARp")

We can perform likelihood ratio tests to compare the fitted objects as follows:

lr.test(fskewlmmAR1, fskewlmm)

## 
## Model selection criteria:
##               logLik      AIC      BIC
## fskewlmmAR1 -850.458 1722.916 1764.089
## fskewlmm    -862.978 1745.957 1783.387
## 
##     Likelihood-ratio Test
## 
## chi-square statistics =  25.04094 
## df =  1 
## p-value =  5.612572e-07 
## 
## The null hypothesis that both models represent the 
## data equally well is rejected at level  0.05

acfresid(fskewlmmAR1, maxLag = 4, calcCI = TRUE) %>% plot()

A possible approach to reduce the effects of outliers when analyzing longitudinal data is to replace usual normality assumptions by heavier-than-normal tailed distributions.

In a previous workSchumacher, F. L., Lachos, V. H., & Matos, L. A. (2021). Scale mixture of skew‐normal linear mixed models with within‐subject serial dependence. Statistics in medicine, 40(7)., we proposed the use the Scale Mixture of Normal (SMN) distributions:

\[ \left( \begin{array}{c} \textbf{b}_i \\ \boldsymbol{\epsilon}_i \end{array}\right) \stackrel {{\rm ind}}{\sim} \textrm{SMN}_{q+n_i} \left( \left(\begin{array}{c} \mathbf{0}\\ \mathbf{0} \end{array} \right), \left(\begin{array}{cc} \mathbf{D} & \mathbf{0} \\ \mathbf{0} & \boldsymbol{\Sigma}_i \end{array} \right); H \right), \]

where

• the tail behavior depends on the distribution function specified by \(H\), including the distributions normal, Student-t, Slash, and contaminated-normal,

• \(\boldsymbol{\Sigma}_i = \sigma_e^2 \mathbf{R}_i\), with \(\mathbf{R}_i = \mathbf{R}_i(\boldsymbol{\phi})\), \(\boldsymbol{\phi}=(\phi_1,\ldots,\phi_p)^\top\), is a structured matrix.

Updating the previous fitted model to consider the \(t\) and CN distributions to miceweight data:

ftAR1 <- update(fskewlmmAR1, distr = "t")
fcnAR1 <- update(fskewlmmAR1, distr = "cn")

criteria(list(normal = fskewlmmAR1, t = ftAR1, cn = fcnAR1)) %>% kable() 

summary(fcnAR1)

## Linear mixed models with distribution cn and dependency structure ARp 
## Call:
## smn.lmm(data = miceweight, formFixed = weight ~ treat * weekt, 
##     groupVar = "mouseID", formRandom = ~weekt, depStruct = "ARp", 
##     distr = "cn")
## 
## Distribution cn with nu = 0.4349053 0.1426085 
## 
## Random effects: 
##   Formula: ~weekt
##   Structure:  
##   Estimated variance (D):
##             (Intercept)     weekt
## (Intercept)   5.2180664 0.7633666
## weekt         0.7633666 0.1430489
## 
## Fixed effects: weight ~ treat * weekt
## with approximate confidence intervals
##                     Value Std.error CI 95% lower CI 95% upper
## (Intercept)   42.60777956 1.2565464   40.1449938   45.0705653
## treatT1        0.60421060 1.4760881   -2.2888688    3.4972900
## treatT2       -0.38441693 1.5243888   -3.3721641    2.6033302
## weekt          1.15530811 0.4107661    0.3502213    1.9603949
## treatT1:weekt  0.46052571 0.4764984   -0.4733940    1.3944455
## treatT2:weekt -0.04879976 0.4680949   -0.9662488    0.8686493
## 
## Dependency structure: ARp
##   Estimate(s):
##    sigma2      phi1 
## 2.8646272 0.7740477 
## 
## Model selection criteria:
##    logLik      AIC      BIC
##  -815.797 1657.594 1706.253
## 
## Number of observations: 312 
## Number of groups: 52

In applications, skewness is often handled by log-transformations. However, interpreting model estimates with respect to log-transformed outcomes is challenging and usually non-intuitive.

Instead, we propose to consider skewed distributions in such situations, which can naturally accommodate such feature.

In this context, the flexible family of distributions Scale Mixture of Skew-Normal can be considered

\[ \left( \begin{array}{c} \textbf{b}_i \\ \boldsymbol{\epsilon}_i \end{array}\right) \stackrel {{\rm ind}}{\sim} \textrm{SMSN}_{q+n_i} \left( \left(\begin{array}{c} c\boldsymbol{\Delta} \\ \mathbf{0} \end{array} \right), \left(\begin{array}{cc} \mathbf{D} & \mathbf{0} \\ \mathbf{0} & \boldsymbol{\Sigma}_i \end{array} \right), \left(\begin{array}{c} \boldsymbol{\lambda} \\ \mathbf{0} \end{array} \right); H \right), \]

• the chosen location parameter for \(\mathbf{b}_i\) ensures that \(E\{\mathbf{b}_i\} = E\{\boldsymbol{\epsilon}_i\} = \mathbf{0}\), so \(E\{\mathbf{Y}_i\} =\textbf{X}_i\boldsymbol{\beta}\),

• \(\boldsymbol{\lambda}\) regulates the skewness,

• when \(\boldsymbol{\lambda}=\mathbf{0}\), the SMN is obtained as a special case.

We can fit a skewed model to the miceweight data with AR(1) correlation using the following syntax:

fscnAR1 <- smsn.lmm(formFixed = weight ~ treat*weekt, formRandom = ~ weekt, 
                   groupVar = "mouseID", data = miceweight, depStruct = "ARp",
                   distr = "scn")

criteria(list(normal = fskewlmmAR1, cn = fcnAR1, scn = fscnAR1)) %>% kable() 

lr.test(fscnAR1, fcnAR1)

## 
## Model selection criteria:
##           logLik      AIC      BIC
## fscnAR1 -809.443 1648.887 1705.032
## fcnAR1  -815.797 1657.594 1706.253
## 
##     Likelihood-ratio Test
## 
## chi-square statistics =  12.70716 
## df =  2 
## p-value =  0.001740502 
## 
## The null hypothesis that both models represent the 
## data equally well is rejected at level  0.05

summary(fscnAR1)

## Linear mixed models with distribution scn and dependency structure ARp 
## Call:
## smsn.lmm(data = miceweight, formFixed = weight ~ treat * weekt, 
##     groupVar = "mouseID", formRandom = ~weekt, depStruct = "ARp", 
##     distr = "scn")
## 
## Distribution scn with nu = 0.4749987 0.1378151 
## 
## Random effects: 
##   Formula: ~weekt
##   Structure:  
##   Estimated variance (D):
##             (Intercept)     weekt
## (Intercept)    12.49990 1.5979104
## weekt           1.59791 0.2805946
## 
## Fixed effects: weight ~ treat * weekt
## with approximate confidence intervals
##                     Value Std.error CI 95% lower CI 95% upper
## (Intercept)   44.21203515 1.6914356   40.8968824   47.5271879
## treatT1       -0.26795376 1.3757840   -2.9644409    2.4285333
## treatT2       -0.48368047 1.3266640   -3.0838941    2.1165332
## weekt          1.38022226 0.4030417    0.5902751    2.1701694
## treatT1:weekt  0.36485313 0.4302992   -0.4785177    1.2082240
## treatT2:weekt -0.09452354 0.4091909   -0.8965230    0.7074759
## 
## Dependency structure: ARp
##   Estimate(s):
##    sigma2      phi1 
## 2.4720220 0.6774895 
## 
## Skewness parameter estimate: 14.55327 4.307805
## 
## Model selection criteria:
##    logLik      AIC      BIC
##  -809.443 1648.887 1705.032
## 
## Number of observations: 312 
## Number of groups: 52

We can assess the goodness of fit using a Healy-type plot:

grid.arrange(healy.plot(fskewlmmAR1, calcCI = TRUE) + geom_point(size = 1.5),
             healy.plot(fscnAR1, calcCI = TRUE) + geom_point(size = 1.5), nrow = 1)

We can check for residual correlation based on ACF:

grid.arrange(plot(acfresid(fskewlmmAR1, calcCI = TRUE, maxLag = 5)),
             plot(acfresid(fscnAR1, calcCI = TRUE, maxLag = 5)), nrow = 1)

We can evaluate atypical observations using the Mahalanobis distance:

grid.arrange(plot(mahalDist(fskewlmmAR1), nlabels = 2),
             plot(mahalDist(fscnAR1), nlabels = 2), nrow = 1)

In a recent workCastro, K.T., Matos, L.A., Schumacher, F.L. (2024+) Diagnostic analysis in scale mixture of skew-normal linear mixed models. Under review., we explored diagnostic analysis in SMSN-LMM by deriving:

• Case-deletion measures;

• Local influence (under several perturbations).

Such tools will be incorporated in the skewlmm package soon.

Censoring can occur when measurement techniques are subject to limits of detection, above or below which the exact measure is not detectable .

Missing data can occur in longitudinal studies for different reasons. E.g.:

• Loss to follow up;

• Dropout (deliberate discontinuation);

• Intermittent missing (missed appointments).

A missing observation can be considered as a censored observation for which the censoring interval is the full real line (\(-\infty,\infty\)).

Assume that some mice (more affected by the diet) dropped out on the study after the 8th week of diet:

• Under a missing at random mechanism, a maximum likelihood-based estimates (e.g., LMM or GLMM) will usually yield unbiased estimates (as opposed to marginal models, such as GEE).

• For improve the performance, a model-based approach based on Expectation-Maximization (EM) algorithm can be used for likelihood-based estimation, where missing values are treated as latent variable

• Full model specification can be found in Matos et al.Matos, L. A., Prates, M. O., Chen, M. H., & Lachos, V. H. (2013). Likelihood-based inference for mixed-effects models with censored response using the multivariate-t distribution. Statistica Sinica, 1323-1345.

• This method is implemented in the context of censored data as part of the R package skewlmm.

smn.clmm(data, formFixed, groupVar, formRandom = ~1, depStruct = "UNC",
         ci, lcl, ucl, timeVar = NULL, distr = "norm",
         nufix = FALSE, pAR = 1, control = lmmControl())

fitCens <- smn.clmm(formFixed = weight ~ treat*weekt, formRandom = ~ weekt, 
                   groupVar = "mouseID", data = miceweight_miss, depStruct = "ARp",
                   ci = "miss", lcl = "lcl")

fitCensT <- update(fitCens, distr = "t")

criteria(list(normal = fitCens, t = fitCensT)) %>% kable() 

summary(fitCensT)$tableFixed %>% kable(digits=4)