Open as a notebook

Working with Missing Data

Motivation

Suppose a longitudinal study is intended to measure some feature of a subject at a number of points in time. To keep track of which measured value corresponds with each time point, the data is stored in an array: index 0 corresponds with the first time, index 1 corresponds with the second, and so on.

import os
os.environ['SCIPY_ARRAY_API'] = '1'
import numpy as np

# Measurements at each of five times
data = np.asarray([1.2, 2.3, 3.4, 4.5, 5.6])

If no measurement is made at a designated time, how should the corresponding element of the array be filled? In general, how do we deal with missing data given the constraints that arrays must be rectanglar - not ragged - and all entries must be filled?

Sentinel Values

In the simplest approach, a “sentinel value”, which does not appear elsewhere in the valid measurement data, is chosen to represent a missing measurement. Frequently, this is a value that lies outside the possible range of measurement, such as a negative value when the possible range is strictly positive.

# A sentinel value, -1., can be used to represent a
# missing measurement of a positive quantity
data = np.asarray([1.1, 1.2, 1.3, -1., 1.5])

This approach satisfies the requirements of arrays, and in principle we can track the fact that the data is missing. But how can we perform valid computations with this data, given than -1 must not be interpreted as a measurement?

Manual Sentinel Value Removal

Suppose we wish to take the harmonic mean of the valid measurements using scipy.stats.hmean. One approach is to manually eliminate the sentinel values, producing a temporary array of a smaller size, and to pass this temporary array to hmean.

from scipy import stats
temp = data[data > 0]
stats.hmean(temp)

np.float64(1.2585258525852587)

But suppose we have more than one subject, with different missing measurements for each subject.

data = np.asarray([[1.1, 1.2, 1.3, -1., 1.5],   # four valid measurements, subject 1
                   [2.9, -1., -1., 2.6, 2.5]])  # three valid measurements, subject 2

We cannot follow the same approach here:

temp = data[data > 0]
temp

array([1.1, 1.2, 1.3, 1.5, 2.9, 2.6, 2.5])

stats.hmean(temp, axis=-1)

np.float64(1.6249727109374565)

data[data > 0] produces a one-dimensional array, so hmean is not be able to produce separate harmonic means for each subject.

One solution is to loop manually over the rows:

res = []
for row in data:
    temp = row[row > 0]
    res.append(stats.hmean(temp))
res = np.asarray(res)
res

array([1.25852585, 2.65617661])

This is valid, but cumbersome and potentially slow for datasets with many subjects. Fortunately, SciPy provides alternatives.

A Common Choice: nan_policy='omit'

If the floating point value NaN (Not a Number) is used as the sentinel:

NaN = np.nan
# data[data < 0] = NaN, or more explicitly:
data = np.asarray([[1.1, 1.2, 1.3, NaN, 1.5],   # four valid measurements, subject 1
                   [2.9, NaN, NaN, 2.6, 2.5]])  # three valid measurements, subject 2

then passing the option nan_policy='omit' instructs SciPy to automatically remove NaNs from each slice of the data while performing the computation.

stats.hmean(data, axis=-1, nan_policy='omit')

array([1.25852585, 2.65617661])

Almost all reducing statistics in scipy.stats support nan_policy='omit'. Coverage is nearly complete because it is implemented in the generic way: looping over the slices, and eliminating the NaNs before performing the operation for each slice. As discussed, this Python-level looping can be slow when there are many slices, so nan_policy='omit' is offered merely for batch calculation convenience, not for speed. Another problem with this approach is that it overloads the meaning of NaN, which is ordinarily used as the result of an invalid calculation, like 0 / 0. Finally, this option is not offered for backends other than NumPy (e.g. CuPy, PyTorch).

Fortunately, there is a more principled approach that can be compatible with alternative backends and faster for large batches.

Masked Arrays

Instead of using sentinel values, fill the space of missing values with arbitrary data, and use a second, boolean array of the same shape - a “mask” - to keep track of which elements are missing.

data = np.asarray([[1.1, 1.2, 1.3, 1.4, 1.5],
                   [2.9, 2.8, 2.7, 2.6, 2.5]])
mask = np.asarray([[False, False, False,  True, False],
                   [False,  True,  True, False, False]])

The Traditional Option: MaskedArray

NumPy offers numpy.ma.MaskedArray for working with masked data, and functions in scipy.stats.mstats were provided to work with these NumPy masked arrays.

