8.1 Dependent Data

In many real-world applications, observations are not independent but exhibit correlations due to shared characteristics. Below are common forms of dependent data:

• Multivariate measurements on individuals: Multiple attributes measured on the same person may be correlated (e.g., blood pressure, cholesterol, and BMI).
• Clustered measurements: Individuals within a shared environment (e.g., families, schools, hospitals) often exhibit correlated responses.
• Repeated measurements: When the same individual is measured multiple times, correlations arise naturally.
  • Example: Tracking cholesterol levels of a person over time.
  • If these repeated measurements follow an experimental design where treatments were applied initially, they are referred to as repeated measures (Schabenberger and Pierce 2001).
• Longitudinal data: Repeated measurements taken over time in an observational study are known as longitudinal data (Schabenberger and Pierce 2001).
• Spatial data: Measurements taken from individuals in nearby locations (e.g., residents of the same neighborhood) often exhibit spatial correlation.

Since standard linear regression assumes independent observations, these correlations violate its assumptions. Thus, Linear Mixed Models (LMMs) provide a framework to account for these dependencies.

A Linear Mixed Model (also called a Mixed Linear Model) consists of two components:

• Fixed effects: Parameters associated with variables that have the same effect across all observations (e.g., gender, age, diet, time).
• Random effects: Individual-specific variations or correlation structures (e.g., subject-specific effects, spatial correlations), leading to dependent (correlated) errors.

A key advantage of LMMs is that they model random subject-specific effects, rather than including individual dummy variables. This provides:

• A reduction in the number of parameters to estimate, avoiding overfitting.
• A framework for inference on a population level, rather than restricting conclusions to observed individuals.

8.1.1 Motivation: A Repeated Measurements Example

Consider a scenario where we analyze repeated measurements of a response variable $Y_{ij}$ for the $i$-th subject at time $j$:

• $i = 1, \dots, N$ (subjects)
• $j = 1, \dots, n_i$ (measurements per subject)

We define the response vector for subject $i$ as:

$$\mathbf{Y}_i = \begin{bmatrix} Y_{i1} \\ Y_{i2} \\ \vdots \\ Y_{in_i} \end{bmatrix}$$

To model this data, we use a two-stage hierarchical approach:

