Polynomial Regression Calculator – Free Curve Fitting Tool

A polynomial regression calculator fits a curved line through your data by adding squared, cubed and higher powers of X to an ordinary least squares model.

Degree 2 gives one bend, degree 3 gives two bends, and so on. The model is still linear in its coefficients, so it is solved with the same least squares maths as straight-line regression. In the worked example on this page, 12 fertiliser plots give the fitted curve below, an R squared of 0.9910, and an optimum dose of 163.35 kg N per hectare.

Rule of thumb: start at degree 2 and only go higher if adjusted R squared rises, BIC falls, and the highest order term has p below 0.05. Above degree 3 you are usually fitting noise, not signal.

Polynomial regression is a regression model in which the relationship between the predictor X and the outcome Y is modelled as a polynomial of degree n. Instead of forcing a straight line through data that clearly bends, you add extra columns to the design matrix, X squared for one bend, X cubed for two bends, and let least squares find the coefficients.

The single most misunderstood point: polynomial regression is a linear model. The word linear in "linear regression" refers to the coefficients, not to the shape of the curve. The model Y = b₀ + b₁X + b₂X² is linear in b₀, b₁ and b₂, so it is fitted with exactly the same normal equations, and the same t tests and F tests apply. Only a model such as Y = a e^bX, where a coefficient sits inside a nonlinear function, needs true nonlinear least squares.

What each degree can do:

• Degree 1, a straight line. No bends. This is simple linear regression.
• Degree 2, quadratic. One bend, a single peak or a single trough. This covers most real research questions: an optimum dose, a saturation point, a performance peak.
• Degree 3, cubic. Two bends. Growth curves that accelerate then flatten, and S shapes.
• Degree 4 and above. Three or more bends. Almost always overfitting outside physics and engineering calibration work.

A worked one-liner: for the 12 fertiliser plots used throughout this page, the degree 2 fit is Y = 2.1544 + 0.042440X − 0.00012991X². Because b₂ is negative the curve opens downward, and the peak sits at X = −b₁ ÷ (2b₂) = 163.35 kg N per hectare, giving a predicted yield of 5.62 t per hectare.

Who uses it: agronomists find optimum fertiliser and irrigation rates, pharmacologists model dose response saturation, engineers calibrate sensors, economists model diminishing returns, sports scientists describe performance across the lifespan, and ecologists fit growth curves. Any time theory says "it rises then it falls", polynomial regression is the first model to try.

If you can write the answer as one equation and a reviewer can check it by hand, polynomial regression is doing its job. If you find yourself at degree 5 trying to chase every wiggle, the honest move is a spline instead.

Is polynomial regression linear or nonlinear?

Polynomial regression is a linear model. Linear describes the coefficients, not the curve on the page. Because Y = b₀ + b₁X + b₂X² is a straight-line combination of its coefficients, it is fitted with ordinary least squares and uses the same t tests and F tests as simple regression. Only a coefficient trapped inside a function, as in Y = a e^bX, makes a model genuinely nonlinear.

What degree should I use in polynomial regression?

Start at degree 2, then check four things before going higher: the p value of the highest order term, the incremental F test against the next lower degree, whether adjusted R squared rises and BIC falls, and whether the residual plot still shows an arc. If all four agree, the extra power has earned its place. Above degree 3 you are almost always fitting noise.

What is a good R squared for polynomial regression?

It depends on the field, not on a universal threshold. Calibration and laboratory work routinely exceeds 0.99, agronomy and ecology sit around 0.80 to 0.95, and 0.30 can be a genuine result in survey research. Always read adjusted R squared beside it, because plain R squared can never fall when you add a term and will reward a model that has simply memorised your data.

Can polynomial regression predict outside the range of my data?

No. Polynomials diverge sharply beyond the fitted range, and the higher the degree the faster it happens. A cubic that fits beautifully from 0 to 200 can predict a negative value at 260. Restrict every prediction to the span of X you actually measured, and if a reviewer asks about higher values, the honest answer is that you would need to collect them.

How many data points does polynomial regression need?

A degree n model needs at least n + 1 distinct X values before it can be fitted at all, but that is the arithmetic minimum, not a usable sample. Aim for roughly 10 observations per estimated coefficient, so about 30 points for a quadratic and 40 for a cubic, spread over at least six distinct X levels so the curve is identifiable.

4.1 Read the equation before you read anything else

The fitted equation is the result. Everything else on the page is evidence about how much you should trust it. Read it left to right. The intercept b₀ is the predicted value of Y when X is zero, which is meaningful only if zero is inside the range you actually measured. The linear coefficient b₁ is the slope at X = 0, not the slope everywhere, because in a curved model the slope changes as X changes. The quadratic coefficient b₂ controls the bend: negative means the curve opens downward and has a maximum, positive means it opens upward and has a minimum.

This is the single most common misreading in applied work. In a straight line model, one coefficient is the whole story. In a polynomial model, no coefficient means anything on its own. The slope at any point X is the derivative, b₁ + 2b₂X + 3b₃X² and so on, so you must state the X value you are talking about. Saying "the effect of fertiliser was 0.045 t per kg" is wrong in a quadratic model. Saying "at 50 kg N per hectare the marginal effect was 0.0294 t per kg, falling to 0.0060 t per kg at 140 kg" is right.

If you centred X, which this calculator does by default, the intercept becomes the predicted Y at the mean of X, which is nearly always more useful than the predicted Y at zero. The curve itself is identical either way. Centring changes the labels on the coefficients, not the fit.

4.2 What R squared tells you, and the trap inside it

R squared is the proportion of the variation in Y that the fitted curve explains. It runs from 0 to 1. An R squared of 0.9910 means the curve accounts for 99.10 per cent of the up and down movement in Y, leaving 0.90 per cent unexplained.

The trap is that R squared can never fall when you add a term. Add X cubed and R squared goes up, even if the cubic term is pure noise, because the extra column gives least squares one more knob to turn. Push the degree to n minus 1 and R squared reaches exactly 1.000, with a curve that passes through every point and predicts nothing. That is why R squared alone must never decide the degree.

