Estimation of a Linear Model in Terms of Intra-Class Correlations of the Residual Error and the Regressors

Juha Lappi

University of Eastern Finland, Joensuu, Finland.

DOI: 10.4236/ojs.2022.122013                

Abstract

Objectives: The objective is to analyze the interaction of the correlation structure and values of the regressor variables in the estimation of a linear model when there is a constant, possibly negative, intra-class correlation of residual errors and the group sizes are equal. Specifically: 1) How does the variance of the generalized least squares (GLS) estimator (GLSE) depend on the regressor values? 2) What is the bias in estimated variances when ordinary least squares (OLS) estimator is used? 3) In what cases are OLS and GLS equivalent. 4) How can the best linear unbiased estimator (BLUE) be constructed when the covariance matrix is singular? The purpose is to make general matrix results understandable. Results: The effects of the regressor values can be expressed in terms of the intra-class correlations of the regressors. If the intra-class correlation of residuals is large, then it is beneficial to have small intra-class correlations of the regressors, and vice versa. The algebraic presentation of GLS shows how the GLSE gives different weight to the between-group effects and the within-group effects, in what cases OLSE is equal to GLSE, and how BLUE can be constructed when the residual covariance matrix is singular. Different situations arise when the intra-class correlations of the regressors get their extreme values or intermediate values. The derivations lead to BLUE combining OLS and GLS weighting in an estimator, which can be obtained also using general matrix theory. It is indicated how the analysis can be generalized to non-equal group sizes. The analysis gives insight to models where between-group effects and within-group effects are used as separate regressors.

Keywords

Best Linear Unbiased Estimator, Ordinary Least-Squares, Generalized Least Squares, Singular Correlation Matrix, Between-Group Effects, Within-Group Effects

1. Introduction

Ordinary least squares (OLS) and generalized least squares (GLS) are basic methods to estimate coefficients of regressor variables in regression equations that are linear with respect to coefficients. Let OLSE and GLSE denote the corresponding estimators. OLS minimizes the sum of squared residuals (deviations from the estimated model), while GLS minimizes the estimation variances of the coefficient estimates if the residual errors (deviations from the true model) have non-equal variances and/or are correlated between different observations. OLS estimates are unbiased also for cases where GLS estimates have smaller variances, but the estimated variances of the coefficient estimates can be biased, leading to biased t- and F-tests. GLS estimation utilizes the inverse of the covariance matrix of residual errors. Thus, GLSE does not exist if the inverse does not exist, i.e., the covariance matrix is singular. If the inverse exists, then GLSE is the best linear unbiased estimator, BLUE. It requires deeper matrix algebra to derive BLUE when GLSE does not exist. Further complications are caused if there are linear dependencies between predictor variables, mathematically if the model matrix, to be presented shortly, does not have full column rank. In practice, such difficulties can be avoided by dropping linearly dependent predictors. It is assumed here that there are no such dependencies.

In the history of statistics, attention has been given to analyze what are the necessary and sufficient conditions that guarantee that OLSE is BLUE. Reference [1] describes the influential role of C. R: Rao in this question.

When analyzing grouped data, a common model for correlated residual errors is to assume that there is a constant correlation between all observations of each group, and members of different groups are uncorrelated. Let ${\rho }_e$ denote this intra-class correlation, and let i. denote group, and j denote an observation number within a group. A non-negative intra-class correlation of residual error is implied by assuming a variance component model $e_{ij}=v_i+{\epsilon }_{ij}$, where $v_i$ is the random group effect with mean zero and variance ${\sigma }_v^2$ and ${\epsilon }_{ij}$ is a random individual effect with variance ${\sigma }_{\epsilon }^2$ and uncorrelated with $v_i$. Thus, ${\rho }_e={\sigma }_v^2/\left({\sigma }_v^2+{\sigma }_{\epsilon }^2\right)$. The intra-class correlation of residual errors can be also negative. If ${\rho }_e$ is positive, then two observations taken from a given group are, on average, more alike than two observations taken from different groups. When ${\rho }_e$ is negative, then two observations taken from a given group are, on average, more different that two observations taken randomly from the whole population. Positive ${\rho }_e$ implies that group averages have a larger variance than a random sample of equal size from the whole population, while negative ${\rho }_e$ implies that group averages vary less than averages of samples from the whole population. On the extreme, all group means are equal, and all the variation is between individuals. While assuming a positive ${\rho }_e$, it is not necessary to specify the group size (of course, when dealing with data, group size plays an important role). The effect of a negative ${\rho }_e$ is always dependent on the group size, and it does not make sense to assume constant negative ${\rho }_e$ for different group sizes. The possibility of a negative ${\rho }_e$ is too often ignored. Even if the overall ${\rho }_e$ is positive, due to natural variation between groups, the competition between individuals from limited resources is decreasing ${\rho }_e$, and this competition should be described when analyzing grouped data. The lower limit of ${\rho }_e$ can be obtained from the condition $\mathrm{var}\left({\displaystyle {\sum }_{j=1}^n{e}_{ij}}\right)=n{\sigma }_e^2+n\left(n-1\right){\rho }_e{\sigma }_e^2\ge 0$, which implies that ${\rho }_e\ge -1/\left(n-1\right)$. The upper limit of ${\rho }_e$ is 1. This paper assumes a constant ${\rho }_e$ and a constant group size, which allows a straightforward treatment of negative ${\rho }_e$.

