Generalized Raking for Survey Weighting
Tags: ruby, python, statistics
Published: 2021-06-29 15:28:44

In the world of surveys, it is very common that our acquired responses need to be weighted in order to achieve a sample that is _representative_ of some target population. This process of _weighting_ simply consists of assigning a _weight_ (a.k.a. _factor_) to each respondent, and calculating all survey results as a weighted sum of respondents.

For example, we might have surveyed 100 male respondents and 150 female respondents but were targeting a male / female ratio of 48% / 52%. In this simple case, we could achieve the target ratio by weighting the male responses by a factor of `0.48 / (100 / (100 + 150)) = 1.2` and weighting the female responses by `0.52 / (150 / (100 + 150) = 0.867`.
The technical term for this method of computing weights is _Post-Stratification_.

However, in a more complex scenario, where we have _many_ different measurable demographic targets, how can we determine weights for all the survey respondents?

## Raking

At Potloc, it is very common that our clients desire survey populations matching a lot of such targets. For example, we might have targets looking like this:
- 42% male
- 58% female
- 20% students
- 80% non-students
- 15% dog owners
- ...

In this setting, weights cannot be calculated using a simple ratio as in the male/female example shown above. Here, we instead need to rely on more involved algorithms, notably a process called **_Raking_**.

### Iterative Proportional Fitting

One common approach to solve the problem of finding good weights that will satisfy our demographic targets is _Iterative Proportional Fitting_. Typically in the industry, when the term "raking" is used it refers to this algorithm. In this method, weights for each respondents are computed **for a single target at a time** using Post-Stratification. By iteratively computing this for each target and repeating a few times, the weights end up converging to values that satisfy our targets.

Great! Problem solved!

...but what if we could do even better? 🤔

### Generalized Raking

Beyond satisfying the demographic targets, the most desirable property for the weights is that they should be as close as possible to _1_. Indeed, weights that are really large mean that those respondents' responses will count for a lot more than the "average" respondent in our survey results. For example, a respondent with a weight of 10 will count for 10 times more than the average respondent, and 100 times more than a respondent with weight 0.1 . Similarly, small weights mean that some responses will have very little impact on the final results.

Unfortunately, Iterative Proportional Fitting does nothing to encourage weights to be close to _1_, which leads to sub-optimal weights. This is where **_Generalized Raking_**, an algorithm introduced by _Deville et al._ (1992), comes into play.

> _Note:_ This is where we get into the more mathematical part of this blog post 🤓. Don't care about this part? No worries! Simply skip to the next section!

The authors of this paper formulated the weighting problem as a constrained optimization method where the objective is that the weights are as close to one as possible and where the constraint is that the targets are matched. Mathematically this looks like this:

{% katex %}
 \argmin_{w} \; G(w) \;\; \text{s.t.} \; X^T w = T \newline
 G(x) = x(\log(x) - 1) + 1
{% endkatex %}

where {% katex inline %}G(x){% endkatex %} is the _raking_ function which encourages weights to be close to _1_, {% katex inline %}w{% endkatex %} is the vector of weights, {% katex inline %}T{% endkatex %} is the vector of targets (in absolute numbers, not percentages) and {% katex inline %}X{% endkatex %} is the {% katex inline %}(\text{numRespondents} \times \text{numTargets}){% endkatex %} matrix of responses. The matrix {% katex inline %}X{% endkatex %} is binary where cells are filled with a '1' if the respondent belongs to the target category and '0' otherwise.

In other words, this is saying that we want to optimize the weights to be as close to 1 as possible while satisfying the target constraints. This is achieved by minimizing the function {% katex inline %}G(x){% endkatex %} which looks like this (notice the global minimum at {% katex inline %}x=1{% endkatex %}!):
![image](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/3nsmhkkonuiknwa4h8dj.png) 

## The Generalized Raking Algorithm

