Lecture 30: Weighted regression

161.251 Regression Modelling

Introduction

So far we have changed the mean model: which predictors enter and how they enter. Weighted least squares changes the error model. Some responses are measured more precisely than others, so their deviations from the fitted mean need not contribute equally.

If observation \(i\) has error variance \(\sigma_i^2\), the usual precision weight is proportional to \(1/\sigma_i^2\). A more precise observation receives a larger weight. We will start with means based on different numbers of people, then use the same rule when individual error variances change across the predictor range.

Weighting by Sample Size

The CPS5 data contain individual wages and characteristics for 528 employees. Here we focus on education (ED) and hourly wage (WG). Goldberger first averages wages within each education level, from 6 to 18 years, then regresses those group means on education. Most employees had 12, 14 or 16 years of education.

Each row here is a mean hourly wage for one education level, not one person. Suppose individual wages within the education groups have a common variance \(\sigma^2\) and are independent. A mean based on \(n_j\) people then has variance \(\sigma^2/n_j\), so inverse-variance weighting gives \(w_j \propto n_j\). This explains why the group containing three people should have less influence than a group containing many more.

CPS5 grouped data

Download CPS5grouped.csv

## CPS5grouped = read.csv("CPS5grouped.csv", header = TRUE)
CPS5grouped

   EDgroup  AvWage  SEMean   StDev   N
1        6  4.4567 0.79805 1.38226   3
2        7  5.7700 0.83690 1.87136   5
3        8  5.9787 0.62779 2.43142  15
4        9  7.3317 1.29904 4.50000  12
5       10  7.3182 0.64444 2.65711  17
6       11  6.5844 0.64084 3.32992  27
7       12  7.8182 0.25206 3.72159 218
8       13  7.8351 0.66765 4.06114  37
9       14 11.0223 0.92157 6.89643  56
10      15 10.6738 1.39739 5.03837  13
11      16 10.8361 0.63981 5.35304  70
12      17 13.6150 1.42417 6.97699  24
13      18 13.5310 1.09328 6.08715  31

Grouped data scatterplot

ggplot(data = CPS5grouped, mapping = aes(x = EDgroup, y = AvWage)) + geom_point(mapping = aes(size = N)) +
    geom_smooth(method = "lm") + labs(x = "Education (years)", y = "Average wage ($)")

Grouped regression summary

The group means lie close to a line, and their unweighted regression has an \(R^2\) above 0.9. This describes variation among 13 group means. The individual-level \(R^2\) also includes variation among people within each education group, so the two values are not directly comparable.

grouped.lm <- lm(AvWage ~ EDgroup, data = CPS5grouped)
summary(grouped.lm)

Call:
lm(formula = AvWage ~ EDgroup, data = CPS5grouped)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.5637 -0.7350  0.1266  0.7158  1.3198 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.01445    0.87462  -0.017    0.987    
EDgroup      0.72410    0.06958  10.406 4.96e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9387 on 11 degrees of freedom
Multiple R-squared:  0.9078,    Adjusted R-squared:  0.8994 
F-statistic: 108.3 on 1 and 11 DF,  p-value: 4.958e-07

The problem with equal weighting

The first point represents only three people, yet the unweighted fit gives it the same influence as the much larger 12-, 14- and 16-year education groups. If we want each employee to contribute equally, we need to weight those group means by their sizes.

The individual data show more variation because they retain the differences between people within each education group. The grouped and ungrouped fits therefore describe different response variation; we should not read their \(R^2\) values as a direct comparison of predictive quality.

Ungrouped CPS5 data

What does the relationship look like when we use individual values in place of group means?

Download CPS5.csv

## CPS5 = read.csv("CPS5.csv", header = TRUE)
head(CPS5)

  ED SO BL HP FE MS EX UN   WG
1 10  0  0  0  0  1 27  0  9.0
2 12  0  0  0  0  1 20  0  5.5
3 12  0  0  0  1  0  4  0  3.8
4 12  0  0  0  1  1 29  0 10.5
5 12  0  0  0  0  1 40  1 15.0
6 16  0  0  0  1  1 27  0  9.0

Ungrouped scatterplot

ggplot(data = CPS5, mapping = aes(x = ED, y = WG)) + geom_point() + geom_smooth(method = "lm") +
    labs(x = "Education (years)", y = "Wage ($)")

`geom_smooth()` using formula = 'y ~ x'

Ungrouped regression summary

cps5.lm = lm(WG ~ ED, data = CPS5)
summary(cps5.lm)

Call:
lm(formula = WG ~ ED, data = CPS5)

