Weighted Gaussian Regression: A Worked Analysis

From a Weighted Model to an Analysis

A weighted Gaussian regression is useful when observations have different known relative precisions. For example, a clinical mean based on many measurements may be more precise than one based on only a few. A weight records this difference in precision: larger weights mean smaller error variance, not larger expected responses. In a complete analysis, the weights affect the fitted coefficients and their uncertainty, while the residuals help estimate the common dispersion scale.

This tutorial follows one analysis from its model specification through estimation and reporting. The design is treated as fixed, the observations are independent, and the precision weights are known. These conditions matter: the standard errors and intervals below depend on the stated variance model.

Model: For observations $i=1,\ldots,n$, suppose $$ Y_i=\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta}+\varepsilon_i,\qquad \varepsilon_i\ \text{independent},\qquad \varepsilon_i\sim N\left(0,\frac{\phi}{w_i}\right), $$ where $w_i>0$ are known precision weights, $\phi>0$ is the unknown dispersion, and $\mathbf{x}_i^\mathsf{T}$ is row $i$ of the design matrix $\mathbf{X}$. Equivalently, $E[Y_i]=\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta}$ and $\operatorname{Var}(Y_i)=\phi/w_i$.

Write $\mathbf{y}=(Y_1,\ldots,Y_n)^\mathsf{T}$, let $\mathbf{W}$ be the diagonal matrix with entries $w_1,\ldots,w_n$, and let $p$ be the number of columns of $\mathbf{X}$. We assume $\mathbf{X}$ has full column rank, so $n\ge p$ and $\mathbf{X}^\mathsf{T}\mathbf{W}\mathbf{X}$ is invertible. The weighted residual sum of squares at $\boldsymbol{\beta}$ is $$ Q(\boldsymbol{\beta})=(\mathbf{y}-\mathbf{X}\boldsymbol{\beta})^\mathsf{T}\mathbf{W}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}). $$

Worked Example: Specifying the Precision-Weight Model

Suppose five independent clinical measurements follow a straight-line mean model, with predictor values $x=(0,1,2,3,4)$ and weights $w=(1,2,1,2,1)$. The weights are relative precisions. In particular, an observation with weight 2 has half the error variance of one with weight 1, all else equal.

Use an intercept and slope, so $\mathbf{x}_i^\mathsf{T}=(1,x_i)$ and $\boldsymbol{\beta}=(\beta_0,\beta_1)^\mathsf{T}$. If, for illustration, $\phi=2$, the five variances are $$ \frac{2}{1}=2,\quad \frac{2}{2}=1,\quad \frac{2}{1}=2,\quad \frac{2}{2}=1,\quad \frac{2}{1}=2. $$ Thus the weights do not change the assumed mean line; they specify how the observations differ in variance around that line. The value $\phi=2$ here only illustrates the variance interpretation. In the analysis below, dispersion will be estimated from the data.

Solving the Weighted Normal Equations

For Gaussian errors with the stated variances, maximizing the likelihood over $\boldsymbol{\beta}$ is equivalent to minimizing $Q(\boldsymbol{\beta})$. The minimizer is the weighted least-squares estimate. The following derivation shows why the normal equations must hold and why their solution is unique under the full-rank condition.

Result: If all $w_i>0$ and $\mathbf{X}$ has full column rank, the unique minimizer of $Q(\boldsymbol{\beta})$ is $$ \hat{\boldsymbol{\beta}}=(\mathbf{X}^\mathsf{T}\mathbf{W}\mathbf{X})^{-1}\mathbf{X}^\mathsf{T}\mathbf{W}\mathbf{y}. $$ It solves the weighted normal equations $\mathbf{X}^\mathsf{T}\mathbf{W}\mathbf{X}\hat{\boldsymbol{\beta}}=\mathbf{X}^\mathsf{T}\mathbf{W}\mathbf{y}$.

Derivation. Expanding the quadratic gives $$ Q(\boldsymbol{\beta})=\mathbf{y}^\mathsf{T}\mathbf{W}\mathbf{y} -2\boldsymbol{\beta}^\mathsf{T}\mathbf{X}^\mathsf{T}\mathbf{W}\mathbf{y} +\boldsymbol{\beta}^\mathsf{T}\mathbf{X}^\mathsf{T}\mathbf{W}\mathbf{X}\boldsymbol{\beta}. $$ Its gradient with respect to $\boldsymbol{\beta}$ is $-2\mathbf{X}^\mathsf{T}\mathbf{W}\mathbf{y}+2\mathbf{X}^\mathsf{T}\mathbf{W}\mathbf{X}\boldsymbol{\beta}$. Setting that gradient to zero gives the normal equations. For any nonzero vector $\mathbf{a}$, $$ \mathbf{a}^\mathsf{T}\mathbf{X}^\mathsf{T}\mathbf{W}\mathbf{X}\mathbf{a} =(\mathbf{X}\mathbf{a})^\mathsf{T}\mathbf{W}(\mathbf{X}\mathbf{a}) =\sum_{i=1}^n w_i(\mathbf{x}_i^\mathsf{T}\mathbf{a})^2>0, $$ because full column rank implies $\mathbf{X}\mathbf{a}\ne\mathbf{0}$ and each $w_i>0$. Hence $\mathbf{X}^\mathsf{T}\mathbf{W}\mathbf{X}$ is positive definite, and $Q$ is strictly convex. The normal-equation solution therefore exists, is unique, and is the minimizer. $\square$

