Machine Learning Regression 2

Import library

library(tidyverse)
library(MASS)
library(gridExtra)
library(psych)
library(mgcv)
library(gratia)
library(rpart)
library(rpart.plot) 
library(vip)
library(randomForest)
library(kernlab)
library(Metrics)

Load dataset

df <- read_csv("data.csv")

Rows: 466 Columns: 4── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
dbl (4): Pysics, Science, Statistics, Math
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(df)

glimpse(df)

Rows: 466
Columns: 4
$ Pysics     <dbl> 64, 74, 60, 84, 80, 75, 66, 77, 70, 89, 73, 78, 67, 80, 72, 68, 81, 64, 79, 70, 76, 71, 81, 85, 61, 74, 85, 60, 83, 85, …
$ Science    <dbl> 67, 74, 59, 88, 88, 68, 59, 71, 71, 84, 76, 74, 63, 81, 75, 63, 89, 66, 81, 76, 80, 64, 77, 83, 67, 71, 82, 60, 88, 80, …
$ Statistics <dbl> 69, 63, 55, 89, 78, 68, 55, 66, 60, 80, 69, 63, 67, 83, 73, 76, 80, 59, 79, 69, 68, 75, 84, 87, 70, 75, 78, 57, 80, 76, …
$ Math       <dbl> 68, 67, 57, 91, 82, 63, 65, 73, 67, 85, 82, 67, 57, 74, 83, 71, 85, 69, 87, 71, 80, 66, 73, 90, 61, 73, 84, 52, 92, 87, …

Exploratory Analysis

Summary Statistics

describe(df)[c(2,3,4,5,8,9,10)]

Data Visualization

Correlation matrix

pairs.panels(df, method = "pearson", hist.col = "#00AFBB")

Data distribution

p <- list()
for (name in colnames(df)){
  p[[name]] <- ggplot(data = df, aes_string(x=name)) + 
                geom_histogram(aes(y=..density..), binwidth = 10, colour="black", fill="white") + 
                geom_density(alpha=.2, fill="#FF6666") 
}
grid.arrange(p[[1]], p[[2]], p[[3]], p[[4]], ncol=2)

apply(X = df, MARGIN = 2, FUN = shapiro.test)

$Pysics

    Shapiro-Wilk normality test

data:  newX[, i]
W = 0.8713, p-value < 2.2e-16


$Science

    Shapiro-Wilk normality test

data:  newX[, i]
W = 0.90791, p-value = 3.374e-16


$Statistics

    Shapiro-Wilk normality test

data:  newX[, i]
W = 0.95389, p-value = 6.929e-11


$Math

    Shapiro-Wilk normality test

data:  newX[, i]
W = 0.96325, p-value = 2.132e-09

Scatter plot

p <- list()
var_list <- names(df)[-4]
for(name in var_list){
  p[[name]] <- ggplot(data=df, aes_string(x=name, y="Math")) + geom_point() + geom_smooth(method=lm, formula = y ~ x)
}

grid.arrange(p[[1]], p[[2]], p[[3]], ncol=2)

Discover nonlinear pattern

grid.arrange(draw(gam(Math ~ s(Pysics, k = 50), data = df, method = "REML"), residuals = TRUE),
             draw(gam(Math ~ s(Science, k = 50), data = df, method = "REML"), residuals = TRUE),
             draw(gam(Math ~ s(Statistics, k = 50), data = df, method = "REML"), residuals = TRUE),
ncol=2)

Modelling

Split train-test

set.seed(123)
index    <- createDataPartition(df$Math, p = 0.8, list = FALSE)
df_train <- df[index,]
df_test  <- df[-index,]

p <- list()
for (name in colnames(df)){
  p[[name]] <- ggplot() + 
               geom_density(data=df_train, aes_string(x=name),alpha=.2, color="blue") +
               geom_density(data=df_test, aes_string(x=name),alpha=.2, color="red")
}
grid.arrange(p[[1]], p[[2]], p[[3]], p[[4]], ncol=2)

Regression Model

m.reg <- lm(Math ~., data = df_train)
summary(m.reg)

Call:
lm(formula = Math ~ ., data = df_train)

Residuals:
    Min      1Q  Median      3Q     Max 
-50.954  -7.878  -0.254   7.956  27.443 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 52.415608   4.820313  10.874  < 2e-16 ***
Pysics      -0.004225   0.036979  -0.114   0.9091    
Science      0.090978   0.040450   2.249   0.0251 *  
Statistics   0.213926   0.049076   4.359 1.69e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 11.18 on 370 degrees of freedom
Multiple R-squared:  0.06903,   Adjusted R-squared:  0.06148 
F-statistic: 9.145 on 3 and 370 DF,  p-value: 7.492e-06

Stepwise Regression (Variable Selection)

m.stepreg <- stepAIC(m.reg, direction = "both", trace = TRUE)

Start:  AIC=1810.02
Math ~ Pysics + Science + Statistics

             Df Sum of Sq   RSS    AIC
- Pysics      1      1.63 46283 1808.0
<none>                    46282 1810.0
- Science     1    632.76 46914 1813.1
- Statistics  1   2376.82 48658 1826.8

Step:  AIC=1808.04
Math ~ Science + Statistics

             Df Sum of Sq   RSS    AIC
<none>                    46283 1808.0
+ Pysics      1      1.63 46282 1810.0
- Science     1    631.37 46915 1811.1
- Statistics  1   2391.08 48674 1824.9

summary(m.stepreg)

