Cubic Regression: Steps in R, Python, SPSS, Excel

Cubic regression fits the model y = b0 + b1x + b2x² + b3x³ to data by least squares. It is a polynomial regression of degree 3: still linear in its four coefficients, but able to bend twice. This guide is for researchers who need to fit one in R, Python or SPSS, report the coefficients correctly, and decide whether the cubic term earns its place over a quadratic or linear model.

Cubic polynomial regression (this page) is not cubic spline regression. A cubic polynomial uses one global equation with four coefficients. A cubic spline joins separate cubic pieces at knots. If you need local flexibility, see SciPy interpolation options instead.

What cubic regression is used for in research

Use it when a scatter plot shows a relationship that curves in one direction and then the other (an S-shape or a single hump and trough within the observed range), and a straight line or a quadratic model leaves a systematic pattern in the residuals. The cubic term lets the fitted curve change direction up to twice. Typical uses are describing a calibration or growth curve over a limited range and testing whether a trend is non-linear (see how to detect a non-linear relationship).

It is a descriptive device for the observed range. A cubic has no mechanistic meaning, and outside the data it heads to plus or minus infinity, so do not extrapolate with it.

Install or access (official sources)

• R: poly() and lm() are in the base stats package, which ships with R from the R Project (r-project.org). No add-on is needed.
• Python: NumPy provides numpy.polynomial.Polynomial.fit (install NumPy from numpy.org or PyPI).
• SPSS Statistics: Curve Estimation is a procedure in IBM SPSS Statistics, under Analyze > Regression.
• Excel: a polynomial trendline of order 3 is available on a chart.

Step-by-step: fit a cubic in R, Python, SPSS and Excel

R

1. Fit with a raw cubic so the coefficients match the equation: fit3 <- lm(y ~ poly(x, 3, raw = TRUE), data = d). The R manual states that poly() returns orthogonal polynomials of degree 1 to degree by default and uses raw polynomials when raw = TRUE.
2. Or fit the orthogonal version for testing: fit3o <- lm(y ~ poly(x, 3), data = d). The fitted values and R² are identical to the raw version; only the coefficients differ, and the orthogonal ones do not map to b1, b2, b3 of the equation above.
3. Fit the quadratic and compare: fit2 <- lm(y ~ poly(x, 2, raw = TRUE), data = d), then anova(fit2, fit3).
4. Read the result with summary(fit3). Check the degree limit: with the default orthogonal polynomials, the manual says degree must be less than the number of unique points.

Python (NumPy)

1. from numpy.polynomial import Polynomial
2. p = Polynomial.fit(x, y, 3). The signature is Polynomial.fit(x, y, deg, domain=None, rcond=None, full=False, w=None, window=None, symbol='x').
3. Recover coefficients in the original units of x: p.convert().coef. Per the NumPy documentation, when domain is None the fit uses the minimal domain covering x, so p.coef alone is in a scaled basis and will not match your equation.

SPSS Statistics

1. Choose Analyze > Regression > Curve Estimation…
2. Move the outcome to Dependent and the predictor to Independent (a variable, or Time).
3. Tick Cubic in the Models list (the procedure offers 11 models, including linear and quadratic, so you can compare them in one run).
4. Keep Include constant in equation (the default), and tick Plot models and Display ANOVA table.

IBM defines the cubic model as Y = b0 + (b1 * t) + (b2 * t**2) + (b3 * t**3).

Excel

Add a scatter chart, add a trendline, choose the Polynomial type and set the order to 3. Show the equation on the chart if you need the coefficients. Microsoft documents the polynomial trendline as a least-squares fit.

Worked example (simulated data)

This is a simulated illustration, not real study data. We generated 40 points at x = 0 to 10 from y = 2 + 1.5x − 0.45x² + 0.03x³ plus normal noise (SD 0.6, NumPy seed 20261107), then fitted degrees 2 and 3 by least squares.