Worked Example: Computing the Weighted Regression Coefficients

Continue with $x=(0,1,2,3,4)$, $w=(1,2,1,2,1)$, and observed responses $y=(2,3,5,4,7)$. For an intercept and slope, the weighted normal equations use the four sums below: $$ \sum_i w_i=7,\qquad \sum_i w_i x_i=14,\qquad \sum_i w_i x_i^2=40,\qquad \sum_i w_i y_i=28,\qquad \sum_i w_i x_i y_i=68. $$ For example, $\sum_i w_ix_i y_i=1(0)(2)+2(1)(3)+1(2)(5)+2(3)(4)+1(4)(7)=68$.

The normal equations are $$ \begin{pmatrix}7&14\\14&40\end{pmatrix} \begin{pmatrix}\hat{\beta}_0\\\hat{\beta}_1\end{pmatrix} = \begin{pmatrix}28\\68\end{pmatrix}. $$ The determinant is $7(40)-14(14)=84$, so the matrix is invertible. Its inverse and the resulting estimate are $$ \begin{pmatrix}7&14\\14&40\end{pmatrix}^{-1} =\frac{1}{84}\begin{pmatrix}40&-14\\-14&7\end{pmatrix}, \qquad \begin{pmatrix}\hat{\beta}_0\\\hat{\beta}_1\end{pmatrix} =\frac{1}{84}\begin{pmatrix}40&-14\\-14&7\end{pmatrix} \begin{pmatrix}28\\68\end{pmatrix} =\begin{pmatrix}2\\1\end{pmatrix}. $$ The fitted mean is therefore $\hat{\mu}(x)=2+x$. Checking directly, the fitted values are $2,3,4,5,6$; the residuals are $0,0,1,-1,1$. Substituting the estimate back into the normal equations also checks the solution: $7(2)+14(1)=28$ and $14(2)+40(1)=68$.

Estimating Dispersion and Quantifying Uncertainty

The weighted residual sum of squares measures the remaining discrepancies on the scale specified by the weights. For the example, it is $0^2+2(0^2)+1(1^2)+2((-1)^2)+1(1^2)=4$. Because two regression coefficients were fitted to five observations, the residual degrees of freedom are $n-p=3$. The usual unbiased estimator of dispersion is therefore $\hat{\phi}=4/3$, rather than the maximum-likelihood estimate $4/5$.

Result: Under the independent Gaussian model with fixed design, known positive weights, and full column rank, the coefficient estimate has covariance $$ \operatorname{Cov}(\hat{\boldsymbol{\beta}})=\phi(\mathbf{X}^\mathsf{T}\mathbf{W}\mathbf{X})^{-1}. $$ Also, the residual sum of squares divided by $\phi$ has a chi-squared distribution with $n-p$ degrees of freedom and is independent of $\hat{\boldsymbol{\beta}}$. Consequently, for any fixed vector $\mathbf{c}$ with $\mathbf{c}^\mathsf{T}\boldsymbol{\beta}$ as the contrast of interest, the statistic $$ \frac{\mathbf{c}^\mathsf{T}\hat{\boldsymbol{\beta}}-\mathbf{c}^\mathsf{T}\boldsymbol{\beta}} {\sqrt{\hat{\phi}\,\mathbf{c}^\mathsf{T}(\mathbf{X}^\mathsf{T}\mathbf{W}\mathbf{X})^{-1}\mathbf{c}}} $$ has a $t$ distribution with $n-p$ degrees of freedom.

Justification. Let $\mathbf{W}^{1/2}$ be the diagonal matrix with entries $\sqrt{w_i}$, and define $\mathbf{Z}=\mathbf{W}^{1/2}\mathbf{X}$. Then $\mathbf{W}^{1/2}\mathbf{y}=\mathbf{Z}\boldsymbol{\beta}+\mathbf{u}$, where $\mathbf{u}=\mathbf{W}^{1/2}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta})\sim N(\mathbf{0},\phi\mathbf{I})$. The matrix $\mathbf{P}=\mathbf{Z}(\mathbf{Z}^\mathsf{T}\mathbf{Z})^{-1}\mathbf{Z}^\mathsf{T}$ is the orthogonal projection onto the $p$-dimensional column space of $\mathbf{Z}$. The fitted coefficient error is a linear function of $\mathbf{P}\mathbf{u}$; the weighted residual vector is $(\mathbf{I}-\mathbf{P})\mathbf{u}$. These are independent because their covariance is $\phi\mathbf{P}(\mathbf{I}-\mathbf{P})=\mathbf{0}$, and jointly Gaussian. The residual sum of squares is $\mathbf{u}^\mathsf{T}(\mathbf{I}-\mathbf{P})\mathbf{u}$. The matrix $\mathbf{I}-\mathbf{P}$ is an orthogonal projection of rank $n-p$, so this quadratic form divided by $\phi$ is chi-squared with $n-p$ degrees of freedom. Finally, a coefficient contrast standardized using the true $\phi$ is standard normal; dividing it by the square root of an independent chi-squared variable divided by its degrees of freedom gives the stated $t$ statistic. $\square$

