0
$\begingroup$

I have a few questions about handling missing data during matrix/array operations. Essentially, I'm doing a vectorized linear regression to perform a bunch of linear regressions in a few matrix/array operations. For a random data X_many and y_many of shapes (2, 50, 3) and (2, 50), I use the following code to estimate the coefficient matrix:

X_op = np.einsum('ikj,ilj->ilk', (np.linalg.inv(np.einsum('ikj,ikl->ijl', X_many, X_many))), X_many)
Beta_hat = np.einsum('ikj,ik->ij', X_op, y_many)

The issue here is that X_many has some missing values. Specifically, in some groups, some variable values are completely unavailable. When I run the code above with those missing values an entire row of coefficients is nan in the matrix. But if I did the regressions one by one, I could omit the missing values, estimate a set of coefficients, and then assume that the missing value coefficients are 0. That way I would still be able to do estimation, but with less data since some coefficients will be 0 (see code below):

for x_i, y_i in zip(X_many, y_many):
    x_i = x_i[:, np.isfinite(x_i.sum(axis=0))]
    beta_hat_i = np.linalg.inv(x_i.T @ x_i) @ x_i.T @ y_i

Is there a way to get the coefficients of individual runs, but still use the array products method? This got me thinking that if I could impute the missing values with data that'd yield a coefficient of 0, then the two answers would be identical. Imputing all 0s didn't work since it resulted in a Singular matrix error, thus I decided to use the second-best alternative and impute with a series that has a very low absolute correlation with y (white noise). I didn't know how to generate that either so I kinda brute forced it and it kinda works?

TLDR: How do I generate a series with exactly 0 correlation with y? If I can't, then is there maybe a better imputation method or just a method that'd allow me to compute these coefficients using array dot products?

Code to reproduce results:

import numpy as np
import statsmodels.api as sm

# data gen
np.random.seed(42)
n = 50
X = np.random.random((n, 2))
X = sm.add_constant(X)
beta = [200, 100, 50]
e = np.random.random(n) * 25
y = np.dot(X, beta) + e

X2 = np.random.random((n, 2))
X2 = sm.add_constant(X2)
beta2 = [500, 10, 5]
e2 = np.random.random(n) * 5
y2 = np.dot(X2, beta2) + e

X[:, 1] = np.NaN # missing values

X_many = np.array([X, X2])
y_many = np.array([y, y2])

print(X_many.shape)
print(y_many.shape)

for x_i, y_i in zip(X_many, y_many):
    x_i = x_i[:, np.isfinite(x_i.sum(axis=0))]
    beta_hat_i = np.linalg.inv(x_i.T @ x_i) @ x_i.T @ y_i
    print(beta_hat_i)

np.random.seed(42)
X_many_new = X_many.copy()
noise = np.random.normal(0, 1e12, (100, np.isnan(X_many_new).sum()))
X_many_new[np.isnan(X_many_new)] = noise[np.argmin(np.abs(np.corrcoef(noise, y)[-1][:-1]))]

X_op = np.einsum('ikj,ilj->ilk', (np.linalg.inv(np.einsum('ikj,ikl->ijl', X_many_new, X_many_new))), X_many_new)
Beta_hat = np.einsum('ikj,ik->ij', X_op, y_many)
print(Beta_hat)
```
$\endgroup$

1 Answer 1

0
$\begingroup$

If anyone is interested in how to impute the missing series with values that would have a corresponding coefficient of 0 and would generate results equivalent to the omitted value regression.

To achieve this a specific condition has to be true. The partial correlation between x and y given the other exogenous variables has to be zero. For this to be true the conditional covariance has to also be 0.

Thus for missing values $x$, target $y$ and other exogenous variables $z$: $$ cov(x,y|z) = 0 $$ This implies that given $Z$: $$ Z = z(z^Tz)^{-1}z^T $$ and standardised (zero mean) ols residuals $r_1$ and $r_2$: $$ r_1 = x - Zx\\ r_2 = y - Zy\\ r_1 r_2 = 0\\ (x - Zx)(y - Zy)=0\\ (I - Z)x(y - Zy)=0\\ ((I - Z^T)(y - Zy))x=0 $$ To conclude to make sure that the above condition holds, the missing imputation values x should be generated so that they are orthogonal to the vector $(I - Z^T)(y - Zy).$

See working example below:

import numpy as np
import statsmodels.api as sm
np.set_printoptions(suppress=False, formatter={'float_kind':'{:f}'.format})

def generate_vector_orthogonal_unit_vector(y):
    # Step 1: Normalize y
    y_hat = y - y.mean()
    y_hat = y_hat / np.linalg.norm(y_hat)
    # Step 2: Generate a random vector v (not parallel to y)
    v = np.random.normal(0, 1, len(y))
    # Step 3: Orthogonalize v with respect to y
    x = v - np.dot(v, y_hat) * y_hat
    x = x - x.mean()  # Normalize x
    return x

# data gen
np.random.seed(42)
n = 50
X = np.random.random((n, 2))
X = sm.add_constant(X)
beta = [200, 100, 50]
e = np.random.random(n) * 25
y = np.dot(X, beta) + e

X2 = np.random.random((n, 2))
X2 = sm.add_constant(X2)
beta2 = [500, 10, 5]
e2 = np.random.random(n) * 5
y2 = np.dot(X2, beta2) + e

X[:, 1] = np.NaN # missing values

X_many = np.array([X, X2])
y_many = np.array([y, y2])

print(X_many.shape)
print(y_many.shape, '\n')

# actual coefs
for x_i, y_i in zip(X_many, y_many):
    x_i = x_i[:, np.isfinite(x_i.sum(axis=0))]
    beta_hat_i = np.linalg.inv(x_i.T @ x_i) @ x_i.T @ y_i
    print(beta_hat_i)
print('')

# fill missing values with effectively missing series
np.random.seed(42)
groups = np.concatenate([X_many, y_many[:, :, np.newaxis]], axis=-1)
for i in range(groups.shape[0]):
    group_finite = groups[i][:, np.isfinite(groups[i].sum(axis=0))]
    y = group_finite[:, -1]
    z = group_finite[:, :-1]
    Z = z @ np.linalg.inv(z.T @ z) @ z.T
    I = np.eye(Z.shape[0], Z.shape[1])
    for j in range(groups.shape[2]):
        if np.isnan(groups[i, :, j].sum()):
            # main logic goes here
            r_leak = (I - Z.T) @ (y - Z @ y)
            v_orthogonal = generate_vector_orthogonal_unit_vector(r_leak)
            groups[i, :, j] = v_orthogonal
X_many_new = groups[:, :, :-1]

# Estimating coefs with array method
X_op = np.einsum('ikj,ilj->ilk', (np.linalg.inv(np.einsum('ikj,ikl->ijl', X_many_new, X_many_new))), X_many_new)
Beta_hat = np.einsum('ikj,ik->ij', X_op, y_many)
print(Beta_hat)
```
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.