Residuals:
   Min     1Q Median     3Q    Max 
-8.068 -3.163 -0.700  2.289 34.709 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.60468    1.10318  -1.455    0.146    
ED           0.81395    0.08281   9.829   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.733 on 526 degrees of freedom
Multiple R-squared:  0.1552,    Adjusted R-squared:  0.1536 
F-statistic:  96.6 on 1 and 526 DF,  p-value: < 2.2e-16

We see there is an outlier, which could be a typo (decimal point in the wrong place?) However the main focus here is that there is a lot more variation among the individual data than there is among the averaged data, so the \(R^2\) is much lower. Also the regression coefficients have different values.

Why Grouped Regression Results Differ

Formulas not examinable but conclusion is important.

To explain the difference between these two regressions, consider the regression of \(Y\) on \(X\) for all individuals, and the regression of group means \(\bar Y\) on \(\bar X\). In the CPS5 example, everyone within an education group has the same ED value, so \(x_{ij}=\bar x_j\). The equivalence below depends on that condition. In Least Squares we choose \(\beta_0\) and \(\beta_1\) to minimise
\[SS_{error} = \sum^n_{i=1}(y_i - \beta_0- \beta_1 x_i)^2 ~~~~ (1) \] However the regression for the averages treats all \(m\) groups as being of equal worth: \[ \sum^m_{j=1}(\bar y_j - \beta_0- \beta_1 \bar x_j)^2 ~~~~~~ (2) \]
That is, the computations will try just as hard to fit a value that arises from a small group of data as a value from a large group.

How are these two sums of squares related?

If we take the sum in (1) and re-order it into Education level groups we get \[SS_{error} = \sum^m_{j=1}\sum^{n_j}_{i=1}(y_{ij} - \beta_0- \beta_1 x_{ij})^2 ~~~~ (3) \]

Now suppose for a moment we pretend that all the \(y_{ij}\) values are the same as the group mean \(\bar y_j\). The sum on the right-hand side reduces to \[SS_{error} = \sum^m_{j=1} n_j(\bar y_{ j} - \beta_0 - \beta_1 \bar x_j)^2 ~~~~ (4) \]

So the difference between this sum and the sum (2) is that the squared residual is multiplied by \(n_j\), the number of people in that group. We call this multiplier the weight \(w_j\).

Essentially we interpret the weight as the amount of information represented by the datapoint. So here in this simple example “amount of information” is easily measured as the sample size \(w_j = n_j\) at each education level.

Error SS decomposition

But in fact we don’t have all the \(y_{ij}\) equal. Instead what we do is to add and subtract \(\bar y_j\) inside the parentheses in equation (3), and then expand the square: \[\begin{aligned}SS_{error} &= \sum^m_{j=1}\sum^{n_j}_{i=1}(y_{ij} -\bar y_j + \bar y_j - \beta_0- \beta_1 x_{ij})^2\\&= \sum^m_{j=1}\sum^{n_j}_{i=1}(y_{ij} -\bar y_j)^2 + \sum^m_{j=1}\sum^{n_j}_{i=1}(\bar y_j - \beta_0- \beta_1 \bar x_j)^2 ~~~~~(*)\\&= \sum^m_{j=1}\sum^{n_j}_{i=1}(y_{ij} -\bar y_j)^2 + \sum^m_{j=1} n_j~(\bar y_j - \beta_0- \beta_1 \bar x_j)^2\end{aligned}\]

(*) Equality happens because the cross-product term \(\sum^m_{j=1} \sum^{n_j}_{i=1} ~2~(y_{ij} -\bar y_j) (\bar y_j - \beta_0- \beta_1 \bar x_j) =0\).

Key insight: SS decomposition

The key fact is that the error SS splits into two parts:

\[SS_{error} = \sum^m_{j=1}\sum^{n_j}_{i=1}(y_{ij} -\bar y_j)^2 + \sum^m_{j=1} n_j~(\bar y_j - \beta_0- \beta_1 \bar x_j)^2\] - The first term is the “pure error” sum of squares. It records how much individual wages vary around their education-group means \(\bar y_j\). - The second term is the weighted sum of squares for fitting the group means, with each squared residual multiplied by the group size \(n_j\).

Because everyone in an education group has the same ED value, the within-group sum of squares does not depend on the regression coefficients. Minimising the individual-level sum of squares therefore gives the same intercept and slope as a regression of the group means weighted by group size. This equality of coefficient estimates depends on the predictor being constant within each group. It does not say that the two analyses have the same residual degrees of freedom, standard errors or R-squared: the grouped analysis has lost the within-group variation.

Weighted regression in R

