1
$\begingroup$

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: 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.

$\endgroup$
1
$\begingroup$

Assuming that "Period" in the plot is the same as your response variable "variable", the data shown in the plot does not seem to match the model that underlies lm. Your data seems to be generated by the following model:

  1. "Period" is uniformly distributed between $P_{min}$ and $200$
  2. The distrubution parameter $P_{min}$ varies with the predictors

Both points contradict the assumptions of lm: 1. contradicts the normality assumption, 2. contradicts the assumption of constant variance.

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.

$\endgroup$
3
  • $\begingroup$ Thank you for the comment. What we want to fit to are the "bottom contour values = 99" in fact. In the script you can see that the model is weighted in an inverse way to the central value (99). So the closer the 'value of a dots is to 99, the more weight they get. What we wanted to model is the plain that described the dots forming the 99% surface. The function of the plain would be used to categorize any point in the period - age uncertainty - window space. Above the surface the values would be considered OK. $\endgroup$ Jun 29 at 8:46
  • $\begingroup$ Weighting with the inverse height is not equivalent to extracting the bottom surface. As your data sems to be approximately discrete and falls within a cell raster, it should be possible to extract the bottom surface by simple rules, e.g. no points below in a certain (unc,window) neighborhood. Alternatively you can estimate the surface normal at each point (the C++ library PCL, e.g., has a routine for this) and only keep points where the z-direction of the surface normal is negative. $\endgroup$
    – cdalitz
    Jun 29 at 13:09
  • $\begingroup$ Thank you for the comments! I applied logistic polynomial multiple regressions with the predictors being Und, MSR and the period explaining the probability of successfully determining an implanted period with a shift < 5% of the period itself. The results seem ok; the quadratic multiple regression gave the most accurate estimation of the success rate according to AIC, BIC and pseudo R2. $\endgroup$ Jul 1 at 11:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.