Rows with Missing Data

Is there any observations with a high number of missing observations?

# include any rows where any column is NA.  Always include id column
sdata_v1 %>% dplyr::filter_all( any_vars(is.na(.))) -> x1  

# only include those columns where some row has NA in that column
# this dataset has no id column / so no primary key
x1 %>% select_if( ~ any(is.na(.) ) ) -> x2

# But the order of the missing observations is same in x1 and x2
# so add back the id and put it on the left side of tibble
x2 %>% mutate( id = x1$id  ) %>% relocate(id)-> missing_train_rows
  
dim(missing_train_rows)

## [1] 529  28

We find there are 529 observations with NA observations. So 20.6% have NA values. We conclude that the NA values are sufficiently well distributed that removal of all observations with NA values is impractical.

In the marginal table below, we investigate if certain observations have a high incidence of NA values amongst their predictors. The answer shows that the incident of multiple

Number of Rows with NA

## [1] "There are 812 cells with NA values out of 82144 equivalent to: 0.99%"

Overall, only 1 percent of cells have missing values.
The concentration of missing values seems to decline at a geometric ratio suggesting independence among the occurrence of missing values.

Scatterplot of Predictor Variables

The scatterplots of the predictor variables against the dependent variable PH indicate no clear linear relationship with any predictor variable. The scatterplots indicate some variable measurements are not continuous and instead fall on specific values on the x-axis.

Density Plots

The density plots of each predictor variable show many predictor variables with non-normal distributions.

Normality Test

The Shapiro-Wilk test of normality is performed on all the numeric predictor variables and the resulting list confirms no p-value above 0.05 and thus no predictor variable follows a normal distribution per evaluation.

## [1] 0

Distribution by Brand Code

The distribution of the Brand Code, the lone non-numeric predictor variable, shows a high count of brand B, comprising almost half of the samples. Also note, some values are missing for the Brand Code variable.

Boxplot of PH by Brand Code

The boxplot of PH values by Brand Code certainly captures overlap of the middle two quartiles, which seems reasonable given the PH range is rather narrow. Of note, the median values of the 4 labeled brands and the unlabeled group do appear distinct besides Brand A and the unlabeled.

VSURF

To find significant variables for the models, the VSURF algorithm is applied to the provided dataset. The algorithm identifies 14 variables as significant to the predictions.

• Mnf Flow
• Usage cont
• Carb Rel
• Filler Level
• Alch Rel
• Bowl Setpoint
• Balling Lvl
• Carb Flow
• Temperature
• Oxygen Filler
• Balling
• Pressure Vacuum
• Filler Speed
• Hyd Pressure1

## Thresholding step
## Estimated computational time (on one core): 12.3 sec.
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |====                                                                  |   5%
  |                                                                            
  |=======                                                               |  10%
  |                                                                            
  |==========                                                            |  15%
  |                                                                            
  |==============                                                        |  20%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |========================                                              |  35%
  |                                                                            
  |============================                                          |  40%
  |                                                                            
  |================================                                      |  45%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |======================================                                |  55%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |==============================================                        |  65%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |========================================================              |  80%
  |                                                                            
  |============================================================          |  85%
  |                                                                            
  |===============================================================       |  90%
  |                                                                            
  |==================================================================    |  95%
  |                                                                            
  |======================================================================| 100%
## Interpretation step (on 26 variables)
## Estimated computational time (on one core): between 799.5 sec. and  2356.1 sec.
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |===                                                                   |   4%
  |                                                                            
  |=====                                                                 |   8%
  |                                                                            
  |========                                                              |  12%
  |                                                                            
  |===========                                                           |  15%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |================                                                      |  23%
  |                                                                            
  |===================                                                   |  27%
  |                                                                            
  |======================                                                |  31%
  |                                                                            
  |========================                                              |  35%
  |                                                                            
  |===========================                                           |  38%
  |                                                                            
  |==============================                                        |  42%
  |                                                                            
  |================================                                      |  46%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |========================================                              |  58%
  |                                                                            
  |===========================================                           |  62%
  |                                                                            
  |==============================================                        |  65%
  |                                                                            
  |================================================                      |  69%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |======================================================                |  77%
  |                                                                            
  |=========================================================             |  81%
  |                                                                            
  |===========================================================           |  85%
  |                                                                            
  |==============================================================        |  88%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |===================================================================   |  96%
  |                                                                            
  |======================================================================| 100%