Even if the matrix formulas of OLSE and GLSE are simple, it may remain intuitively unclear how the correlation structure of residual errors and the values of the regressor variables interact. It is shown that this interaction can be described in terms of the empirical intra-class correlations of the regressor variables. Closed form algebraic formulas are first derived for a special case where a simple linear model is estimated. These derivations are then generalized to a model with several regressors. They give a concrete and easily understood explanation and illustration for the matrix theory results which go beyond the mathematical expertise of most scientists who are just interested to apply statistical methods in research

This paper discusses the following questions: 1) How do the variances of the GLS estimates depend on the intra-class correlation of residuals and on the values of the x-variable in a simple linear model? 2) What are the ratios of the GLS variances and OLS variances? 3) What is the bias in the estimated variances of the coefficient estimates when OLS is used? 4) How does GLS utilize the group means of the regressors and the deviations from the group means? 5) How this formulation shows why OLS is sometimes equal to GLS? 6) How can a BLUE be constructed when the correlation matrix of residual errors is singular, i.e., ${\rho }_e=1$ or ${\rho }_e=-1/\left(n-1\right)$ and the intra-call correlations of regressors gets any values? Similarly, as with ${\rho }_e$, the extreme values of the regressor correlations need to be considered separately. This last BLUE part may not be of interest to practicing scientists, but it may serve as an introduction to the BLUE problem for researchers, who cannot fluently read matrix theory. For matrix theory experts this part may be trivial or self-evident.

Reference [2] compared the GLS estimation of a simple regression line with OLS estimation using differences in the predictor variable. This work is here put into a general context. The results reported in [3] regarding the bias of F-tests when OLS is used, are exemplified and generalized. Algebraic formulas show how OLSE, GLSE, BLUE and singular correlation matrices are related, and thus illustrate the matrix results presented in [4] [5], and [6]. Reference [5] provides the matrix conditions that are used to show that a suggested estimator is BLUE.

Let us then move to more formal definitions. Let us assume a model $y=Xb+e$, where $y$ is a N dimensional random vector, $X$ is $N\times p$ model matrix, $b$ is p dimensional fixed coefficient vector, $E\left(e\right)=0$, and $\mathrm{var}\left(e\right)={\sigma }^{2}V$. If $b$ is estimated with the GLS estimator ${\hat{b}}_{\text{GLS}}={\left(X^{\prime }{V}^{-1}X\right)}^{-1}{X}^{\prime }{V}^{-1}y$, then $\mathrm{var}\left({\hat{b}}_{\text{GLS}}\right)={\sigma }^2{\left(X^{\prime }{V}^{-1}X\right)}^{-1}$. If $b$ is estimated with the OLS estimator ${\hat{b}}_{\text{OLS}}={\left(X^{\prime }X\right)}^{-1}{X}^{\prime }y$ then $E\left({\hat{b}}_{\text{OLS}}\right)=b$, and $\mathrm{var}\left({\hat{b}}_{\text{OLS}}\right)={\sigma }^2{\left(X^{\prime }X\right)}^{-1}{X}^{\prime }VX{\left(X^{\prime }X\right)}^{-1}$. In the OLS regression, ${\sigma }^2$ is estimated with ${\hat{\sigma }}_{\text{OLS}}^2={\left(N-p\right)}^{-1}{u}^{\prime }u$, where $u$ is the residual vector $y-X{\hat{b}}_{\text{OLS}}$. Then, $u^{\prime }u=e^{\prime }Me$ where $M=I-X{\left(X^{\prime }X\right)}^{-1}{X}^{\prime }$. The expected value of $e^{\prime }Me$ is ${\sigma }^2\mathrm{trace}\left(MV\right)$. Thus, $E\left({\hat{\sigma }}_{\text{OLS}}^2\right)={\left(N-p\right)}^{-1}\mathrm{trace}\left(MV\right){\sigma }^2$.

It is now assumed that that data consist of $K\ge 2$ groups, indexed with i, each group has $n\ge 2$ observations, indexed with j. It is denoted that $N=Kn$. $V$ is block diagonal so that each block $V_i$ has a compound symmetry correlation structure where the non-diagonal elements are equal to the intra-class correlation ${\rho }_e$. The term “intra-class correlation” is used for historical reasons, but otherwise the term “group” is used. The intercept is the first parameter in $b$, i.e., the first column in $X$ is vector of ones, $1_N$. It is assumed further that other columns are centered, i.e., the sum of elements of the column is zero. This assumption makes the analysis simpler. Total and group averages are denoted as $\bar{x}$, ${\bar{x}}_i$, $\bar{y}$ and ${\bar{y}}_i$.

Closed form algebraic equations are initially derived for the simple linear model $y_{ij}=a+b{x}_{ij}+e_{ij}$. A model with several regressors is presented later. Centering means that the working model is $y_{ij}=a^{\prime }+b\left(x-\bar{x}\right)+e_{ij}$. Estimates of $a^{\prime }$ and b are uncorrelated both in OLS and GLS. After estimating $a^{\prime }$ and b, an estimate of a is obtained from $\hat{a}={\hat{a}}^{\prime }-\hat{b}\bar{x}$. Thereafter $\mathrm{var}\left(\hat{a}\right)=\mathrm{var}\left({\hat{a}}^{\prime }\right)+{\bar{x}}^2\mathrm{var}\left(\hat{b}\right)$, and $\mathrm{cov}\left(\hat{a},\hat{b}\right)=\mathrm{var}\left({\hat{a}}^{\prime }\right)$.

2. GLS in the Simple Linear Model

First, let us assume that $-1/\left(n-1\right)<{\rho }_e<1$, which means that $V$ is non-singular. Then $V_i^{-1}$ is a matrix with diagonal elements:

$\alpha =\left(1+\left(n-2\right){\rho }_e\right)/\left(\left(1-{\rho }_e\right)\left(1+\left(n-1\right){\rho }_e\right)\right)$ (1)

