data(mtcars)
model_full <- lm(mpg ~ wt + hp + disp, data = mtcars)

# optional AIC-based stepwise selection
# direction: "both" (default), "forward", or "backward"
model <- MASS::stepAIC(model_full, direction = "both", trace = FALSE)
AIC(model_full); AIC(model)

summary(model)$coefficients          # estimate, std. error, t value, p-value
confint(model, level = 0.95)         # 95% confidence intervals
effectsize::standardize_parameters(model, method = "refit")  # standardised betas
anova(model)                         # sequential (Type I) ANOVA

sm  <- summary(model)
res <- residuals(model)
sm$r.squared; sm$adj.r.squared; sm$sigma   # R2, adjusted R2, residual std. error
sm$fstatistic                              # F, df1, df2 (overall model test)
sqrt(mean(res^2))                          # RMSE
mean(abs(res))                             # MAE
AIC(model); BIC(model)
DAAG::press(model)                         # leave-one-out PRESS statistic

car::vif(model)                 # multicollinearity (VIF / GVIF)
car::ncvTest(model)             # non-constant error variance (heteroscedasticity)

# normality of residuals: Shapiro-Wilk up to n = 5000, Jarque-Bera beyond
r <- residuals(model)
if (length(r) > 5000) moments::jarque.test(as.numeric(r)) else shapiro.test(r)

lmtest::dwtest(model, alternative = "two.sided")  # autocorrelation (Durbin-Watson)
acf(residuals(model), plot = FALSE)               # ACF of residuals (n >= 10)

v <- car::vif(model)
if (is.matrix(v)) v[, ncol(v)]^2 else v   # values on the VIF scale

n <- nobs(model); p <- length(coef(model))
hatvalues(model)        # leverage        (flag > 2p/n)
cooks.distance(model)   # Cook's distance (flag > 4/n)
rstandard(model)        # standardised residuals
rstudent(model)         # studentised residuals (|.| > 3 extreme)
dffits(model)           # DFFITS

# a transformation is applied on the model's left-hand side, e.g. log(mpg) ~ ...
model_log <- lm(log(mpg) ~ wt + hp, data = mtcars)

# predictions are returned to the original response scale for reporting;
# note that exp() of a log-scale prediction estimates the MEDIAN response
pred_log <- predict(model_log, newdata = mtcars)
exp(pred_log)

newdata <- data.frame(wt = 3.0, hp = 120, disp = 160)
predict(model, newdata)                                  # fitted value
predict(model, newdata, interval = "confidence", level = 0.95)  # mean response
# the regression-line plot displays the 95% confidence band around the fit

df <- na.omit(mtcars)
n  <- nrow(df)

set.seed(123)                              # reproducible split
train_pct  <- 80                           # user-set, clamped to 50-95
idx        <- sample(seq_len(n), size = floor((train_pct / 100) * n))
trainSet   <- df[idx, ]
testSet    <- df[-idx, ]

fit   <- lm(mpg ~ wt + hp + disp, data = trainSet)
preds <- predict(fit, newdata = testSet)   # back-transformed if Y was transformed

actual <- testSet$mpg
ss_res <- sum((actual - preds)^2)
ss_tot <- sum((actual - mean(actual))^2)

sqrt(mean((actual - preds)^2))   # RMSE
mean((actual - preds)^2)         # MSE
ss_res                           # test-set SSE (not leave-one-out PRESS)
1 - ss_res / ss_tot              # test R-squared