## Prediction step (on 17 variables)
## Maximum estimated computational time (on one core): 970.2 sec.
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |====                                                                  |   6%
  |                                                                            
  |========                                                              |  12%
  |                                                                            
  |============                                                          |  18%
  |                                                                            
  |================                                                      |  24%
  |                                                                            
  |=====================                                                 |  29%
  |                                                                            
  |=========================                                             |  35%
  |                                                                            
  |=============================                                         |  41%
  |                                                                            
  |=================================                                     |  47%
  |                                                                            
  |=====================================                                 |  53%
  |                                                                            
  |=========================================                             |  59%
  |                                                                            
  |=============================================                         |  65%
  |                                                                            
  |=================================================                     |  71%
  |                                                                            
  |======================================================                |  76%
  |                                                                            
  |==========================================================            |  82%
  |                                                                            
  |==============================================================        |  88%
  |                                                                            
  |==================================================================    |  94%
  |                                                                            
  |======================================================================| 100%

## ** VSURF results **
## The results object is a list of length around 20.
## Most interesting components are the following:
## 
##   name                description                                       
## 1 "$varselect.thres"  "variables selected after thesholding step"       
## 2 "$varselect.interp" "variables selected after interpretation step"    
## 3 "$varselect.pred"   "variables selected after prediction step"        
## 4 "$ord.imp$x"        "mean VI in decreasing order for all variables"   
## 5 "$ord.imp$ix"       "indices of the ordering of all variables VI mean"
## 6 "mean.perf"         "mean OOB rate for RF build on all variables"     
## 
##  For more information about VSURF outputs see the VSURF help page

## 
##  VSURF computation time: 36.2 mins 
## 
##  VSURF selected: 
##  26 variables at thresholding step (in 7.4 secs)
##  17 variables at interpretation step (in 23.4 mins)
##  16 variables at prediction step (in 12.6 mins)

##  [1]  9 19 30 16 29 26 31 20 18 23 25 24 27 17 12 28

Data Partition into Training and Test Sets

Next we partition the student evaluation data set by training and test sets. We will retain the terminology of test set to refer to the validation data set used for model performance assessment outside of the cross validation process.

set.seed(19372)
train_indices = createDataPartition( sdata_v1$PH , times = 1, list = FALSE, p = 0.80 )

sdata_test    = sdata_v1[-train_indices , ]
sdata_train   = sdata_v1[ train_indices , ]
sdataY_test   = sdata_v1$PH[-train_indices ]
sdataY_train  = sdata_v1$PH[ train_indices ]

Drop Near Zero Variance Predictors

We conclude there are no near zero variance predictors or constant predictors to be dropped from the training data set.

No Near Zero Variance Predictors in Training

Outliers and Zeros

Next we identify and fix zeros and outliers in the predictors before imputing missing values.
The reason is that KNN algorithm may be affected by outliers. The table below uses skim to report the minimum, maximum, median, mean and standard deviation for each predictor. In addition, we construct Z-score metrics that tells us how many standard deviations is the minimum and maximum.

\[ \text{zscore_min}(X_i) = \frac{ min(X_i) - \mu(X_i) }{\sigma(X_i)}\] \[\text{zscore_max}(X_i) = \frac{ max(X_i) - \mu(X_i) }{\sigma(X_i )}\]

Predictors with Outliers by Z-Score

