Multivariate Normal and Linear Regression

Simple linear regression least squares

Simple linear regression models each response as a linear function of one regressor, Y_i=\beta_0+\beta_1x_i+\epsilon_i, where the errors have mean zero and common variance \sigma^2. The ordinary least squares (OLS) estimators minimize \sum_i(y_i-\beta_0-\beta_1x_i)^2. With S_{xx}=\sum_i(x_i-\bar x)^2 and S_{xy}=\sum_i(x_i-\bar x)(y_i-\bar y), the solution is \hat\beta_1=S_{xy}/S_{xx} and \hat\beta_0=\bar y-\hat\beta_1\bar x. Since \hat\beta_1=\sum_i c_iY_i with c_i=(x_i-\bar x)/S_{xx}, the weights satisfy \sum_i c_i=0, and this fact drives several results below.

The standard variances and covariance are \operatorname{Var}(\hat\beta_1)=\sigma^2/S_{xx}, \operatorname{Var}(\hat\beta_0)=\sigma^2\left(\frac1n+\frac{\bar x^2}{S_{xx}}\right), and \operatorname{Cov}(\hat\beta_0,\hat\beta_1)=-\bar x\sigma^2/S_{xx}. The covariance contains \sigma^2, so it is not free of the error variance. Adding a known constant c to an estimator shifts its mean but leaves its variance unchanged, so \operatorname{Var}(\hat\beta_0+c)=\operatorname{Var}(\hat\beta_0). Because the weights c_i sum to zero and the errors are uncorrelated with common variance, \operatorname{Cov}(\bar Y,\hat\beta_1)=\frac{\sigma^2}{n}\sum_i c_i=0.

When the points lie exactly on a line, the OLS fit reproduces that line, so the slope estimate is the common slope and the intercept is \bar y-\hat\beta_1\bar x. The variance of the fitted mean response at a point x_0 is \sigma^2\left[\frac1n+\frac{(x_0-\bar x)^2}{S_{xx}}\right], which grows as x_0 moves away from \bar x. For regression through the origin, y_i=\beta x_i+\epsilon_i, the estimator is \hat\beta=\sum x_iy_i/\sum x_i^2; for example, the points (1,3) and (2,5) give \hat\beta=13/5=2.6. If the response is transformed linearly, the slope estimate is transformed accordingly, and the estimate can be recovered by writing the new sums in terms of the old ones.

Unbiasedness is an important property. Under the model, E(\text{SSE})=(n-2)\sigma^2, so \text{SSE}/(n-2) is unbiased for \sigma^2. The regression sum of squares equals \hat\beta_1^2S_{xx}, and its expectation is \sigma^2+\beta_1^2S_{xx}, which equals \sigma^2 only when \beta_1=0. For testing H_0:\beta_1=0, the F statistic is F=\text{MSR}/\text{MSE} with \text{MSR}=S_{xy}^2/S_{xx} and \text{MSE}=\text{SSE}/(n-2), and it follows F(1,n-2) under H_0.

Under normal errors, the joint confidence region for (\beta_0,\beta_1) is the set of \beta for which a quadratic form in \hat\beta-\beta is bounded by a constant times an F(2,n-2) quantile, and its boundary is an ellipse.

Finally, the OLS estimators are not always the maximum likelihood estimators. If the errors are Laplace distributed, the MLE minimizes the sum of absolute deviations, and it always exists, while the least squares estimators exist and are unique when the design matrix has full column rank. For n\ge 2 with at least two distinct values of \log k, the columns \mathbf 1 and (\log k)_k are linearly independent.

Quadratic forms of normal vectors

A quadratic form in a random vector Y has the form Y^TAY, where A is a symmetric matrix. Its distribution is especially simple when Y\sim N_n(\mathbf 0,I_n). If A is idempotent, then Y^TAY\sim\chi^2_r with r=\operatorname{rank}(A)=\operatorname{tr}(A). For example, if A=\frac13\mathbf 1\mathbf 1^T in three dimensions, then A is idempotent of rank 1, so Y^TAY=\frac13\left(\sum_iY_i\right)^2\sim\chi^2_1.

For a general nonsingular covariance, if W\sim N_p(\mathbf 0,V) with V positive definite, then W^TV^{-1}W\sim\chi^2_p. This is why X-Y\sim N_p(\mathbf 0,2\Sigma) for i.i.d. X,Y\sim N_p(\mu,\Sigma) gives the scaled form (X-Y)^T(2\Sigma)^{-1}(X-Y)\sim\chi^2_p. The sample mean satisfies \bar X\sim N_p(\mu,\Sigma/n), so n(\bar X-\mu)^T\Sigma^{-1}(\bar X-\mu)\sim\chi^2_p. Note that the centering and the scale matrix must match the distribution exactly; a squared expression such as (X_1+X_2)^2 is not chi-square merely because the variance is known, unless the mean is subtracted.

For a \chi^2_2 variable W, the tail probability has the closed form P(W\ge w)=e^{-w/2}, so a probability event written as W\ge 10 gives e^{-5}. For X\sim N_{n}(\mathbf 0,I), X^TX\sim\chi^2_n and E(X^TX)=n.

Any linear map of a normal vector is normal: if Y\sim N_n(\mu,\Sigma), then AY\sim N(A\mu,A\Sigma A^T). A matrix such as (X^TX+\lambda I_p)^{-1}X^TY is linear in Y, so it is multivariate normal, but it is generally biased and its covariance is a sandwich form, not the inverse of the coefficient matrix.