Adjusted R squared fixes this by charging a penalty for every coefficient you spend. It can fall. When adjusted R squared stops rising as you add a power, you have found the honest stopping point. In the worked example, adjusted R squared goes 0.7208 at degree 1, 0.9890 at degree 2, 0.9879 at degree 3, then 0.9864 at degree 4. The huge jump is from 1 to 2. Everything after that is noise dressed up as improvement.

4.3 The overall F test: does the curve beat a flat line?

The F test in the ANOVA table compares your fitted polynomial against a model with nothing in it but the mean of Y. A significant F says the model explains more than chance would. In the worked example, F(2, 9) = 493.58 with p < 0.001, so the model is doing real work.

Be clear about what this test does not tell you. It does not say the degree is right. A degree 5 model fitted to quadratic data will also produce a wildly significant F, because two of its terms are real and three are junk riding along. The overall F is a screening test. Degree selection is done in section 4.4.

4.4 Choosing the degree, the part that actually matters

Four pieces of evidence point at the right degree, and you should look at all four before deciding.

1. The p value of the highest order term. If the X cubed term has p = 0.22, that cubic bend is not distinguishable from noise. Drop back to degree 2. This is the fastest check and it is what reviewers look at first.
2. The incremental F test. This asks whether the extra sum of squares bought by one more power is large relative to the residual variance. It is algebraically the square of the t test on that coefficient when only one term is added, so the two agree, but the F framing extends naturally to adding several terms at once.
3. Adjusted R squared and BIC. Adjusted R squared should rise. BIC should fall. BIC penalises extra terms harder than AIC does, so when they disagree, BIC is the conservative choice and the one to quote in a paper.
4. The residual plot. If residuals versus fitted values still show a clear arc or wave, curvature remains and the degree is too low. If they look like a random cloud, stop.

Never choose the degree by eyeballing which curve looks prettiest on the scatter plot. A higher degree always looks like a tighter fit, which is exactly the illusion overfitting creates.

4.5 The coefficient table: standard errors, t values and p values

Each row of the coefficient table gives an estimate, its standard error, a t value and a p value. The t value is simply the estimate divided by its standard error, and it asks whether that coefficient could plausibly be zero. Degrees of freedom are n minus the number of coefficients, so a degree 2 model on 12 points has 12 minus 3 = 9 degrees of freedom.

Two warnings. First, do not delete a lower order term because its p value is large while a higher order term is significant. If X squared stays, X stays with it. Removing the linear term forces the vertex of the parabola to sit at X = 0, which is a strong and almost always false assumption. This rule is called hierarchy or marginality, and violating it is a common referee complaint.

Second, in an uncentred fit the coefficients are correlated with each other by construction, because X and X squared move together. That inflates standard errors and can make individually significant terms look non-significant. Centring X, which is the default here, largely removes the problem.

4.6 Confidence intervals and what they are really saying

Each coefficient comes with a confidence interval, by default 95 per cent. Read it as the range of coefficient values that are compatible with your data at that confidence level. If the interval excludes zero, the coefficient is significant at the matching alpha, so the interval and the p value always agree.

Intervals are more informative than p values because they carry units and magnitude. An interval of [−0.000149, −0.000111] for b₂ says the downward curvature is small but very precisely estimated. An interval of [−0.004, 0.002] would say you cannot even agree on the direction of the bend.

4.7 The turning point, and why it is often the real answer

For a quadratic, the vertex sits at X = −b₁ ÷ (2b₂). This is usually the number the research question was really asking for: the optimum dose, the most efficient speed, the age of peak performance. The calculator reports it whenever the fit is quadratic and the vertex falls inside a sensible range.

Two cautions. If the vertex lies outside the range of X you actually measured, do not report it as an optimum. All you have shown is that the curve is still bending at the edge of your data. And the vertex is a ratio of two estimated coefficients, so its own uncertainty is wider than either coefficient's interval suggests. For a formal confidence interval on a turning point, use the delta method or a bootstrap rather than quoting the point estimate alone.

4.8 Residual diagnostics, four questions in plain words

Is the shape right? Residuals versus fitted values should look like a shapeless cloud. An arc means the degree is too low. This is the most important plot on the page.

Is the spread constant? A funnel that widens to the right means heteroscedasticity: the model is more wrong for large values. The Breusch-Pagan test puts a p value on it. Coefficients stay unbiased but standard errors, and therefore p values, become untrustworthy. Fix it with a log transform of Y or with robust standard errors.

Are the residuals normal? The Q-Q plot should show points close to the diagonal. The Jarque-Bera test checks skewness and kurtosis together. Mild departures do not matter much at moderate sample sizes, because the central limit theorem is doing quiet work in the background. Severe skew or a few huge outliers do matter.

Are the observations independent? The Durbin-Watson statistic runs from 0 to 4, and 2 means no autocorrelation. Below about 1.5 suggests positive autocorrelation, which is common in time ordered data and makes p values far too small. Above about 2.5 suggests negative autocorrelation. If your X is time, take this seriously and consider a time series model instead.

4.9 The condition number: your overfitting alarm

Powers of X are highly correlated with each other, and the higher the degree, the worse it gets. The condition number of the design matrix measures how close the columns are to being redundant. Below 30 is comfortable. Between 30 and 100 is worth noting. Above 1000 means the coefficients are numerically unstable, and a tiny change in one data value can swing them wildly.

Centring X, which this calculator does by default, usually drops the condition number by orders of magnitude at no cost to the fit. If it is still enormous after centring, the honest reading is that your degree is too high for the amount of data you have.

4.10 Prediction, interpolation and the extrapolation cliff

Inside the range of X you measured, a well fitted polynomial predicts well. Outside it, polynomials do something no other model does quite so violently: they take off. A cubic that fits your data beautifully from 0 to 200 can predict a negative yield at 260 or a physically impossible one at 400. The higher the degree, the sooner and the harder this happens.

The rule is simple and there are no exceptions worth learning. Never predict outside the range of X in your data. If a reviewer asks what happens at higher X, the correct answer is that you would need to measure it.

4.11 Effect size, and how much the curve really matters

R squared is itself the effect size for a regression model. Cohen's benchmarks, adapted for f squared = R squared ÷ (1 − R squared), put 0.02 as small, 0.15 as medium and 0.35 as large, which correspond to R squared values of about 0.02, 0.13 and 0.26. Field data routinely beat these, and calibration data in engineering can sit above 0.999.