The non-diagonal elements are equal to:

$\beta =-{\rho }_e/\left(\left(1-{\rho }_e\right)\left(1+\left(n-1\right){\rho }_e\right)\right)$ (2)

Thereafter $X^{\prime }{V}^{-1}X$ is

$X^{\prime }{V}^{-1}X=\left(\begin{array}{cc}N\left(\alpha +\left(n-1\right)\beta \right)& 0\\ 0& \alpha {\displaystyle \underset{i=1}{\overset{K}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\left(x_{ij}-\bar{x}\right)}^2}}+\beta {\displaystyle \underset{i=1}{\overset{K}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\displaystyle \underset{j^{\prime }\ne j}{\overset{}{\sum }}\left(x_{ij}-\bar{x}\right)\left(x_{i{j}^{\prime }}-\bar{x}\right)}}}\end{array}\right)$ (3)

The second diagonal element in (3) can be expressed as ${\displaystyle \underset{i=1}{\overset{K}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\left(x_{ij}-\bar{x}\right)}^2}}\left(\alpha +\beta R_x\right)$ using the empirical intra-class correlation $R_x$, which can be presented in two equivalent forms:

$R_x={\displaystyle \underset{i=1}{\overset{K}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\displaystyle \underset{j^{\prime }\ne j}{\overset{}{\sum }}\left(x_{ij}-\bar{x}\right)\left(x_{i{j}^{\prime }}-\bar{x}\right)}}}/\left(n-1\right){\displaystyle \underset{i=1}{\overset{K}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\left(x_{ij}-\bar{x}\right)}^2}}$ (4)

$R_x=n^2{\displaystyle \underset{i}{\overset{K}{\sum }}{\left({\bar{x}}_i-\bar{x}\right)}^2}/\left(n-1\right){\displaystyle \underset{i=1}{\overset{K}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\left(x_{ij}-\bar{x}\right)}^2}}-\frac{1}{n-1}$ (5)

Reference [7] explains how the intra-class correlation, dating back to Fisher, generalizes the correlation idea in order to measure the similarity of group members. Equation (4) appears as covariance divided by variance, and (5) is almost equal to the sample variance of group means divided by the total sample variance, which resembles the definition of ${\rho }_e$. Equation (5) indicates that $R_x\ge -1/\left(n-1\right)$. The lower limit is obtained when all group means ${\bar{x}}_i$ are equal. The upper limit $R_x=1$ is obtained when $x_{ij}={\bar{x}}_i$ for all i and j. If $R_x=0$, then ${\displaystyle \underset{i}{\overset{K}{\sum }}{\left({\bar{x}}_i-\bar{x}\right)}^2}/K$, the population variance of ${\bar{x}}_i$, is equal to the population variance ${\displaystyle \underset{i=1}{\overset{K}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\left(x_{ij}-\bar{x}\right)}^2}}/N$ divided by n, resembling the variance of the mean of uncorrelated random variables.

Matrix ${\left(X^{\prime }{V}^{-1}X\right)}^{-1}$ is obtained by taking the reciprocals of the diagonal elements of (3). The variances of the estimates of $a^{\prime }$ considered here are proportional to $C_{a^{\prime }}={\sigma }^2/N$, and the variances of the estimate of b are proportional

to $C_b={\sigma }^2/{\displaystyle \underset{i=1}{\overset{K}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\left(x_{ij}-\bar{x}\right)}^2}}$. When ${\rho }_e=0$, $\mathrm{var}\left({\hat{a}}^{\prime }\right)=C_a$, and $\mathrm{var}\left(\hat{b}\right)=C_b$. First, we note that ${{\hat{a}}^{\prime }}_{\text{GLS}}={{\hat{a}}^{\prime }}_{\text{OLS}}=\bar{y}$. Thus, we directly derive that $\mathrm{var}\left({\hat{a}}^{\prime }\right)=C_a\left(1+\left(n-1\right)r_e\right)$. If $r_e=-1/\left(n-1\right)$, then $\mathrm{var}\left({\hat{a}}^{\prime }\right)=0$. For b, we obtain:

$\mathrm{var}\left({\hat{b}}_{\text{GLS}}\right)=C_b\frac{\left(1-{\rho }_e\right)\left(1+\left(n-1\right){\rho }_e\right)}{1+\left(n-2\right){\rho }_e-\left(n-1\right)R_x{\rho }_e}$ (6)

If ${\rho }_e=0$ or ${\rho }_e=R_x$, then $\mathrm{var}\left({\hat{b}}_{\text{GLS}}\right)=C_b$. It holds that $\mathrm{var}\left({\hat{b}}_{\text{GLS}}\right)>C_b$, if ${\rho }_e$ is between 0 and $R_x$. Equation (6) is an increasing function of $R_x{\rho }_e$, an increasing function of $R_x$ when ${\rho }_e>0$, and a decreasing function of $R_x$ when ${\rho }_e<0$. When $R_x=1$, then $\mathrm{var}\left({\hat{b}}_{\text{GLS}}\right)=C_b\left(1+\left(n-1\right){\rho }_e\right)$. When $R_x=-1/\left(n-1\right)$, then ${\mathrm{var}}_{\text{GLS}}\left(\hat{b}\right)=C_b\left(1-{\rho }_e\right)$. The dependency of $\mathrm{var}\left({\hat{b}}_{\text{GLS}}\right)/C_b$ on ${\rho }_e$ for different values of $R_x$ is shown in Figure 1.

3. OLS in Simple Linear Model When ${\rho }_e\ne 0$

