Team 5 - Data 624 - Project 2
Ohannes (Hovig) Ohannessian, Niteen Kumar, Gurpreet Singh, Peter Goodridge
4/20/2019
student_evaluation_gh_file <- "https://github.com/hovig/Team5-Data624-Project2/raw/master/StudentEvaluation.xlsx"
student_data_gh_file <- "https://github.com/hovig/Team5-Data624-Project2/raw/master/StudentData.xlsx"
se_temp_file <- tempfile(fileext = ".xlsx")
sd_temp_file <- tempfile(fileext = ".xlsx")
se_data <- GET(student_evaluation_gh_file, authenticate(Sys.getenv("GITHUB_PAT"), ""), write_disk(path = se_temp_file))
df_data <- GET(student_data_gh_file, authenticate(Sys.getenv("GITHUB_PAT"), ""), write_disk(path = sd_temp_file))
se_data <- readxl::read_excel(se_temp_file)
sd_data <- readxl::read_excel(sd_temp_file)hist.data.frame(se_data)
hist.data.frame(sd_data)
##Explore and Process Data
bev_raw <- read_csv('https://raw.githubusercontent.com/hovig/Team5-Data624-Project2/master/StudentData.csv') #%>%
#drop_na(PH)bev_raw %>%
ggplot(aes(PH, fill=PH > 9)) +
geom_histogram(bins=30) +
theme_bw() +
theme(legend.position='center') +
labs(y='Count', title='PH Levels in Dataset')
bev_raw <- bev_raw %>% filter(!is.na(bev_raw$PH), bev_raw$PH < 9) dim(bev_raw)## [1] 2566 33str(bev_raw)## Classes 'spec_tbl_df', 'tbl_df', 'tbl' and 'data.frame': 2566 obs. of 33 variables:
## $ Brand Code : chr "B" "A" "B" "A" ...
## $ Carb Volume : num 5.34 5.43 5.29 5.44 5.49 ...
## $ Fill Ounces : num 24 24 24.1 24 24.3 ...
## $ PC Volume : num 0.263 0.239 0.263 0.293 0.111 ...
## $ Carb Pressure : num 68.2 68.4 70.8 63 67.2 66.6 64.2 67.6 64.2 72 ...
## $ Carb Temp : num 141 140 145 133 137 ...
## $ PSC : num 0.104 0.124 0.09 NA 0.026 0.09 0.128 0.154 0.132 0.014 ...
## $ PSC Fill : num 0.26 0.22 0.34 0.42 0.16 0.24 0.4 0.34 0.12 0.24 ...
## $ PSC CO2 : num 0.04 0.04 0.16 0.04 0.12 0.04 0.04 0.04 0.14 0.06 ...
## $ Mnf Flow : num -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 ...
## $ Carb Pressure1 : num 119 122 120 115 118 ...
## $ Fill Pressure : num 46 46 46 46.4 45.8 45.6 51.8 46.8 46 45.2 ...
## $ Hyd Pressure1 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ Hyd Pressure2 : num NA NA NA 0 0 0 0 0 0 0 ...
## $ Hyd Pressure3 : num NA NA NA 0 0 0 0 0 0 0 ...
## $ Hyd Pressure4 : num 118 106 82 92 92 116 124 132 90 108 ...
## $ Filler Level : num 121 119 120 118 119 ...
## $ Filler Speed : num 4002 3986 4020 4012 4010 ...
## $ Temperature : num 66 67.6 67 65.6 65.6 66.2 65.8 65.2 65.4 66.6 ...
## $ Usage cont : num 16.2 19.9 17.8 17.4 17.7 ...
## $ Carb Flow : num 2932 3144 2914 3062 3054 ...
## $ Density : num 0.88 0.92 1.58 1.54 1.54 1.52 0.84 0.84 0.9 0.9 ...
## $ MFR : num 725 727 735 731 723 ...
## $ Balling : num 1.4 1.5 3.14 3.04 3.04 ...
## $ Pressure Vacuum : num -4 -4 -3.8 -4.4 -4.4 -4.4 -4.4 -4.4 -4.4 -4.4 ...
## $ PH : num 8.36 8.26 8.94 8.24 8.26 8.32 8.4 8.38 8.38 8.5 ...
## $ Oxygen Filler : num 0.022 0.026 0.024 0.03 0.03 0.024 0.066 0.046 0.064 0.022 ...
## $ Bowl Setpoint : num 120 120 120 120 120 120 120 120 120 120 ...
## $ Pressure Setpoint: num 46.4 46.8 46.6 46 46 46 46 46 46 46 ...
## $ Air Pressurer : num 143 143 142 146 146 ...
## $ Alch Rel : num 6.58 6.56 7.66 7.14 7.14 7.16 6.54 6.52 6.52 6.54 ...
## $ Carb Rel : num 5.32 5.3 5.84 5.42 5.44 5.44 5.38 5.34 5.34 5.34 ...
## $ Balling Lvl : num 1.48 1.56 3.28 3.04 3.04 3.02 1.44 1.44 1.44 1.38 ...
## - attr(*, "spec")=
## .. cols(
## .. `Brand Code` = col_character(),
## .. `Carb Volume` = col_double(),
## .. `Fill Ounces` = col_double(),
## .. `PC Volume` = col_double(),
## .. `Carb Pressure` = col_double(),
## .. `Carb Temp` = col_double(),
## .. PSC = col_double(),
## .. `PSC Fill` = col_double(),
## .. `PSC CO2` = col_double(),
## .. `Mnf Flow` = col_double(),
## .. `Carb Pressure1` = col_double(),
## .. `Fill Pressure` = col_double(),
## .. `Hyd Pressure1` = col_double(),
## .. `Hyd Pressure2` = col_double(),
## .. `Hyd Pressure3` = col_double(),
## .. `Hyd Pressure4` = col_double(),
## .. `Filler Level` = col_double(),
## .. `Filler Speed` = col_double(),
## .. Temperature = col_double(),
## .. `Usage cont` = col_double(),
## .. `Carb Flow` = col_double(),
## .. Density = col_double(),
## .. MFR = col_double(),
## .. Balling = col_double(),
## .. `Pressure Vacuum` = col_double(),
## .. PH = col_double(),
## .. `Oxygen Filler` = col_double(),
## .. `Bowl Setpoint` = col_double(),
## .. `Pressure Setpoint` = col_double(),
## .. `Air Pressurer` = col_double(),
## .. `Alch Rel` = col_double(),
## .. `Carb Rel` = col_double(),
## .. `Balling Lvl` = col_double()
## .. )hist.data.frame(bev_raw)
table(bev_raw$`Brand Code`)##
## A B C D
## 293 1235 303 615summary(bev_raw)## Brand Code Carb Volume Fill Ounces PC Volume
## Length:2566 Min. :5.040 Min. :23.63 Min. :0.07933
## Class :character 1st Qu.:5.293 1st Qu.:23.92 1st Qu.:0.23933
## Mode :character Median :5.347 Median :23.97 Median :0.27133
## Mean :5.370 Mean :23.97 Mean :0.27719
## 3rd Qu.:5.453 3rd Qu.:24.03 3rd Qu.:0.31200
## Max. :5.700 Max. :24.32 Max. :0.47800
## NA's :10 NA's :38 NA's :39
## Carb Pressure Carb Temp PSC PSC Fill
## Min. :57.00 Min. :128.6 Min. :0.00200 Min. :0.0000
## 1st Qu.:65.60 1st Qu.:138.4 1st Qu.:0.04800 1st Qu.:0.1000
## Median :68.20 Median :140.8 Median :0.07600 Median :0.1800
## Mean :68.19 Mean :141.1 Mean :0.08468 Mean :0.1953
## 3rd Qu.:70.60 3rd Qu.:143.8 3rd Qu.:0.11200 3rd Qu.:0.2600
## Max. :79.40 Max. :154.0 Max. :0.27000 Max. :0.6200
## NA's :27 NA's :26 NA's :33 NA's :23
## PSC CO2 Mnf Flow Carb Pressure1 Fill Pressure
## Min. :0.00000 Min. :-100.20 Min. :105.6 Min. :34.60
## 1st Qu.:0.02000 1st Qu.:-100.00 1st Qu.:119.0 1st Qu.:46.00
## Median :0.04000 Median : 70.70 Median :123.2 Median :46.40
## Mean :0.05645 Mean : 24.68 Mean :122.6 Mean :47.92
## 3rd Qu.:0.08000 3rd Qu.: 140.80 3rd Qu.:125.4 3rd Qu.:50.00
## Max. :0.24000 Max. : 229.40 Max. :140.2 Max. :60.40
## NA's :39 NA's :32 NA's :17
## Hyd Pressure1 Hyd Pressure2 Hyd Pressure3 Hyd Pressure4
## Min. :-0.80 Min. : 0.0 Min. :-1.20 Min. : 62.00
## 1st Qu.: 0.00 1st Qu.: 0.0 1st Qu.: 0.00 1st Qu.: 86.00
## Median :11.40 Median :28.6 Median :27.60 Median : 96.00
## Mean :12.46 Mean :21.0 Mean :20.49 Mean : 96.31
## 3rd Qu.:20.20 3rd Qu.:34.6 3rd Qu.:33.40 3rd Qu.:102.00
## Max. :58.00 Max. :59.4 Max. :50.00 Max. :142.00
## NA's :11 NA's :15 NA's :15 NA's :27
## Filler Level Filler Speed Temperature Usage cont
## Min. : 55.8 Min. : 998 Min. :63.60 Min. :12.08
## 1st Qu.: 98.3 1st Qu.:3888 1st Qu.:65.20 1st Qu.:18.38
## Median :118.4 Median :3982 Median :65.60 Median :21.80
## Mean :109.3 Mean :3688 Mean :65.96 Mean :21.00
## 3rd Qu.:120.0 3rd Qu.:3998 3rd Qu.:66.40 3rd Qu.:23.76
## Max. :161.2 Max. :4030 Max. :76.20 Max. :25.90
## NA's :15 NA's :53 NA's :11 NA's :5
## Carb Flow Density MFR Balling
## Min. : 26 Min. :0.300 Min. : 31.4 Min. :0.346
## 1st Qu.:1170 1st Qu.:0.900 1st Qu.:706.3 1st Qu.:1.496
## Median :3030 Median :0.980 Median :724.0 Median :1.648
## Mean :2473 Mean :1.175 Mean :704.0 Mean :2.201
## 3rd Qu.:3188 3rd Qu.:1.620 3rd Qu.:731.0 3rd Qu.:3.292
## Max. :5104 Max. :1.920 Max. :868.6 Max. :4.012
## NA's :2 NA's :207
## Pressure Vacuum PH Oxygen Filler Bowl Setpoint
## Min. :-6.600 Min. :7.880 Min. :0.00240 Min. : 70.0
## 1st Qu.:-5.600 1st Qu.:8.440 1st Qu.:0.02200 1st Qu.:100.0
## Median :-5.400 Median :8.540 Median :0.03340 Median :120.0
## Mean :-5.216 Mean :8.545 Mean :0.04635 Mean :109.3
## 3rd Qu.:-5.000 3rd Qu.:8.680 3rd Qu.:0.05960 3rd Qu.:120.0
## Max. :-3.600 Max. :8.940 Max. :0.40000 Max. :140.0
## NA's :11 NA's :2
## Pressure Setpoint Air Pressurer Alch Rel Carb Rel
## Min. :44.00 Min. :140.8 Min. :5.280 Min. :4.960
## 1st Qu.:46.00 1st Qu.:142.2 1st Qu.:6.540 1st Qu.:5.340
## Median :46.00 Median :142.6 Median :6.560 Median :5.400
## Mean :47.61 Mean :142.8 Mean :6.898 Mean :5.437
## 3rd Qu.:50.00 3rd Qu.:143.0 3rd Qu.:7.240 3rd Qu.:5.550
## Max. :52.00 Max. :148.2 Max. :8.620 Max. :6.060
## NA's :12 NA's :6 NA's :7
## Balling Lvl
## Min. :0.000
## 1st Qu.:1.380
## Median :1.480
## Mean :2.052
## 3rd Qu.:3.140
## Max. :3.660
## NA's :1describe(bev_raw %>% select(-`Brand Code`))## vars n mean sd median trimmed mad min
## Carb Volume 1 2556 5.37 0.11 5.35 5.37 0.11 5.04
## Fill Ounces 2 2528 23.97 0.09 23.97 23.98 0.08 23.63
## PC Volume 3 2527 0.28 0.06 0.27 0.27 0.05 0.08
## Carb Pressure 4 2539 68.19 3.54 68.20 68.12 3.56 57.00
## Carb Temp 5 2540 141.09 4.03 140.80 140.99 3.85 128.60
## PSC 6 2533 0.08 0.05 0.08 0.08 0.05 0.00
## PSC Fill 7 2543 0.20 0.12 0.18 0.18 0.12 0.00
## PSC CO2 8 2527 0.06 0.04 0.04 0.05 0.03 0.00
## Mnf Flow 9 2566 24.68 119.50 70.70 21.20 160.57 -100.20
## Carb Pressure1 10 2534 122.56 4.72 123.20 122.52 4.45 105.60
## Fill Pressure 11 2549 47.92 3.18 46.40 47.71 2.37 34.60
## Hyd Pressure1 12 2555 12.46 12.43 11.40 10.87 16.90 -0.80
## Hyd Pressure2 13 2551 21.00 16.38 28.60 21.10 13.34 0.00
## Hyd Pressure3 14 2551 20.49 15.97 27.60 20.54 13.64 -1.20
## Hyd Pressure4 15 2539 96.31 13.10 96.00 95.46 11.86 62.00
## Filler Level 16 2551 109.25 15.70 118.40 111.04 9.19 55.80
## Filler Speed 17 2513 3688.11 769.63 3982.00 3920.25 47.44 998.00
## Temperature 18 2555 65.96 1.38 65.60 65.80 0.89 63.60
## Usage cont 19 2561 21.00 2.98 21.80 21.25 3.17 12.08
## Carb Flow 20 2564 2473.00 1069.56 3030.00 2604.91 323.21 26.00
## Density 21 2566 1.17 0.38 0.98 1.15 0.15 0.30
## MFR 22 2359 704.05 73.90 724.00 718.16 15.42 31.40
## Balling 23 2566 2.20 0.93 1.65 2.13 0.37 0.35
## Pressure Vacuum 24 2566 -5.22 0.57 -5.40 -5.25 0.59 -6.60
## PH 25 2566 8.55 0.17 8.54 8.55 0.18 7.88
## Oxygen Filler 26 2555 0.05 0.04 0.03 0.04 0.02 0.00
## Bowl Setpoint 27 2564 109.34 15.29 120.00 111.36 0.00 70.00
## Pressure Setpoint 28 2554 47.61 2.04 46.00 47.60 0.00 44.00
## Air Pressurer 29 2566 142.83 1.21 142.60 142.58 0.59 140.80
## Alch Rel 30 2560 6.90 0.51 6.56 6.84 0.06 5.28
## Carb Rel 31 2559 5.44 0.13 5.40 5.43 0.12 4.96
## Balling Lvl 32 2565 2.05 0.87 1.48 1.98 0.21 0.00
## max range skew kurtosis se
## Carb Volume 5.70 0.66 0.39 -0.47 0.00
## Fill Ounces 24.32 0.69 -0.02 0.87 0.00
## PC Volume 0.48 0.40 0.35 0.67 0.00
## Carb Pressure 79.40 22.40 0.18 -0.01 0.07
## Carb Temp 154.00 25.40 0.24 0.23 0.08
## PSC 0.27 0.27 0.85 0.65 0.00
## PSC Fill 0.62 0.62 0.93 0.77 0.00
## PSC CO2 0.24 0.24 1.73 3.71 0.00
## Mnf Flow 229.40 329.60 0.00 -1.87 2.36
## Carb Pressure1 140.20 34.60 0.03 0.10 0.09
## Fill Pressure 60.40 25.80 0.55 1.41 0.06
## Hyd Pressure1 58.00 58.80 0.78 -0.15 0.25
## Hyd Pressure2 59.40 59.40 -0.31 -1.56 0.32
## Hyd Pressure3 50.00 51.20 -0.32 -1.57 0.32
## Hyd Pressure4 142.00 80.00 0.56 0.61 0.26
## Filler Level 161.20 105.40 -0.85 0.05 0.31
## Filler Speed 4030.00 3032.00 -2.88 6.75 15.35
## Temperature 76.20 12.60 2.39 10.25 0.03
## Usage cont 25.90 13.82 -0.54 -1.02 0.06
## Carb Flow 5104.00 5078.00 -0.99 -0.57 21.12
## Density 1.92 1.62 0.54 -1.22 0.01
## MFR 868.60 837.20 -5.09 30.46 1.52
## Balling 4.01 3.67 0.60 -1.41 0.02
## Pressure Vacuum -3.60 3.00 0.53 -0.03 0.01
## PH 8.94 1.06 -0.33 -0.08 0.00
## Oxygen Filler 0.40 0.40 2.41 8.89 0.00
## Bowl Setpoint 140.00 70.00 -0.97 -0.06 0.30
## Pressure Setpoint 52.00 8.00 0.20 -1.60 0.04
## Air Pressurer 148.20 7.40 2.25 4.72 0.02
## Alch Rel 8.62 3.34 0.88 -0.85 0.01
## Carb Rel 6.06 1.10 0.50 -0.30 0.00
## Balling Lvl 3.66 3.66 0.60 -1.51 0.02###Zero variance
df1 <- bev_raw %>% select(-`Brand Code`) %>% mutate_each(funs(as.numeric(.)))%>%complete.cases()%>%
as.data.frame()
names(df1)[nearZeroVar(df1)]## character(0)nzv <- nearZeroVar(bev_raw,saveMetrics= TRUE)
nzv[nzv$nzv,]## freqRatio percentUnique zeroVar nzv
## Hyd Pressure1 31 9.547935 FALSE TRUE###Box plot
df.m <- melt(bev_raw %>% select(-MFR, -`Filler Speed`, -`Carb Flow`,-`Bowl Setpoint`,`Carb Pressure1`,
`Hyd Pressure4`, `Air Pressurer`, `Carb Temp`, `Filler Level`, `Mnf Flow`), id.var = "Brand Code")
p <-ggplot(data = df.m, aes(x=variable, y=value)) +
geom_boxplot(outlier.colour="red", outlier.shape=8, outlier.size=4,aes(fill=variable)) +
scale_y_continuous(name = "Predictors for PH", breaks = seq(0, 2, 0.5)) +
coord_flip()
p
###Normality
bev_raw%>%
select(-`Brand Code`) %>%
select(2:20) %>%
gather() %>%
ggplot(aes(value)) +
facet_wrap(~ key, scales = "free") +
geom_density() 
bev_raw%>%
select(-`Brand Code`) %>%
select(`Carb Flow`, Density, MFR, Balling, `Pressure Vacuum`, PH, `Oxygen Filler`, `Bowl Setpoint`, `Pressure Setpoint`, `Air Pressurer`, `Alch Rel`,`Carb Rel`,`Balling Lvl`) %>%
gather() %>%
ggplot(aes(value)) +
facet_wrap(~ key, scales = "free") +
geom_density() 
##Prep Data
We perform 3 data preperatoin steps - Remove the near-zero variance variables we previously mentioned - Impute missing values with a Random Forest Regression with the MICE package - Create dummy variables for the categorical variable, brand
Random Forest was chosen as the regression method in imputation because it requires very little tuning, allowing and there is no interface to tune the imputation model in MICE.
bev <- select(bev_raw, -PH) %>%
union_all(se_data) #since we're not using the labels (y values), combining the sets for data clearning is fine
cols <- str_replace(colnames(bev), '\\s', '')
colnames(bev) <- cols
bev <- mutate(bev, BrandCode = ifelse(is.na(BrandCode), 'Unknown', BrandCode))
y <- bev_raw$PH
zero_vars <- nearZeroVar(bev[1:2566, ])
bev_new <- bev[, -zero_vars]
pred <- mice(data = bev_new, m = 3, method = 'rf', maxit = 3)
bev_imputed <- complete(pred)
form <- dummyVars(~ ., data = bev_imputed)
bev_imputed <- predict(form, bev_imputed) %>% data.frame() %>% as.matrix()###Test Train Split
We will split the data into train/test sets before running cross validaiton. This will allow us to make sure that predictions are consistent at all levels of PH.
bev_eval <- bev_imputed[2567:2833,]
bev_train <- bev_imputed[1:2566,]
samples <- createDataPartition(y, p = .75, list = F)
x_train <- bev_train[samples, ]
x_test<- bev_train[-samples, ]
y_train <- y[samples]
y_test <- y[-samples]##Build Models
For state-of-the-art prediction quality, we will use a model stack. This will consist of tuning models separately and then combining the candidate models into a metamodel that will formulate predictions with a linear combination of our tuned models.
From the perspective of understanding the manufacturing process, the model stack will also provide benefits. The stack is like a panel of experts, each looking at the data through slightly different lenses to form their diagnoses. By looking at the predictors each model uses, we can gather assemble a complete picture of the factors that affect our manufacturing process.
#xgb_grid <- expand.grid(eta = c(.01,.025), nrounds = c(750,1000), max_depth = c(5,6,7), gamma = c(.1,0), colsample_bytree = c(.8), min_child_weight = c(.4,.8), subsample = c(.8))
#cub_grid <- expand.grid(.committees = c(1,3,5), .neighbors = c(1,3,5,7,9))
#mars_grid <- expand.grid(.degree = c(1,2), .nprune = seq(16,36,4))
#dart_grid <- expand.grid(eta = c(.025), nrounds = c(1000), gamma = c(.001,.1), skip_drop = c(.4,.8), rate_drop = c(.4,.8), #max_depth = c(6,7,8), colsample_bytree = c(.8), min_child_weight = c(.8), subsample = c(.8))
#rf_grid <- expand.grid(mtry = c(5,10,15,20,25,30))
xgb_grid <- expand.grid(eta = c(.01), nrounds = c(1000), max_depth = c(6), gamma = c(0), colsample_bytree = c(.8), min_child_weight = c(.8), subsample = c(.8))
cub_grid <- expand.grid(.committees = c(5), .neighbors = c(7))
mars_grid <- expand.grid(.degree = c(2), .nprune = 24)
dart_grid <- expand.grid(eta = c(.01), nrounds = c(1000), gamma = c(.1), skip_drop = c(.6), rate_drop = c(.4), max_depth = c(6), colsample_bytree = c(.6), min_child_weight = c(.6), subsample = c(.6))
rf_grid <- expand.grid(mtry = 25)
tuning_list <-list(
caretModelSpec(method="xgbTree", tuneGrid = xgb_grid),
caretModelSpec(method="xgbDART", tuneGrid = dart_grid),
caretModelSpec(method="cubist", tuneGrid = cub_grid),
caretModelSpec(method="rf", tuneGrid = rf_grid, importance = TRUE),
caretModelSpec(method="bagEarth", tuneGrid = mars_grid)
)
cl <- makeCluster(detectCores())
registerDoParallel(cl)
set.seed(42)
my_control <- trainControl(method = 'cv',
number = 5,
savePredictions = 'final',
index = createFolds(y_train, 5),
allowParallel = TRUE)
mod_list <- caretList(
x = x_train,
y = y_train,
preProcess = c('center', 'scale'),
trControl = my_control,
tuneList = tuning_list
)###MARS
#helper function
show_results <- function(model){
rslts <- model$results %>%
arrange(RMSE)
head(rslts, 10) %>% kable () %>% kable_styling(bootstrap_options = "striped", full_width = F)
}show_results(mod_list$bagEarth)| degree | nprune | RMSE | Rsquared | MAE | RMSESD | RsquaredSD | MAESD |
|---|---|---|---|---|---|---|---|
| 2 | 24 | 0.1330667 | 0.4156957 | 0.0970321 | 0.006257 | 0.0419238 | 0.0020659 |
#plot(mod_list$mars, metric = "RMSE")predict(mod_list$bagEarth, x_test) %>%
postResample(y_test)## RMSE Rsquared MAE
## 0.1176685 0.5226173 0.0894521###RandomForest
show_results(mod_list$rf)| mtry | RMSE | Rsquared | MAE | RMSESD | RsquaredSD | MAESD |
|---|---|---|---|---|---|---|
| 25 | 0.1219214 | 0.511375 | 0.0929628 | 0.0034692 | 0.0260702 | 0.002122 |
predict(mod_list$rf, x_test) %>%
postResample(y_test)## RMSE Rsquared MAE
## 0.09458247 0.69663264 0.06935676###Cubist
show_results(mod_list$cubist)| committees | neighbors | RMSE | Rsquared | MAE | RMSESD | RsquaredSD | MAESD |
|---|---|---|---|---|---|---|---|
| 5 | 7 | 0.1277949 | 0.4602262 | 0.0963549 | 0.0023973 | 0.0183168 | 0.0019012 |
#plot(mod_list$cub, metric = "RMSE")predict(mod_list$cubist, x_test) %>%
postResample(y_test)## RMSE Rsquared MAE
## 0.10141406 0.64766552 0.07322895###XGB Trees
show_results(mod_list$xgbTree)| eta | nrounds | max_depth | gamma | colsample_bytree | min_child_weight | subsample | RMSE | Rsquared | MAE | RMSESD | RsquaredSD | MAESD |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.01 | 1000 | 6 | 0 | 0.8 | 0.8 | 0.8 | 0.1241533 | 0.4842761 | 0.0936952 | 0.0026617 | 0.0214108 | 0.0013011 |
#plot(mod_list$xgbt, metric = "RMSE")predict(mod_list$xgbTree, x_test) %>%
postResample(y_test)## RMSE Rsquared MAE
## 0.10021691 0.65562489 0.07580665###XGB Dart
show_results(mod_list$xgbDART)| eta | nrounds | gamma | skip_drop | rate_drop | max_depth | colsample_bytree | min_child_weight | subsample | RMSE | Rsquared | MAE | RMSESD | RsquaredSD | MAESD |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.01 | 1000 | 0.1 | 0.6 | 0.4 | 6 | 0.6 | 0.6 | 0.6 | 0.1319881 | 0.4478713 | 0.1048579 | 0.0023489 | 0.0175528 | 0.001896 |
#plot(mod_list$xgbd, metric = "RMSE")predict(mod_list$xgbDART, x_test) %>%
postResample(y_test)## RMSE Rsquared MAE
## 0.11195350 0.59448452 0.08782905Assesment
While the XGBTee model offers the lowest RMSE, it still only explains about 50% of the variance in PH. There is plenty of room for the other models to offer value. We will use the hyperparameters yielding the best RMSE in our final ensemble that will be built later.
###Variable Importances
p1 <- varImp(mod_list$cubist) %>%
plot(top = 10, main = 'Cubist')
p2 <- varImp(mod_list$xgbTree) %>%
plot(top = 10, main = 'XGBTrees')
p3 <- varImp(mod_list$xgbDART) %>%
plot(top = 10, main = 'XGBDart')
p4 <- varImp(mod_list$bagEarth) %>%
plot(top = 10, main = 'MARS')
p5 <- varImp(mod_list$rf) %>%
plot(top = 10, main = 'RF')
grid.arrange(p1,p2,p3,p4,p5, nrow = 3, ncol = 2)
Observations:
- The top predictor, MnFlow, is consistent throughout the set of models
- The XGBoost models are similar in variables, but less similar in their distributions of importance
- After the top 5, divergence increases
Based on what we see, there are more factors contributing to PH with the ensemble than with a single model. We’ll confirm this by looking at the count by predictor among each top 10:
Variable Counts among the Models
top_10 <- function(model) {
varImp(model)$importance %>%
rownames_to_column('var') %>%
arrange(desc(Overall)) %>%
head(10) %>%
select(var)
}
rbind(top_10(mod_list$xgbDART), top_10(mod_list$xgbTree), top_10(mod_list$cubist), top_10(mod_list$bagEarth), top_10(mod_list$rf)) %>%
group_by(var) %>%
summarise(ModelCount = n()) %>%
arrange(desc(ModelCount)) %>%
kable () %>% kable_styling(bootstrap_options = "striped", full_width = F)| var | ModelCount |
|---|---|
| AlchRel | 5 |
| MnfFlow | 5 |
| PressureVacuum | 5 |
| BrandCodeC | 4 |
| CarbPressure1 | 4 |
| OxygenFiller | 4 |
| Temperature | 4 |
| FillerSpeed | 3 |
| Usagecont | 3 |
| AirPressurer | 2 |
| Balling | 2 |
| BallingLvl | 2 |
| CarbRel | 2 |
| HydPressure3 | 2 |
| BowlSetpoint | 1 |
| CarbFlow | 1 |
| Density | 1 |
Best Single Tree
Because 4 of the 5 models in the ensemble are tree-based, looking a single decision tree modle will give us an idea of not just what predictors are important, but how they are being used to form predictions.
rp_mod <- rpart::rpart(ph ~ ., data = data.frame(bev_train, ph = y), method = 'anova')
rpart.plot(rp_mod, uniform=TRUE)
####Test Model Stack
Now that our models are optimized and understood, we will form the ensemble and test the results. We first look at their prediction correlations to confirm that the stack, intuitivley, be effective. The lower the correlation between the models, the more effective the stack will be.
resamples <- resamples(mod_list)
modelCor(resamples)## xgbTree xgbDART cubist rf bagEarth
## xgbTree 1.0000000 0.9952425 0.9484601 0.8493977 0.9537161
## xgbDART 0.9952425 1.0000000 0.9293228 0.8291363 0.9531146
## cubist 0.9484601 0.9293228 1.0000000 0.7179804 0.8370685
## rf 0.8493977 0.8291363 0.7179804 1.0000000 0.9517525
## bagEarth 0.9537161 0.9531146 0.8370685 0.9517525 1.0000000Our model list has negative correlations, indicating that it will be effective.
ensemble1 <- caretStack(
mod_list,
method="glmnet",
metric="RMSE",
trControl=trainControl(
method="cv",
number=5,
savePredictions="final"
)
)
ensemble1## A glmnet ensemble of 2 base models: xgbTree, xgbDART, cubist, rf, bagEarth
##
## Ensemble results:
## glmnet
##
## 7704 samples
## 5 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 6164, 6164, 6162, 6163, 6163
## Resampling results across tuning parameters:
##
## alpha lambda RMSE Rsquared MAE
## 0.10 0.0002459543 0.1188060 0.5255875 0.08876547
## 0.10 0.0024595432 0.1189930 0.5241183 0.08912062
## 0.10 0.0245954322 0.1202218 0.5156488 0.09072844
## 0.55 0.0002459543 0.1188111 0.5255451 0.08879073
## 0.55 0.0024595432 0.1194029 0.5208463 0.08978722
## 0.55 0.0245954322 0.1208517 0.5178956 0.09168230
## 1.00 0.0002459543 0.1188191 0.5254764 0.08882444
## 1.00 0.0024595432 0.1197060 0.5184685 0.09027092
## 1.00 0.0245954322 0.1223586 0.5177983 0.09355144
##
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were alpha = 0.1 and lambda
## = 0.0002459543.We see a several percent improvement over the best single model.
preds <- predict(ensemble1, x_test)
preds %>% postResample(y_test)## RMSE Rsquared MAE
## 0.09070157 0.71474571 0.06441936Predicted Vs Actual
preds_df <- data.frame(predicted = preds, actual = y_test)
ggplot(preds_df) + geom_point(aes(x = predicted, y = actual)) + geom_smooth(aes(x = predicted, y = actual))
On the test set, the ensemble’s predictoins are consistent throughout the range of PH. Consitent quality is, theoretically, a feature of ensemble models because the models are weighted such that predictions are maximized. If one model predicts better at lower PH, it will be upweighted over that range. There are several outliers, but when not all of the variance is explainable, this is unavoidable.
##Predict on Eval Set
With our hyperparameters of the individual models optimized, we want to provide each model the entire training set and build the meta-model based on these optimized single models.
#xgb_grid <- expand.grid(eta = c(.01,.025), nrounds = c(750,1000), max_depth = c(5,6,7), gamma = c(.1,0), colsample_bytree = c(.8), min_child_weight = c(.4,.8), subsample = c(.8))
#cub_grid <- expand.grid(.committees = c(1,3,5), .neighbors = c(1,3,5,7,9))
#mars_grid <- expand.grid(.degree = c(1,2), .nprune = seq(16,36,4))
#dart_grid <- expand.grid(eta = c(.025), nrounds = c(1000), gamma = c(.001,.1), skip_drop = c(.4,.8), rate_drop = c(.4,.8), #max_depth = c(6,7,8), colsample_bytree = c(.8), min_child_weight = c(.8), subsample = c(.8))
#rf_grid <- expand.grid(mtry = c(5,10,15,20,25,30))
xgb_grid <- expand.grid(eta = c(.01), nrounds = c(1000), max_depth = c(6), gamma = c(0), colsample_bytree = c(.8), min_child_weight = c(.8), subsample = c(.8))
cub_grid <- expand.grid(.committees = c(5), .neighbors = c(7))
mars_grid <- expand.grid(.degree = c(2), .nprune = 24)
dart_grid <- expand.grid(eta = c(.01), nrounds = c(1000), gamma = c(.1), skip_drop = c(.6), rate_drop = c(.4), max_depth = c(6), colsample_bytree = c(.6), min_child_weight = c(.6), subsample = c(.6))
rf_grid <- expand.grid(mtry = 25)
tuning_list <-list(
caretModelSpec(method="xgbTree", tuneGrid = xgb_grid),
caretModelSpec(method="xgbDART", tuneGrid = dart_grid),
caretModelSpec(method="cubist", tuneGrid = cub_grid),
caretModelSpec(method="rf", tuneGrid = rf_grid, importance = TRUE),
caretModelSpec(method="bagEarth", tuneGrid = mars_grid)
)
my_control <- trainControl(method = 'cv',
number = 5,
savePredictions = 'final',
index = createFolds(y, 5),
allowParallel = TRUE)
mod_list <- caretList(
x = bev_train,
y = y,
preProcess = c('center', 'scale'),
trControl = my_control,
tuneList = tuning_list
)
final_ensemble <- caretStack(
mod_list,
method="glmnet",
metric="RMSE",
trControl=trainControl(
method="cv",
number=5,
savePredictions="final"
)
)
final_ensemble## A glmnet ensemble of 2 base models: xgbTree, xgbDART, cubist, rf, bagEarth
##
## Ensemble results:
## glmnet
##
## 10264 samples
## 5 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 8212, 8210, 8210, 8211, 8213
## Resampling results across tuning parameters:
##
## alpha lambda RMSE Rsquared MAE
## 0.10 0.0002558189 0.1109015 0.5831504 0.08262287
## 0.10 0.0025581892 0.1111623 0.5812696 0.08299034
## 0.10 0.0255818923 0.1129869 0.5692609 0.08520505
## 0.55 0.0002558189 0.1109057 0.5831216 0.08263288
## 0.55 0.0025581892 0.1116436 0.5777214 0.08358466
## 0.55 0.0255818923 0.1137781 0.5711415 0.08621973
## 1.00 0.0002558189 0.1109142 0.5830606 0.08264615
## 1.00 0.0025581892 0.1124799 0.5714072 0.08454082
## 1.00 0.0255818923 0.1155139 0.5710741 0.08809889
##
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were alpha = 0.1 and lambda
## = 0.0002558189.With more data available, the ensemble improve around 6%. Our model is considerably more predictive than a “baseline” single model.
Ensemble Model Importances
plot(varImp(final_ensemble$ens_model))
The MARS model is excluded, but all four tree models contribute the the final prediction
Export Predictions
predictions <- predict(final_ensemble, bev_eval)
predictions %>% tibble::enframe(name = NULL) %>% datatable()predictions %>% tibble::enframe(name = NULL) %>% rownames_to_column() %>% rename(PH = value, row = rowname)## # A tibble: 267 x 2
## row PH
## <chr> <dbl>
## 1 1 8.50
## 2 2 8.37
## 3 3 8.49
## 4 4 8.56
## 5 5 8.40
## 6 6 8.50
## 7 7 8.46
## 8 8 8.55
## 9 9 8.51
## 10 10 8.64
## # ... with 257 more rowswrite.csv(predictions, file = "predictions.csv")y <- randomForest(x = x_train, y = y_train, importance=TRUE, ntree=1000)
predictions <- predict(y, newdata=x_test)
predictions %>% tibble::enframe(name = NULL) %>% datatable()
predictions %>% tibble::enframe(name = NULL) %>% rownames_to_column() %>% rename(PH = value, row = rowname)
write.csv(predictions, file = "predictions.csv")