To do a weighted regression we have to specify weights in the lm() command. Below we compare the weighted least squares (WLS) estimates from the averaged data to the ungrouped ordinary least squares (OLS) estimates from the raw data.

cps.wls = lm(AvWage ~ EDgroup, weights = N, data = CPS5grouped)
summary(cps.wls)

Call:
lm(formula = AvWage ~ EDgroup, data = CPS5grouped, weights = N)

Weighted Residuals:
   Min     1Q Median     3Q    Max 
-6.944 -3.971  2.699  4.151  9.217 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.60468    1.27442  -1.259    0.234    
EDgroup      0.81395    0.09567   8.508 3.62e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.467 on 11 degrees of freedom
Multiple R-squared:  0.8681,    Adjusted R-squared:  0.8561 
F-statistic: 72.39 on 1 and 11 DF,  p-value: 3.621e-06

OLS comparison

summary(cps5.lm)

Call:
lm(formula = WG ~ ED, data = CPS5)

Residuals:
   Min     1Q Median     3Q    Max 
-8.068 -3.163 -0.700  2.289 34.709 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.60468    1.10318  -1.455    0.146    
ED           0.81395    0.08281   9.829   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.733 on 526 degrees of freedom
Multiple R-squared:  0.1552,    Adjusted R-squared:  0.1536 
F-statistic:  96.6 on 1 and 526 DF,  p-value: < 2.2e-16

The weighted group-mean fit gives the same intercept and slope as the individual-level fit, apart from rounding in the supplied means. Education is constant within each group. The standard errors and \(R^2\) need not match: the grouped fit has only 13 rows and does not include the within-group variation in wages.

Comparing WLS vs OLS plots

The following plot shows the averaged data with the unweighted (dashed) and weighted (solid red) lines. The weighted line is closer to where more data are.

ggplot(data = CPS5grouped, mapping = aes(x = EDgroup, y = AvWage)) + geom_point(mapping = aes(size = N)) +
    geom_smooth(method = "lm", mapping = aes(weight = N), col = "red", se = FALSE) +
    geom_smooth(method = "lm", linetype = "dashed", se = FALSE) + labs(x = "Education (years)",
    y = "Average wage ($)")

Residuals for Weighted Regression

Clearly a plot of ordinary residuals could be misleading in the context of weighted least squares. The low weight points on the left are badly fitted and this could make them look like influential outliers.

augment(cps.wls, CPS5grouped) |>
    ggplot(mapping = aes(x = EDgroup, y = .resid)) + geom_point() + geom_hline(yintercept = 0,
    linetype = "dotted") + labs(y = "Simple residuals", x = "Education group (years)")

Pearson residuals

Instead we use weighted, or Pearson, residuals \(r_j = \sqrt{w_j}(y_j - \hat y_j)\).

augment(cps.wls, CPS5grouped) |>
    mutate(pearson = .resid * sqrt(N)) |>
    ggplot(mapping = aes(x = EDgroup, y = pearson)) + geom_point() + geom_hline(yintercept = 0,
    linetype = "dotted") + labs(y = "Pearson (weighted) residuals", x = "Education group (years)")

The weighted residuals are measured on the scale used by the weighted fit. A small weighted residual means the observation is close to its fitted value relative to the assumed variance; it is not a complete influence analysis and does not make a point scientifically unimportant. The ordinary and weighted residual panels have different scales, so compare their patterns rather than their vertical magnitudes.

Application: Weighting For Non-Constant Variance

Suppose we have data for which the variance is non-constant. So we want to give more weight to those points which have small variance (or are known more precisely) and give less weight to those points that have bigger variance (or are known less precisely).

This still fits with the idea that Weight represents the amount of information carried by each datapoint.

Supervisors and Workers

The following data is on the number of supervisors \(Y\) (e.g. team leaders) and the number of workers \(X\) at a number of industrial establishments.
The data are from the book by Chatterjee and Price “Regression Analysis by Example”

Download supervisors.csv

ggplot(supervisors, mapping = aes(x = X, y = Y)) + geom_point() + geom_smooth(method = "lm") +
    labs(x = "Number of workers (X)", y = "Number of supervisors (Y)")

OLS residual diagnostics

Note that the line is overshooting the cluster of points on the left.
The residual plots show curvature and heteroscedasticity.
The fourth graph shows a simple linear model is highly influenced by the points at the right, which have very large residuals.

sup.lm <- lm(Y ~ X, data = supervisors)
par(mfrow = c(1, 2))
plot(sup.lm, which = c(1, 5))

Heteroscedasticity pattern

