Iterative proportional fitting

Given a two-way (I × J)-table of counts ${\displaystyle (x_{ij})}$, where the cell values are assumed to be Poisson or multinomially distributed, we wish to estimate a decomposition ${\displaystyle {\hat {m}}_{ij}=a_{i}b_{j}}$ for all i and j such that ${\displaystyle ({\hat {m}}_{ij})}$ is the maximum likelihood estimate (MLE) of the expected values ${\displaystyle (m_{ij})}$ leaving the marginals ${\displaystyle \textstyle x_{i+}=\sum _{j}x_{ij}\,}$ and ${\displaystyle \textstyle x_{+j}=\sum _{i}x_{ij}\,}$ fixed. The assumption that the table factorizes in such a manner is known as the model of independence (I-model). Written in terms of a log-linear model, we can write this assumption as ${\displaystyle \log \ m_{ij}=u+v_{i}+w_{j}+z_{ij}}$, where ${\displaystyle m_{ij}:=\mathbb {E} (x_{ij})}$, ${\displaystyle \sum _{i}v_{i}=\sum _{j}w_{j}=0}$ and the interaction term vanishes, that is ${\displaystyle z_{ij}=0}$ for all i and j.

Choose initial values ${\displaystyle {\hat {m}}_{ij}^{(0)}:=1}$ (different choices of initial values may lead to changes in convergence behavior), and for ${\displaystyle \eta \geq 1}$ set

${\displaystyle {\hat {m}}_{ij}^{(2\eta -1)}={\frac {{\hat {m}}_{ij}^{(2\eta -2)}x_{i+}}{\sum _{k=1}^{J}{\hat {m}}_{ik}^{(2\eta -2)}}}}$

${\displaystyle {\hat {m}}_{ij}^{(2\eta )}={\frac {{\hat {m}}_{ij}^{(2\eta -1)}x_{+j}}{\sum _{k=1}^{I}{\hat {m}}_{kj}^{(2\eta -1)}}}.}$

Notes:

• Convergence does not depend on the actual distribution. Distributional assumptions are necessary for inferring that the limit ${\displaystyle ({\hat {m}}_{ij}):=\lim _{\eta \rightarrow \infty }({\hat {m}}_{ij}^{(\eta )})}$ is an MLE indeed.
• IPFP can be manipulated to generate any positive marginals by replacing ${\displaystyle x_{i+}}$ by the desired row marginal ${\displaystyle u_{i}}$ (analogously for the column marginals).
• IPFP can be extended to fit the model of quasi-independence (Q-model), where ${\displaystyle m_{ij}=0}$ is known a priori for ${\displaystyle (i,j)\in S}$. Only the initial values have to be changed: Set ${\displaystyle {\hat {m}}_{ij}^{(0)}=0}$ if ${\displaystyle (i,j)\in S}$ and 1 otherwise.

Assume the same setting as in the classical IPFP. Alternatively, we can estimate the row and column factors separately: Choose initial values ${\displaystyle {\hat {b}}_{j}^{(0)}:=1}$, and for ${\displaystyle \eta \geq 1}$ set

${\displaystyle {\hat {a}}_{i}^{(\eta )}={\frac {x_{i+}}{\sum _{j}{\hat {b}}_{j}^{(\eta -1)}}},}$

${\displaystyle {\hat {b}}_{j}^{(\eta )}={\frac {x_{+j}}{\sum _{i}{\hat {a}}_{i}^{(\eta )}}}}$

Setting ${\displaystyle {\hat {m}}_{ij}^{(2\eta )}={\hat {a}}_{i}^{(\eta )}{\hat {b}}_{j}^{(\eta )}}$, the two variants of the algorithm are mathematically equivalent (can be seen by formal induction).

Notes:

• In matrix notation, we can write ${\displaystyle ({\hat {m}}_{ij})={\hat {a}}{\hat {b}}^{T}}$, where ${\displaystyle {\hat {a}}=({\hat {a}}_{1},\ldots ,{\hat {a}}_{I})^{T}=\lim _{\eta \rightarrow \infty }{\hat {a}}^{(\eta )}}$ and ${\displaystyle {\hat {b}}=({\hat {b}}_{1},\ldots ,{\hat {b}}_{J})^{T}=\lim _{\eta \rightarrow \infty }{\hat {b}}^{(\eta )}}$.
• The factorization is not unique, since it is ${\displaystyle m_{ij}=a_{i}b_{j}=(\gamma a_{i})({\frac {1}{\gamma }}b_{j})}$ for all ${\displaystyle \gamma >0}$.
• The factor totals remain constant, i.e. ${\displaystyle \sum _{i}{\hat {a}}_{i}^{(\eta )}=\sum _{i}{\hat {a}}_{i}^{(1)}}$ for all ${\displaystyle \eta \geq 1}$ and ${\displaystyle \sum _{j}{\hat {b}}_{j}^{(\eta )}=\sum _{j}{\hat {b}}_{j}^{(0)}}$ for all ${\displaystyle \eta \geq 0}$.
• To fit the Q-model, where ${\displaystyle m_{ij}=0}$ a priori for ${\displaystyle (i,j)\in S}$, set ${\displaystyle \delta _{ij}=0}$ if (${\displaystyle i,j)\in S}$ and ${\displaystyle \delta _{ij}=1}$ otherwise. Then

${\displaystyle {\hat {a}}_{i}^{(\eta )}={\frac {x_{i+}}{\sum _{j}\delta _{ij}{\hat {b}}_{j}^{(\eta -1)}}},}$

${\displaystyle {\hat {b}}_{j}^{(\eta )}={\frac {x_{+j}}{\sum _{i}\delta _{ij}{\hat {a}}_{i}^{(\eta )}}}}$

${\displaystyle {\hat {m}}_{ij}^{(2\eta )}=\delta _{ij}{\hat {a}}_{i}^{(\eta )}{\hat {b}}_{j}^{(\eta )}}$

Obviously, the I-model is a particular case of the Q-model.

The Problem: Let ${\displaystyle M=(m_{ij}^{(0)})\in \mathbb {R} ^{I\times J}}$ be the initial matrix with nonnegative entries, ${\displaystyle u\in \mathbb {R} ^{I}}$ a vector of specified row marginals (e.i. row sums) and ${\displaystyle v\in \mathbb {R} ^{J}}$ a vector of column marginals. We wish to compute a matrix ${\displaystyle {\hat {M}}=({\hat {m}}_{ij})\in \mathbb {R} ^{I\times J}}$ similar to M and consistent with the predefined marginals, meaning

${\displaystyle {\hat {m}}_{i+}=\sum _{j=1}^{n}{\hat {m}}_{ij}=u_{i}}$

and

${\displaystyle {\hat {m}}_{+j}=\sum _{i=1}^{m}{\hat {m}}_{ij}=v_{j}}$

Define the diagonalization operator ${\displaystyle diag:\mathbb {R} ^{k}\longrightarrow \mathbb {R} ^{k\times k}}$, which produces a (diagonal) matrix with its input vector on the main diagonal and zero elsewhere. Then, for ${\displaystyle \eta \geq 0}$, set

${\displaystyle M^{(2\eta +1)}={\text{diag}}(r^{(\eta +1)})M^{(2\eta )}}$

${\displaystyle M^{(2\eta +2)}=M^{(2\eta +1)}{\text{diag}}(s^{(\eta +1)})}$

where

${\displaystyle r_{i}^{\eta +1}={\frac {u_{i}}{\sum _{j}m_{ij}^{(2\eta )}}}}$

${\displaystyle s_{j}^{\eta +1}={\frac {v_{j}}{\sum _{i}m_{ij}^{(2\eta +1)}}}}$