What happens if OLS is used when ${\rho }_e\ne 0$ ? Using $\mathrm{var}\left({\hat{b}}_{\text{OLS}}\right)$, provided in the introduction section, we obtain:

$\mathrm{var}\left({\hat{b}}_{\text{OLS}}\right)=C_b\left(1+\left(n-1\right)R_x{\rho }_e\right)$ (7)

When ${\rho }_e\ne 0$ and OLS are used to estimate the variances of the parameter estimates, the bias combines the bias of the estimate of ${\sigma }^2$ and the errors of ${\mathrm{var}}_{\text{OLS}}\left({{\hat{a}}^{\prime }}_{\text{OLS}}\right)/{\sigma }^2$ and ${\mathrm{var}}_{\text{OLS}}\left({{\hat{b}}^{\prime }}_{\text{OLS}}\right)/{\sigma }^2$ where ${\mathrm{var}}_{\text{OLS}}$ indicates the variances implied by OLS assumptions (these are not biases, as they are not related to the expected values). In our case:

$E\left({\hat{\sigma }}_{\text{OLS}}^2\right)=\left(1-\frac{\left(n-1\right)\left(1+R_x\right){\rho }_e}{N-2}\right){\sigma }^2$ (8)

It holds that $E\left({\hat{\sigma }}_{\text{OLS}}^2\right)<{\sigma }^2\iff {\rho }_e>0$. Now,

Figure 1. The dependency of $\mathrm{var}\left({\hat{b}}_{\text{GLS}}\right)/C_b$ in (6) on ${\rho }_e$, the intra-class correlation of residual errors, for different values of $R_x$, the intra-class correlation of x (denoted on the curves) when n = 3. Blue lines are for negative values of $R_x$ and red lines for positive values. The thick line is for $R_x=0$. Note that the same variance is obtained for ${\rho }_e=0$ and ${\rho }_e=R_x$.

$\begin{array}{c}\mathrm{v}\hat{a}{r}_{\text{OLS}}\left({{\hat{a}}^{\prime }}_{\text{OLS}}\right)={\hat{\sigma }}_{\text{OLS}}^2/N={\hat{\sigma }}_{\text{OLS}}^2\mathrm{var}\left({{\hat{a}}^{\prime }}_{\text{OLS}}\right)/\left(N\mathrm{var}\left({{\hat{a}}^{\prime }}_{\text{OLS}}\right)\right)\\ ={\hat{\sigma }}_{\text{OLS}}^2\mathrm{var}\left({{\hat{a}}^{\prime }}_{\text{OLS}}\right)/\left(1+\left(n-1\right){\rho }_e\right)\end{array}$ (9)

In (9), $\mathrm{v}\hat{a}{r}_{\text{OLS}}$ denotes the variance estimate we obtain in OLS. If it is combined with (8), we can obtain $E\left(\mathrm{v}\hat{a}{r}_{\text{OLS}}\left({\hat{a}}^{\prime }\right)\right)$ in terms of $\mathrm{var}\left({a^{\prime }}_{\text{OLS}}\right)$. The estimated variance is biased downwards if ${\rho }_e>0$. If we similarly combine (8) and (7) we obtain:

$E\left(\mathrm{v}\hat{a}{r}_{\text{OLS}}\left({\hat{b}}_{\text{OLS}}\right)\right)=\frac{N-2-\left(n-1\right)\left(1+R_x\right){\rho }_e}{\left(N-2\right)\left(1+\left(n-1\right)R_x{\rho }_e\right)}\mathrm{var}\left({\hat{b}}_{\text{OLS}}\right)$ (10)

When K increases, the bias coefficient in (10) approaches $1/\left(1+\left(n-1\right){\rho }_e{R}_x\right)$, which again demonstrates the adverse effect of ${\rho }_e{R}_x$. The direction of the bias can be seen using the difference between the numerator and denominator of (10). Thus

$E\left(\mathrm{v}\hat{a}{r}_{\text{OLS}}\left({\hat{b}}_{\text{OLS}}\right)\right)<\mathrm{var}\left({\hat{b}}_{\text{OLS}}\right)\iff \left(N-2\right)R_x{\rho }_e+\left(1+R_x\right){\rho }_e>0$ (11)

When K increases (and thus N also increases), then the direction of bias is dependent on the sign of $R_x{\rho }_e$. The condition can be stated in a simpler form as $\left(N-1\right)R_x{\rho }_e+{\rho }_e>0$, but (11) is needed for comparison with [3] in which two inequalities are derived which taken together, imply that the F-test in OLS leads to p-values that are too small. In our case, conditions in [3] are $\left(n-1\right)\left(1+R_x\right){\rho }_e>0$ and $\left(n-1\right)R_x{\rho }_e>0$, which together imply the validity of the inequality in (11). Moreover, (11) can solve the direction of the bias in the case where sub-conditions indicate opposing directions. When $\mathrm{v}\hat{a}{r}_{\text{OLS}}\left({\hat{b}}_{\text{OLS}}\right)$ is used in a t-test, its bias also produces bias in the computed p-values. Reversing the inequality in (11), we obtain a condition for obtaining too large p-values.

4. BLUE for Several Regressors, Nonsingular V

Let us now assume a general linear model $y=Xb+e$, the first column of $X$ being $1_N$ and the other columns centered. It is here assumed that $V$ is nonsingular, i.e., $-1/\left(n-1\right)<{\rho }_e<1$. Let x and z denote two regressors. Both $\alpha$ and $\beta$ are inversely proportional to $\left(1-{\rho }_e\right)\left(1+\left(n-1\right){\rho }_e\right)$. Thus, $\alpha$ and $\beta$ increase when ${\rho }_e$ approaches $-1/\left(n-1\right)$ or 1, leading to singular $V_i^{}$ when either of these limits are obtained. Let us define $W_i=\left(1-{\rho }_e\right)\left(1+\left(n-1\right){\rho }_e\right)V_i^{-1}$. The diagonal elements of $W_i$ are equal to $\left(1+\left(n-2\right){\rho }_e\right)$, the nondiagonal elements are equal to $-{\rho }_e$. Let $W$ combine all blocks. Using $W$, we obtain a linear estimator:

${\hat{b}}_W={\left(X^{\prime }WX\right)}^{-1}{X}^{\prime }Wy$ (12)

Element $ij$ of $X^{\prime }W$ is for x:

$\begin{array}{c}{u}_{ij}=\left(1+\left(n-2\right){\rho }_e\right)\left(x_{ij}-\bar{x}\right)-{\rho }_e{\displaystyle \underset{j^{\prime }\ne j}{\overset{n}{\sum }}\left(x_{i{j}^{\prime }}-\bar{x}\right)}\\ =\left(1+\left(n-1\right){\rho }_e\right)\left(x_{ij}-{\bar{x}}_i\right)+\left(1-{\rho }_e\right)\left({\bar{x}}_i-\bar{x}\right)\end{array}$ (13)

The first diagonal element of $X^{\prime }WX$, corresponding to the intercept, is $N\left(1-{\rho }_e\right)$, and an element for x is:

$X^{\prime }W{X}_{xx}={\displaystyle \underset{i=1}{\overset{K}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}\left[\left(1+\left(n-1\right){\rho }_e\right){\left(x_{ij}-{\bar{x}}_i\right)}^2+\left(1-{\rho }_e\right){\left({\bar{x}}_i-\bar{x}\right)}^2\right]}}$ (14)

Equation (14) explains why (6) is an increasing function of $R_x{\rho }_e$ : with a large $R_x{\rho }_e$ GLS places large weighting on the component with small variation. The non-diagonal elements are zero for the first row and column, and generally:

$\begin{array}{c}{X}^{\prime }W{X}_{xz}={\displaystyle \underset{i=1}{\overset{K}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}\left(1+\left(n-1\right){\rho }_e\right)\left(x_{ij}-{\bar{x}}_i\right)\left(z_{ij}-{\bar{z}}_i\right)}}\\ +{\displaystyle \underset{i=1}{\overset{K}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}\left(1-{\rho }_e\right)\left({\bar{x}}_i-\bar{x}\right)\left({\bar{z}}_i-\bar{z}\right)}}\end{array}$ (15)

Note that on the second-row summation over n could be dropped by multiplying with n.

The first element of $X^{\prime }Wy$ is $N\left(1-{\rho }_e\right)\bar{y}$, and the others are:

$\begin{array}{c}{X}^{\prime }W{y}_x={\displaystyle \underset{i=1}{\overset{K}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}\left(1+\left(n-1\right){\rho }_e\right)\left(x_{ij}-{\bar{x}}_i\right)\left(y_{ij}-{\bar{y}}_i\right)}}\\ +{\displaystyle \underset{i=1}{\overset{K}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}\left(1-{\rho }_e\right)\left({\bar{x}}_i-\bar{x}\right)\left({\bar{y}}_i-\bar{y}\right)}}\end{array}$ (16)

The averages ${\bar{y}}_i$ and $\bar{y}$ do not contribute, but they may increase (or decrease) understanding. They do not indicate centering of y because they are not involved in the first element of $X^{\prime }Wy$. In GLS, the same estimates are obtained with any scaling of $V^{-1}$, thus ${\hat{b}}_W={\hat{b}}_{\text{GLS}}$ for a non-singular $V$, but may be computationally more stable. Variances are then computed using

$\mathrm{var}\left({\hat{b}}_W\right)=\left(\left(1-{\rho }_e\right)\left(1+1/\left(n-1\right)\right)\right){\left(X^{\prime }WX\right)}^{-1}$ (17).

5. BLUE When V is Singular

5.1. BLUE in the Simple Linear Regression with Singular V

Let us introduce the estimation with a singular $V$ using the simple linear model. $y_{ij}=a+b{x}_{ij}+e_{ij}$. If ${\rho }_e=1$ and $R_x\ne 1$ we can estimate b in the model with zero variance with $\hat{b}=\left(y_{ij}-y_{i{j}^{\prime }}\right)/\left(x_{ij}-x_{i{j}^{\prime }}\right)$ using two observations such that $x_{ij}\ne x_{i{j}^{\prime }}$. As ${\rho }_e=1$ implies that $e_{ij}={\bar{e}}_i$ for all observations, a can be estimated by computing the arithmetic mean (which is OLS estimator) of $y_{ij}-\hat{b}{x}_{ij}$ or $y_i-\hat{b}{\bar{x}}_i$, the latter computation implying correct variance estimator ${\sigma }^2/K$.

If. $R_x\ne 1$, we can find at least two observations in the same group having different x values. If there are more observations deviating from the group means, there is an infinite number of estimating equations which all produce this same estimate with zero variance.

If ${\rho }_e=1$ and $R_x=1$, then all y and x values are equal to the group means in all groups (for y with probability one), and b needs to be estimated with OLS using group means. This means that $\mathrm{var}\left(\hat{b}\right)$ is n times larger than we would get for ${\rho }_e=0$, as could be anticipated from Figure 1 for $n=3$. The same estimates are produced using observation level OLS regression, which produces biased estimates for the variances $\mathrm{var}\left(\hat{b}\right)$ and $\mathrm{var}\left(\hat{a}\right)$.

