Iterative proportional fitting


The iterative proportional fitting procedure is an iterative algorithm for estimating cell values of a contingency table such that the marginal totals remain fixed and the estimated table decomposes into an outer product.
IPF is a method that has been "re-invented" many times, e.g. G.U. Yule in 1912 in relation to standardizing cross-tabulations and Kruithof in 1937
in relation to telephone traffic, and expanded upon by Deming and Stephan in 1940, it has seen various extensions and related research. A rigorous proof of convergence by means of differential geometry is due to Fienberg. He interpreted the family of contingency tables of constant crossproduct ratios as a particular -dimensional manifold of constant interaction and showed that the IPFP is a fixed-point iteration on that manifold. Nevertheless, he assumed strictly positive observations. Generalization to tables with zero entries is still considered a hard and only partly solved problem.
An exhaustive treatment of the algorithm and its mathematical foundations can be found in the book of Bishop et al.. The first general proof of convergence, built on non-trivial measure theoretic theorems and entropy minimization, is due to Csiszár.
Relatively new results on convergence and error behavior have been published by Pukelsheim and Simeone
. They proved simple necessary and sufficient conditions for the convergence of the IPFP for arbitrary two-way tables by analysing an -error function.
Other general algorithms can be modified to yield the same limit as the IPFP, for instance the Newton–Raphson method and
the EM algorithm. In most cases, IPFP is preferred due to its computational speed, numerical stability and algebraic simplicity.

Algorithm 1 (classical IPF)

Given a two-way -table of counts, where the cell values are assumed to be Poisson or multinomially distributed, we wish to estimate a decomposition for all i and j such that is the maximum likelihood estimate of the expected values leaving the marginals and fixed. The assumption that the table factorizes in such a manner is known as the model of independence. Written in terms of a log-linear model, we can write this assumption as, where, and the interaction term vanishes, that is for all i and j.
Choose initial values , and for set
Notes:
Assume the same setting as in the classical IPFP.
Alternatively, we can estimate the row and column factors separately: Choose initial values, and for set
Setting, the two variants of the algorithm are mathematically equivalent.
Notes:
Obviously, the I-model is a particular case of the Q-model.

Algorithm 3 (RAS)

The Problem: Let be the initial matrix with nonnegative entries, a vector of specified
row marginals and a vector of column marginals. We wish to compute a matrix similar to M and consistent with the predefined marginals, meaning
and
Define the diagonalization operator, which produces a matrix with its input vector on the main diagonal and zero elsewhere. Then, for, set
where
Finally, we obtain

Discussion and comparison of the algorithms

Although RAS seems to be the solution of an entirely different problem, it is indeed identical to the classical IPFP. In practice,
one would not implement actual matrix multiplication, since diagonal matrices are involved. Reducing the operations to the necessary ones,
it can easily be seen that RAS does the same as IPFP. The vaguely demanded 'similarity' can be explained as follows: IPFP
maintains the crossproduct ratios, e.i.
since
This property is sometimes called structure conservation and directly leads to the geometrical interpretation of contingency tables and the proof of convergence in the seminal paper of Fienberg.
Nevertheless, direct factor estimation is under all circumstances the best way to deal with IPF: Whereas classical IPFP needs
elementary operations in each iteration step, factor estimation needs only
operations being at least one order in magnitude faster than classical IPFP.

Existence and uniqueness of MLEs

Necessary and sufficient conditions for the existence and uniqueness of MLEs are complicated in the general case, but sufficient conditions for 2-dimensional tables are simple:
If unique MLEs exist, IPFP exhibits linear convergence in the worst case, but exponential convergence has also been observed. If a direct estimator exists, IPFP converges after 2 iterations. If unique MLEs do not exist, IPFP converges toward the so-called extended MLEs by design, but convergence may be arbitrarily slow and often computationally infeasible.
If all observed values are strictly positive, existence and uniqueness of MLEs and therefore convergence is ensured.

Goodness of fit

Checking if the assumption of independence is adequate, one uses the Pearson X-squared statistic
or alternatively the likelihood-ratio test statistic
Both statistics are asymptotically -distributed, where is the number of degrees of freedom.
That is, if the p-values and are not too small, there is no indication to discard the hypothesis of independence.

Interpretation

If the rows correspond to different values of property A, and the columns correspond to different values of property B, and the hypothesis of independence is not discarded, the properties A and B are considered independent.

Example

Consider a table of observations :

Right-handedLeft-handedTOTAL
Male43952
Female44448
TOTAL8713100

For executing the classical IPFP, we first initialize the matrix with ones, leaving the marginals untouched:

Right-handedLeft-handedTOTAL
Male1152
Female1148
TOTAL8713100

Of course, the marginal sums do not correspond to the matrix anymore, but this is fixed in the next two iterations of IPFP. The first iteration deals with the row sums:

Right-handedLeft-handedTOTAL
Male262652
Female242448
TOTAL8713100

Note that, by definition, the row sums always constitute a perfect match after odd iterations, as do the column sums for even ones. The subsequent iteration updates the matrix column-wise:

Right-handedLeft-handedTOTAL
Male45.246.7652
Female41.766.2448
TOTAL8713100

Now, both row and column sums of the matrix match the given marginals again.
The variables are now independent, meaning the odds ratio is 1. This can be checked in either dimension: for both male and female, the odds of right-handed vs. left-handed are, since. Similarly, for both right-handed and left-handed, the odds of being male vs. female are, since.
For a 2×2 table, an exact solution is possible and iteration converges in a single pair of steps, and in fact a closed-form solution is to just take the outer product of the frequencies and divide by the population size, which yields the same values as above:

Right-handedLeft-handedTOTAL
Male87·52/10013·52/10052
Female87·48/10013·48/10048
TOTAL8713100

However for larger tables an exact solution is not always possible, and multiple iteration steps are necessary.
The p-value of this matrix approximates to, meaning: gender and left-handedness/right-handedness can be considered independent.

Implementation

The R package mipfp provides a multi-dimensional implementation of the traditional iterative proportional fitting procedure. The package allows the updating of a N-dimensional array with respect to given target marginal distributions.
Python has an equivalent package, ipfn that can be installed via pip. The package supports numpy and pandas input objects.