Using scatterplots of each of the 7 identified predictors below, we see no obvious data errors exist in the distribution of the predictors with outlier values.
Isolated points exist but lie within a reasonable proximity to neighboring points in the scatter plot. Therefore no data corrections for outliers will be applied to the training data.

Handle Special Cases

The code below will transform the special case situations of data pre-processing identified earlier:

• Replace the MFR missing values with its median value from the training set
• Replace the missing BrandCode with a fictional brand code E.

These rules are applied consistently to be the training, validation and test data sets prior to subsequent pre processing.

KNN Imputation

We apply the caret pre processing function to that training, validation and test data sets.

# Build the caret function to preprocess the Chemical data and impute missing values.
# There is a bug in caret which causes tibbles to be rejected.  They need to be cast as data frames.
# ---------------------------------------------------------------------------------
preProcFunc = preProcess(as.data.frame(sdata_v2_train[,3:33]) , method = c("BoxCox", "center", "scale",  "knnImpute") )

# Becomes the source data for the model building
sdata_train_pp = predict( preProcFunc,  as.data.frame(sdata_v2_train ) )

# Becomes the final version of test data for validation
sdata_test_pp  = predict( preProcFunc,  as.data.frame(sdata_v2_test) )

# Need to generate predictions based on test data without known response
edata_test_pp  = predict( preProcFunc,  as.data.frame(edata_v2) )


sdata_trainX_pp  = sdata_train_pp %>% select( -PH)
sdata_testX_pp   = sdata_test_pp %>% select(-PH)
edata_testX_pp   = edata_test_pp %>% select( -PH)

Variable type: character

Variable type: numeric

Variable type: character

Variable type: numeric

Variable type: character

Variable type: logical

Variable type: numeric

Exporting the Post-Processed Data Sets

Our last step of preprocessing is to export the data sets. This need only but done once.

readr::write_csv(sdata_train_pp, "sdata_train_pp.csv", append = FALSE )

readr::write_csv(sdata_test_pp, "sdata_test_pp.csv", append = FALSE )

readr::write_csv(edata_test_pp, "edata_test_pp.csv", append = FALSE )

The above data files have the following properties and relationships to the original raw files:

Pre-Processed Data Files

An important note about the id column. The id in the file edata_test_pp.csv is the row number of the same observation in StudentEvaluation.xlsx. But the same is not true for sdata_train_pp.csv and sdata_test_pp.csv because they are sampled from the source file StudentData.xlsx. The 10th row of sdata_test_pp.csv is not necessarily the 10th row of StudentData.xlsx rather one should use the id field as the corresponding row number of the original file.

Importing the Pre-Processed Data

We load the pre-processed datasets for training and testing in this section. We load the csv file for the training data set below and display the first few column characteristics from the file.

sdata_train_pp = read_csv("sdata_train_pp.csv", 
                          col_types = cols( .default = "?", id = "i" , BrandCode = "f")) %>% 
  as_tibble() %>% remove_rownames()

sdata_test_pp  = read_csv("sdata_test_pp.csv", col_types = cols( .default = "?", id = "i" , BrandCode = "f")) %>% 
  as_tibble() %>% remove_rownames()

sdata_trainX_pp  = sdata_train_pp %>% select( -PH, -id )
sdata_testX_pp   = sdata_test_pp %>% select(-PH , -id)

tuningFull = TRUE   # Set to TRUE if running the lengthy tuning process for all models

Gradient Boosted Trees

We consider the same dataset and variable importance problem using gradient boosted trees. The gbm library and variable importance measures are run within caret training procedure. The gradient boosting approach relies on a weak learner in an additive manner that minimizes the given loss function.

set.seed(1027)

gbmControlFull <- trainControl(method = "repeatedcv", 
                           number = 10,  # for debug mode, set this to a low number (1)
                           repeats = 10 ,
                           selectionFunction = "best", 
                           verboseIter = FALSE)