If ${\rho }_e=-1/\left(n-1\right)$ and $R_x\ne -1/\left(n-1\right)$, then we can find at least two groups i and $i^{\prime }$ so that ${\bar{x}}_i\ne {\bar{x}}_{i^{\prime }}$, and b can be estimated with $\left({\bar{y}}_i-{\bar{y}}_{i^{\prime }}\right)/\left({\bar{x}}_i-{\bar{x}}_{i^{\prime }}\right)$ with zero variance. If there are more than two differing group means, there is an infinite number of estimating equations producing same estimates. After estimating b, a can be estimated with the arithmetic mean of ${\bar{y}}_i-b{\bar{x}}_i$ with variance ${\sigma }^2/K$. The same estimate is obtained with the arithmetic mean of $y_{ij}-b{x}_{ij}$, but standard OLS computations would produce a biased estimate for the estimation variance.

If ${\rho }_e=-1/\left(n-1\right)$ and $R_x=-1/\left(n-1\right)$, i.e. all group means of ${\bar{e}}_i$ and ${\bar{x}}_i$ are equal (for ${\bar{e}}_i$ with probability one), then BLUE of b and a can be obtained by dropping one observation from each group and doing OLS regression in the remaining data. This means that $\mathrm{var}\left(\hat{b}\right)$ is equal to $n/\left(n-1\right)$ times the OLS variance obtained for ${\rho }_e=0$. For $n=3$, $n/\left(n-1\right)=1.5$ as could be anticipated from Figure 1.

5.2. BLUE in the Multiple Regression When ${\rho }_e=1$

When ${\rho }_e=1$, then $V$ is singular. Also $W$ is singular, all of its diagonal elements are equal to $n-1$, and nondiagonal elements are equal to −1. The estimator cannot be applied directly because the first diagonal element of $X^{\prime }WX$ is $N\left(1-{\rho }_e\right)=0$ thus producing singular $X^{\prime }WX$. According to (14) the diagonal elements are zero also for such predictors x having $R_x=1$, i.e. no variation in $x_{ij}-{\bar{x}}_i$. If all predictors have $R_x=1$ then OLSE is BLUE, because $C\left(VX\right)\subset C\left(X\right)$, see [4].

A general BLUE can be obtained using the decomposition $X=\left(\begin{array}{cc}{X}_1& X_2\end{array}\right)$ where $X_1$ contains $1_N$ and all columns having $R_x=1$ and $X_2$ contains other columns. It is now required that K ≥ 3. Let $b={\left(\begin{array}{cc}{b^{\prime }}_1& {b^{\prime }}_2\end{array}\right)}^{\prime }$ decompose $b$ similarly. If we weight columns of $X_1$ as done in OLS, i.e. replace columns of ${X^{\prime }}_{1}W$ in Equation (12) with $X_1$ we get an estimator

${\hat{b}}_{\text{GOLS}}=By={\left(\begin{array}{cc}{X^{\prime }}_1{X}_1& {X^{\prime }}_1{X}_2\\ {X^{\prime }}_{2}W{X}_1& {X^{\prime }}_{2}W{X}_2\end{array}\right)}^{-1}\left(\begin{array}{c}{X^{\prime }}_1\\ {X^{\prime }}_{2}W\end{array}\right)y$ (18)

I suggest that the estimator is called a GOLS estimator as it combines OLS and GLS estimation principles in the same estimator. Using the standard formula for computing the inverse of a partitioned matrix, the estimator can be presented as