In the residuals-versus-fitted plot, the horizontal axis shows fitted numbers of supervisors. Compare the vertical spread around fitted \(Y\) values of 50, 100 and 150. The spread increases as the fitted value rises. Since fitted \(Y\) also rises with \(X\), a reasonable working model is \(\mbox{sd}(\varepsilon_i) \propto x_i\), or \(\sigma_i=kx_i\) for some constant \(k\).

Now let us suppose we want to fit a straight-line model \[y_i = \beta_0 + \beta_1 x_i + \varepsilon_i ~ \] We would ordinarily solve this by minimising the \[SS_{error} = \sum( y_i - \beta_0 - \beta_1 x_i)^2\] where the components of the sum (the residuals \(y_i - \beta_0 - \beta_1 x_i\)) have constant standard deviation. Except now the standard deviation is not constant.

But what we could do is divide through by \(x_i\). The model becomes \[\begin{aligned}\frac{y_i}{x_i} &= \frac{\beta_0}{x_i} + \beta_1\frac{x_i}{x_i} + \frac{\varepsilon_i}{x_i}\\Y_i &=\beta_0 X_i + \beta_1 + \eta_i\end{aligned}\] where \(Y_i = \frac{y_i}{x_i}\), \(X_i = \frac{1}{x_i}\) and \(\eta_i = \frac{\varepsilon_i}{x_i}\).

Constant variance after transformation

The error term \(\eta_i\) now has constant standard deviation as

\[\mbox{sd}\left( \frac{\varepsilon_i}{x_i}\right) = \frac{ kx_i}{x_i} = k.\]

We now have a linear model with constant variance so can solve by least squares in the usual way.

We could do the transformation ourselves and then just use lm() normally, but another way to look at it is we’re minimising

\[\sum \left(\frac{y_i}{x_i} - \frac{\beta_0}{x_i} - \beta_1 \right)^2 ~~= ~~ \sum \frac{1}{x_i^2}\left(y_i - \beta_0 -\beta_1 x_i \right)^2\]

which is a weighted regression with \(w_i = 1/ x_i^2\). So the task with this transformation is the same as weighted least squares with these weights.

WLS for supervisors

sup.wls = lm(Y ~ X, data = supervisors, weights = 1/X^2)
summary(sup.wls)

Call:
lm(formula = Y ~ X, data = supervisors, weights = 1/X^2)

Weighted Residuals:
      Min        1Q    Median        3Q       Max 
-0.041477 -0.013852 -0.004998  0.024671  0.035427 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 3.803296   4.569745   0.832    0.413    
X           0.120990   0.008999  13.445 6.04e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.02266 on 25 degrees of freedom
Multiple R-squared:  0.8785,    Adjusted R-squared:  0.8737 
F-statistic: 180.8 on 1 and 25 DF,  p-value: 6.044e-13

WLS vs OLS supervisors plot

ggplot(data = supervisors, mapping = aes(x = X, y = Y)) + geom_point() + geom_smooth(method = "lm",
    mapping = aes(weight = 1/X^2), col = "red", se = FALSE) + geom_smooth(method = "lm",
    linetype = "dotted", se = FALSE) + labs(x = "Number of workers (X)", y = "Number of supervisors (Y)")

Weighted residuals plot

The weighted line sits closer to the low-\(X\) observations, where the assumed error variance is smaller and the weights are larger. The OLS line gives these points the same per-observation weight as the higher-\(X\) observations, despite their greater residual spread.

Let’s check the weighted residuals below. Their spread is more similar across the range of \(X\), although a little curvature may remain. This is more consistent with the proposed variance model; the plot cannot establish that the variance is exactly constant.

augment(sup.wls) |>
    mutate(pearson = .resid/X) |>
    ggplot(mapping = aes(x = X, y = pearson)) + geom_point() + geom_hline(yintercept = 0,
    linetype = "dotted") + labs(y = "Pearson (weighted) residuals", x = "Number of workers (X)")

Basic rule for Weighting

The common rule is inverse variance: \(w_i \propto 1/\operatorname{Var}(Y_i \mid X_i)\). In the CPS5 example each row was a group mean; for supervisors each row was one establishment. Those different rows led to different expressions for the weight.

Application: Meta Analysis CPS5grouped data

Meta-analysis combines estimates from several studies. In a simple inverse-variance meta-analysis, each study contributes an estimate with its own standard error. The estimate is what we weight, so its weight is proportional to \(1/\mathrm{SE}^2\). An estimate with a smaller standard error contributes more.

