# How to avoid overfitting in multiple linear regression (lm) in R, as the polynomial degree increases so do the fit stats

I am fitting multiple linear regressions to a data set in which the fitted plane has to approximate the values equal to 99; reddish-orangeish dots in the figure below. I am testing fits with a polynomial degree of 1 to 7.

My problem is that using the common performance metrics of a linear model, the higher the polynomial degree the better the fit (obviously), but it (i) does not make sense scientifically, and (ii) the model (fitted surface) cannot be used afterwards because of its shape (it curves back in places). To be able to use the function of the plane afterward I would need a smooth monotonic surface tilting in a given direction like this one: .

The first degree model provides just this, but has the worst fit regarding AIC, BIC, adj. R2 etc.

The data is here (mdata_nem_interp2) and the model code I use in R is this one.


library(performance)
library(stats)

centrum=99

Weight_nem_interp=1/(centrum-mdata_nem_interp2\$value)^2
Weight_nem_interp[mapply(is.infinite, Weight_nem_interp)] <- 1

per_lm1 <- lm(variable ~ polym(Unc, Window, degree=1, raw=TRUE),data = mdata_nem_interp2, weights = Weight_nem_interp)
per_lm2 <- lm(variable ~ polym(Unc, Window, degree=2, raw=TRUE),data = mdata_nem_interp2, weights = Weight_nem_interp)
per_lm3 <- lm(variable ~ polym(Unc, Window, degree=3, raw=TRUE),data = mdata_nem_interp2, weights = Weight_nem_interp)
per_lm4 <- lm(variable ~ polym(Unc, Window, degree=4, raw=TRUE),data = mdata_nem_interp2, weights = Weight_nem_interp)
per_lm5 <- lm(variable ~ polym(Unc, Window, degree=5, raw=TRUE),data = mdata_nem_interp2, weights = Weight_nem_interp)
per_lm6 <- lm(variable ~ polym(Unc, Window, degree=6, raw=TRUE),data = mdata_nem_interp2, weights = Weight_nem_interp)
per_lm7 <- lm(variable ~ polym(Unc, Window, degree=7, raw=TRUE),data = mdata_nem_interp2, weights = Weight_nem_interp)

#Evaluation

#Get the performance metrics
model_performance(per_lm1)
# Indices of model performance

AIC       |       BIC |    R2 | R2 (adj.) |   RMSE |  Sigma
-----------------------------------------------------------
56882.588 | 56907.695 | 0.282 |     0.282 | 76.943 | 18.289

> model_performance(per_lm2)
# Indices of model performance

AIC       |       BIC |    R2 | R2 (adj.) |   RMSE |  Sigma
-----------------------------------------------------------
56824.623 | 56868.562 | 0.293 |     0.293 | 56.366 | 18.147

> model_performance(per_lm3)
# Indices of model performance

AIC       |       BIC |    R2 | R2 (adj.) |   RMSE |  Sigma
-----------------------------------------------------------
56801.164 | 56870.210 | 0.299 |     0.298 | 52.652 | 18.084

> model_performance(per_lm4)
# Indices of model performance

AIC       |       BIC |    R2 | R2 (adj.) |   RMSE |  Sigma
-----------------------------------------------------------
56790.610 | 56891.041 | 0.303 |     0.300 | 45.092 | 18.048

> model_performance(per_lm5)
# Indices of model performance

AIC       |       BIC |    R2 | R2 (adj.) |   RMSE |  Sigma
-----------------------------------------------------------
56786.532 | 56924.623 | 0.306 |     0.302 | 47.961 | 18.025

> model_performance(per_lm6)
# Indices of model performance

AIC       |       BIC |    R2 | R2 (adj.) |   RMSE |  Sigma
-----------------------------------------------------------
56775.245 | 56957.276 | 0.310 |     0.305 | 59.946 | 17.984

> model_performance(per_lm7)
# Indices of model performance

AIC       |       BIC |    R2 | R2 (adj.) |   RMSE |  Sigma
-----------------------------------------------------------
56755.145 | 56987.390 | 0.316 |     0.310 | 62.748 | 17.920



My question(s),

• (i) which alternative metric(s) could I use to evaluate the performance of the models from a different aspect, which gives a penalty to overfitting,
• (ii) should I try to use a different model fit?
• (iii) is there any literature that I could use to still argue why I have to use the linear fit?

Thanks for any suggestions.

1. "Period" is uniformly distributed between $$P_{min}$$ and $$200$$
2. The distrubution parameter $$P_{min}$$ varies with the predictors
I would thus recommend to build a hand-crafted model according to 1. and 2. and to predict $$P_{min}$$ as a function of the two other variables. As the ML estimator for the lower bound of a uniform distribution is just the smallest observed value, I would guess (but have not proved!) that an equivalent approach would be to only use the bottom contour values of "Period" and do a linear fit with only these values as response.