library(glmnet)
library(varSelRF)
library(FSelector)
library(mlbench)
library(plot3D)
library(knitr)
library(rgl)
knit_hooks$set(webgl = hook_webgl)

x <- sort(rnorm(1000))
y <- rnorm(1000)
z <- rnorm(1000) + atan2(x,y)
plot3d(x, y, z, col=rainbow(1000))

Your browser does not support the HTML5 canvas element.

You must enable Javascript to view this page properly.

We chose to have 3 categories.
A Gain in rank -> Black – Improves rank
No Change in rank -> Red
A Loss in rank -> Green – Falls behind
Load the data

dat <- read.csv("C:/Users/Prashan/Dropbox (MIT)/MIT/Predictive Analytics/code/data/NASCAR_5f/allf_3c_phoenix2_2014_prototype_sel.csv")
plot3d(dat[,1], dat[,2], dat[,3], col=rainbow(1000))

Your browser does not support the HTML5 canvas element.

You must enable Javascript to view this page properly.

Start with 43 features, create the x-attributes, y-label matrices

features <- 3:45 
complete <- which(rowSums(is.na(dat[, features]))==0)

datc <- dat[complete, ]
x <- as.matrix(datc[, features])
y <- datc[,2]
cl <- factor(y)

#column one contains the names of the datapoints (i.e cereal xxx)
nascar_data_with_label=datc[,(-1)]
selected_features<-c()

Your browser does not support the HTML5 canvas element.

You must enable Javascript to view this page properly.

Variable selection from random forests using OOB error

knit_hooks$set(webgl = hook_webgl)
rf.vs1 <- varSelRF(x, cl, ntree = 200, ntreeIterat = 100,
vars.drop.frac = 0.2)
rf.vs1$selected.vars

## [1] "X157..before.pit..sqrt.of..Binned.avg.pre.slope..prev.epoch.final.rank.127"
## [2] "X34..before.pit..25th.percentile.rank.upto.bef.pit"

  if (length(rf.vs1$selected.vars)==3){
        plot3d(x[,rf.vs1$selected.vars[1]],x[,rf.vs1$selected.vars[2]],x[,rf.vs1$selected.vars[3]],col=y,xlab=rf.vs1$selected.vars[1],ylab=rf.vs1$selected.vars[2],zlab=rf.vs1$selected.vars[3])
  }else if (length(rf.vs1$selected.vars)==2){
        plot(x[,rf.vs1$selected.vars[1]],x[,rf.vs1$selected.vars[2]],col=y,xlab=rf.vs1$selected.vars[1],ylab=rf.vs1$selected.vars[2])
  }

# plot(rf.vs1)

Your browser does not support the HTML5 canvas element.

You must enable Javascript to view this page properly.