In the CPS5 grouped data, each row is a mean wage for N employees. The StDev column estimates the standard deviation of individual wages within that group. The variance of the group mean is therefore approximately \(\mathrm{StDev}_j^2/n_j\), giving the weight \(w_j=n_j/\mathrm{StDev}_j^2\). This uses both group size and the variation among employees in each group.

CPS5weighted <- CPS5grouped |> mutate(Wt = N/ StDev^2)
cps.wls2 <- lm(AvWage ~ EDgroup, weights=Wt, data=CPS5weighted)
summary(cps.wls2)

Call:
lm(formula = AvWage ~ EDgroup, data = CPS5weighted, weights = Wt)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-1.4367  0.1153  0.7399  1.0895  1.7058 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.26605    0.93427   0.285    0.781    
EDgroup      0.65602    0.07828   8.381 4.19e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.21 on 11 degrees of freedom
Multiple R-squared:  0.8646,    Adjusted R-squared:  0.8523 
F-statistic: 70.24 on 1 and 11 DF,  p-value: 4.186e-06

CPS5 re-weighted plot

ggplot(data = CPS5weighted, mapping = aes(x = EDgroup, y = AvWage)) + geom_point(mapping = aes(size = N)) +
    geom_smooth(method = "lm", mapping = aes(weight = Wt), col = "red", se = FALSE) +
    geom_smooth(method = "lm", mapping = aes(weight = N), linetype = "dotted", se = FALSE) +
    labs(x = "Education group (years)", y = "Average Wage ($)")

The standard deviations are much higher for the highly-educated groups, so that counteracts some of the effect of higher sample sizes.

Running Example: Weighted Summaries for Seoul Bike Data

The Seoul bike data also show why grouped summaries need care. If we average demand within temperature bands, some groups represent many more hourly observations than others.

seoul_hourly = read_csv("../data/seoul_bike_hourly.csv")

seoul_grouped = seoul_hourly |>
    mutate(temp_band = cut(temperature, breaks = c(-20, -10, 0, 10, 20, 30, 40),
        include.lowest = TRUE)) |>
    group_by(temp_band) |>
    summarise(avg_bike_count = mean(rented_bike_count), avg_temperature = mean(temperature),
        N = n(), .groups = "drop")

seoul_grouped

# A tibble: 6 × 4
  temp_band avg_bike_count avg_temperature     N
  <fct>              <dbl>           <dbl> <int>
1 [-20,-10]           141.          -12.7    228
2 (-10,0]             210.           -4.19  1226
3 (0,10]              425.            5.04  2143
4 (10,20]             765.           15.3   2264
5 (20,30]            1120.           24.5   2387
6 (30,40]            1106.           33.0    512

ggplot(seoul_grouped, aes(x = avg_temperature, y = avg_bike_count)) + geom_point(aes(size = N)) +
    geom_smooth(method = "lm", se = FALSE, linetype = "dashed") + geom_smooth(method = "lm",
    aes(weight = N), se = FALSE, colour = "red") + labs(x = "Average temperature within band",
    y = "Average bike demand", title = "Aggregated groups with different support should not be treated equally")

These temperature bands summarise unequal numbers of hours. Weighting a band’s mean by its count changes how much that summary influences the line. Here temperature also varies within each band, so the line through weighted band means need not be identical to a line fitted to all the hourly observations. The exact equality in the education example relied on everybody in a group having the same predictor value.

Putting the weights together

We have now used four forms of weight. In every case, start with the variance of the response represented by one row, then take its reciprocal.

What one row represents	Variance of that row	Relative weight	Example in this lecture
A mean of \(n_j\) independent people, with a common individual variance \(\sigma^2\)	\(\sigma^2/n_j\)	\(n_j\)	First CPS5 grouped fit, weights=N
An observation or estimate with variance \(v_i\)	\(v_i\)	\(1/v_i\)	General rule; study estimates in meta-analysis
A mean of \(n_j\) independent people, with individual variance \(\sigma_j^2\) in group \(j\)	\(\sigma_j^2/n_j\)	\(n_j/\sigma_j^2\)	CPS5 re-weighted fit, weights=N/StDev^2
An observation whose error standard deviation is proportional to \(x_i\)	\(k^2x_i^2\)	\(1/x_i^2\)	Supervisors and workers, weights=1/X^2

The extra \(n_j\) in \(n_j/\sigma_j^2\) comes from averaging: the variance of a mean of \(n_j\) independent observations is the individual variance divided by \(n_j\). Likewise, \(1/x_i^2\) follows from inverse variance when the error standard deviation is \(kx_i\). A constant shared by every row, such as \(\sigma^2\) or \(k^2\), can be left out of the relative weights.