More useful than any benchmark is the practical question: how much does Y actually change across the range of X you studied? Multiply that difference by whatever unit matters to the reader, money, kilograms, minutes, and say it out loud. A statistically significant curvature that moves yield by 0.02 t per hectare is not a finding worth a paragraph, however small the p value.

4.12 Common mistakes this calculator prevents

• Choosing degree by R squared alone. The degree table shows adjusted R squared, AIC and BIC beside it, so the penalty is always visible.
• Dropping the linear term. The model always keeps the full hierarchy, so the parabola is free to sit where the data put it.
• Unstable high-degree arithmetic. X is centred and scaled internally, then the coefficients are converted back to the original scale, so degree 5 and 6 fits do not collapse into rounding error.
• Silent missing values. Rows where either X or Y is missing or non-numeric are dropped and counted, and the count is shown, so n is never a surprise.
• Unreported diagnostics. Normality, constant variance, independence and conditioning are all tested and shown, whether or not you go looking for them.
• Extrapolation by accident. The prediction card tells you plainly when the X you asked about is outside the range of the data.

▶ Run the analysis above to auto-fill all five examples with your own results.

When reporting polynomial regression results, four things must appear somewhere in the sentence: the degree of the polynomial, the overall model test with its degrees of freedom, R squared or adjusted R squared, and the test of the highest order term. Everything else is optional detail. A reader who has those four numbers can reconstruct your conclusion; a reader who has only R squared cannot.

Ten steps, start to finish, using the fertiliser dataset that is loaded by default.

1. Load or paste your data. Pick a sample dataset from the dropdown, or paste your own numbers into the X and Y boxes as comma-separated values. Spaces, tabs and new lines work too.
2. Name the group. The group name box is editable. Type something meaningful such as Field trial, north block. That name appears on the chart legend, in the results table and in the exported report.
3. Add more groups if you want to compare curves. Press Add group and paste a second dataset. Each group is fitted separately and drawn on the same chart.
4. Or upload a file instead. Open the Upload tab, choose a CSV or Excel file, pick the X column, then click each Y column you want. Every selected column is loaded as its own group.
5. Choose the degree. Start at 2. If you are not sure, choose Auto select and the calculator will compare degrees 1 to 6 and show you the whole table anyway.
6. Set alpha, the confidence level and decimal places. Alpha 0.05 and 95 per cent intervals are the defaults, and four decimal places suits most data. Small coefficients need five or six.
7. Leave centring on. Centring X makes the arithmetic stable at higher degrees and gives the intercept a useful meaning. The fitted curve is identical either way.
8. Press Run Polynomial Regression. Nothing is calculated before you press it, so you can set everything up in peace.
9. Read the results in order: the fitted equation, then R squared and adjusted R squared, then the p value of the highest order term, then the degree table, then the residual plots.
10. Copy your sentence or export the report. Section 5 fills in with your numbers. Use Download Doc for a plain text report or Download PDF for a printable version.

Excel can fit a polynomial three different ways, and they do not all give you the same information. A chart trendline gives you the equation and R squared in about ten seconds but no p values. LINEST gives you the full statistical output in one array formula. The Analysis ToolPak gives you a familiar regression table but only after you build the X squared column yourself. This section shows all three, with the spreadsheet on screen at every step, using the same 12 fertiliser plots used everywhere else on this page.

8.1 The functions you need

8.2 Step by step

8.3 The awkward cases

8.4 Why Excel, R and Python can disagree

All three fit the same least squares model, so with clean data they agree to many decimal places. Where they part company is in defaults and in numerical conditioning.

• Raw versus orthogonal powers. R's poly(x, 2) uses orthogonal polynomials by default, which give different looking coefficients and different individual p values, though identical fitted values, R squared and overall F. To match Excel exactly you must write poly(x, 2, raw = TRUE) or I(x^2). This is the single most common source of "my numbers do not match" emails.
• Conditioning at high degrees. Excel's LINEST is known to lose precision on high order fits with large X values, where NumPy and R stay accurate because they solve by QR decomposition rather than by inverting X transpose X. At degree 2 or 3 on centred data you will not see a difference. At degree 6 on raw years you certainly will.
• Trendline display precision. Excel's chart trendline rounds the equation for display. The underlying fit is fine; the label lies to you. Always take coefficients from LINEST, not from a chart label.
• Missing values. R and pandas drop incomplete rows for you if you ask them to. Excel refuses to run at all. This is arguably the safer behaviour.
• p values. R and statsmodels print them. Excel makes you compute them from the t or F statistic yourself, which is why so many spreadsheet regressions get reported without any.

For this page's dataset, all three routes give b₂ = −0.00012991, R² = 0.9910 and F(2, 9) = 493.58.

8.5 Charting the curve properly

1. Select A2:B13 and insert a Scatter with only markers chart. Never use a line chart, because it joins points in row order and will draw a zigzag if your X values are unsorted.
2. Add the polynomial trendline at the order you fitted, and extend the forecast zero units forward and backward. Extrapolating a polynomial on a chart is how misleading figures are born.
3. Label both axes with units. A yield axis with no units is not publishable.
4. Add a second chart, a scatter of residuals against fitted values, with a horizontal line at zero. Put it in the supplementary material if not in the paper.
5. For 300 dpi output, right click the chart, choose Save as Picture, and pick PNG. Excel's copy-paste route drops resolution.

8.6 Error messages and what they mean

R fits a polynomial regression in one line, lm(y ~ poly(x, 2, raw = TRUE)), but a bare summary() is not enough to defend a degree choice. The script below is complete and runnable. It fits the model, prints the coefficient table with confidence intervals, compares degrees 1 to 6 with proper incremental F tests, runs four assumption checks, prints a plain-English verdict, and saves one publication-ready figure at 300 dpi. It uses base R only, so there is nothing to install.

9.1 The complete script

R, base only, saves polynomial_regression_plot.png

# ============================================================
# Polynomial regression with full statistics and one figure
# StatsUnlock - https://statsunlock.com/polynomial-regression-calculator/
# Base R only. Nothing to install.
# ============================================================

set.seed(42)

