Tutorial

Consider the following model:

\[ y_{ij} = \beta_{0} + \beta_{1} x_{i} + \gamma_i + \epsilon_{ij}, \quad \epsilon_{ij} \sim \text{iid } N(0, \sigma^2) \]

A concrete example is where:

• \(x_{i} =\) age of \(i\)-th participant in a study

• \(y_{ij} =\) systolic blood pressure measured at \(j\)-th occasion \((\text{say } j = 1, 2, 3)\)

Here we are primarily interested in \((\beta_{0}, \beta_{1})\), which determines the expected \(y\) given \(x\). The individual-specific effect \(\gamma_i\) is included to acknowledge that repeated measurements for the same individuals are likely to be closer to each other than measurements from different individuals.

A dataset is available at https://www.isid.ac.in/~statcomp2026/datasets/nhanes-bp-repeated.csv, which we can read in as follows:

bpsys <- read.csv("http://www.isid.ac.in/~statcomp2026/datasets/nhanes-bp-repeated.csv")
names(bpsys) <- c("SEQN", "age", "gender", "ethnicity", "time", "bpsys")
bpsys$SEQN <- factor(bpsys$SEQN)
xtabs(~ ethnicity + gender, bpsys)

                                     gender
ethnicity                             Female Male
  Mexican American                      1942 1864
  Non-Hispanic Black                    4272 3992
  Non-Hispanic White                    5258 5383
  Other Hispanic                        1549 1483
  Other Race - Including Multi-Racial   2578 2644

bpsys_wm <- subset(bpsys, age >= 20 & gender == "Male" & ethnicity == "Non-Hispanic White")
bpsys_mf <- subset(bpsys, age >= 20 & gender == "Female" & ethnicity == "Mexican American")
bpsys_sub <- bpsys_mf

A wrong option is to fit a linear model:

fm1 <- lm(bpsys ~ 1 + age, data = bpsys_sub)
fm1

Call:
lm(formula = bpsys ~ 1 + age, data = bpsys_sub)

Coefficients:
(Intercept)          age  
    82.9327       0.7576  

A natural and basic question that we may be interested in asking is whether systolic blood pressure is independent of age. For the linear model above, this is equivalent to testing \(H_0: \beta_{1} = 0\), but a more useful summary is to look at a confidence interval for \(\beta_{1}\). Confidence intervals for the model parameters in fm1 can be obtained as follows.

confint(fm1, level = 0.95)

                 2.5 %     97.5 %
(Intercept) 80.1202540 85.7450647
age          0.6997983  0.8153212

A simple and mostly correct option is to take the average of the repeated measurements:

\[ \bar{y}_{i\cdot} = \beta_{0} + \beta_{1} x_{i} + \gamma_i + \bar{\epsilon}_{i\cdot} = \beta_{0} + \beta_{1} x_{i} + \delta_i \]

where \(\delta_i = \gamma_i + \bar{\epsilon}_{i\cdot}\) combines the individual variations \(\gamma_i\) and the measurement errors \(\epsilon_{ij}\). This is reasonable as long as we consider \(\delta_i\) to be mean 0 errors as well. Specifically, if we assume \(\gamma_i \sim N(0, \tau^2)\), independent of \(\epsilon_{ij}\), then

\[ E[\delta_i] = 0, \quad V[\delta_i] = V[\gamma_i + \bar{\epsilon}_{i\cdot}] = \tau^2 + \frac{\sigma^2}{n_i}, \]

where \(n_i\) is the number of replications for the \(i\)th individual. However, if the number of replications \(n_i\) are different for different individuals, then the problem is not straightforward.

Task 1: Compute the average \(\bar{y}_{i\cdot}\) values and fit a suitable linear regression model assuming all \(n_i = 3\). Obtain the corresponding confidence interval for \(\beta_{1}\), and compare it with the interval obtained from fm1 above.

NoteSolution

bpsys_avg <-
    tapply(bpsys_sub, ~ SEQN,
           with, data.frame(age = mean(age),
                            bpsys = mean(bpsys))) |>
        array2DF()
lm(bpsys ~ 1 + age, data = bpsys_avg) |> confint(level = 0.95)

                 2.5 %     97.5 %
