Idea: If \(\sigma_j^2\) is large, then \(\left|y_{i, j}-\bar{y}_j\right|=\left|\hat{\epsilon}_{i, j}\right|\) should be large.

- Let \(z_{i, j}=\left|\hat{\epsilon}_{i j}\right|\)

- Use the ANOVA F-test for across-group differences in the \(z_{i, j}\) ’s

fit.nels<-lm(y.nels~as.factor(g.nels))
z.nels<-abs( fit.nels$res )
anova(lm(z.nels~as.factor(g.nels)) )

Analysis of Variance Table

Response: z.nels
                     Df Sum Sq Mean Sq F value    Pr(>F)    
as.factor(g.nels)   683  27078  39.645  1.6092 < 2.2e-16 ***
Residuals         12290 302776  24.636                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

\[\begin{gathered}\mathbf{y}_j=\mathbf{X}_j \boldsymbol{\beta}+\mathbf{Z}_j \mathbf{a}_j+\epsilon_j \\ \mathbf{E}\left[\begin{array}{l}\mathbf{a}_j \\ \epsilon_j\end{array}\right]=\left[\begin{array}{l}0 \\ 0\end{array}\right] \text { and } \operatorname{Cov}\left[\begin{array}{l}\mathbf{a}_j \\ \epsilon_j\end{array}\right]=\left[\begin{array}{ll}\Psi & 0 \\ 0 & \Sigma\end{array}\right]\end{gathered}\]

• Across-group heterogeneity: \(\Psi\) is the variance-covariance in \(\mathbf{a}_1, \ldots, \mathbf{a}_m\).

• Within-group heterogeneity: \(\Sigma\) is the variance-covariance of \(y_{1_j}, \ldots, y_{n_j j}\).

\[ \begin{aligned} y_{i, j} & =\mu+a_j+\epsilon_{i, j} \\ \left\{a_j\right\} & \sim \text { iid } N\left(0, \tau^2\right) \\ \left\{\epsilon_{i, j}\right\} & \sim \text { iid } N\left(0, \sigma^2\right) \end{aligned} \]

Express this model as \(\mathbf{y}_j=\mathbf{X}_j \boldsymbol{\beta}+\mathbf{Z}_j \mathbf{a}_j+\boldsymbol{\epsilon}_{\boldsymbol{j}}\)

• Regression parameters:

\[ \beta=\mu, a_j=a_j \]

• Design matrices:

\[ \mathbf{X}_j=\mathbf{Z}_j=\left[\begin{array}{c} 1 \\ \vdots \\ 1 \end{array}\right] \quad \text { for each } j \in\{1, \ldots, m\} \]

\[ \begin{aligned} y_{i, j} & =\beta^T \mathbf{x}_{i, j}+\mathbf{a}_j^T \mathbf{x}_{i, j}+\epsilon_{i, j} \\ \left\{\mathbf{a}_j\right\} & \sim \text { iid } N(0, \Psi) \\ \left\{\epsilon_{i, j}\right\} & \sim \text { iid } N\left(0, \sigma^2\right) \end{aligned} \]

Express this model as \(\mathbf{y}_j=\mathbf{X}_j \boldsymbol{\beta}+\mathbf{Z}_j \mathbf{a}_j+\boldsymbol{\epsilon}_j\)

• Design matrices: \[ \mathbf{X}_j=\mathbf{Z}_j=\left[\begin{array}{c} \mathbf{x}_{1, j} \rightarrow \\ \vdots \\ \mathbf{x}_{n_j j} \rightarrow \end{array}\right] \quad \text { fir each } j \in\{1, \ldots, m\} \]

• Regression parameters: \(\boldsymbol{\beta}=\boldsymbol{\beta}, \mathbf{a}_j=\mathbf{a}_j\)

• Covariance terms: \(\Psi=\operatorname{Cov}\left[\mathbf{a}_j\right], \Sigma=\sigma^2 \mathbf{I}\)

This is just a special case where \(\mathbf{X}_j=\mathbf{Z}_j\).

1. OLS is still unbiased, but not BLUE

• Unbiasedness: The OLS estimate for \(\boldsymbol{\beta}\) remains unbiased if the random effects \(\mathbf{a}_j\) are independent of the predictors \(\mathbf{X}_j\). This means that \(\mathbb{E}\left[\mathbf{y}_j \mid \mathbf{X}_j\right]=\mathbf{X}_j \boldsymbol{\beta}\), so OLS can recover the fixed-effect coefficients without bias.

• Not BLUE

2. Variance of OLS estimate is no longer \(\sigma^2\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1}\)

• In the presence of random effects \(\left(\mathbf{a}_j\right)\), the true variance of \(\mathbf{y}_j\) includes both:
  • Within-group error variance ( \(\sigma^2\), from \(\epsilon_j\) ).
  • Between-group variance (from \(\mathbf{Z}_j \mathbf{a}_j\) ).

3. BLUE is (approximately) the MLE returned by lmer

• The Best Linear Unbiased Estimator (BLUE) for \(\boldsymbol{\beta}\) accounts for both fixed effects ( \(\mathbf{X}_j\) ) and random effects \(\left(\mathbf{Z}_j \mathbf{a}_j\right)\) in the data. This is achieved by using the mixed-effects model, where the estimates are computed by maximizing the marginal likelihood of the data: \[\mathbf{y}_j \sim N\left(\mathbf{X}_j \boldsymbol{\beta}, \mathbf{V}\right)\] where \(\mathbf{V}=\mathbf{Z}_j \mathbf{\Psi} \mathbf{Z}_j^{\top}+\sigma^2 \mathbf{I}\) is the covariance matrix incorporating random effects ( \(\left.\mathbf{\Psi}\right)\) and residual variance ( \(\sigma^2\) ).

• The lmer function in R uses the maximum likelihood (ML) methods to estimate \(\boldsymbol{\beta}, \mathrm{\Psi}\), and \(\sigma^2\), which are optimal under the assumed model.

Why Use lmer?

The lmer function (from the lme4 package) provides:

1. Efficient Estimates: It calculates fixed-effect coefficients \((\boldsymbol{\beta})\) and random-effect variances (\(\Psi\)) simultaneously, producing estimates that are optimal under the model.
2. Correct Standard Errors: It computes standard errors that account for the hierarchical structure of the data, reflecting both within-group and between-group variability.
3. Flexibility: Allows modeling of complex group structures (e.g., random intercepts, random slopes, or crossed random effects).