Variable selection from CFS filter (https://cran.r-project.org/web/packages/FSelector/FSelector.pdf)

knit_hooks$set(webgl = hook_webgl)
subset <- cfs(rank_change_label~., datc[,(-1)])
plot3d(x[,subset[1]],x[,subset[2]],x[,subset[3]],col=y,xlab=rf.vs1$selected.vars[1],ylab=rf.vs1$selected.vars[2],zlab=rf.vs1$selected.vars[3])
subset

## [1] "X76..after.pit..indicator.2.tire.change..ref.19"                                                       
## [2] "X86..before.pit..indicator.front..lte.8..of.the.pack.pre.pit"                                          
## [3] "X89..before.pit..Indicator.2.tire.AND.bottom..gte.23..of.pack.at.pit.entry...features.76...feature.87."
## [4] "X140..after.pit..Time.in.pits.before.outing"                                                           
## [5] "X160..before.pit..sqrt.of.Binned.avg.prev.slope.HI.RES.133"

dat <- read.csv("C:/Users/Prashan/Dropbox (MIT)/MIT/Predictive Analytics/code/data/NASCAR_5f/fold1/train/fold_1_train.csv")

plot3d(dat[,1],dat[,2],dat[,3])

Your browser does not support the HTML5 canvas element.

You must enable Javascript to view this page properly.

Variable selection from Consistency-based filter (https://cran.r-project.org/web/packages/FSelector/FSelector.pdf)

knit_hooks$set(webgl = hook_webgl)
subset <- consistency(rank_change_label~., datc[,(-1)])
plot3d(x[,subset[1]],x[,subset[2]],x[,subset[3]],col=y,xlab=subset[1],ylab=subset[2],zlab=subset[4])
subset

## [1] "X26..before.pit..incoming.rank.before.pit.bef.leg" 
## [2] "X34..before.pit..25th.percentile.rank.upto.bef.pit"
## [3] "X35..before.pit..75th.percentile.rank.upto.bef.pit"
## [4] "X82..before.pit..starting.position.of.the.car"

selected_features<-c(selected_features,subset)

Your browser does not support the HTML5 canvas element.

You must enable Javascript to view this page properly.

Variable selection from RReliefF filter (https://cran.r-project.org/web/packages/FSelector/FSelector.pdf)

knit_hooks$set(webgl = hook_webgl)
weights <- relief(rank_change_label~., datc[,(-1)], neighbours.count = 5, sample.size = 20)
subset <- cutoff.k(weights,5)
plot3d(x[,subset[1]],x[,subset[2]],x[,subset[3]],col=y,xlab=subset[1],ylab=subset[2],zlab=subset[3])
subset

## [1] "X77..after.pit..indicator.4.tire.change..ref.19"                        
## [2] "X26..before.pit..incoming.rank.before.pit.bef.leg"                      
## [3] "X139..after.pit...Binned.avg.pre.slope..pit.exit.rank.HI.RES"           
## [4] "X34..before.pit..25th.percentile.rank.upto.bef.pit"                     
## [5] "X61..before.pit..NBHD.prev.epoch..number.of.4.pit.types.of.neighborhood"

selected_features<-c(selected_features,subset)
#RandomForest filter
#1=mean decrease in accuracy, 2=mean decrease in node impurity
# weights <- random.forest.importance(rank_change_label~., datc[,(-1)], importance.type = 2)
# subset <- cutoff.biggest.diff(weights)
# subset

Your browser does not support the HTML5 canvas element.

You must enable Javascript to view this page properly.

Variable selection from Information Gain (https://cran.r-project.org/web/packages/FSelector/FSelector.pdf)

knit_hooks$set(webgl = hook_webgl)
weights <- information.gain(rank_change_label~., nascar_data_with_label)
subset <- cutoff.biggest.diff(weights)
plot3d(x[,subset[1]],x[,subset[2]],x[,subset[3]],col=y,xlab=subset[1],ylab=subset[2],zlab=subset[3])
subset

##  [1] "X26..before.pit..incoming.rank.before.pit.bef.leg"                                  
##  [2] "X82..before.pit..starting.position.of.the.car"                                      
##  [3] "X34..before.pit..25th.percentile.rank.upto.bef.pit"                                 
##  [4] "X35..before.pit..75th.percentile.rank.upto.bef.pit"                                 
##  [5] "X137..before.pit...Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES"              
##  [6] "X157..before.pit..sqrt.of..Binned.avg.pre.slope..prev.epoch.final.rank.127"         
##  [7] "X158..before.pit..square.of..Binned.avg.pre.slope..prev.epoch.final.rank.127"       
##  [8] "X163..before.pit..sqrt.of..Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES.137"  
##  [9] "X164..before.pit..square.of..Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES.137"
## [10] "X165..before.pit..exp.of..Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES.137"

selected_features<-c(selected_features,subset)
unique_f<-unique(selected_features)
unique_f

##  [1] "X26..before.pit..incoming.rank.before.pit.bef.leg"                                  
##  [2] "X34..before.pit..25th.percentile.rank.upto.bef.pit"                                 
##  [3] "X35..before.pit..75th.percentile.rank.upto.bef.pit"                                 
##  [4] "X82..before.pit..starting.position.of.the.car"                                      
##  [5] "X77..after.pit..indicator.4.tire.change..ref.19"                                    
##  [6] "X139..after.pit...Binned.avg.pre.slope..pit.exit.rank.HI.RES"                       
##  [7] "X61..before.pit..NBHD.prev.epoch..number.of.4.pit.types.of.neighborhood"            
##  [8] "X137..before.pit...Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES"              
##  [9] "X157..before.pit..sqrt.of..Binned.avg.pre.slope..prev.epoch.final.rank.127"         
## [10] "X158..before.pit..square.of..Binned.avg.pre.slope..prev.epoch.final.rank.127"       
## [11] "X163..before.pit..sqrt.of..Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES.137"  
## [12] "X164..before.pit..square.of..Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES.137"
## [13] "X165..before.pit..exp.of..Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES.137"

subset_x=x[,unique_f]

Your browser does not support the HTML5 canvas element.

You must enable Javascript to view this page properly.

lasso with a multinomial loss. Here we train on the full set of features and see that even with maxit=1000000 there is no convergence.

train <- sample(nrow(x[]), round(nrow(x) * 0.7))
cvfit <- cv.glmnet(x[train, ], y[train], family="multinomial",maxit=1000000)

## Warning: from glmnet Fortran code (error code -70); Convergence for 70th
## lambda value not reached after maxit=1000000 iterations; solutions for
## larger lambdas returned

## Warning: from glmnet Fortran code (error code -79); Convergence for 79th
## lambda value not reached after maxit=1000000 iterations; solutions for
## larger lambdas returned

## Warning: from glmnet Fortran code (error code -68); Convergence for 68th
## lambda value not reached after maxit=1000000 iterations; solutions for
## larger lambdas returned

plot(cvfit) # let's look at the cross validated error to choose lambda

Your browser does not support the HTML5 canvas element.

You must enable Javascript to view this page properly.

We train on the data with a subset of features

train <- sample(nrow(subset_x), round(nrow(subset_x) * 0.7))
cvfit <- cv.glmnet(subset_x[train, ], y[train], family="multinomial",maxit=1000000)
plot(cvfit) # let's look at the cross validated error to choose lambda

Your browser does not support the HTML5 canvas element.

You must enable Javascript to view this page properly.

coef(cvfit, s="lambda.1se")

## $`1`
## 14 x 1 sparse Matrix of class "dgCMatrix"
##                                                                                             1
## (Intercept)                                                                         0.4742207
## X26..before.pit..incoming.rank.before.pit.bef.leg                                   .        
## X34..before.pit..25th.percentile.rank.upto.bef.pit                                  .        
## X35..before.pit..75th.percentile.rank.upto.bef.pit                                  .        
## X82..before.pit..starting.position.of.the.car                                       .        
## X77..after.pit..indicator.4.tire.change..ref.19                                     .        
## X139..after.pit...Binned.avg.pre.slope..pit.exit.rank.HI.RES                        .        
## X61..before.pit..NBHD.prev.epoch..number.of.4.pit.types.of.neighborhood             .        
## X137..before.pit...Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES               .        
## X157..before.pit..sqrt.of..Binned.avg.pre.slope..prev.epoch.final.rank.127          .        
## X158..before.pit..square.of..Binned.avg.pre.slope..prev.epoch.final.rank.127        .        
## X163..before.pit..sqrt.of..Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES.137   .        
## X164..before.pit..square.of..Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES.137 .        
## X165..before.pit..exp.of..Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES.137    .        
## 
## $`2`
## 14 x 1 sparse Matrix of class "dgCMatrix"
##                                                                                              1
## (Intercept)                                                                         -0.5066086
## X26..before.pit..incoming.rank.before.pit.bef.leg                                    .        
## X34..before.pit..25th.percentile.rank.upto.bef.pit                                   .        
## X35..before.pit..75th.percentile.rank.upto.bef.pit                                   .        
## X82..before.pit..starting.position.of.the.car                                        .        
## X77..after.pit..indicator.4.tire.change..ref.19                                      .        
## X139..after.pit...Binned.avg.pre.slope..pit.exit.rank.HI.RES                         .        
## X61..before.pit..NBHD.prev.epoch..number.of.4.pit.types.of.neighborhood              .        
## X137..before.pit...Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES                .        
## X157..before.pit..sqrt.of..Binned.avg.pre.slope..prev.epoch.final.rank.127           .        
## X158..before.pit..square.of..Binned.avg.pre.slope..prev.epoch.final.rank.127         .        
## X163..before.pit..sqrt.of..Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES.137    .        
## X164..before.pit..square.of..Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES.137  .        
## X165..before.pit..exp.of..Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES.137     .        
## 
## $`3`
## 14 x 1 sparse Matrix of class "dgCMatrix"
##                                                                                              1
## (Intercept)                                                                         0.03238792
## X26..before.pit..incoming.rank.before.pit.bef.leg                                   .         
## X34..before.pit..25th.percentile.rank.upto.bef.pit                                  .         
## X35..before.pit..75th.percentile.rank.upto.bef.pit                                  .         
## X82..before.pit..starting.position.of.the.car                                       .         
## X77..after.pit..indicator.4.tire.change..ref.19                                     .         
## X139..after.pit...Binned.avg.pre.slope..pit.exit.rank.HI.RES                        .         
## X61..before.pit..NBHD.prev.epoch..number.of.4.pit.types.of.neighborhood             .         
## X137..before.pit...Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES               .         
## X157..before.pit..sqrt.of..Binned.avg.pre.slope..prev.epoch.final.rank.127          .         
## X158..before.pit..square.of..Binned.avg.pre.slope..prev.epoch.final.rank.127        .         
## X163..before.pit..sqrt.of..Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES.137   .         
## X164..before.pit..square.of..Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES.137 .         
## X165..before.pit..exp.of..Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES.137    .

yhat <- predict(cvfit, newx=subset_x[-train, ], s="lambda.1se", type="class")
mean(y[-train] != yhat) # classification error

## [1] 0.5208333

table(y[-train], yhat) # confusion matrix

##    yhat
##      1
##   1 23
##   2  9
##   3 16

We look at the pairwise plots of the features selected by the various algorithms.

# pairs(x, pch="o", col=y, lower.panel=NULL)
pairs(subset_x, pch="o", col=y, lower.panel=NULL)

Using only the selected subset of unique features

pca <- prcomp(subset_x, scale=TRUE)
plot(pca$x, col=y, pch=20)

The first two principal component directions represent the following linear combinations of nutrients:

pca$rotation[,1:2]

##                                                                                            PC1
## X26..before.pit..incoming.rank.before.pit.bef.leg                                   0.31454322
## X34..before.pit..25th.percentile.rank.upto.bef.pit                                  0.31871934
## X35..before.pit..75th.percentile.rank.upto.bef.pit                                  0.31791246
## X82..before.pit..starting.position.of.the.car                                       0.30629681
## X77..after.pit..indicator.4.tire.change..ref.19                                     0.11170985
## X139..after.pit...Binned.avg.pre.slope..pit.exit.rank.HI.RES                        0.28167293
## X61..before.pit..NBHD.prev.epoch..number.of.4.pit.types.of.neighborhood             0.10354919
## X137..before.pit...Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES               0.32586853
## X157..before.pit..sqrt.of..Binned.avg.pre.slope..prev.epoch.final.rank.127          0.31373333
## X158..before.pit..square.of..Binned.avg.pre.slope..prev.epoch.final.rank.127        0.29601771
## X163..before.pit..sqrt.of..Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES.137   0.32028511
## X164..before.pit..square.of..Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES.137 0.31600091
## X165..before.pit..exp.of..Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES.137    0.08584174
##                                                                                             PC2
## X26..before.pit..incoming.rank.before.pit.bef.leg                                   -0.06449678
## X34..before.pit..25th.percentile.rank.upto.bef.pit                                  -0.06482819
## X35..before.pit..75th.percentile.rank.upto.bef.pit                                  -0.06754430
## X82..before.pit..starting.position.of.the.car                                       -0.06171546
## X77..after.pit..indicator.4.tire.change..ref.19                                      0.47702897
## X139..after.pit...Binned.avg.pre.slope..pit.exit.rank.HI.RES                        -0.08189330
## X61..before.pit..NBHD.prev.epoch..number.of.4.pit.types.of.neighborhood             -0.70435628
## X137..before.pit...Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES                0.01570011
## X157..before.pit..sqrt.of..Binned.avg.pre.slope..prev.epoch.final.rank.127           0.03736663
## X158..before.pit..square.of..Binned.avg.pre.slope..prev.epoch.final.rank.127         0.17699471
## X163..before.pit..sqrt.of..Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES.137   -0.01798382
## X164..before.pit..square.of..Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES.137  0.06562220
## X165..before.pit..exp.of..Binned.avg.pre.slope..prev.epoch.final.rank.HI.RES.137     0.46398917