# Multivariate regression with weighted least squares in python?

I have a multivariate regression problem that I need to solve using the weighted least squares method. In particular, I have a dataset X which is a 2D array. It consists of a number of observations, n, and each observation is represented by one row. Each observation also consists of a number of features, m. So that means each row has m columns. Therefore my dataset X is a n×m array.

Given a test data observation, multivariate regression should produce a function that predicts the response vector y, which is a 2D array as well. This function will consist of m coefficients, i.e. one coefficient/parameter for each of the m features of the test input.

This solution is already available here: Statsmodels' WLS, except that they don't support a 2D response vector yet. In other words, when I fit the data, I have to provide my dataset X, but can only provide a 1D array as the response y.

In addition, I also need a 2D weights vector, similar in dimension to the response vector y.

Is there a Python implementation of WLS multivariate regression where y and the weights can be 2D vectors?

Or if not a direct implementation, can any of the existing packages be used as an implementation somehow, by a small amount of adjustment?

## Edit

To make my question clearer, these are the parameters I would give in and the results I would need to get out:

Input:

• X : a 2D dataset, like 10x3, which is 10 observations with 3 features each.

• y : which is also a 2D vector, in this case 10x2. In other words, a 2-value response vector for each observation. (I'm doing classification and there are two possible classes).

• weights: a 2D response vector which is also 10x2, like y.

The 10 above is an arbitrary number of rows. Ultimately, however many observations I have, that's just what the number of rows is going to be, for all vectors above.

Output I need:

• the coefficients of the regression. Given that my response and weight vectors are 2D, I believe the coefficients would also be a 2D array, probably 3x2 or 2x3.
• Can you clarify what you mean by multivariate regression? Commented May 1, 2014 at 16:22
• I just edited my question to clarify it. Commented May 1, 2014 at 16:37
• This is probably a question for stackoverflow, btw. Commented May 1, 2014 at 17:14
• In what way is your response a 2-D vector? The description sounds like your response is one variable and that all your data X is presented by 2-D. Commented May 1, 2014 at 17:15
• I say my response is a 2D vector because the parameters that I supply to the multivariate regression function are: (1) X (a 2D dataset, like 10x3, which is 10 observations with 3 features each). (2) y, which is also a 2D vector, in this case 10x2 (a 2-value response vector for each observation). (3) weights, a 2D response vector which is also 10x2, like y. Commented May 1, 2014 at 19:51

It's still not entirely clear to me what you want to do, but if your weights are 1d, you can (ab)use sm.WLS to do this.

import numpy as np
import statsmodels.api as sm
np.random.seed(12345)

N = 30

X = np.random.uniform(-20, 20, size=(N,10))
beta = np.random.randn(11)

weights = np.random.uniform(1, 20, size=(N,))
weights = weights/weights.sum()

y = np.dot(X, beta) + weights*np.random.uniform(-100, 100, size=(N,))

Y = np.c_[y,y,y]

mod = sm.WLS(Y, X, weights=1/weights).fit()


If your weights are not 1d, WLS will indeed break, because it's not designed for this case. You can use a loop over WLS or just roll your own solution depending on what exactly you want to do.

weights = np.random.uniform(1, 20, size=(N,3))
weights = weights/weights.sum(0)
y = np.dot(X, beta)[:,None] + weights*np.random.uniform(-100, 100, size=(N,3))


This is the entirety of the WLS solution for each equation, assuming this is what you want to do

beta_hat = np.array([np.linalg.pinv(1/weights[:,i,None]**.5 * X).dot(y[:,i]) for i in range(3)])

• Wow that is a wonderful one-liner! However, I've just edited my question - I've now mentioned about the weights. I do in fact need 2D weights, of the same shape as the response vector y. I'm still new at python, but in the second block of code that you gave, I understand you've used 2D weights? Commented May 1, 2014 at 19:53
• Yes, the weights are 2d but they're applied equation by equation like sm.WLS(y[:i], X, weights = weights[:,i].fit(), if that's not what you want, maybe you can get there from what's posted. Commented May 2, 2014 at 2:39

In module sklearn, linear_model provides many regression functions, which will satisfy your demand.

For example, lasso.fit(X,y) where y has shape = (n_samples,) or (n_samples, n_targets). In your situation, n_targets = 2.

• I need a 2D array of weights too, though. Does lasso fit include weights? And, importantly (since I'm not familiar with this), can lasso fit be used for weighted least squares regression? I'm using this for gentleboosting in order to classify images. The features I'm using are HoG features. So each observation (or row) in $X$ will consist of many columns, i.e. many HoG features. Commented May 5, 2014 at 3:44
• Take also a look at sklearn.linear_model.lassocv Commented Jul 21, 2021 at 13:42
• @user961627: if you have sample_size smaller than qty of features, then statsmodels is not a good choice, but you can weighten your data in neural network (coded with tensorflow or keras) & create any dimensionality output Commented Aug 26, 2023 at 16:41