(Intercept) 77.9249455 87.4417149
age          0.6666168  0.8618804

Approach A: Mixed Effects Model

A natural model is the following mixed effects model:

\[ y_{ij} = \beta_{0} + \beta_{1} x_{i} + \gamma_i + \epsilon_{ij}, \quad \gamma_i \sim N(0, \tau^2), \quad \epsilon_{ij} \sim N(0, \sigma^2) \]

We can formulate, in principle, the corresponding likelihood \(L(\beta_{0}, \beta_{1}, \tau^2, \sigma^2)\) and maximize it, although in practice the details can be complicated.

library(lme4)

Loading required package: Matrix

fm2 <- lmer(bpsys ~ 1 + age + (1|SEQN), data = bpsys_sub, REML = FALSE)
summary(fm2)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: bpsys ~ 1 + age + (1 | SEQN)
   Data: bpsys_sub

      AIC       BIC    logLik -2*log(L)  df.resid 
   8975.0    8995.6   -4483.5    8967.0      1257 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-9.2848 -0.4166 -0.0177  0.4207  4.9087 

Random effects:
 Groups   Name        Variance Std.Dev.
 SEQN     (Intercept) 239.77   15.48   
 Residual              22.28    4.72   
Number of obs: 1261, groups:  SEQN, 421

Fixed effects:
            Estimate Std. Error t value
(Intercept) 82.70441    2.41367   34.27
age          0.76368    0.04953   15.42

Correlation of Fixed Effects:
    (Intr)
age -0.948

anova(fm2, fm1) # test for tau = 0

Data: bpsys_sub
Models:
fm1: bpsys ~ 1 + age
fm2: bpsys ~ 1 + age + (1 | SEQN)
    npar   AIC     BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)    
fm1    3 10599 10614.5 -5296.5     10593                         
fm2    4  8975  8995.6 -4483.5      8967 1626.1  1  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

As before, confidence intervals for the parameters can be obtained using confint().

confint(fm2, level = 0.95)

Computing profile confidence intervals ...

                 2.5 %    97.5 %
.sig01      14.4618666 16.626420
.sigma       4.5033716  4.955578
(Intercept) 77.9626970 87.445757
age          0.6663939  0.860982

Notice that the confidence intervals are considerably wider than the ones obtained from the wrong model fm1, which gains false confidence by ignoring the correlation in the errors within subjects. On the other hand, they should closely match the confidence interval you get after averaging.

Approach B: Regularization

Another option is to take a regularization approach by saying that we want both \(\gamma_i^2\) and \(\epsilon_{ij}^2\) to be small. Of course we need to give some relative weights to these errors, so we may minimize, given \(\lambda > 0\):

\[ q(\boldsymbol{\beta}, \boldsymbol{\gamma}) = \sum_{i}\sum_{j} \epsilon_{ij}^2 + \lambda \sum_{i} \gamma_i^2 = \sum_{i}\sum_{j} (y_{ij} - \beta_{0} - \beta_{1} x_{i} - \gamma_i)^2 + \lambda \sum_{i} \gamma_i^2. \]

In terms of a suitable design matrix \(X\), this may be written as:

\[ q(\boldsymbol{\beta}, \boldsymbol{\gamma}) = (\mathbf{y} - X\boldsymbol{\theta})'( \mathbf{y} - X\boldsymbol{\theta}) + \lambda \boldsymbol{\gamma}'\boldsymbol{\gamma}, \quad \text{where } \boldsymbol{\theta} = \begin{pmatrix} \boldsymbol{\beta} \\ \boldsymbol{\gamma} \end{pmatrix} \]

Here \(X = [X_1 \mid X_2]\), \(X\boldsymbol{\theta} = X_1 \boldsymbol{\beta} + X_2 \boldsymbol{\gamma}\), and:

\[q(\boldsymbol{\beta}, \boldsymbol{\gamma}) = (\mathbf{y} - X_1\boldsymbol{\beta} - X_2\boldsymbol{\gamma})'(\mathbf{y} - X_1\boldsymbol{\beta} - X_2\boldsymbol{\gamma}) + \lambda \boldsymbol{\gamma}'\boldsymbol{\gamma}\]

Differentiating w.r.t. \(\boldsymbol{\theta}\), we get the following score equations:

• Objective function:

\[ q(\boldsymbol{\beta}, \boldsymbol{\gamma}) = \mathbf{y}'\mathbf{y} + \boldsymbol{\beta}'X_1'X_1\boldsymbol{\beta} + \boldsymbol{\gamma}'X_2'X_2\boldsymbol{\gamma} + 2\boldsymbol{\beta}'X_1'X_2\boldsymbol{\gamma} - 2\mathbf{y}'X_1\boldsymbol{\beta} - 2\mathbf{y}'X_2\boldsymbol{\gamma} + \lambda \boldsymbol{\gamma}'\boldsymbol{\gamma} \]

• Derivative with respect to \(\boldsymbol{\beta}\):

\[ \nabla_{\boldsymbol{\beta}} q(\boldsymbol{\beta}, \boldsymbol{\gamma}) = 2X_1'X_1\boldsymbol{\beta} + 2X_1'X_2\boldsymbol{\gamma} - 2X_1'\mathbf{y} = 0 \implies X_1'X_1\boldsymbol{\beta} + X_1'X_2\boldsymbol{\gamma} = X_1'\mathbf{y} \]

• Derivative with respect to \(\boldsymbol{\gamma}\):

\[ \nabla_{\boldsymbol{\gamma}} q(\boldsymbol{\beta}, \boldsymbol{\gamma}) = 2X_2'X_2\boldsymbol{\gamma} + 2X_2'X_1\boldsymbol{\beta} - 2X_2'\mathbf{y} + 2\lambda\boldsymbol{\gamma} = 0 \implies (X_2'X_2 + \lambda I)\boldsymbol{\gamma} + X_2'X_1\boldsymbol{\beta} = X_2'\mathbf{y} \]

Together, this gives the system:

\[\begin{bmatrix} X_1'X_1 & X_1'X_2 \\ X_2'X_1 & X_2'X_2 + \lambda I \end{bmatrix} \begin{pmatrix} \boldsymbol{\beta} \\ \boldsymbol{\gamma} \end{pmatrix} = \begin{pmatrix} X_1'\mathbf{y} \\ X_2'\mathbf{y} \end{pmatrix}\]

\[\left(X'X + \lambda \begin{bmatrix} 0 & 0 \\ 0 & I \end{bmatrix}\right) \begin{pmatrix} \boldsymbol{\beta} \\ \boldsymbol{\gamma} \end{pmatrix} = X'\mathbf{y}\]

We could also use a \(L_1\) penalty on the \(\gamma_i\)-s, and try to minimize

\[ q(\boldsymbol{\beta}, \boldsymbol{\gamma}) = \sum_{i}\sum_{j} \epsilon_{ij}^2 + \lambda \sum_{i} \gamma_i^2 = \sum_{i}\sum_{j} (y_{ij} - \beta_{0} - \beta_{1} x_{i} - \gamma_i)^2 + \lambda \sum_{i} \lvert \gamma_i \rvert. \]

Task 2: Use glmnet() to fit both Ridge (alpha=0) and lasso (alpha=1) to solve this problem. Use CV to choose \(\lambda\).

Approach C: Bayesian Approach

Another approach is a Bayesian one. If we assume that \(\tau^2\) and \(\sigma^2\) are known, then:

\[\boldsymbol{\gamma} \sim N_{p_2}(\mathbf{0}, \tau^2 I_{p_2})\]

\[\mathbf{y} \mid \boldsymbol{\gamma} \sim N_n(X_1\boldsymbol{\beta} + X_2\boldsymbol{\gamma}, \sigma^2 I_n)\]

which gives a posterior distribution of \(\boldsymbol{\gamma}\):

\[ \boldsymbol{\gamma} \mid \mathbf{y} \sim N_{p_2}(\dots, \dots) \]

If we complete the calculations, we will see that we the posterior mean only depends on the ratio \(\tau^2 / \sigma^2\).

\(\boldsymbol{\beta}\) can be handled similarly if we take an improper prior of the form \(\boldsymbol{\beta} \sim N_{p_1}(\mathbf{0}, \infty I)\).

In a “full Bayesian” approach, we can additionally postulate a prior for \(\sigma^2\) and \(\tau^2\) as well.