# ---- 1. Your data ------------------------------------------
# Replace with your own values, or read a CSV:
#   d <- read.csv("data.csv")
#   d <- d[complete.cases(d$x, d$y), ]
#   x <- d$x ; y <- d$y
x <- c(0, 20, 40, 60, 80, 100, 120, 140, 160, 180, 200, 220)
y <- c(2.18, 2.90, 3.60, 4.23, 4.78, 5.34, 5.23, 5.36, 5.70, 5.61, 5.42, 5.23)

xlab <- "Fertiliser N (kg/ha)"
ylab <- "Yield (t/ha)"
degree <- 2
alpha  <- 0.05

keep <- is.finite(x) & is.finite(y)
x <- x[keep]; y <- y[keep]
n <- length(y)
cat(sprintf("n = %d usable pairs\n\n", n))

# ---- 2. Fit the polynomial ---------------------------------
# raw = TRUE gives the ordinary powers of x, which is what Excel,
# NumPy and every textbook formula use. Leave it out and R fits
# ORTHOGONAL polynomials instead: same curve, different coefficients.
fit <- lm(y ~ poly(x, degree, raw = TRUE))
s   <- summary(fit)
b   <- unname(coef(fit))

print(s)

# ---- 3. The fitted equation --------------------------------
eq <- sprintf("%.6g", b[1])
for (k in 1:degree) {
  eq <- paste0(eq, sprintf(" %+.6g*x^%d", b[k + 1], k))
}
cat("\nEquation:  y =", eq, "\n")

# ---- 4. Confidence intervals -------------------------------
cat("\nConfidence intervals for the coefficients:\n")
print(round(confint(fit, level = 1 - alpha), 8))

# ---- 5. Fit statistics -------------------------------------
e    <- resid(fit)
sse  <- sum(e^2)
sst  <- sum((y - mean(y))^2)
p    <- degree + 1
dfe  <- n - p
fs   <- s$fstatistic
cat(sprintf("\nR squared          = %.4f\n", s$r.squared))
cat(sprintf("Adjusted R squared = %.4f\n", s$adj.r.squared))
cat(sprintf("Residual SE (s)    = %.4f on %d df\n", s$sigma, dfe))
cat(sprintf("RMSE               = %.4f\n", sqrt(sse / n)))
cat(sprintf("F(%d, %d)           = %.2f,  p = %.4g\n",
            fs[2], fs[3], fs[1], pf(fs[1], fs[2], fs[3], lower.tail = FALSE)))
cat(sprintf("AIC = %.2f   BIC = %.2f\n",
            n * log(sse / n) + 2 * p, n * log(sse / n) + p * log(n)))

# ---- 6. Degree comparison and incremental F ----------------
cat("\nDeg     R2   AdjR2      AIC      BIC    F_inc    p_inc\n")
prev <- NULL
maxd <- min(6, n - 2)
for (d in 1:maxd) {
  md  <- lm(y ~ poly(x, d, raw = TRUE))
  sm  <- summary(md)
  ss  <- sum(resid(md)^2)
  pd  <- d + 1
  fin <- NA; pin <- NA
  if (!is.null(prev)) {
    a   <- anova(prev, md)
    fin <- a$F[2]; pin <- a$`Pr(>F)`[2]
  }
  cat(sprintf("%3d %7.4f %7.4f %8.2f %8.2f %8.3f %8.4g\n",
              d, sm$r.squared, sm$adj.r.squared,
              n * log(ss / n) + 2 * pd, n * log(ss / n) + pd * log(n), fin, pin))
  prev <- md
}

# ---- 7. Assumption checks ----------------------------------
z  <- (e - mean(e)) / sqrt(mean((e - mean(e))^2))
sk <- mean(z^3); ku <- mean(z^4) - 3
jb <- n / 6 * (sk^2 + ku^2 / 4)
dw <- sum(diff(e)^2) / sum(e^2)
g  <- e^2 / (sse / n)
bp <- 0.5 * sum((fitted(lm(g ~ poly(x, degree, raw = TRUE))) - mean(g))^2)
cn <- kappa(cbind(1, poly(scale(x), degree, raw = TRUE)), exact = TRUE)
cat(sprintf("\nJarque-Bera normality : JB = %.3f, p = %.4f\n", jb, pchisq(jb, 2, lower.tail = FALSE)))
cat(sprintf("Breusch-Pagan variance: BP = %.3f, p = %.4f\n", bp, pchisq(bp, degree, lower.tail = FALSE)))
cat(sprintf("Durbin-Watson         : DW = %.3f\n", dw))
cat(sprintf("Condition number      : %.1f (on centred and scaled x)\n", cn))

# ---- 8. Plain-English verdict ------------------------------
cat("\nVERDICT\n")
pmodel <- pf(fs[1], fs[2], fs[3], lower.tail = FALSE)
if (pmodel < alpha) {
  cat(sprintf("  The degree %d model is significant, F(%d, %d) = %.2f, p = %.4g.\n",
              degree, fs[2], fs[3], fs[1], pmodel))
  cat(sprintf("  It explains %.2f%% of the variance in %s.\n", 100 * s$r.squared, ylab))
} else {
  cat("  The model is not significant. A flat line describes these data just as well.\n")
}
phigh <- coef(s)[degree + 1, 4]
if (phigh < alpha) {
  cat(sprintf("  The highest order term is significant (p = %.4g), so degree %d is justified.\n",
              phigh, degree))
} else {
  cat(sprintf("  The highest order term is not significant (p = %.4g). Try degree %d.\n",
              phigh, degree - 1))
}
if (degree == 2) {
  vx   <- -b[2] / (2 * b[3])
  vy   <- b[1] + b[2] * vx + b[3] * vx^2
  kind <- if (b[3] < 0) "maximum" else "minimum"
  cat(sprintf("  Turning point (%s) at x = %.2f, predicted y = %.4f%s\n", kind, vx, vy,
              if (vx >= min(x) && vx <= max(x)) ""
              else "  [OUTSIDE the data range, do not report as an optimum]"))
}

# ---- 9. One publication-ready figure -----------------------
png("polynomial_regression_plot.png", width = 3300, height = 1320, res = 300)
par(mfrow = c(1, 2), mar = c(4.4, 4.4, 3, 1))
grid_x <- seq(min(x), max(x), length.out = 300)
plot(x, y, pch = 19, col = "#0f766e", cex = 1.2, xlab = xlab, ylab = ylab,
     main = sprintf("Polynomial fit, R2 = %.4f", s$r.squared))