8.1.1.1 Stage 1: Regression Model (Within-Subject Variation) {#sec-stage-1-regression-model-(within-subject-variation}

We first model how the response changes over time for each subject:

$$\mathbf{Y}_i = Z_i \beta_i + \epsilon_i$$

where:

• $Z_i$ is an $n_i \times q$ matrix of known covariates (e.g., time, treatment).
• $\beta_i$ is a $q \times 1$ vector of subject-specific regression coefficients.
• $\epsilon_i$ represents random errors, often assumed to follow $\epsilon_i \sim N(0, \sigma^2 I)$.

At this stage, each individual has their own unique regression coefficients $\beta_i$. However, estimating a separate $\beta_i$ for each subject is impractical when $N$ is large. Thus, we introduce a [second stage](#sec-stage-2-parameter-model-(between-subject-variation) to impose structure on $\beta_i$.

8.1.1.2 Stage 2: Parameter Model (Between-Subject Variation) {#sec-stage-2-parameter-model-(between-subject-variation}

We assume that the subject-specific coefficients $\beta_i$ arise from a common population distribution:

$$\beta_i = K_i \beta + b_i$$

where:

• $K_i$ is a $q \times p$ matrix of known covariates.
• $\beta$ is a $p \times 1$ vector of global parameters (fixed effects).
• $\mathbf{b}_i$ are random effects, assumed to follow $\mathbf{b}_i \sim N(0, D)$.

Thus, each individual’s regression coefficients $\beta_i$ are modeled as deviations from a population-level mean $\beta$, with subject-specific variations $b_i$. This hierarchical structure enables:

• Borrowing of strength: Individual estimates $\beta_i$ are informed by the overall population distribution.
• Improved generalization: The model captures both fixed and random variability efficiently.

The full Linear Mixed Model can then be written as:

$$\mathbf{Y}_i = Z_i K_i \beta + Z_i b_i + \epsilon_i$$

where:

• $Z_i K_i \beta$ represents the fixed effects component.
• $Z_i b_i$ represents the random effects component.
• $\epsilon_i$ represents the residual errors.

This formulation accounts for both within-subject correlations (via random effects) and between-subject variability (via fixed effects), making it a powerful tool for analyzing dependent data.

---

8.1.2 Example: Linear Mixed Model for Repeated Measurements

8.1.2.1 Stage 1: Subject-Specific Model

The first stage models the response variable $Y_{ij}$ for subject $i$ at time $t_{ij}$:

$$Y_{ij} = \beta_{1i} + \beta_{2i} t_{ij} + \epsilon_{ij}$$

where:

• $j = 1, \dots, n_i$ represents different time points for subject $i$.
• $\beta_{1i}$ is the subject-specific intercept (baseline response for subject $i$).
• $\beta_{2i}$ is the subject-specific slope (rate of change over time).
• $\epsilon_{ij} \sim N(0, \sigma^2)$ are independent errors.

In matrix notation, the model is written as:

$$\mathbf{Y}_i = \begin{bmatrix} Y_{i1} \\ Y_{i2} \\ \vdots \\ Y_{in_i} \end{bmatrix}, \quad \mathbf{Z}_i = \begin{bmatrix} 1 & t_{i1} \\ 1 & t_{i2} \\ \vdots & \vdots \\ 1 & t_{in_i} \end{bmatrix},$$

$$\beta_i = \begin{bmatrix} \beta_{1i} \\ \beta_{2i} \end{bmatrix}, \quad \epsilon_i = \begin{bmatrix} \epsilon_{i1} \\ \epsilon_{i2} \\ \vdots \\ \epsilon_{in_i} \end{bmatrix}.$$

Thus, the model can be rewritten as:

$$\mathbf{Y_i = Z_i \beta_i + \epsilon_i}.$$

---

8.1.2.2 Stage 2: Population-Level Model

Since estimating a separate $\beta_{1i}$ and $\beta_{2i}$ for each subject is impractical, we assume they follow a population distribution:

$$\begin{aligned} \beta_{1i} &= \beta_0 + b_{1i}, \\ \beta_{2i} &= \beta_1 L_i + \beta_2 H_i + \beta_3 C_i + b_{2i}. \end{aligned}$$

where:

• $L_i, H_i, C_i$ are indicator variables for treatment groups:
  • $L_i = 1$ if the subject belongs to treatment group 1, else 0.
  • $H_i = 1$ if the subject belongs to treatment group 2, else 0.
  • $C_i = 1$ if the subject belongs to treatment group 3, else 0.
• $\beta_0$ represents the average baseline response across all subjects.
• $\beta_1, \beta_2, \beta_3$ are the average time effects for the respective treatment groups.
• $b_{1i}, b_{2i}$ are random effects representing subject-specific deviations.

This structure implies that while the intercept $\beta_{1i}$ varies randomly across subjects, the slope $\beta_{2i}$ depends both on treatment and random subject-specific deviations.

In matrix form, this is:

$$\begin{aligned} \mathbf{K}_i &= \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & L_i & H_i & C_i \end{bmatrix}, \\ \beta &= \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \beta_3 \end{bmatrix}, \\ \mathbf{b}_i &= \begin{bmatrix} b_{1i} \\ b_{2i} \end{bmatrix}, \\ \beta_i &= \mathbf{K_i \beta + b_i}. \end{aligned}$$

---

8.1.2.3 Final Mixed Model Formulation

Substituting Stage 2 into Stage 1, we obtain the full mixed model:

$$\mathbf{Y}_i = \mathbf{Z}_i (\mathbf{K}_i \beta + \mathbf{b}_i) + \epsilon_i.$$

Expanding:

$$\mathbf{Y}_i = \mathbf{Z}_i \mathbf{K}_i \beta + \mathbf{Z}_i \mathbf{b}_i + \epsilon_i.$$

Interpretation:

• $\mathbf{Z}_i \mathbf{K}_i \beta$ represents the fixed effects (population-level trends).
• $\mathbf{Z}_i \mathbf{b}_i$ represents the random effects (subject-specific variations).
• $\epsilon_i$ represents the residual errors.

Assumptions:

• Random effects: $\mathbf{b}_i \sim N(0, D)$, where $D$ is the variance-covariance matrix.
• Residual errors: $\epsilon_i \sim N(0, \sigma^2 I)$.
• Independence: $\mathbf{b}_i$ and $\epsilon_i$ are independent.

---

To estimate $\hat{\beta}$, one might consider a sequential approach:

1. Estimate $\hat{\beta}_i$ for each subject in the first stage.
2. Estimate $\hat{\beta}$ in the second stage by replacing $\beta_i$ with $\hat{\beta}_i$.

However, this method introduces several problems:

• Loss of information: Summarizing $\mathbf{Y}_i$ solely by $\hat{\beta}_i$ discards valuable within-subject variability.
• Ignoring uncertainty: Treating estimated values $\hat{\beta}_i$ as known values leads to underestimated variability.
• Unequal sample sizes: Subjects may have different numbers of observations ($n_i$), which affects variance estimation.

To address these issues, we adopt the Linear Mixed Model (LMM) framework (Laird and Ware 1982).

---

8.1.2.4 Reformulating the Linear Mixed Model

Substituting [Stage 2](#sec-stage-2-parameter-model-(between-subject-variation) into [Stage 1](#sec-stage-2-parameter-model-(between-subject-variation), we obtain:

$$\mathbf{Y}_i = \mathbf{Z}_i \mathbf{K}_i \beta + \mathbf{Z}_i \mathbf{b}_i + \mathbf{\epsilon}_i.$$

Defining $\mathbf{X}_i = \mathbf{Z}_i \mathbf{K}_i$ as an $n_i \times p$ matrix of covariates, we rewrite the model as:

$$\mathbf{Y}_i = \mathbf{X}_i \beta + \mathbf{Z}_i \mathbf{b}_i + \mathbf{\epsilon}_i.$$

where:

• $i = 1, \dots, N$ (subjects).

• $\beta$ are fixed effects, common to all subjects.

• $\mathbf{b}_i$ are subject-specific random effects, assumed to follow:

  $$\mathbf{b}_i \sim N_q(\mathbf{0,D}).$$

• $\mathbf{\epsilon}_i$ represents residual errors:

  $$\mathbf{\epsilon}_i \sim N_{n_i}(\mathbf{0,\Sigma_i}).$$

• Independence assumption: $\mathbf{b}_i$ and $\mathbf{\epsilon}_i$ are independent.

• Dimension notation:

  • $\mathbf{Z}_i$ is an $n_i \times q$ matrix of known covariates for random effects.
  • $\mathbf{X}_i$ is an $n_i \times p$ matrix of known covariates for fixed effects.

---

8.1.2.5 Hierarchical (Conditional) Formulation

Rewriting the LMM in hierarchical form:

$$\begin{aligned} \mathbf{Y}_i | \mathbf{b}_i &\sim N(\mathbf{X}_i \beta+ \mathbf{Z}_i \mathbf{b}_i, \mathbf{\Sigma}_i), \\ \mathbf{b}_i &\sim N(\mathbf{0,D}). \end{aligned}$$

where:

• The first equation states that, given the subject-specific random effects $\mathbf{b}_i$, the response follows a normal distribution with mean $\mathbf{X}_i \beta+ \mathbf{Z}_i \mathbf{b}_i$ and covariance $\mathbf{\Sigma}_i$.
• The second equation states that the random effects $\mathbf{b}_i$ follow a multivariate normal distribution with mean zero and covariance $D$.

We denote the respective probability density functions as:

$$f(\mathbf{Y}_i |\mathbf{b}_i) \quad \text{and} \quad f(\mathbf{b}_i).$$

---

Using the general marginalization formula:

$$\begin{aligned} f(A,B) &= f(A|B)f(B) \\ f(A) &= \int f(A,B)dB = \int f(A|B) f(B) dB \end{aligned}$$

we obtain the marginal density of $\mathbf{Y}_i$:

$$f(\mathbf{Y}_i) = \int f(\mathbf{Y}_i | \mathbf{b}_i) f(\mathbf{b}_i) d\mathbf{b}_i.$$

Solving this integral, we obtain:

$$\mathbf{Y}_i \sim N(\mathbf{X_i \beta}, \mathbf{Z_i D Z_i'} + \mathbf{\Sigma_i}).$$

Interpretation:

• Mean structure: The expectation remains $\mathbf{X}_i \beta$, the fixed effects.

• Variance structure: The covariance now incorporates random effect variability:

  $$\mathbf{Z_i D Z_i'} + \mathbf{\Sigma_i}.$$

  This shows that the random effects contribute additional correlation between observations.

🔹 Key Takeaway:
The marginal formulation of LMM no longer includes $Z_i b_i$ in the mean, but instead incorporates it into the covariance structure, introducing marginal dependence in $\mathbf{Y}_i$.

Technically, rather than “averaging out” the random effect $b_i$, we add its contribution to the variance-covariance matrix.

---

Continue with the example:

$$Y_{ij} = (\beta_0 + b_{1i}) + (\beta_1L_i + \beta_2 H_i + \beta_3 C_i + b_{2i})t_{ij} + \epsilon_{ij}.$$

For each treatment group:

$$Y_{ij} = \begin{cases} \beta_0 + b_{1i} + (\beta_1 + b_{2i}) t_{ij} + \epsilon_{ij}, & L_i = 1 \\ \beta_0 + b_{1i} + (\beta_2 + b_{2i}) t_{ij} + \epsilon_{ij}, & H_i = 1 \\ \beta_0 + b_{1i} + (\beta_3 + b_{2i}) t_{ij} + \epsilon_{ij}, & C_i = 1. \end{cases}$$

Interpretation:

• Intercepts and slopes are subject-specific: Each subject has their own baseline response and rate of change.
• Treatment groups affect slopes, but not intercepts:
  • $\beta_0$ is the common intercept across groups.
  • Slopes differ by treatment: $\beta_1$ for $L$, $\beta_2$ for $H$, $\beta_3$ for $C$.
• Random effects introduce within-subject correlation:
  • $b_{1i}$ allows individual variation in baseline response.
  • $b_{2i}$ allows subject-specific deviations in slopes.

---

In the hierarchical model form, we again express the LMM as:

$$\begin{aligned} \mathbf{Y}_i | \mathbf{b}_i &\sim N(\mathbf{X}_i \beta + \mathbf{Z}_i \mathbf{b}_i, \mathbf{\Sigma}_i)\\ \mathbf{b}_i &\sim N(\mathbf{0,D}) \end{aligned}$$

where:

• $\mathbf{X}_i$ and $\mathbf{Z}_i$ are design matrices for fixed and random effects, respectively.
• $\mathbf{b}_i$ represents subject-specific random effects.
• $\mathbf{\Sigma}_i$ is the residual error covariance matrix.
• $\mathbf{D}$ is the random effects covariance matrix.

The fixed-effects parameter vector is:

$$\beta = (\beta_0, \beta_1, \beta_2, \beta_3)'.$$

From the model structure, we define:

$$\mathbf{X}_i = \mathbf{Z}_i \mathbf{K}_i.$$

Expanding,

$$\mathbf{Z}_i = \begin{bmatrix} 1 & t_{i1} \\ 1 & t_{i2} \\ \vdots & \vdots \\ 1 & t_{in_i} \end{bmatrix}, \quad \mathbf{K}_i = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & L_i & H_i & C_i \end{bmatrix}.$$

Multiplying:

$$\mathbf{X}_i = \begin{bmatrix} 1 & t_{i1}L_i & t_{i1}H_i & t_{i1}C_i \\ 1 & t_{i2}L_i & t_{i2}H_i & t_{i2}C_i \\ \vdots & \vdots & \vdots & \vdots \\ 1 & t_{in_i}L_i & t_{in_i}H_i & t_{in_i}C_i \end{bmatrix}.$$

The random effects vector is:

$$\mathbf{b}_i = \begin{bmatrix} b_{1i} \\ b_{2i} \end{bmatrix}.$$

The covariance matrix $\mathbf{D}$ for random effects is:

$$\mathbf{D} = \begin{bmatrix} d_{11} & d_{12} \\ d_{12} & d_{22} \end{bmatrix}.$$

We assume conditional independence:

$$\mathbf{\Sigma}_i = \sigma^2 \mathbf{I}_{n_i}.$$

This means that, given the random effects $\mathbf{b}_i$ and $\beta$, the responses for subject $i$ are independent.

---

To derive the marginal model, we integrate out the random effects $\mathbf{b}_i$. This leads to:

$$Y_{ij} = \beta_0 + \beta_1 L_i t_{ij} + \beta_2 H_i t_{ij} + \beta_3 C_i t_{ij} + \eta_{ij},$$

where:

$$\eta_i \sim N(\mathbf{0}, \mathbf{Z}_i \mathbf{D} \mathbf{Z}_i' + \mathbf{\Sigma}_i).$$

Thus, the full marginal model is:

$$\mathbf{Y_i} \sim N(\mathbf{X}_i \beta, \mathbf{Z}_i \mathbf{D} \mathbf{Z}_i' + \mathbf{\Sigma}_i).$$

---

Example: Case Where $n_i = 2$

For a subject with two observations, we compute:

$$\mathbf{Z}_i = \begin{bmatrix} 1 & t_{i1} \\ 1 & t_{i2} \end{bmatrix}.$$

Then:

$$\mathbf{Z}_i \mathbf{D} \mathbf{Z}_i' = \begin{bmatrix} 1 & t_{i1} \\ 1 & t_{i2} \end{bmatrix} \begin{bmatrix} d_{11} & d_{12} \\ d_{12} & d_{22} \end{bmatrix} \begin{bmatrix} 1 & 1 \\ t_{i1} & t_{i2} \end{bmatrix}.$$

Expanding the multiplication:

$$\mathbf{Z}_i \mathbf{D} \mathbf{Z}_i' = \begin{bmatrix} d_{11} + 2d_{12}t_{i1} + d_{22} t_{i1}^2 & d_{11} + d_{12}(t_{i1} + t_{i2}) + d_{22}t_{i1}t_{i2} \\ d_{11} + d_{12}(t_{i1} + t_{i2}) + d_{22}t_{i1}t_{i2} & d_{11} + 2d_{12}t_{i2} + d_{22} t_{i2}^2 \end{bmatrix}.$$

Finally, incorporating the residual variance:

$$\text{Var}(Y_{i1}) = d_{11} + 2d_{12}t_{i1} + d_{22} t_{i1}^2 + \sigma^2.$$

Interpretation of the Marginal Model

1. Correlation in the Errors:
   • The off-diagonal elements of $\mathbf{Z}_i \mathbf{D} \mathbf{Z}_i'$ introduce correlation between observations within the same subject.
   • This accounts for the fact that repeated measurements on the same individual are not independent.
2. Variance Structure:

   • The variance function of $Y_{ij}$ is quadratic in time:

     $$\text{Var}(Y_{ij}) = d_{11} + 2d_{12}t_{ij} + d_{22}t_{ij}^2 + \sigma^2.$$

   • The curvature of this variance function is determined by $d_{22}$.

   • If $d_{22} > 0$, variance increases over time.

🔹 Key Takeaway:
In the hierarchical model, random effects contribute to the mean structure.
In the marginal model, they affect the covariance structure, leading to heteroskedasticity (changing variance over time) and correlation between repeated measures.

---

8.1.3 Random-Intercepts Model

The Random-Intercepts Model is obtained by removing random slopes, meaning:

• All subject-specific variability in slopes is attributed only to treatment differences.
• The model allows each subject to have their own intercept, but within each treatment group, all subjects share the same slope.

The hierarchical (conditional) model is:

$$\begin{aligned} \mathbf{Y}_i | b_i &\sim N(\mathbf{X}_i \beta + 1 b_i, \mathbf{\Sigma}_i), \\ b_i &\sim N(0, d_{11}). \end{aligned}$$

where:

• $b_i$ is the random intercept for subject $i$, assumed to follow $N(0, d_{11})$.
• $\mathbf{X}_i$ contains the fixed effects, which include treatment and time.
• $\mathbf{\Sigma}_i$ represents residual variance, typically assumed to be $\sigma^2 \mathbf{I}$ (independent errors).

Since there are no random slopes, the only source of subject-specific variability is the intercept.

Integrating out $b_i$ and assuming conditional independence $\mathbf{\Sigma}_i = \sigma^2 \mathbf{I}_{n_i}$, the marginal distribution of $\mathbf{Y}_i$ is:

$$\mathbf{Y}_i \sim N(\mathbf{X}_i \beta, 11'd_{11} + \sigma^2 \mathbf{I}).$$

The marginal covariance matrix is:

$$\begin{aligned} \text{Cov}(\mathbf{Y}_i) &= 11'd_{11} + \sigma^2 I \\ &= \begin{bmatrix} d_{11} + \sigma^2 & d_{11} & d_{11} & \dots & d_{11} \\ d_{11} & d_{11} + \sigma^2 & d_{11} & \dots & d_{11} \\ d_{11} & d_{11} & d_{11} + \sigma^2 & \dots & d_{11} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ d_{11} & d_{11} & d_{11} & \dots & d_{11} + \sigma^2 \end{bmatrix}. \end{aligned}$$

The correlation matrix is obtained by standardizing the covariance matrix:

$$\text{Corr}(\mathbf{Y}_i) = \begin{bmatrix} 1 & \rho & \rho & \dots & \rho \\ \rho & 1 & \rho & \dots & \rho \\ \rho & \rho & 1 & \dots & \rho \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \rho & \rho & \rho & \dots & 1 \end{bmatrix},$$

where:

$$\rho = \frac{d_{11}}{d_{11} + \sigma^2}.$$

This correlation structure is known as compound symmetry, meaning:

• Constant variance: $\text{Var}(Y_{ij}) = d_{11} + \sigma^2$ for all $j$.
• Equal correlation: $\text{Corr}(Y_{ij}, Y_{ik}) = \rho$ for all $j \neq k$.
• Positive correlation: Any two observations from the same subject are equally correlated.

Interpretation of $\rho$ (Intra-Class Correlation)

• $\rho$ is called the intra-class correlation coefficient (ICC).

• It measures the proportion of total variability that is due to between-subject variability:

  $$\rho = \frac{\text{Between-Subject Variance}}{\text{Total Variance}} = \frac{d_{11}}{d_{11} + \sigma^2}.$$

• If $\rho$ is large:

  • The inter-subject variability ($d_{11}$) is large relative to the intra-subject variability ($\sigma^2$).
  • This means subjects differ substantially in their intercepts.
  • Responses from the same subject are highly correlated.

• If $\rho$ is small:

  • The residual error variance dominates, meaning individual differences in intercepts are small.
  • Measurements from the same subject are weakly correlated.

---

Summary of the Random-Intercepts Model

Feature	Explanation
Intercepts	Random, subject-specific ($b_i$)
Slopes	Fixed within each treatment group
Covariance Structure	Compound symmetry (constant variance, equal correlation)
Intra-Class Correlation ($\rho$)	Measures between-subject variability relative to total variability
Interpretation	Large $\rho$ → Strong subject-level differences, Small $\rho$ → Mostly residual noise

---

8.1.4 Covariance Models in Linear Mixed Models

Previously, we assumed that the within-subject errors are conditionally independent, meaning:

$$\mathbf{\Sigma}_i = \sigma^2 \mathbf{I}_{n_i}.$$

However, real-world data often exhibit correlated errors, particularly in longitudinal studies where observations over time are influenced by underlying stochastic processes.

To model this dependence, we decompose the error term into two components:

$$\epsilon_i = \epsilon_{(1)i} + \epsilon_{(2)i},$$

where:

• $\epsilon_{(1)i}$ (Serial Correlation Component):
  • Represents subject-specific response to time-varying stochastic processes.
  • Captures dependency across observations for the same subject.
• $\epsilon_{(2)i}$ (Measurement Error Component):
  • Represents pure measurement error, assumed independent of $\epsilon_{(1)i}$.

Thus, the full LMM formulation becomes:

$$\mathbf{Y_i} = \mathbf{X_i \beta} + \mathbf{Z_i b_i} + \mathbf{\epsilon_{(1)i}} + \mathbf{\epsilon_{(2)i}}.$$

where:

• Random effects: $\mathbf{b_i} \sim N(\mathbf{0, D})$.
• Measurement errors: $\epsilon_{(2)i} \sim N(\mathbf{0, \sigma^2 I_{n_i}})$.
• Serial correlation errors: $\epsilon_{(1)i} \sim N(\mathbf{0, \tau^2 H_i})$.
• Independence assumption: $\mathbf{b}_i$ and $\epsilon_i$ are mutually independent.

---

To model the correlation structure of the serial correlation component $\epsilon_{(1)i}$, we define the $n_i \times n_i$ correlation (or covariance) matrix $\mathbf{H}_i$.

The $(j,k)$th element of $\mathbf{H}_i$ is denoted as:

$$h_{ijk} = g(t_{ij}, t_{ik}),$$

where:

• $h_{ijk}$ represents the covariance (or correlation) between time points $t_{ij}$ and $t_{ik}$.

• $g(t_{ij}, t_{ik})$ is a decreasing function that defines the correlation between time points $t_{ij}$ and $t_{ik}$.

• Typically, we assume that this function depends only on the time difference (stationarity assumption):

  $$h_{ijk} = g(|t_{ij} - t_{ik}|)$$

  meaning that the correlation only depends on the absolute time lag.

• To be a valid correlation matrix, we require:

  $$g(0) = 1.$$

  This ensures that each observation has a perfect correlation with itself.

---

Common Choices for $g(|t_{ij} - t_{ik}|)$

Several functions can be used to define the decay in correlation as time differences increase.

1. Exponential Correlation (Continuous-Time AR(1))

$$g(|t_{ij} - t_{ik}|) = \exp(-\phi |t_{ij} - t_{ik}|)$$

• Decay rate: Controlled by $\phi > 0$.
• Interpretation:
  • Observations closer in time are more correlated.
  • The correlation decreases exponentially as time separation increases.
• Use case: Often used in biological and economic time-series models.

---

2. Gaussian Correlation (Squared Exponential Kernel)

$$g(|t_{ij} - t_{ik}|) = \exp(-\phi (t_{ij} - t_{ik})^2)$$

• Decay rate: Faster than exponential.
• Interpretation:
  • If $\phi$ is large, correlation drops sharply as time separation increases.
  • Produces smooth correlation structures.
• Use case: Used in spatial statistics and machine learning (Gaussian processes).

---

3. Autoregressive (AR(1)) Correlation

A First-Order Autoregressive Model (AR(1)) assumes that the value at time $t$ depends on its previous value:

$$\alpha_t = \phi \alpha_{t-1} + \eta_t,$$

where

• $\eta_t \sim \text{i.i.d. } N(0, \sigma^2_\eta)$ (white noise process).

• $\phi$ is the autocorrelation coefficient.

Then, the covariance between two observations at different times is:

$$\text{Cov}(\alpha_t, \alpha_{t+h}) = \frac{\sigma^2_\eta \phi^{|h|}}{1 - \phi^2}.$$

Thus, the correlation between time points $t$ and $t+h$ is:

$$\text{Corr}(\alpha_t, \alpha_{t+h}) = \phi^{|h|}.$$

For a sequence of $T$ time points, the resulting Toeplitz correlation matrix is:

$$\text{Corr}(\alpha_T) = \begin{bmatrix} 1 & \phi^1 & \phi^2 & \dots & \phi^{T-1} \\ \phi^1 & 1 & \phi^1 & \dots & \phi^{T-2} \\ \phi^2 & \phi^1 & 1 & \dots & \phi^{T-3} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \phi^{T-1} & \phi^{T-2} & \phi^{T-3} & \dots & 1 \end{bmatrix}.$$

• Decay rate: $\phi$ controls how fast correlation decreases.
• Use case:
  • Common in time series models.
  • Appropriate when observations are equally spaced in time.

Properties of the AR(1) Structure:

• Correlation decreases with time lag:
  • Observations closer in time are more correlated.
  • Decay rate is controlled by $\phi$.
• Toeplitz structure:
  • The covariance matrix exhibits a banded diagonal pattern.
  • Useful in longitudinal and time-series data.
• Applicability:
  • Small $\phi$ ($\approx 0$) → Weak autocorrelation (errors are mostly independent).
  • Large $\phi$ ($\approx 1$) → Strong autocorrelation (highly dependent observations).

---

8.1.5 Covariance Structures in Mixed Models

Structure	Covariance Function $g(|t_{ij} - t_{ik}|)$	Decay Behavior	Correlation Matrix Pattern	Key Properties	Common Use Cases
Compound Symmetry	Constant: $g(|t_{ij} - t_{ik}|) = \rho$	No decay	Constant off-diagonal	Equal correlation across all time points (homogeneous)	Random-intercepts model, repeated measures
AR(1) (Autoregressive)	$\phi^{|t_{ij} - t_{ik}|}$	Exponential decay	Toeplitz (banded diagonal)	Strong correlation for nearby observations, weak for distant ones	Time series, longitudinal data
Exponential	$\exp(-\phi |t_{ij} - t_{ik}|)$	Smooth decay	Spatially continuous	Models continuous correlation decline over time/space	Growth curves, ecological models
Gaussian (Squared Exponential)	$\exp(-\phi (t_{ij} - t_{ik})^2)$	Rapid decay	Spatially continuous	Smooth decay but stronger locality effect	Spatial processes, geostatistics
Toeplitz	Varies by lag, $g(|t_{ij} - t_{ik}|) = c_h$	General pattern	General symmetric matrix	Arbitrary structure, allows irregular spacing	Irregular time points, flexible spatial models
Unstructured	Fully parameterized, no constraints	No fixed pattern	General symmetric matrix	Allows any correlation pattern, but many parameters	Small datasets with unknown correlation
Banded	Zero correlation beyond a certain lag	Abrupt cutoff	Banded diagonal	Assumes only close observations are correlated	Large datasets, high-dimensional time series
Spatial Power	$\phi^{|s_{ij} - s_{ik}|}$	Exponential decay	Distance-based structure	Used when correlation depends on spatial distance	Spatial statistics, environmental data
Matérn Covariance	Function of $|t_{ij} - t_{ik}|^\nu$	Flexible decay	Spatially continuous	Generalization of Gaussian and exponential structures	Geostatistics, spatiotemporal models

---

Choosing the Right Covariance Structure

Scenario	Recommended Structure
Repeated measures with equal correlation	Compound Symmetry
Longitudinal data with time dependence	AR(1), Toeplitz
Continuous-time process	Exponential, Gaussian
Spatial correlation	Spatial Power, Matérn
Irregularly spaced time points	Toeplitz, Unstructured
High-dimensional data	Banded, AR(1)

---

References

Laird, Nan M, and James H Ware. 1982. “Random-Effects Models for Longitudinal Data.” Biometrics, 963–74.

Schabenberger, Oliver, and Francis J Pierce. 2001. Contemporary Statistical Models for the Plant and Soil Sciences. CRC press.