Weighted Standard Errors in R

Weighted Averages

Weighted averages show up often in data analysis. There are two reasons for this, one conceptual and one practical. The conceptual reason is that just about every estimand of interest is itself an average. Why is this? That pesky pesky fundamental problem of causal inference. Since causal effects are defined by comparing potential outcomes, and only one outcome can ever be observed for an individual, the problem of causal inference is a missing data problem.¹ As a consequence, we require multiple units to make any inference, and have to make do with an average instead of individual causal effects. ATE? Average! CATE? Conditional Average! LATE? Local Average! Everything is an average.

Normally, the set up for causal inference tends to make an implicit assumption that the causal effect for every unit is constant. Practically, we probably should never believe that to be true for anything interesting. The treatment effect should be heterogeneous, which means that the estimand of interest is going to be a weighted average.

The practical reason that weighted averages show up often is that the definition of an average itself is a definition of a weighted sum. Normally it is introduced like this:

\[ \mu = \frac{1}{n}\sum_i^n x_i \]

We can generalize that definition to deal with weights.

\[ \mu = \frac{\sum_i^nw_ix_i}{\sum_i^nw_i} \]

where the \(w_i\) indicates the weight for each unit. If the weights are all equal \(w_i = \frac{1}{n}\), and that formula reduces to the usual definition of an average.

In R by hand

weighted.average <- function(x, w){
    ## Sum of the weights 
    sum.w <- sum(w, na.rm = T)
    ## Sum of the weighted $x_i$ 
    xw <- sum(w*x, na.rm = T)
    
    ## Return the weighted average 
    return(xw/sum.w)
}

In action, we can see that with equal weights this reduces to the usual mean function.

x <- c(1,2,3,4)
w <- c(1,1,1,1)

## test against R's built in mean function
weighted.average(x, w) == mean(x)

## [1] TRUE

When the weights are unequal, this returns what we expect.

x2 <- c(1,2,3,4)
w2 <- c(2,1,2,3)

## Test against R's built in weighted.mean function
weighted.average(x2, w2) == weighted.mean(x2, w2)

## [1] TRUE

Weighted Variance

So far so good, and at least in R there are default functions to do these calculations for us.

The trouble arises when we are interested in the uncertainty of our calculation, as would occur in a sample of data. We would like to know what is the variance of a weighted average. To find our formula we can derive this from the definition of variance.

Define \(W = \sum_i^n w_i\) and \(\bar{x} = \frac{1}{W}\sum_i^nw_ix_i\) which is the weighted average. Then:

\[\begin{aligned} V[\bar{x}] &= V[\frac{1}{W}\sum_i^nw_ix_i]\\ V[\bar{x}] &= \sum_i^nV[\frac{w_i}{W}x_i] \\ V[\bar{x}] &= \sum_i^n\left(\frac{w_i}{W}\right)^2V[x_i] \\ V[\bar{x}] &= V[x]\sum_i^n\left(\frac{w_i}{W}\right)^2 \end{aligned}\]

What’s going on in the derivation. Line 1 is simply the definition of variance. In Line 2 we use a property that the variance of independent sums is equal to the sum of the variances. Line 3 is an application of the fact that \(V[aX] = a^2V[X]\). Line 4 is due to the fact that sum of the variances of \(x_i\) is just \(V[x]\) which can pulled out of the summation. Note that if all the weights are equal, the weights fraction reduces to \(\frac{1}{n}\). To see this, consider the example above with equal weights. Each \(\frac{w_i}{W} = \frac{1}{n}\) so the total sum is \(n\frac{1}{n^2}\) which reduces to \(\frac{1}{n}\) after cancelling terms.

Normally when we consider a weighted average, we are thinking of the case where we want to give more weight to some points rather than others. For example, if we take a survey of 100 people and weight by self-identified race, some points will get larger weight.² In my work, this is the dominant weighting scenario to worry about. If we use the usual standard error of the mean formula described in footnote 2 in this case we will end up with a very unstable estimate of the standard error of the mean as weights become uneven.

## Function to calculate standard error with weights 
## Incorrect in case where weights imply meaning and not
## uncertainty
bad_se <- function(x, w){
    n = length(x)
    numerator = sum(w*(x-weighted.average(x,w))^2)
    denominator = sum(w)
    correction = n/(n-1)
    var_x = (correction * (numerator/denominator))/n
    return(sqrt(var_x))
}

To see why, consider a case where we place a massive amount of weight on one observation in a sequence. Here we generate every integer from 1 to 100 and see what happens when we use the usual standard error.

x4 <- seq(1,100,1)
w4 <- c(rep(1,99),10000)
avg = round(weighted.average(x4,w4),2)
se_avg = round(bad_se(x4,w4),2)

The average is close to but not quite 100, which makes sense because the last value of x has a massive weight compared to every other value. The standard error though is functionally zero because the last weight is dominating every other value. What is happening here is that the effective degrees of freedom is very small. The last point has such a large weight that it dominates all the other points, and we are not correcting for this switch.³ Instead of doing what we expect, become larger because the weights are uneven, this estimate of the standard error of the mean goes to zero.

We should thus correct for this loss of degrees of freedom by getting the effective sample size and use it in our standard error calculation.⁴

\[ n_{eff} = \frac{(\sum_i^nw_i)^2}{\sum_i^nw_i^2} \]

Note that when all the weights are equal, the effective sample size is identical to the nominal sample size because the numerator is \(n^2\) and the denominator is \(n\). However, if the weights are dominated by a single point, then the effective sample size will converge to 1. To see this imagine that we have two points, and the second is 100 times as important as the first.

\[ n = \frac{(1+100)^2}{(1+10000)} \\ n = \frac{101^2}{10001} \\ n = 1.02 \]

By using the effective n, we will get an unbiased estimate of the standard error. Here’s how to write this function in R.

weighted.se.mean <- function(x, w, na.rm = T){
    ## Remove NAs 
    if (na.rm) {
      i <- !is.na(x)
        w <- w[i]
        x <- x[i]
    }
    
    ## Calculate effective N and correction factor
    n_eff <- (sum(w))^2/(sum(w^2))
    correction = n_eff/(n_eff-1)
    
    ## Get weighted variance 
    numerator = sum(w*(x-weighted.average(x,w))^2)
    denominator = sum(w)
    
    ## get weighted standard error of the mean 
    se_x = sqrt((correction * (numerator/denominator))/n_eff)
    return(se_x)
}

We can now apply it to our example case.

x4 <- seq(1,100,1)
w4 <- c(rep(1,99),10000)
correct_avg = round(weighted.average(x4,w4),3)
correct_se = round(weighted.se.mean(x4,w4),3)

Now we get an estimate of the standard error that is far more in line with our expectation.

As a coda, in case you do not want to roll your own function there is a weighted variance function in the Hmisc R package. This could be used with the appropriate corrections to get the same answer.