How to interpret the coefficients returned from a multivariate cubic regression (polynomial degree 3) when using linearRegression().coef_?

Question

I am trying to fit a hyperplane to a dataset which includes 2 features and 1 target variable. I processed the features using PolynomialFeatures.fit_transform() and PolynomialFeature(degree = 3), and then fitted those features and target variable into a LinearRegression() model. When I use LinearRegression().coef_ to get the coefficients in order to write out a function for the hyperplane (I want the written out function itself), 10 coefficients are returned and I don't know how to interpret them into a function. I know that for a PolynomialFeature(degree = 2) model, 6 coefficients are returned and the function looks like m[0] + x1*m[1] + x2*m[2] + (x1**2)*m[3] + (x2**2)*m[4] + x1*x2*m[5] where m is the list of coefficients returned in that order. How would I interpret the cubic one?

Here is what my code for thee cubic model looks like:

poly = polyF(degree = 3)
x_poly = poly.fit_transform(x)
model = linR()
model.fit(x_poly, y)

model.coef_
(returns):
array([ 0.00000000e+00, -1.50603348e+01,  2.33283686e+00,  6.73172519e-01,
      -1.93686431e-01, -7.30930307e-02, -9.31687047e-03,  3.48729458e-03,
       1.63718406e-04,  2.26682333e-03])

Answer 1

So if (X1,X2) transforms to (1,X1,X2,X1^2,X1X2,X2^2)

Then (X1,X2,X3) should transform to (1, X1, X2, X3, X1X2, X1X3, X2X3, X1^2 * X2, X2^2 * X3, X3^2 * X1)

Answer 2

I know I am very late for this, but I had the same question and I found out on sklearn's documentation that you can easily have access to the exponents for each coefficient using PolynomialFeatures.fit(x).powers_, Although I guess they're not shown in the most intuitive way. Because of this, I quickly wrapped together a code that can display it in a more "readable" manner, so to say:

def contrasts(poly):
  c_ = ['1']
  for i in range(poly.powers_.shape[0]-1):
    app = str()
    for j in range(poly.powers_.shape[1]):
      if poly.powers_[i+1][j] > 1:
        app += f'x{j+1}**{poly.powers_[i+1][j]}'
      elif poly.powers_[i+1][j] == 1:
        app += f'x{j+1}'
    c_.append(app)
  return c_
poly_ = polyF(degree = 3)
exponents = contrasts(poly_.fit(x))

Hopefully this runs just fine for you! It did for me and can be used with at least second and third degree polynomials (I think this is good to use for n-degree polys, but I haven't tested so take this with a grain of salt).

Answer 3

I was facing the same question and developed the following code block to print the fit equation. To do so, it was necessary to include_bias=True in PolynomialFeatures and to set fit_intercept=False in LinearRegression, as opposed to conventional use:

import pandas as pd
from sklearn import linear_model
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import PolynomialFeatures

def polyReg():

  seed=12341

  df=pd.read_csv("input.txt", delimiter=', ', engine='python')

  X=df[["x1","x2","x3"]]
  y=df["y"]    

  poly=PolynomialFeatures(degree=2,include_bias=True)
  poly_X=poly.fit_transform(X)

  X_train,X_test,y_train,y_test=train_test_split(poly_X,y,test_size=0.5,random_state=seed)

  regression=linear_model.LinearRegression(fit_intercept=False)

  fit=regression.fit(X_train,y_train)    

  variable_names=poly.get_feature_names_out(X.columns)
  variable_names=np.core.defchararray.replace(variable_names.astype(str),' ','*')    

  fit_coeffs=["{:0.5g}".format(x) for x in fit.coef_]

  arr_list=[fit_coeffs,variable_names]
  fit_equation=np.apply_along_axis(join_txt, 0, arr_list)
  fit_equation='+'.join(fit_equation)
  fit_equation=fit_equation.replace("*1+","+")
  fit_equation=fit_equation.replace("+-","-")

  print("Fit equation:")    
  print(fit_equation)

def join_txt(text,delim='*'):

   return np.asarray(delim.join(text),dtype=object)