lines(grid_x, predict(fit, data.frame(x = grid_x)), col = "#c2410c", lwd = 2.6)
abline(lm(y ~ x), col = "#94a3b8", lwd = 1.8, lty = 2)
legend("bottomright", bty = "n", cex = .85,
       legend = c("Observed", sprintf("Degree %d fit", degree), "Straight line"),
       col = c("#0f766e", "#c2410c", "#94a3b8"), pch = c(19, NA, NA), lty = c(NA, 1, 2))
plot(fitted(fit), e, pch = 19, col = "#7c3aed", cex = 1.2,
     xlab = "Fitted values", ylab = "Residual", main = "Residuals vs fitted")
abline(h = 0, col = "#64748b", lwd = 1.4)
dev.off()

cat("\nFigure saved as polynomial_regression_plot.png (300 dpi)\n")
cat("Generated by STATS UNLOCK - statsunlock.com\n")

9.2 What each part does

9.3 What the figure shows

Two panels side by side. On the left, your data as points with the fitted polynomial curve in orange and the ordinary straight line as a dashed grey reference, so the reader can see exactly what the extra term bought. On the right, residuals against fitted values with a horizontal line at zero. That right panel is the one a reviewer looks at first: a shapeless cloud means the degree is right, and any visible arc means a term is missing.

9.4 Problems and fixes

9.5 One-liners worth remembering

Python gives you two routes. numpy.polyfit is one line and returns the coefficients, but no standard errors and no p values. statsmodels gives you the full regression table in a form a journal will accept. The script below deliberately builds the statistics from NumPy and SciPy directly, so you can see exactly where every number comes from, and so it runs anywhere without statsmodels installed. It prints the same numbers as the R script on the same data.

Install with pip install numpy scipy matplotlib if you do not already have them. For the shorter route, pip install statsmodels and use the one-liner in section 10.5.

10.1 The complete script

Python 3, saves polynomial_regression_plot.png

# ============================================================
# Polynomial regression with full statistics and one figure
# StatsUnlock - https://statsunlock.com/polynomial-regression-calculator/
# Needs numpy, scipy and matplotlib only.
# ============================================================
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from scipy import stats

np.random.seed(42)

# ---- 1. Your data ------------------------------------------
# Replace with your own values, or read a CSV:
#   import pandas as pd
#   df = pd.read_csv("data.csv").dropna(subset=["x", "y"])
#   x, y = df["x"].to_numpy(float), df["y"].to_numpy(float)
x = np.array([0, 20, 40, 60, 80, 100, 120, 140, 160, 180, 200, 220], dtype=float)
y = np.array([2.18, 2.90, 3.60, 4.23, 4.78, 5.34, 5.23, 5.36, 5.70, 5.61, 5.42, 5.23])

xlab, ylab = "Fertiliser N (kg/ha)", "Yield (t/ha)"
degree = 2
alpha = 0.05

ok = np.isfinite(x) & np.isfinite(y)
x, y = x[ok], y[ok]
n = len(y)

def fit_poly(x, y, d):
    X = np.vander(x, d + 1, increasing=True)
    beta, *_ = np.linalg.lstsq(X, y, rcond=None)
    fitted = X @ beta
    resid = y - fitted
    p = d + 1
    dfe = n - p
    sse = float(resid @ resid)
    sst = float(((y - y.mean()) ** 2).sum())
    ssr = sst - sse
    mse = sse / dfe
    XtXinv = np.linalg.inv(X.T @ X)
    se = np.sqrt(np.diag(mse * XtXinv))
    tval = beta / se
    pval = 2 * stats.t.sf(np.abs(tval), dfe)
    r2 = 1 - sse / sst
    adj = 1 - (1 - r2) * (n - 1) / dfe
    fstat = (ssr / (p - 1)) / mse
    fp = stats.f.sf(fstat, p - 1, dfe)
    return dict(beta=beta, se=se, t=tval, p=pval, fitted=fitted, resid=resid,
                sse=sse, ssr=ssr, sst=sst, dfe=dfe, mse=mse, s=np.sqrt(mse),
                r2=r2, adj=adj, f=fstat, fp=fp, p_n=p,
                aic=n * np.log(sse / n) + 2 * p,
                bic=n * np.log(sse / n) + p * np.log(n))

m = fit_poly(x, y, degree)

# ---- 2. The fitted equation --------------------------------
terms = [f"{m['beta'][0]:.6g}"] + [
    f"{m['beta'][k]:+.6g}*x^{k}" if k > 1 else f"{m['beta'][k]:+.6g}*x"
    for k in range(1, degree + 1)]
print("=" * 62)
print(f"POLYNOMIAL REGRESSION, degree {degree}, n = {n}")
print("=" * 62)
print("Equation:  y =", " ".join(terms))

# ---- 3. Coefficient table ----------------------------------
tcrit = stats.t.ppf(1 - alpha / 2, m["dfe"])
print(f"\n{'Term':<10}{'Estimate':>14}{'SE':>13}{'t':>9}{'p':>11}"
      f"{'CI low':>14}{'CI high':>14}")
for k in range(degree + 1):
    name = "Intercept" if k == 0 else ("x" if k == 1 else f"x^{k}")
    b, s = m["beta"][k], m["se"][k]
    print(f"{name:<10}{b:>14.6g}{s:>13.6g}{m['t'][k]:>9.3f}"
          f"{m['p'][k]:>11.4g}{b - tcrit * s:>14.6g}{b + tcrit * s:>14.6g}")

# ---- 4. Fit statistics -------------------------------------
print(f"\nR squared          = {m['r2']:.4f}")
print(f"Adjusted R squared = {m['adj']:.4f}")
print(f"Residual SE (s)    = {m['s']:.4f} on {m['dfe']} df")
print(f"RMSE               = {np.sqrt(m['sse'] / n):.4f}")
print(f"F({degree}, {m['dfe']})           = {m['f']:.2f},  p = {m['fp']:.4g}")
print(f"AIC = {m['aic']:.2f}   BIC = {m['bic']:.2f}")

# ---- 5. Degree comparison and incremental F ----------------
print(f"\n{'Deg':>4}{'R2':>10}{'AdjR2':>10}{'AIC':>10}{'BIC':>10}"
      f"{'F_inc':>10}{'p_inc':>11}")