Finally, we obtain ${\displaystyle {\hat {M}}=\lim _{\eta \rightarrow \infty }M^{(\eta )}.}$

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 (and thus RAS) maintains the crossproduct ratios, e.i.

${\displaystyle {\frac {m_{ij}^{(0)}m_{hk}^{(0)}}{m_{ik}^{(0)}m_{hj}^{(0)}}}={\frac {m_{ij}^{(\eta )}m_{hk}^{(\eta )}}{m_{ik}^{(\eta )}m_{hj}^{(\eta )}}}\ \forall \ \eta \geq 0{\text{ and }}i\neq h,\quad j\neq k}$

since ${\displaystyle m_{ij}^{(\eta )}=a_{i}^{(\eta )}b_{j}^{(\eta )}.}$

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 (1970).

Nevertheless, direct factor estimation (algorithm 2) is under all circumstances the best way to deal with IPF: Whereas classical IPFP needs

${\displaystyle IJ(2+J)+IJ(2+I)=I^{2}J+IJ^{2}+4IJ\,}$

elementary operations in each iteration step (including a row and a column fitting step), factor estimation needs only

${\displaystyle I(1+J)+J(1+I)=2IJ+I+J\,}$

operations being at least one order in magnitude faster than classical IPFP.

Necessary and sufficient conditions for the existence and uniqueness of MLEs are complicated in the general case (see[11]), but sufficient conditions for 2-dimensional tables are simple:

• the marginals of the observed table do not vanish (that is, ${\displaystyle x_{i+}>0,\ x_{+j}>0}$) and
• the observed table is inseparable (i.e. the table does not permute to a block-diagonal shape).

If unique MLEs exist, IPFP exhibits linear convergence in the worst case (Fienberg 1970), but exponential convergence has also been observed (Pukelsheim and Simeone 2009). If a direct estimator (i.e. a closed form of ${\displaystyle ({\hat {m}}_{ij})}$) exists, IPFP converges after 2 iterations. If unique MLEs do not exist, IPFP converges toward the so-called extended MLEs by design (Haberman 1974), 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.

Checking if the assumption of independence is adequate, one uses the Pearson X-squared statistic

${\displaystyle X^{2}=\sum _{i,j}{\frac {(x_{ij}-{\hat {m_{ij}}})^{2}}{\hat {m_{ij}}}}}$

or alternatively the likelihood-ratio test (G-test) statistic

${\displaystyle G=2\sum _{i,j}x_{ij}\log \ {\frac {x_{ij}}{{\hat {m}}_{ij}}}}$.

Both statistics are asymptotically ${\displaystyle \mathrm {X} _{r}^{2}}$-distributed, where ${\displaystyle r=(I-1)(J-1)}$ is the number of degrees of freedom. That is, if the p-values ${\displaystyle 1-\mathrm {X} _{r}^{2}(X^{2})}$ and ${\displaystyle 1-\mathrm {X} _{r}^{2}(G)}$ are not too small (> 0.05 for instance), there is no indication to discard the hypothesis of independence.

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.

Consider a table of observations (taken from the entry on contingency tables):

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

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:

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:

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 ${\displaystyle 87/13\approx 6.69}$, since ${\displaystyle 45.24/6.76=41.76/6.24=87/13}$. Similarly, for both right-handed and left-handed, the odds of being male vs. female are ${\displaystyle 52/48=13/12\approx 1.08}$, since ${\displaystyle 45.24/41.76=6.76/6.24=13/12}$.

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:

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 ${\displaystyle p(X^{2})\approx 0.1824671}$, meaning: gender and left-handedness/right-handedness can be considered independent.

The R package mipfp (currently in version 3.1) provides a multi-dimensional implementation of the traditional iterative proportional fitting procedure.[12] The package allows the updating of a N-dimensional array with respect to given target marginal distributions (which, in turn can be multi-dimensional).

Python has an equivalent package, ipfn[13][14] that can be installed via pip. The package supports numpy and pandas input objects.