# Multiple linear regression and model build in light of regression diagnostics

I have a dataset of approx. 200 observations, consisting of Profit which is my dependent variable and is continuous, and the independent variables are Turnover (also continuous), and 3 additional nominal variables, let's call them Ownership (who is dichotomous), Item (who takes 6 possible values), and Mark (with 3 possible values).

My goal is to identify the combinations of these variables that give rise to relatively "higher" values of Profit (i.e. which constellations provide better profits). As 3 out of 4 explanatory variables are categorical, I chose the multiple linear regression approach with interactions (but in parallel I'm also exploring the option of clustering of mixed type data). Using the following commands in R (notice that I'm using the square of the turnover as it helped with the linear relationship part of the model)

reg_6 <- lm(Profit ~ Turnover_sq + Ownership + Item + Mark
+ Turnover_sq*Ownership + Turnover_sq*Item + Turnover_sq*Mark, data = myData)
summary(reg_6)
plot(reg_6, main = "Reg_6")


the output I get back is:

Call:
lm(formula = Profit ~ Turnover_sq + Ownership + Item + Mark +
Turnover_sq * Ownership + Turnover_sq * Item + Turnover_sq *
Mark, data = myData)

Residuals:
Min       1Q   Median       3Q      Max
-1279149   -28941   -12803    13844  1014036

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)                         8.330e+04  8.692e+04   0.958  0.33916
Turnover_sq                         2.551e-10  5.321e-09   0.048  0.96182
OwnershipP                         -4.822e+04  8.015e+04  -0.602  0.54816
ItemLok.Ind.Log.Aff                 1.224e+04  3.851e+04   0.318  0.75096
ItemMix P,F,G                       3.092e+05  1.536e+05   2.014  0.04554 *
ItemPark                            2.497e+05  1.276e+05   1.956  0.05197 .
ItemSk., Offent. Fast.             -1.871e+04  8.333e+04  -0.224  0.82263
ItemVäg & Gata                      5.721e+03  1.314e+05   0.044  0.96532
MarkM&V                             2.073e+04  5.038e+04   0.411  0.68120
MarkV                               6.236e+00  4.343e+04   0.000  0.99989
Turnover_sq:OwnershipP              7.323e-08  1.527e-08   4.796 3.37e-06 ***
Turnover_sq:ItemLok.Ind.Log.Aff    -4.994e-08  1.733e-08  -2.882  0.00443 **
Turnover_sq:ItemMix P,F,G          -1.955e-08  1.291e-09 -15.145  < 2e-16 ***
Turnover_sq:ItemPark               -7.282e-09  8.969e-10  -8.119 7.03e-14 ***
Turnover_sq:ItemSk., Offent. Fast. -3.749e-09  4.104e-09  -0.913  0.36220
Turnover_sq:ItemVäg & Gata         -1.697e-08  9.789e-09  -1.734  0.08467 .
Turnover_sq:MarkM&V                 1.471e-08  5.295e-09   2.777  0.00606 **
Turnover_sq:MarkV                   2.292e-08  1.065e-08   2.153  0.03266 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 213100 on 181 degrees of freedom
Multiple R-squared:  0.9061,    Adjusted R-squared:  0.8973
F-statistic: 102.7 on 17 and 181 DF,  p-value: < 2.2e-16


and the output of the regression diagnostics plot(reg_6) is

Given this setting, I see a plethora of problems with the setting I have chosen, and my question(s) is on how to potentially address the issues listed here:

1. For starters, the linear relationship part of the model seems to be sketchy at best (although maybe that's just my subjective interpretation of the first graph presented, so what other methods or data transformations might be suggested for this?

2. Additionally, there seems to be a wacky structure present in the standardized residuals, indicating heteroscedasticity, and last but not least,

3. The residuals are decidedly heavy-tailed (I have come across various questions here on CV where this may or may not constitute a problem, but the previous two items certainly do)

A final reason I am inclined to significantly doubt the potential predictive abilities of this model is that a 10-fold cross-validation

set.seed(123)
train.control <- trainControl(method = "cv", number = 10)
# Train the model
model <- train(Profit ~ Turnover_sq + Ownership + Item + Mark
+ Turnover_sq*Ownership + Turnover_sq*Item + Turnover_sq*Mark,
data = myData, method = "lm", trControl = train.control)
# Summarize results
print(model)


results in quite a drop in $$R^2$$ (I have omitted the " prediction from a rank-deficient fit may be misleading" error):

Linear Regression

199 samples
4 predictor

No pre-processing
Resampling: Cross-Validated (10 fold)
Summary of sample sizes: 179, 179, 179, 179, 180, 179, ...
Resampling results:

RMSE     Rsquared   MAE
4188959  0.5618473  1076385

Tuning parameter 'intercept' was held constant at a value of TRUE


Any input, guidance, or suggestion as to how (or even if) this modelling approach can be improved is welcome.

• When using the squared value you should also retain the linear term otherwise you are making a strong functional assumption. Centering prior to squaring is also a good idea and perhaps rescale it too. – Robert Long Jul 14 '20 at 11:07