prev = None
maxd = min(6, n - 2)
for d in range(1, maxd + 1):
    md = fit_poly(x, y, d)
    if prev is None:
        finc = pinc = float("nan")
    else:
        finc = ((prev["sse"] - md["sse"]) / 1) / md["mse"]
        pinc = stats.f.sf(finc, 1, md["dfe"])
    print(f"{d:>4}{md['r2']:>10.4f}{md['adj']:>10.4f}{md['aic']:>10.2f}"
          f"{md['bic']:>10.2f}{finc:>10.3f}{pinc:>11.4g}")
    prev = md

# ---- 6. Assumption checks ----------------------------------
r = m["resid"]
skew, exkurt = stats.skew(r), stats.kurtosis(r)
jb = n / 6 * (skew ** 2 + exkurt ** 2 / 4)
jb_p = stats.chi2.sf(jb, 2)
dw = float(((np.diff(r)) ** 2).sum() / (r @ r))
g = r ** 2 / (m["sse"] / n)
Xb = np.vander(x, degree + 1, increasing=True)
gh = Xb @ np.linalg.lstsq(Xb, g, rcond=None)[0]
bp = 0.5 * float(((gh - g.mean()) ** 2).sum())
bp_p = stats.chi2.sf(bp, degree)
cond = np.linalg.cond(np.vander((x - x.mean()) / x.std(ddof=1), degree + 1, increasing=True))
print(f"\nJarque-Bera normality : JB = {jb:.3f}, p = {jb_p:.4f}")
print(f"Breusch-Pagan variance: BP = {bp:.3f}, p = {bp_p:.4f}")
print(f"Durbin-Watson         : DW = {dw:.3f}")
print(f"Condition number      : {cond:.1f} (on centred and scaled x)")

# ---- 7. Plain-English verdict ------------------------------
print("\nVERDICT")
if m["fp"] < alpha:
    print(f"  The degree {degree} model is significant, F({degree}, {m['dfe']}) = "
          f"{m['f']:.2f}, p = {m['fp']:.4g}. It explains {100 * m['r2']:.2f}% of the variance.")
else:
    print("  The model is not significant. A flat line describes these data just as well.")
if m["p"][degree] < alpha:
    print(f"  The highest order term is significant (p = {m['p'][degree]:.4g}), "
          f"so degree {degree} is justified.")
else:
    print(f"  The highest order term is not significant (p = {m['p'][degree]:.4g}). "
          f"Try degree {degree - 1}.")
if degree == 2:
    vertex = -m["beta"][1] / (2 * m["beta"][2])
    peak = np.polyval(m["beta"][::-1], vertex)
    inside = x.min() <= vertex <= x.max()
    kind = "maximum" if m["beta"][2] < 0 else "minimum"
    print(f"  Turning point ({kind}) at x = {vertex:.2f}, predicted y = {peak:.4f}"
          f"{'' if inside else '  [OUTSIDE the data range, do not report as an optimum]'}")
print("  Normality OK." if jb_p > .05 else "  Residuals depart from normality.")
print("  Constant variance OK." if bp_p > .05 else "  Variance is not constant.")
print("  Independence OK." if 1.5 < dw < 2.5 else "  Autocorrelation is present.")

# ---- 8. One publication-ready figure -----------------------
grid = np.linspace(x.min(), x.max(), 300)
fig, ax = plt.subplots(1, 2, figsize=(11, 4.4))
ax[0].scatter(x, y, s=52, color="#0f766e", zorder=3, label="Observed")
ax[0].plot(grid, np.polyval(m["beta"][::-1], grid), color="#c2410c", lw=2.4,
           label=f"Degree {degree} fit")
b1 = fit_poly(x, y, 1)["beta"]
ax[0].plot(grid, np.polyval(b1[::-1], grid), color="#94a3b8", lw=1.6, ls="--",
           label="Straight line")
ax[0].set_xlabel(xlab); ax[0].set_ylabel(ylab)
ax[0].set_title(f"Polynomial fit, R$^2$ = {m['r2']:.4f}")
ax[0].legend(frameon=False, fontsize=9); ax[0].grid(alpha=.25)
ax[1].axhline(0, color="#64748b", lw=1.2)
ax[1].scatter(m["fitted"], r, s=52, color="#7c3aed", zorder=3)
ax[1].set_xlabel("Fitted values"); ax[1].set_ylabel("Residual")
ax[1].set_title("Residuals vs fitted"); ax[1].grid(alpha=.25)
fig.suptitle("Polynomial Regression - StatsUnlock", fontsize=11)
fig.tight_layout()
fig.savefig("polynomial_regression_plot.png", dpi=300)
print("\nFigure saved as polynomial_regression_plot.png (300 dpi)")
print("Generated by STATS UNLOCK - statsunlock.com")

10.2 What each part does

10.3 What the figure shows

The same two panels as the R figure, drawn to the same 300 dpi standard: the fitted curve over your data with the straight line for comparison on the left, and residuals against fitted values on the right. Both scripts save the file under the same name, so you can run them one after the other and compare the images pixel by pixel if you want reassurance that the two languages agree.

10.4 Problems and fixes

10.5 One-liners worth remembering

Seven lookup tables you will actually reach for while reading a polynomial regression output. All critical values are two tailed unless stated.

Table 1. Critical t values for coefficient tests

Residual degrees of freedom are n minus (degree + 1). A degree 2 fit on 12 points has 9 df.

What to conclude: if the absolute t value for a coefficient exceeds the critical value in your alpha column, that coefficient differs from zero. The worked example has t = −15.61 on 9 df against a critical 2.262, so the quadratic term is very clearly real.

Table 2. Critical F values at alpha = .05

Numerator df is the number of terms in the model excluding the intercept, which equals the degree. Denominator df is n minus (degree + 1).

What to conclude: the worked example gives F(2, 9) = 493.58 against a critical 4.26, so the curve explains far more than chance. Use the df1 = 1 column for an incremental F test, because stepping up one degree adds exactly one term.

Table 3. Choosing the degree

What to conclude: minimum n is the point at which the model can be fitted at all, not the point at which it can be trusted. A rule that serves well in practice is at least 10 observations per estimated coefficient, so a quadratic wants about 30 points before you would call it well estimated.