gbmControlLight <- trainControl(method = "repeatedcv", 
                           number = 5,  # for debug mode, set this to a low number (1)
                           repeats = 1 ,
                           selectionFunction = "best", 
                           verboseIter = FALSE)

if(tuningFull == TRUE)
{
   gbmControl = gbmControlFull
   gbmTune = expand.grid( n.trees = seq(250, 1000, by = 250),
                                                  shrinkage = c( 0.02, 0.05,  0.1) ,
                                                  interaction.depth = c(  3, 7) ,
                                                  n.minobsinnode = 15
                                                  )
} else  {
   gbmControl = gbmControlLight
   gbmTune =  expand.grid( n.trees = c( 2000), 
                                                  shrinkage = c( 0.02, 0.05 , 0.1 ) ,
                                                  interaction.depth = c(  3, 7) ,
                                                  n.minobsinnode = 4
                                                  )
}

(gbmTune = caret::train( x = sdata_trainX_pp, 
                          y = sdata_train_pp$PH , 
                          method = "gbm",
                          tuneGrid =  gbmTune ,
                          verbose = FALSE,
                          metric = "RMSE" ,
                          trControl = gbmControl ) )

## Stochastic Gradient Boosting 
## 
## 2055 samples
##   32 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 10 times) 
## Summary of sample sizes: 1850, 1849, 1849, 1851, 1850, 1849, ... 
## Resampling results across tuning parameters:
## 
##   shrinkage  interaction.depth  n.trees  RMSE       Rsquared   MAE      
##   0.02       3                   250     0.1434366  0.3445322  0.1130661
##   0.02       3                   500     0.1399034  0.3781615  0.1093662
##   0.02       3                   750     0.1381810  0.3948955  0.1075668
##   0.02       3                  1000     0.1371433  0.4053153  0.1065292
##   0.02       7                   250     0.1375662  0.3998266  0.1072462
##   0.02       7                   500     0.1339808  0.4329093  0.1037283
##   0.02       7                   750     0.1326817  0.4453548  0.1025127
##   0.02       7                  1000     0.1320709  0.4515327  0.1019375
##   0.05       3                   250     0.1391827  0.3846115  0.1085754
##   0.05       3                   500     0.1370326  0.4073776  0.1063487
##   0.05       3                   750     0.1364019  0.4163811  0.1056864
##   0.05       3                  1000     0.1360931  0.4213869  0.1053766
##   0.05       7                   250     0.1337931  0.4330732  0.1032574
##   0.05       7                   500     0.1324935  0.4472241  0.1018857
##   0.05       7                   750     0.1319768  0.4531566  0.1013501
##   0.05       7                  1000     0.1317495  0.4562800  0.1010483
##   0.10       3                   250     0.1379832  0.3995082  0.1072565
##   0.10       3                   500     0.1371777  0.4135910  0.1063966
##   0.10       3                   750     0.1368283  0.4201240  0.1059677
##   0.10       3                  1000     0.1366484  0.4244918  0.1056588
##   0.10       7                   250     0.1343586  0.4341016  0.1034324
##   0.10       7                   500     0.1335822  0.4437819  0.1027288
##   0.10       7                   750     0.1332925  0.4472134  0.1023291
##   0.10       7                  1000     0.1332832  0.4482130  0.1022717
## 
## Tuning parameter 'n.minobsinnode' was held constant at a value of 15
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were n.trees = 1000, interaction.depth =
##  7, shrinkage = 0.05 and n.minobsinnode = 15.

As the model output indicates the greater tree depth outperforms the lower counts and the shrinkage value of 0.05 produces the best RMSE value. The shrinkage, or learning rate, typically works better with smaller values, yet require longer processing time.

## Saving 7 x 5 in image

The plot of the tuned model clearly shows the better performance with the increased max tree depth but also indicates the lower shrinkage value doesn’t guarantee a lower RMSE. For the max tree depth of 3 and 7, the lowest shrinkage value does not produce the lowest RMSE for that tree depth.

## Saving 7 x 5 in image