In principle, the masked array approach is advantageous because it has the potential to avoid conflating missing NaN values with invalid NaN values. It can also be faster in batch calculations with many slices, because batched masked array calculations can be implemented to ignore masked values without introducing Python for loops.

However, the mstats.hmean function is now deprecated along with the mstats namespace and all other uses of MaskedArray in scipy.stats.

x = np.ma.MaskedArray(data, mask=mask)
stats.mstats.hmean(x, axis=-1)

/tmp/ipykernel_1401/2803316482.py:2: DeprecationWarning: `scipy.stats.mstats` is deprecated as of SciPy 2.0.0 and will be removed in SciPy 2.4.0. See function documentation for alternatives.
  stats.mstats.hmean(x, axis=-1)
/tmp/ipykernel_1401/2803316482.py:2: DeprecationWarning: `scipy.stats.mstats.hmean` is deprecated as of SciPy 2.0.0 and will be removed, along with the `scipy.stats.mstats` namespace, in SciPy 2.4.0. For similar functionality, use `scipy.stats.hmean` with MArray(s) instead of NumPy masked array(s). 
  stats.mstats.hmean(x, axis=-1)
/tmp/ipykernel_1401/2803316482.py:2: DeprecationWarning: Support for NumPy masked arrays is deprecated as of SciPy 2.0.0 and will be removed in SciPy 2.4.0. See function documentation for alternatives.
  stats.mstats.hmean(x, axis=-1)

array([1.25852585, 2.65617661])

There are several reasons.

The first is that NumPy masked arrays do, in fact, conflate invalid and missing values. Consider the following example:

x = np.asarray([0, 1, 2, 3, 4])
np.sum(x / x)

/tmp/ipykernel_1401/2028185835.py:2: RuntimeWarning: invalid value encountered in divide
  np.sum(x / x)

np.float64(nan)

Ordinary NumPy arrays warn that 0 / 0 produces NaN, and this invalid value propagates in the sum. It is impossible to get a valid numerical result when a NaN is involved in arithmetic.

Yet NumPy masked arrays seem to provide a number.

y = np.ma.MaskedArray(x)
np.sum(y / y)

np.float64(4.0)

This occurs because NaNs arising from invalid numerical calculations involving NumPy masked arrays are masked without warning (and subsequently ignored).

y / y

masked_array(data=[--, 1.0, 1.0, 1.0, 1.0],
             mask=[ True, False, False, False, False],
       fill_value=1e+20)

This can lead to invalid calculations producing apparently valid but actually bogus numerical results, which is clearly unsafe in scientific computing. Rather than alerting the user to the invalid result so it can be fixed, NumPy masked arrays hide the problem and produce erroneous values.

The second reason for deprecating mstats is that its function interfaces and implementations were entirely separate from those of stats, and often neglected in terms of maintenance and enhancements. Consider, for instance, the stark difference in documentation thoroughness and feature completeness between stats.mannwhitneyu and mstats.mannwhitneyu. Eliminating mstats in favor of adding missing data support to stats allows SciPy maintainers to provide users with one, complete implementation with a single interface.

The final reason for deprecating mstats and support for NumPy masked arrays is the rise in support for array API standard compatible arrays throughout SciPy. NumPy masked arrays themselves are mostly unmaintained and do not conform to the Array API Standard, so when adding support for high priority libraries like CuPy, JAX, and PyTorch, which are compatible with the standard, it is difficult to also preserve support for the legacy MaskedArray type.

The Modern Option: MArray

Fortunately, as support for the Array API Standard closes this window, it opens the door toward a new way of supporting masked data. Specifically, MArray is an array API standard compatible array type that wraps the functionality of other array backends and endows them with support for masks.

from marray import numpy as xp  # or:
# from marray import torch as xp
# from marray import cupy as xp
x = xp.asarray(data, mask=mask)
x

MArray(
    array([[1.1, 1.2, 1.3,   _, 1.5],
           [2.9,   _,   _, 2.6, 2.5]]),
    array([[False, False, False,  True, False],
           [False,  True,  True, False, False]])
)

stats.hmean(x, axis=-1)

MArray(array([1.25852585, 2.65617661]), array([False, False]))

Consequently, existing users of mstats and NumPy masked arrays are advised to begin using the corresponding stats functions with MArrays where possible (see function documentation) and nan_policy='omit' otherwise.