${\hat{b}}_{\text{GOLS}}=By=\left(\begin{array}{c}{\left({X^{\prime }}_1{X}_1\right)}^{-1}{X^{\prime }}_1-{\left({X^{\prime }}_1{X}_1\right)}^{-1}{X^{\prime }}_1{X}_2{\left(X_2\text{'}W{X}_2\right)}^{-1}{X^{\prime }}_{2}W\\ {\left({X^{\prime }}_{2}W{X}_2\right)}^{-1}{X^{\prime }}_{2}W\end{array}\right)y$ (19)

Thus

${\hat{b}}_2={\left({X^{\prime }}_{2}W{X}_2\right)}^{-1}{X^{\prime }}_{2}Wy$ (20)

Thus, ${\hat{b}}_2$ is the GLSE of $b_2$ as estimated solely using $X_2$. Now $\mathrm{var}\left({\hat{b}}_2\right)={\left({X^{\prime }}_{2}W{X}_2\right)}^{-1}{X^{\prime }}_{2}WVW{X}_2{\left({X^{\prime }}_{2}W{X}_2\right)}^{-1}=0$. Noting that the second row of $B$ is part of the first row, we obtain:

${\hat{b}}_1={\left({X^{\prime }}_1{X}_1\right)}^{-1}{X^{\prime }}_{1}y-{\left({X^{\prime }}_1{X}_1\right)}^{-1}{X^{\prime }}_1{\hat{b}}_2$ (21)

${\hat{b}}_1={\hat{b}}_{1,\text{OLS}}-{\left({X^{\prime }}_1{X}_1\right)}^{-1}{X^{\prime }}_1{X}_2{\hat{b}}_2$ (22)

where ${\hat{b}}_{1,\text{OLS}}$ is the OLS estimator of $b_1$ obtained ignoring $X_2$.

Using (22), we get $\mathrm{var}\left({\hat{b}}_1\right)={\sigma }^{2}n{\left({X^{\prime }}_1{X}_1\right)}^{-1}$, i.e. the same variance what we get when $b_2=0$ and we regress group means of y on group means of all variables in $X_1$. This is natural because all y values and x values are equal in all observations in each group.

Denote that $M=I-X{\left(X^{\prime }X\right)}^{-1}{X}^{\prime }$. An estimator $Ay$ is BLUE for $Xb$, if and only if $A\left(\begin{array}{cc}X& VM\end{array}\right)=\left(\begin{array}{cc}X& 0\end{array}\right)$, see [5]. If $X$ has full column rank (as we have), then $b$ is estimable and $By$ is BLUE for $b$ if and only if $B\left(\begin{array}{cc}X& VM\end{array}\right)=\left(\begin{array}{cc}I& 0\end{array}\right)$. Then, that $BX=I$ can be seen directly from taking into account that ${X^{\prime }}_{2}W{X}_1=0$. Noting that $WV=0$ and ${X^{\prime }}_{1}V=n{X^{\prime }}_1$ we obtain that:

$BVM={\left(n{\left({X^{\prime }}_1{X}_1\right)}^{-1}{X^{\prime }}_1-n{\left({X^{\prime }}_1{X}_1\right)}^{-1}{X^{\prime }}_{1}X{\left(X^{\prime }X\right)}^{-1}{X}^{\prime },0\right)}^{\prime }$ (23)

Matrix $X{\left(X^{\prime }X\right)}^{-1}{X}^{\prime }$ is a projector to the space spanned by $X$. Now $C\left(n{\left({X^{\prime }}_1{X}_1\right)}^{-1}{X^{\prime }}_1\right)\subset C\left(X\right)$, thus $n{\left({X^{\prime }}_1{X}_1\right)}^{-1}{X^{\prime }}_{1}X{\left(X^{\prime }X\right)}^{-1}{X}^{\prime }=n{\left({X^{\prime }}_1{X}_1\right)}^{-1}{X^{\prime }}_1$, which completes the demonstration that $BVM=0$.

The GOLS estimator can be put into general matrix theory context using [5]. Now $C\left(X_1\right)$ and $C\left(X_2\right)$ are disjoint and $C\left(X_1\right)\subset C\left(V\right)$ and $C\left(X_2\right)\not\subset C\left(V\right)$. In such a case (20) provides BLUE for $b_2$. The estimator $X_1{\left({X^{\prime }}_1{X}_1\right)}^{-1}{X^{\prime }}_{1}y$ is an unbiased for $X_1{b}_1$. This unbiased estimator is updated into BLUE in (22), being an application of proposition 10.5 on p. 228 of [5].

5.3. BLUE When ${\rho }_e=-1/\left(n-1\right)$

If ${\rho }_e=-1/\left(n-1\right)$ then all elements of $W_i$ are equal to $1/\left(n-1\right)$. Let us organize $X$ into $X=\left(\begin{array}{cc}{X}_1& X_2\end{array}\right)$, where $X_1$ contains all such predictors for which $R_x=-1/\left(n-1\right)$ and $X_2$ contains all other columns. Equation provides again a BLUE with $\mathrm{var}\left({\hat{b}}_2\right)=0$ and the same arguments can be used to prove it also for this case. Now $\mathrm{var}\left({\hat{b}}_1\right)={\sigma }^2\frac{n}{n-1}{\left({X^{\prime }}_1{X}_1\right)}^{-1}$, which is the same variance we obtain if $b_2=0$ and $b_1$ is estimated from data where one redundant observation is dropped from each group.

It may a useful exercise to get confidence in multiple regressor derivations to show how the estimators derived in Section 5.1 for the simple linear model are produced also with the matrix formulas. Matrix formulas use in computations of all observations while formulas in section 5.1 utilize only nonredundant information.

6. Discussion

When $y_{ij}$ is regressed on $x_{ij}$, it is implicitly assumed that the regression of $y_{ij}$ on $x_{ij}-{\bar{x}}_i$, and the regression of ${\bar{y}}_i$ on ${\bar{x}}_i-\bar{x}$, have the same slope. GLS estimation provides different weights to $x_{ij}-{\bar{x}}_i$ and ${\bar{x}}_i-\bar{x}$ when attempting to utilize the correlation structure, although $\left(1+\left(n-1\right){\rho }_e\right)\left(x_{ij}-{\bar{x}}_i\right)$ and $\left(1-{\rho }_e\right)\left({\bar{x}}_i-\bar{x}\right)$ are put into the same column of $X$ as shown in (14). When developing mixed models for grouped data, it is often necessary to consider both $x_{ij}-{\bar{x}}_i$ and ${\bar{x}}_i$ as separate predictors, see [8]. It is natural to assume that $y_{ij}$ is related to $x_{ij}-{\bar{X}}_i$ and $x_{ij}-{\bar{X}}_i$, where ${\bar{X}}_i$ is the mean of x in the whole group. Reference [8] suggests a solution, based on a multivariate mixed model, which solves the measurement error bias problem that occurs when ${\bar{x}}_i$ computed from a sample from group i is used to “measure” ${\bar{X}}_i$.

The mean of variable $x_{ij}-{\bar{x}}_i$ is zero in each group, thereby indicating a negative within-group correlation. If $y_{ij}$ is linearly related to $x_{ij}-{\bar{x}}_i$, the negative within-group correlation of $x_{ij}-{\bar{x}}_i$ is also transmitted to $y_{ij}$. Hence, the addition of predictor $x_{ij}-{\bar{x}}_i$ to the models makes the estimated variance of the random intercept larger, as the “negative variance” is subtracted from the residual variance, see [8] and [9].

The mixed model formulation always leads to non-negative intra-class correlations. However, negative intra-class correlations are needed in situations where the group means have a smaller variance than that implied by the assumption of uncorrelated individuals. In the marginal interpretation of mixed model equations, a negative definite variance matrix of random effects that maximize the likelihood is allowed, leading to negative correlations and negative variances, provided the marginal variance of the y-vector is positive semi-definite, see [10]. The significance of a negative intra-class coefficient is dependent on the group size.

When the x-variable is a random variable with a theoretical intra-class correlation ${\rho }_x$, the empirical intra-class correlation $R_x$ is a random variable. The expected value of $R_x$ approaches ${\rho }_x$ when K increases. Simulations with normally distributed variables show that, for given values of n and K, $\mathrm{var}\left(R_x\right)$ can be well described with function $c{\left({\rho }_x+1/\left(n-1\right)\right)}^s{\left(1-{\rho }_x\right)}^q$. Negatively correlated variables can be simulated using principal components.

Reference [2] discussed the estimation of the slope in a simple linear model assuming that $n=2$, $x_{i1}=0$ and $x_{i2}=t_i$ for all i. To compare the formulas presented above with the formulas in [2], denote that $p=K{\bar{t}}^2/{\displaystyle {\sum }_{i=1}^K{t}_i^2}$ where $\bar{t}$ is the average of $t_i$, then $0<p\le 1$, and $p=1$ when all $t_i$ ’s values are equal. It holds that $R_x=-p/\left(2-p\right)$. They compare the variance of GLSE to the variance of the OLS estimate obtained by regressing $y_{i2}-y_{i1}$ on $t_i$. The standard OLSE of the slope has, however, a smaller variance than their OLSE based on differences when ${\rho }_e>0.5$. For ${\rho }_e=0$, the ratio of the variance of the standard OLSE to the variance of their OLSE is $2-p$. When all $t_i$ ’s are equal, then $p=1$, $R_x=-1$, and all three methods provide the same estimate.

If the group sizes $n_i$ vary in a data set, it is not reasonable to assume the same negative intra-class coefficient ${\rho }_e$ for all groups. Thus ${\rho }_e$ needs to be made group specific ${\rho }_{ei}$. Reference [7] generalizes the intra-class correlation to non-equal group sizes, but this generalization does not provide meaningful analysis for the interaction of the intra-class correlation of the residual and the values of regressors. However, it is possible to define a measure for the between-group variation of x which gets value of zero if there is no variation between groups, and use this measure to analyze how the intra-class correlation of residual errors and the between-group and within-group variation of predictors interact. This analysis can be extended to auto-correlation structures., and is under development.

In experiments, (6) informs us that $R_x$ should be $-1/\left(n-1\right)$ (all ${\bar{x}}_i$ should be equal) if it is anticipated that ${\rho }_e>0$. If it is anticipated that ${\rho }_e<0$, then the same treatment level should be applied to the whole group.

The derivations of this paper are fully covered with the general matrix theory. The purpose of the paper is to make general matrix formula understandable using algebraic derivations which show how between-group variation and within-group variations of residual errors are connected to the between-group and within-group variations of regressors in the estimation of a general linear model.

Acknowledgements

I thank Ronald Christensen for taking seriously the first version dealing with the simple linear model, Simo Puntanen for patient guidance to the conditions which prove that an estimator is BLUE, and Jarkko Isotalo for indicating how the GOLS estimator is related to the general BLUE theory. I take full responsibility if I have not understood their explanations correctly. They fly smoothly in linear spaces and subspaces where I need to walk watching each step.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Markiewics, A., Puntanen, S. and Styan, P.H. (2022) The Legend of the Equality of OLSE and BLUE: Highlighted by C. R. Rao in 1967. In: Arnold, B.C., Balakrishnan, N. and Coelho, C.A., Eds., Methodology and Applications of Statistics, Springer, Cham.
[2]	Jeske, D.R. and Myhre, J.M. (2018) Regression Using Pairs vs. Regression on Differences: A Real-life Case Study for a Master’s Level Methods Class, The American Statistician, 72, 163-168. https://doi.org/10.1080/00031305.2017.1292956
[3]	Christensen, R. (2003) Significantly insignificant F tests. The American Statistician. 57, 27-32. https://doi.org/10.1198/0003130031108
[4]	Puntanen, S. and Styan, G.P.H. (1989) The Equality of the Ordinary Least Squares Estimator and the Best Linear Unbiased Estimator. The American Statistician, 43, 153-161. https://doi.org/10.1080/00031305.1989.10475644
[5]	Puntanen, S., Styan, G.P.H. and Isotalo, J. (2011) Matrix Tricks for Linear Statistical Models: Our Personal Top Twenty. Springer, Heidelberg. https://doi.org/10.1007/978-3-642-10473-2
[6]	Christensen, R. (2011) Plane Answers to Complex Questions. The Theory of Linear Models. 4th Edition, Springer, New York. https://doi.org/10.1007/978-1-4419-9816-3
[7]	Kendall, M. and Stuart, A. (1979) The Advanced Theory of Statistics: Volume 2: Inference and Relationship. 4th Edition, Charles Griffin, London.
[8]	Mehtätalo, L. and Lappi, J. (2020) Biometry for Forestry and Environmental Data with Examples in R. Chapman and Hall/CRC, New York. https://doi.org/10.1201/9780429173462
[9]	Snijders, T. and Bosker, R.J. (1994) Modeled Variance in Two-Level Models. Sociological Methods Research, 22, 342-363. https://doi.org/10.1177/0049124194022003004
[10]	Pryseley, A., Tchonlafi, C., Verbeke, G. and Molenberghs, G. (2011) Estimating Negative Variance Components from Gaussian and Non-Gaussian Data: A Mixed Models Approach. Computational Statistics and Data Analysis, 55, 1071-1085. https://doi.org/10.1016/j.csda.2010.09.002