While it is possible to solve this optimization problem using general methods such as [Sequential Least Squares Programming](https://docs.scipy.org/doc/scipy/reference/optimize.minimize-slsqp.html), the authors of Generalized Raking have devised a more efficient and robust algorithm for this specific problem:

1. Initialize variables
  - A {% katex inline %}(\text{numRespondents} \times 1){% endkatex %} vector {% katex inline %}w{% endkatex %} to ones
  - A {% katex inline %}(\text{numTargets} \times 1){% endkatex %} vector {% katex inline %}\lambda{% endkatex %} to zeros
  - A {% katex inline %}(\text{numRespondents} \times \text{numRespondents}){% endkatex %} square matrix {% katex inline %}H{% endkatex %} to the Identity matrix
2. While the weights have not converged, repeat:
  1. {% katex inline %}\lambda = \lambda + (X^T H X)^{-1} (T - X^T w){% endkatex %}
  2. {% katex inline %}w = G^{-1}(X \lambda){% endkatex %}
  3. {% katex inline %}H = \text{diag}({G^{-1}}'(X \lambda)){% endkatex %}

Here, {% katex inline %}G^{-1}(x){% endkatex %} is the inverse of the derivative of the raking function, i.e. {% katex inline %}e^x{% endkatex %}.
{% katex inline %}{G^{-1}}'{% endkatex %} is its derivative, in this case also {% katex inline %}e^x{% endkatex %}.

The final value of {% katex inline %}w{% endkatex %} corresponds to the weighting factors we are looking for!

## Implementation

While there are many implementations of this algorithm in R, we were not able to find one in Ruby that could play well with our codebase and be easily maintainable.
We therefore decided to make our own and to share it here for anyone looking for something similar. We started by making an implementation in python with the popular `numpy` library:

```python
import numpy as np

def raking_inverse(x):
  return np.exp(x)

def d_raking_inverse(x):
  return np.exp(x)

def graking(X, T, max_steps=500, tolerance=1e-6):
  # Based on algo in (Deville et al., 1992) explained in detail on page 37 in
  # https://orca.cf.ac.uk/109727/1/2018daviesgpphd.pdf

  # Initialize variables - Step 1
  n, m = X.shape
  L = np.zeros(m) # Lagrange multipliers (lambda)
  w = np.ones(n) # Our weights (will get progressively updated)
  H = np.eye(n)
  success = False
  
  for step in range(max_steps):
    L += np.dot(np.linalg.pinv(np.dot(np.dot(X.T, H), X)), (T - np.dot(X.T, w))) # Step 2.1
    w = raking_inverse(np.dot(X, L)) # Step 2.2
    H = np.diag(d_raking_inverse(np.dot(X, L))) # Step 2.3
    
    # Termination condition:
    loss = np.max(np.abs(np.dot(X.T, w) - T) / T)
    if loss < tolerance:
        success = True
        break
          
  if not success: raise Exception("Did not converge")
  return w
```

### Ruby Implementation

After validating the algorithm in python, we then proceeded to replicate it in Ruby. For this, we had to find an equivalent to `numpy` which we found in [Numo](https://github.com/ruby-numo). Numo is an awesome library for vector and matrix operations, and its [`linalg`](https://github.com/ruby-numo/numo-linalg) sub-library was perfect for us as we needed to compute a matrix pseudo-inverse. This allowed us to translate the code to Ruby almost line by line:

```ruby
require "numo/narray"
require "numo/linalg"

def raking_inverse(x)
  Numo::NMath.exp(x)
end

def d_raking_inverse(x)
  Numo::NMath.exp(x)
end

def graking(X, T, max_steps: 500, tolerance: 1e-6)
  # Based on algo in (Deville et al., 1992) explained in detail on page 37 in
  # https://orca.cf.ac.uk/109727/1/2018daviesgpphd.pdf

  # Initialize variables - Step 1
  n, m = X.shape
  L = Numo::DFloat.zeros(m)
  w = Numo::DFloat.ones(n)
  H_diag = Numo::DFloat.ones(n)

  success = false

  max_steps.times do
    L += Numo::Linalg.pinv((X.transpose * H_diag).dot(X)).dot(T - X.transpose.dot(w)) # Step 2.1
    w = raking_inverse(X.dot(L)) # Step 2.2
    H_diag = d_raking_inverse(X.dot(L)) # Step 2.3

    # Termination condition:
    loss = ((T - X.transpose.dot(w)).abs / T).max
    if loss < tolerance
      success = true
      break
    end
  end

  raise StandardError, "Did not converged" unless success
  w
end
```

You may have noticed that the code doesn't quite match exactly the algorithm described above, notably steps 2.1 and 2.3. This is because we have found it to be vastly faster with Numo to store the sparse matrix {% katex inline %}H{% endkatex %} as a flat vector `h_matrix_diagonal` since it only contains values on the diagonal. As a result, the step of taking the product {% katex inline %}X^T H{% endkatex %} can be rewritten as `X.Transpose * h_matrix_diagonal`, making use of Numo's implicit broadcasting.

In practice, we optimize this code a bit further by exiting early whenever possible (for example if our loss becomes `NaN`) and by allowing to pass as input an initial value for the vector `lambdas` if we believe to have an initialisation value better than the default.

With these few lines of code, we are now able to support complex survey weighting scenarios while having all of our code in our beautiful Ruby monolith 🎉

_Interested in what we do at Potloc? Come join us! [We are hiring](https://jobs.lever.co/Potloc?team=Engineering) 🚀_

## References

- [Gareth Davies' excellent PhD Thesis on the subject](https://orca.cf.ac.uk/109727/1/2018daviesgpphd.pdf)
- [Deville et al.'s 1992 _Generalized Raking Procedures_ paper](https://www.jstor.org/stable/2290793)