Several common misconceptions need care. Each marginal component of a normal vector is normal, and every linear combination is normal. Components are independent only when they are uncorrelated jointly normal variables, which requires a diagonal covariance matrix. The sum of squares of components is chi-square only when the covariance is the identity; otherwise it is a weighted sum of independent chi-square variables.

Independence of quadratic forms and linear forms is governed by matrix products. For Y\sim N_n(\mathbf 0,I), Y^TAY and c^TY are independent when Ac=\mathbf 0, and two quadratic forms with idempotent matrices are independent when A_1A_2=0. Ratios of normal variables can have heavy tails: the ratio of two independent centred normals is standard Cauchy, and a standard Cauchy variable has E|C|=\infty, so an expectation of its absolute value does not exist.

Multiple linear regression

In multiple linear regression, Y=X\beta+\epsilon with p regressors, and the OLS estimator is \hat\beta=(X^TX)^{-1}X^TY. The coefficient of determination is R^2=\text{SSR}/\text{SST}=1-\text{SSE}/\text{SST}, where \text{SST}=\text{SSR}+\text{SSE}. Because R^2 never decreases when regressors are added, the adjusted coefficient penalizes the number of parameters: \bar R^2=1-\dfrac{\text{SSE}/(n-p-1)}{\text{SST}/(n-1)}=1-(1-R^2)\dfrac{n-1}{n-p-1}. For example, with n=25, p=4 and R^2=5/6, the adjusted value is 1-\frac16\cdot\frac{24}{20}=0.80.

The overall F test for H_0:\beta_1=\cdots=\beta_p=0 uses F=\dfrac{\text{MSR}}{\text{MSE}}=\dfrac{\text{SSR}/p}{\text{SSE}/(n-p-1)}, which in terms of R^2 is F=\dfrac{(n-p-1)}{p}\cdot\dfrac{R^2}{1-R^2}. Rejection occurs when F exceeds the upper \alpha point of F(p,n-p-1), and since F is increasing in R^2, the rejection set is an interval R^2>R^2_{\text{crit}}. For instance, with p=2 and error degrees of freedom 22, the boundary satisfies 11R^2/(1-R^2)=F_{0.05,2,22}. Using the same formula, a value of R^2 that gives F below the critical value cannot lead to rejection at that level.

The one-way ANOVA decomposition of the regression is the same test: \text{SST}=\text{SSR}+\text{SSE} with \text{SSR} having p degrees of freedom and \text{SSE} having n-p-1. The test statistic is F=\text{MSR}/\text{MSE}. For example, with n=52 and p=3, the error degrees of freedom are 48, and F=\text{MSR}/\text{MSE} is computed from the two sums of squares.

An important property concerns identifiability and unbiasedness. In a polynomial regression with an intercept, the coefficients of non-constant terms such as t and t^2 are identifiable and have unbiased least squares estimators, whatever the error mean is. The intercept and an additive constant in the errors enter the mean only through their sum, so neither can be estimated unbiasedly on its own.

Conditional distribution of multivariate normal

If \mathbf X=(\mathbf X_1,\mathbf X_2) is multivariate normal with mean (\mu_1,\mu_2) and covariance partitioned as \Sigma_{11},\Sigma_{12},\Sigma_{21},\Sigma_{22}, then the conditional distribution of \mathbf X_1 given \mathbf X_2=\mathbf x_2 is normal with mean \mu_1+\Sigma_{12}\Sigma_{22}^{-1}(\mathbf x_2-\mu_2) and covariance \Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}. The conditional covariance does not depend on the conditioning value. For example, if \Sigma_{12}=(1\;\,1) and \Sigma_{22} is a 2\times 2 matrix, the conditional variance is \Sigma_{11} minus a quadratic form in \Sigma_{12} that uses \Sigma_{22}^{-1}. Note that the conditional mean is affine in the conditioning value, and the coefficient is the matrix \Sigma_{12}\Sigma_{22}^{-1}.

In the bivariate case, E(Y\mid X=x)=\mu_Y+\dfrac{\operatorname{Cov}(X,Y)}{\operatorname{Var}(X)}(x-\mu_X) and \operatorname{Var}(Y\mid X)=\sigma_Y^2(1-\rho^2). The slope of E(Y\mid X=x) in x is \rho\sigma_Y/\sigma_X. Consequently, if both means are zero, the ratio E(Y\mid X=x)/E(X\mid Y=y) for suitable x,y gives \sigma_Y^2/\sigma_X^2, and the ratio of the two conditional variances is the same quantity \sigma_Y^2/\sigma_X^2. Given E(Y\mid X=x) at one point, the correlation can be solved for: if the conditional mean rises by 2\rho per unit of x, then the value at a given x determines \rho.

To compute probabilities, note that any linear combination a^T\mathbf X is normal with mean a^T\mu and variance a^T\Sigma a. Standardizing gives probabilities through \Phi; for example, P(|Z|<c)=1-2\Phi(-c) for standard normal Z.

Independence is also characterized by covariance. Two linear combinations c^T\mathbf X and d^T\mathbf X are independent if and only if c^T\Sigma d=0, since their joint distribution is normal. A residual-type statistic X_2-a^T\mathbf X is therefore independent of X_2 exactly when the covariance of the two terms is zero, which is a linear condition on the coefficient vector a.