Table 4. Interpreting R squared and adjusted R squared

What to conclude: f squared = R² ÷ (1 − R²). Judge R squared against the norms of your field, not against these labels alone. A 0.30 in survey research is impressive; a 0.30 in a calibration curve means something is broken.

Table 5. Durbin-Watson bands

What to conclude: the worked example gives DW = 2.085, comfortably inside the acceptable band. Durbin-Watson only means something when the rows have a natural order such as time or distance. On unordered data it is noise.

Table 6. Condition number bands

What to conclude: centring x before fitting typically drops the condition number by two or three orders of magnitude at no cost to the fitted curve. The worked example sits at 2.7 after centring.

Table 7. Turning it into words

What to conclude: never write p = .000. Write p < .001. And never write "approaching significance" for p = .07. Say the estimate and its confidence interval, and let the reader judge.

Eight worked examples, each with the real data, the real output, and the sentence you would write. Between them they cover the textbook case, a clear result, a borderline one, a null one, a small sample, a significant but trivial effect, an assumption failure, and a genuine cubic.

Design statement. Polynomial regression needs one continuous predictor and one continuous outcome, measured on the same units, with enough distinct X values spread across the range to reveal the shape. The curve you can detect is limited by where you sampled, not by how many observations you took. Ten points clustered at three X values cannot show a bend, however precise each measurement is.

1. State the expected shape before you collect anything. Write down whether theory predicts a peak, a plateau, or an S shape, and therefore what degree you expect. Choosing the degree after seeing the data makes every p value optimistic.
2. Choose the range of X deliberately. The range must cover the region where the bend is expected to happen. If you think yield peaks near 150 kg, sampling 0 to 100 will only ever show you a straight line.
3. Use at least six distinct X levels for a quadratic, and eight for a cubic. More levels beat more replicates at the same levels when the goal is to identify a shape.
4. Space the levels evenly, or slightly denser near the expected bend. Even spacing keeps the design matrix well conditioned. Heavy clustering at one end inflates the standard errors on the higher powers.
5. Randomise the assignment order. If plots, animals or samples are processed in X order, any drift over time gets absorbed into the curve and you will never know.
6. Replicate at each level where you can. Replication gives you a pure error estimate, which lets you run a formal lack of fit test rather than only eyeballing the residuals.
7. Measure X as precisely as Y. Ordinary least squares assumes X is known without error. Noisy X biases polynomial coefficients toward zero and flattens the curve you are trying to find.
8. Record X in units that keep the numbers modest. Years as 2020, 2021, 2022 raised to the sixth power is a numerical disaster. Record year minus 2020 instead.
9. Log every exclusion at the moment you make it. Removed points change a polynomial fit far more than they change a straight line fit, because a curve can pivot around a single influential point.
10. Keep the raw file untouched. Do all cleaning in a script or a separate sheet, so the fit can be reproduced from the original numbers.
11. Aim for about 10 observations per estimated coefficient. That is 30 for a quadratic and 40 for a cubic as a working target.
12. Pre-register the degree and the exclusion rules if the study is confirmatory. If it is exploratory, say so in the abstract.

How to lay the data out

Minimum sample. A quadratic can be fitted with 4 points and a cubic with 5, but neither can be trusted there. Treat 15 as the floor for a quadratic in real research and 25 for a cubic.

Pre-registration line you can copy: "We will fit a degree 2 polynomial regression of Y on X by ordinary least squares. The quadratic term will be tested at alpha = .05. If it is not significant, the linear model will be reported as the primary analysis."

Common mistakes at the collection stage

• Sampling only two or three X levels, which makes any curve unidentifiable.
• Sampling a range that stops before the expected turning point.
• Recording X as a category such as "low, medium, high", which cannot be squared.
• Processing samples in X order, so time drift masquerades as curvature.
• Rounding X heavily, for example to the nearest 50, which throws away the resolution the higher powers need.

This free polynomial regression calculator is designed for one continuous predictor and one continuous outcome, where a scatter plot shows a bend that a straight line cannot follow.

Use it when

• The scatter plot of Y against X shows one clear bend, a peak, or a plateau.
• The residuals from a straight line fit form a visible arc.
• You need a single equation you can quote, differentiate, and hand to a reader.
• You want the turning point: an optimum dose, the most efficient speed, the age of peak performance.
• Theory says the effect of X should weaken, stop, or reverse as X grows.
• You are building an interpolation or calibration curve inside a known range.

Do not use it when

• The outcome is binary, a count, or a proportion. Use logistic, Poisson or beta regression instead, with polynomial terms inside them if you need curvature.
• You want to predict beyond the range of your data. Polynomials diverge violently outside the fitted range.
• The shape has three or more bends. A spline or a generalised additive model will fit it better and more stably.
• Theory gives you an exact functional form, such as exponential decay or a logistic growth curve. Fit that form directly with nonlinear least squares.
• The observations are not independent, for example repeated measures on the same subject. Use a mixed model with polynomial fixed effects.
• X is time and the series is autocorrelated. Use a time series model. A polynomial trend on autocorrelated data produces p values that are far too small.

Four quick examples

Decision guide

1. Plot Y against X first. Always. If there is no visible bend, fit degree 1 and stop.
2. If there is one bend, fit degree 2 and check the quadratic term's p value.
3. If there are two bends, try degree 3 and compare adjusted R squared and BIC against degree 2.
4. If you need three or more bends, stop using polynomials and fit a spline.
5. If the outcome is not continuous, switch to the matching generalised linear model.
6. If the observations are grouped, clustered or repeated, switch to a mixed model.

The five assumptions, how each is checked, and what happens when it fails

Limitations you should state in a paper

• Polynomials extrapolate terribly. Beyond the range of X the curve diverges. If prediction outside the range matters, use a model whose shape is bounded by theory, such as a logistic or asymptotic regression.
• Coefficients are not individually interpretable. Only the whole curve means something. Quote the marginal effect at named X values, or the turning point, instead.
• Degree selection inflates significance. If the degree was chosen after seeing the data, the reported p values are optimistic. Cross validation or a hold-out sample gives an honest error estimate.
• The fit is global, not local. One coefficient set has to describe the whole range, so a feature at one end distorts the fit at the other. Splines and GAMs do not have this problem.
• High degrees are numerically fragile. Above degree 4 on uncentred data, rounding error can dominate. Centring helps, orthogonal polynomials help more.
• Influential points matter more than in linear regression. A single point at the edge of the X range can pivot the whole curve. Always refit without it as a sensitivity check.
• A significant curve is not a causal claim. Polynomial regression describes shape, nothing more. Causal language needs a design that supports it.