The variable importance plot shows the 100% importance of the MnfFlow followed by BrandCode as the only other variable with an importance above 50%. Usagecont, OxygenFiller, AlchRel, and Temperature round out the top 6 variables based on importance with scores greater than 25%.

Overall, the final Gradient Boosted model results in an RMSE of approximately 0.132 and an \(R^2\) less than 0.5.

Cubist

Now we train a Cubist tree using the best selectionFunction on RMSE. The Cubist library and variable importance measures are run within caret training procedure.
The Cubist model is a rules-based approach in which the terminal leaves contains linear regression models. The models are based on predictions defined in the previous nodes of the tree.

set.seed(1027)

cubistControlFull <- trainControl(method = "cv" ,  selectionFunction = "best")
tuneGridFull  <- expand.grid( committees = c( 10, 50, 100 ) ,
                              neighbors = c( 0, 1, 5, 9 )
                                                  ) 

cubistControlLight <- trainControl(method = "cv" ,  selectionFunction = "best")
tuneGridLight <- expand.grid( committees = c( 10, 20 ) , 
                              neighbors = c( 0, 5 )  )

if(tuningFull == TRUE)
{
   cubistControl = cubistControlFull
   cubistGrid = tuneGridFull
} else  {
   cubistControl = cubistControlLight
   cubistGrid = tuneGridLight
}

(cubistTune = caret::train( x = sdata_trainX_pp, 
                          y = sdata_train_pp$PH , 
                          method = "cubist",
                          tuneGrid = cubistGrid ,
                          verbose = FALSE,
                          metric = "RMSE" ,
                          trControl = cubistControl ) )

## Cubist 
## 
## 2055 samples
##   32 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 1850, 1849, 1849, 1851, 1850, 1849, ... 
## Resampling results across tuning parameters:
## 
##   committees  neighbors  RMSE        Rsquared   MAE       
##    10         0          0.10326126  0.6469501  0.07521436
##    10         1          0.10698108  0.6394174  0.07408043
##    10         5          0.09661241  0.6865169  0.06923965
##    10         9          0.09715807  0.6838488  0.06954861
##    50         0          0.10213076  0.6574062  0.07362497
##    50         1          0.10583349  0.6449193  0.07208654
##    50         5          0.09539318  0.6943587  0.06757449
##    50         9          0.09598916  0.6917813  0.06811918
##   100         0          0.10218813  0.6576706  0.07366847
##   100         1          0.10587327  0.6449503  0.07201655
##   100         5          0.09547621  0.6941947  0.06752479
##   100         9          0.09606335  0.6917983  0.06810842
## 
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were committees = 50 and neighbors = 5.

The Cubist model results provide a strong RMSE value and \(R^2\) value with generally a higher number of committees and neighbors.

The line plot confirms the better model per RMSE with neighbors at 5 and committees set to 50.

## Saving 7 x 5 in image

The variable importance plot identifies MnfFlow as the most valuable variable followed by Balling, BallingLvl, AlchRel, and PressureVacuum to round out the top 5 variables. Overall, 7 variables resulted in over 50% importance indicating those variables are contained in over 50% of the rules defining the Cubist model.

## Saving 7 x 5 in image

The final Cubist model confirms the high usage rate of the variable MnfFlow.

Cubist Predictor Usage

The splits plot indicates the usage of each variable by committee/rule. The values are normalized and highlight the high usage rate of MnfFlow along with the use of greater than or less than inequalities for MnfFlow.

## quartz_off_screen 
##                 2

The dotplot allows us to see the distribution of the coefficient parameters for each linear model associated with a committee/rule.
When a predictor has more dots, it is used in more models which suggests the variable is influential. The color density is proportional to the variable importance. The sign of the coefficient predicts the marginal effect of the predictor on the response PH. However, since Cubist uses many models, the concept of the sign of a predictor coefficient can be unclear.
The x-axis value of the dot in the dotplot indicates the true coefficient in the original units unlike the above plot which normalizes the x-axis.