Worked Example: Reporting a Contrast and a Mean-Response Interval

For the fitted example, $n-p=3$, the weighted residual sum of squares is 4, and hence $\hat{\phi}=4/3$. The inverse information matrix is $$ (\mathbf{X}^\mathsf{T}\mathbf{W}\mathbf{X})^{-1} =\frac{1}{84}\begin{pmatrix}40&-14\\-14&7\end{pmatrix}. $$ Consider first the coefficient contrast $\beta_0+2\beta_1$, represented by $\mathbf{c}=(1,2)^\mathsf{T}$. Its estimate is $2+2(1)=4$. Its estimated variance is $$ \hat{\phi}\,\mathbf{c}^\mathsf{T}(\mathbf{X}^\mathsf{T}\mathbf{W}\mathbf{X})^{-1}\mathbf{c} =\frac{4}{3}\cdot\frac{40-56+28}{84} =\frac{4}{21}, $$ so its standard error is $\sqrt{4/21}\approx0.436$. With $t_{0.975,3}\approx3.182$, the 95% confidence interval is $$ 4\ \pm\ 3.182\sqrt{\frac{4}{21}} \quad\Longrightarrow\quad (2.61,\ 5.39). $$

Now consider the mean response at $x_0=2.5$, whose design vector is $\mathbf{x}_0=(1,2.5)^\mathsf{T}$. The estimated mean is $\mathbf{x}_0^\mathsf{T}\hat{\boldsymbol{\beta}}=2+2.5=4.5$. The estimated variance of this fitted mean is $$ \hat{\phi}\,\mathbf{x}_0^\mathsf{T}(\mathbf{X}^\mathsf{T}\mathbf{W}\mathbf{X})^{-1}\mathbf{x}_0 =\frac{4}{3}\cdot\frac{40-70+43.75}{84} =\frac{55}{252}, $$ giving standard error $\sqrt{55/252}\approx0.467$. Its 95% confidence interval is $$ 4.5\ \pm\ 3.182\sqrt{\frac{55}{252}} \quad\Longrightarrow\quad (3.01,\ 5.99). $$ This interval estimates the mean response at $x_0$, not the range of an individual future observation. A future response also includes its own Gaussian error, so its prediction interval has an additional variance contribution.

How to Report the Analysis

A concise report should identify the response, the predictor, and the weight interpretation; give the fitted equation; and state how dispersion and uncertainty were calculated. In this example, the weighted fit is $\hat{\mu}(x)=2+x$, with residual dispersion estimate $4/3$ on 3 residual degrees of freedom. The estimated contrast $\beta_0+2\beta_1$ is 4, with a 95% interval from 2.61 to 5.39. The estimated mean at $x=2.5$ is 4.5, with a 95% interval from 3.01 to 5.99. The intervals use the $t$ reference distribution and the stated Gaussian model.

A useful reporting habit is to keep three quantities distinct: the estimated mean, the estimated uncertainty in that mean, and the residual dispersion. The fitted mean is determined by the weighted regression; its standard error depends on the design and weights as well as $\hat{\phi}$. Dispersion describes the residual scale after accounting for the fitted coefficients.

Common Mistake: Treating Weights as Replication Counts

Wrong reasoning: “A weight of 2 means this response should count as two independent observations, so I will duplicate the row.”
Correction: In this model, a weight of 2 specifies variance $\phi/2$, relative to $\phi$ at weight 1. It does not by itself say that two independent responses were observed. Weighting changes the contribution to the fit and the model-based variance through $\mathbf{W}$; it does not automatically change the number of observations or residual degrees of freedom.
Another common error: using the unweighted residual sum of squares or dividing the weighted residual sum of squares by $n$ to estimate dispersion. Use $\sum_i w_i(Y_i-\mathbf{x}_i^\mathsf{T}\hat{\boldsymbol{\beta}})^2$, and for the unbiased estimate divide by $n-p$.

Key takeaway: In a weighted Gaussian regression with known relative precisions, solve the weighted normal equations, estimate dispersion from weighted residuals using $n-p$ degrees of freedom, and use the resulting covariance formula for coefficient contrasts and mean responses. The interval must match the target: a mean response is not a future individual response.

Check Your Understanding

Try these questions using the model assumptions and calculations above.

If two observations have weights 1 and 4, what is the ratio of their error variances?
For a design with six observations and three fitted coefficients, what residual degrees of freedom are used in the unbiased dispersion estimate?
Why does full column rank of the design matrix make the weighted least-squares solution unique when all weights are positive?
What is the estimated dispersion in the worked data, and which residual sum of squares is used to calculate it?
Explain why an interval for a mean response at a predictor value is generally narrower than a prediction interval for one future response at that value.