Model	R²	Note
Quadratic	0.444	Misses the second bend
Cubic	0.838 (adjusted 0.824)	Nested F test for the cubic term: F(1, 36) = 87.3, p < 0.001

Fitted cubic coefficients in the original x units: b0 = 1.67, b1 = 1.97, b2 = −0.570, b3 = 0.0381 (standard error of b3 = 0.0041). With Polynomial.fit and no conversion, p.coef returned 2.03, −4.39, 0.02, 4.76 for the same curve, because x was rescaled from [0, 10] to [−1, 1]. The converted coefficients matched numpy.polyfit exactly. Scaling also changes the design matrix condition number: about 1,850 for raw x versus about 7.8 for x centred and divided by 5 in this example.

Output and how to read it

• Test the cubic term, not the individual lower-order terms. b1, b2 and b3 are strongly correlated because x, x² and x³ are, so none of the lower-order coefficients is interpretable alone. Keep lower-order terms in the model whenever the cubic term is included.
• Compare nested models. The F test comparing the quadratic with the cubic (anova(fit2, fit3) in R) asks whether the extra bend improves fit beyond chance. A higher R² alone always follows from adding a term; see how to interpret R-squared and use adjusted R².
• Plot the fitted curve over the data and inspect residuals. A cubic that fits well numerically but swings sharply at the ends signals overfitting; see how to read a residual plot and what overfitting is.
• Report the equation with coefficients in original units and state whether predictors were centred or scaled.

Troubleshooting

Symptom	Likely cause and fix
NumPy coefficients do not match your hand-written equation	You read p.coef (scaled basis). Use p.convert().coef.
R coefficients differ from the equation	poly(x, 3) is orthogonal by default; use raw = TRUE to match the equation.
R error about degree	With orthogonal polynomials, degree must be less than the number of unique x values.
Huge or unstable coefficients, large standard errors	Collinear powers of x. Centre or scale x (see variance inflation factor) or use orthogonal polynomials for testing.
Cubic term significant but curve is implausible outside the data	Do not extrapolate; consider a different functional form.

Free alternatives

Tool	Cubic fit	Source
R (stats)	lm(y ~ poly(x, 3, raw = TRUE))	R manual, poly
NumPy	Polynomial.fit(x, y, 3), then .convert().coef	NumPy documentation

For a MATLAB route, see polyfit, fit and fittype in MATLAB.

Frequently asked questions

What is the formula for cubic regression?

y = b0 + b1x + b2x² + b3x³, fitted by least squares. SPSS documents the same form as Y = b0 + (b1 * t) + (b2 * t**2) + (b3 * t**3).

How is cubic regression different from quadratic and linear regression?

Linear has one slope, quadratic adds one bend (x²), and cubic adds a second possible bend (x³). Each is nested in the next, so an F test can compare them. See quadratic regression.

Is cubic regression the same as a cubic spline?

No. Cubic regression is one polynomial over the whole range; a cubic spline joins separate cubic pieces at knots.

Can I do cubic regression in SPSS?

Yes: Analyze > Regression > Curve Estimation…, then tick Cubic.

Is cubic regression a linear model?

Yes in the coefficients, which is why ordinary least squares works, although the fitted curve is non-linear in x.

Sources

• R Core Team, poly: Compute Orthogonal Polynomials, R documentation (stat.ethz.ch/R-manual/R-devel/library/stats/html/poly.html), read 2026-11-07.
• NumPy documentation, numpy.polynomial.polynomial.Polynomial.fit (numpy.org/doc/stable/reference/generated/numpy.polynomial.polynomial.Polynomial.fit.html), read 2026-11-07.
• IBM Documentation, Curve Estimation and Curve Estimation Models, IBM SPSS Statistics (ibm.com/docs), read 2026-11-07.
• Microsoft Learn, XlTrendlineType enumeration (Excel) (learn.microsoft.com/en-us/office/vba/api/excel.xltrendlinetype), read 2026-11-07.