## quartz_off_screen 
##                 2

The boxplot below which we devised constructs a distribution of coefficients for each predictor. We interpret the sign of the predictor \(X\) marginal effect to be the sign of the median of the distribution of coefficients pf \(X\) over those models which include \(X\). The predictors have been sorted by their median.

The plot provides clarity regarding the impact of a given variable as the median is a good estimate of the sign of the variable coefficient. The AlchRel variable box indicates a positive relationship with the dependent variable PH. The most important variable MnfFlow box results in a negative coefficient, and thus an inverse effect on PH. For all the rules in the final Cubist model, the understanding of the coefficient values for each variable provides a clear relationship to the dependent variable despite the obtuse nature of the Cubist model construction.

## Saving 7 x 5 in image

Overall, the final Cubist model results in an RMSE under 0.1 and a \(R^2\) just below of 0.7.

CART

Using the rpart library and its training methods, we consider the best CART based on the RMSE metric. The CART approach follows a classification and regression tree methodology.

set.seed(10393)

library(rpart)
cartGrid = expand.grid( maxdepth = seq( 1, 20, by = 1 )  )

cartControl <- trainControl(method = "boot", 
                            number = 30,  
                            selectionFunction = "best")

(rpartTune = train( x = sdata_trainX_pp, 
                    y = sdata_train_pp$PH ,
                   method = "rpart2" ,
                   metric = "RMSE" ,
                   tuneGrid = cartGrid ,
                   trControl = cartControl ) )

## CART 
## 
## 2055 samples
##   32 predictor
## 
## No pre-processing
## Resampling: Bootstrapped (30 reps) 
## Summary of sample sizes: 2055, 2055, 2055, 2055, 2055, 2055, ... 
## Resampling results across tuning parameters:
## 
##   maxdepth  RMSE       Rsquared   MAE      
##    1        0.1540934  0.2150028  0.1207013
##    2        0.1485972  0.2697600  0.1163125
##    3        0.1439126  0.3150026  0.1130389
##    4        0.1418976  0.3346403  0.1106321
##    5        0.1408007  0.3455111  0.1095038
##    6        0.1398046  0.3557469  0.1085050
##    7        0.1389994  0.3635874  0.1074293
##    8        0.1381955  0.3717156  0.1065335
##    9        0.1373228  0.3802125  0.1057161
##   10        0.1362384  0.3899605  0.1046780
##   11        0.1353910  0.3979319  0.1038336
##   12        0.1345705  0.4055253  0.1028965
##   13        0.1340562  0.4106769  0.1022712
##   14        0.1335638  0.4154690  0.1017284
##   15        0.1330891  0.4202569  0.1012114
##   16        0.1327996  0.4230057  0.1008957
##   17        0.1327006  0.4239743  0.1008136
##   18        0.1326532  0.4245372  0.1007342
##   19        0.1326119  0.4250877  0.1006977
##   20        0.1326119  0.4250877  0.1006977
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was maxdepth = 19.

The model output shows increasing tree depth leads to improved RMSE and \(R^2\) values.

The below plot displays the improving RMSE value with the increasing tree depth with an optimal tree depth of 19.

## Saving 7 x 5 in image

Once again, the variable MnfFlow is the most important variable with an importance rating of 100%. Six variables achieve an importance score of greater than 50% with BrandCode, Temperature, and CarbRel each achieving greater than 75% importance.

## Saving 7 x 5 in image

The CART tree plot of the final model reaffirms the importance of the variable MnfFlow as the initial node splits on MnfFlow. One of the second level nodes is BrandCode, another highly important variable to the CART model. Eight ff the top 10 most important variables in the Cubist model are found in the top 15 variables of the CART model. This overlap of variable importance shows a consistency across the different model approaches.

Overall, the final CART model results in an RMSE of approximately 0.132 and an \(R^2\) less than 0.5.