A polynomial regression calculator earns its place the moment a scatter plot stops looking straight. Adding powers of X to an ordinary least squares model lets you follow a curve while keeping everything that makes linear regression easy to work with: a closed form solution, exact standard errors, familiar t and F tests, and one equation you can write on a single line of a paper.

The whole skill is degree selection. Degree 2 answers most real research questions, because the questions people actually ask, where is the optimum, when do returns start to diminish, at what age does performance peak, are all questions about a single bend. Degree 3 is right when growth accelerates and then flattens. Beyond that you are almost always fitting noise, and the honest move is a spline. Judge the degree with adjusted R squared, BIC, the incremental F test, and the residual plot, never with R squared alone, because R squared cannot fall when a term is added and will happily reward a model that has memorised your data.

Once the degree is settled, read the model as a shape rather than as a list of coefficients. No single coefficient in a polynomial means anything by itself. Quote the marginal effect at a named X value, or quote the turning point and say plainly whether it falls inside the range you measured. Check the four diagnostics, report them even when they pass, and never predict beyond the edge of your data.

Every number this polynomial regression calculator produces can be reproduced in Excel, R and Python using the scripts in sections 8, 9 and 10, and the reference tables in section 11 let you check any critical value by hand. That is the point: a result you can defend, in a form you can hand to a reviewer.

If this polynomial regression calculator produced numbers that appear in your work, cite it the way you would cite any other piece of analysis software.

Run the analysis and the Methods paragraph fills in with your own degree, n, R squared and test statistics.

Eight Stats Unlock calculators that pair naturally with polynomial regression, from the straight line you should always fit first to the models you move to when the assumptions fail.

Browse the full set under Statistical Models on Stats Unlock.

Twenty peer-reviewed journal articles underpinning this polynomial regression calculator, covering least squares curve fitting, degree and model selection, fractional polynomials, residual diagnostics, and the correct reporting of regression results.

1. Box, G. E. P., & Wilson, K. B. (1951). On the experimental attainment of optimum conditions. Journal of the Royal Statistical Society: Series B, 13(1), 1–45. https://doi.org/10.1111/j.2517-6161.1951.tb00067.x
2. Royston, P., & Altman, D. G. (1994). Regression using fractional polynomials of continuous covariates: Parsimonious parametric modelling. Journal of the Royal Statistical Society: Series C, 43(3), 429–467. https://doi.org/10.2307/2986270
3. Royston, P., Ambler, G., & Sauerbrei, W. (1999). The use of fractional polynomials to model continuous risk variables in epidemiology. International Journal of Epidemiology, 28(5), 964–974. https://doi.org/10.1093/ije/28.5.964
4. Altman, D. G., & Royston, P. (2006). The cost of dichotomising continuous variables. BMJ, 332(7549), 1080. https://doi.org/10.1136/bmj.332.7549.1080
5. Breusch, T. S., & Pagan, A. R. (1979). A simple test for heteroscedasticity and random coefficient variation. Econometrica, 47(5), 1287–1294. https://doi.org/10.2307/1911963
6. Durbin, J., & Watson, G. S. (1951). Testing for serial correlation in least squares regression. II. Biometrika, 38(1/2), 159–177. https://doi.org/10.2307/2332325
7. Jarque, C. M., & Bera, A. K. (1980). Efficient tests for normality, homoscedasticity and serial independence of regression residuals. Economics Letters, 6(3), 255–259. https://doi.org/10.1016/0165-1765(80)90024-5
8. Akaike, H. (1974). A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19(6), 716–723. https://doi.org/10.1109/TAC.1974.1100705
9. Schwarz, G. (1978). Estimating the dimension of a model. The Annals of Statistics, 6(2), 461–464. https://doi.org/10.1214/aos/1176344136
10. Burnham, K. P., & Anderson, D. R. (2004). Multimodel inference: Understanding AIC and BIC in model selection. Sociological Methods & Research, 33(2), 261–304. https://doi.org/10.1177/0049124104268644
11. Stone, M. (1974). Cross-validatory choice and assessment of statistical predictions. Journal of the Royal Statistical Society: Series B, 36(2), 111–147. https://doi.org/10.1111/j.2517-6161.1974.tb00994.x
12. Ramsey, J. B. (1969). Tests for specification errors in classical linear least-squares regression analysis. Journal of the Royal Statistical Society: Series B, 31(2), 350–371. https://doi.org/10.1111/j.2517-6161.1969.tb00796.x
13. Shapiro, S. S., & Wilk, M. B. (1965). An analysis of variance test for normality (complete samples). Biometrika, 52(3/4), 591–611. https://doi.org/10.2307/2333709
14. White, H. (1980). A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity. Econometrica, 48(4), 817–838. https://doi.org/10.2307/1912934
15. Long, J. S., & Ervin, L. H. (2000). Using heteroscedasticity consistent standard errors in the linear regression model. The American Statistician, 54(3), 217–224. https://doi.org/10.1080/00031305.2000.10474549
16. Cook, R. D. (1977). Detection of influential observation in linear regression. Technometrics, 19(1), 15–18. https://doi.org/10.2307/1268249
17. Marquardt, D. W. (1970). Generalized inverses, ridge regression, biased linear estimation, and nonlinear estimation. Technometrics, 12(3), 591–612. https://doi.org/10.2307/1267205
18. Narula, S. C. (1979). Orthogonal polynomial regression. International Statistical Review, 47(1), 31–36. https://doi.org/10.2307/1403201
19. Eilers, P. H. C., & Marx, B. D. (1996). Flexible smoothing with B-splines and penalties. Statistical Science, 11(2), 89–121. https://doi.org/10.1214/ss/1038425655
20. Anscombe, F. J. (1973). Graphs in statistical analysis. The American Statistician, 27(1), 17–21. https://doi.org/10.1080/00031305.1973.10478966