# Vectorized Linear Regression with missing values

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))
beta = [200, 100, 50]
e = np.random.random(n) * 25
y = np.dot(X, beta) + e

X2 = np.random.random((n, 2))
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)
$$$$


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))
beta = [200, 100, 50]
e = np.random.random(n) * 25
y = np.dot(X, beta) + e

X2 = np.random.random((n, 2))
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)
$$$$