Help performing polynomial features and linear regression.

I have a dataframe with columns X, Y, Z. I need to perform polynomial features and linear regression to fit a polynomial line to the scatter plot.

Some example data points

X: [-3.5, 8.8,-20,-7.9,-14, 16, -14.5, 12, -4]

Y:[-2, 8, -19.5, -7.5, -13, 14, -14, 12.5, -3]

Z:[-960, 110, 4867, -149, 1493, 2411, 1938, 2201, 1596, -870]

I know how to do two dimensional data with polynomial features but not 3dimensional.

from sklearn.preprocessing import PolynomialFeatures

#X data polynomial features

polf=PolynomialFeatures(degree=2,include_bias=False) Xpolf=polf.fit_transform(X.reshape(-1,1))

and then do linear regression

from sklearn.linear_model import LinearRegression

linreg = LinearRegression() linreg=linreg.fit(Xpolf,y) linreg.intercept_,linreg.coef_

However, I can't figure out what to do for 3 dimensional data.

I tried doing polynomial features on x and y and then passing a dataframe into linear regression. However, when I try to do a test set x, y, z to fit the scatterplot, I can't get it to work.