Call:
lm(formula = Math ~ Science + Statistics, data = df_train)

Residuals:
    Min      1Q  Median      3Q     Max 
-50.952  -7.830  -0.204   7.954  27.437 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 52.18010    4.35144  11.991  < 2e-16 ***
Science      0.09065    0.04030   2.250   0.0251 *  
Statistics   0.21332    0.04873   4.378 1.56e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 11.17 on 371 degrees of freedom
Multiple R-squared:  0.069, Adjusted R-squared:  0.06398 
F-statistic: 13.75 on 2 and 371 DF,  p-value: 1.74e-06

Generalized Additive Model (GAM)

m.gam <- gam(Math ~ s(Pysics, k = 50) + s(Science, k = 50) + s(Statistics, k = 50), 
          data = df, method = "REML")
summary(m.gam)

Family: gaussian 
Link function: identity 

Formula:
Math ~ s(Pysics, k = 50) + s(Science, k = 50) + s(Statistics, 
    k = 50)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  74.5408     0.4985   149.5   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                edf Ref.df     F  p-value    
s(Pysics)     1.699  2.125 0.652   0.5786    
s(Science)    2.227  2.820 3.467   0.0158 *  
s(Statistics) 3.497  4.391 8.064 1.86e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.113   Deviance explained = 12.7%
-REML = 1769.5  Scale est. = 115.82    n = 466

draw(m.gam, residuals = TRUE)

Decision Tree

m.tree <- rpart(formula = Math ~.,
                data = df_train,
                method = "anova")
rpart.plot(m.tree)

Random Forest

m.rf <- randomForest(Math ~ ., data = df_train)
m.rf

Call:
 randomForest(formula = Math ~ ., data = df_train) 
               Type of random forest: regression
                     Number of trees: 500
No. of variables tried at each split: 1

          Mean of squared residuals: 125.8585
                    % Var explained: 5.31

Support Vector Machine (SVM)

m.svm <- ksvm(data = df_train, Math ~., kernel = "rbfdot")
m.svm

Support Vector Machine object of class "ksvm" 

SV type: eps-svr  (regression) 
 parameter : epsilon = 0.1  cost C = 1 

Gaussian Radial Basis kernel function. 
 Hyperparameter : sigma =  0.738963400401823 

Number of Support Vectors : 349 

Objective Function Value : -204.2063 
Training error : 0.735014 

Variable Importance

grid.arrange(arrangeGrob(vip(m.reg), top ="Linear Regression"), 
             arrangeGrob(vip(m.stepreg), top ="Stepwise Regression"), 
             arrangeGrob(vip(m.tree), top ="Decision Tree"),
             arrangeGrob(vip(m.rf), top ="Random Forest"),
             ncol=2)

Prediction

x_test <- df_test[,-4]
y_test <- as.numeric(as.matrix(df_test[,4]))
predict.reg     <- predict(m.reg,x_test)
predict.stepreg <- predict(m.stepreg,x_test)
predict.gam     <- predict(m.gam,x_test)
predict.tree    <- predict(m.tree,x_test)
predict.rf      <- predict(m.rf,x_test)
predict.svm     <- predict(m.svm,x_test)

Model Evaluation

df_pred <- data.frame(actual=y_test, reg=predict.reg, stepreg=predict.stepreg,
                      gam=predict.gam, tree=predict.tree, rf=predict.rf,
                      svm=predict.svm)
df_pred

RMSE

apply(X = df_pred[,-1], MARGIN = 2, FUN = rmse, actual = df_pred$actual) %>% sort()

      rf      gam  stepreg      reg     tree      svm 
10.16479 10.31009 10.65697 10.66888 10.68761 10.74176 

Actual vs Predicted plots

p1 <- ggplot(data=df_pred) +
        geom_point(aes(x=reg, y=actual)) +
        geom_abline(intercept=0, slope=1, color="red", size=2) +
        xlab("predicted") + ggtitle("Linear Regression") +
        xlim(25,100) + ylim(25,100)
p2 <- ggplot(data=df_pred) +
        geom_point(aes(x=stepreg, y=actual)) +
        geom_abline(intercept=0, slope=1, color="red", size=2) +
        xlab("predicted") + ggtitle("Stepwise Regression") +
        xlim(25,100) + ylim(25,100)
p3 <- ggplot(data=df_pred) +
        geom_point(aes(x=gam, y=actual)) +
        geom_abline(intercept=0, slope=1, color="red", size=2) +
        xlab("predicted") + ggtitle("GAM") +
        xlim(25,100) + ylim(25,100)
p4 <- ggplot(data=df_pred) +
        geom_point(aes(x=tree, y=actual)) +
        geom_abline(intercept=0, slope=1, color="red", size=2) +
        xlab("predicted") + ggtitle("Decision Tree") +
        xlim(25,100) + ylim(25,100)
p5 <- ggplot(data=df_pred) +
        geom_point(aes(x=rf, y=actual)) +
        geom_abline(intercept=0, slope=1, color="red", size=2) +
        xlab("predicted") + ggtitle("Random Forest") +
        xlim(25,100) + ylim(25,100)
p6 <- ggplot(data=df_pred) +
        geom_point(aes(x=svm, y=actual)) +
        geom_abline(intercept=0, slope=1, color="red", size=2) +
        xlab("predicted") + ggtitle("SVM") +
        xlim(25,100) + ylim(25,100)
grid.arrange(p1,p2,p3,p4,p5,p6,ncol=3)