ST558_Project-2

ST558 Project 2

Mu-Tien, Lee

• Require package
• Read in data
• Summarize the training data
• Training Model
  • Tree-based model
  • Boosted tree model
• Predicting using the best tree-base model
• Predicting using the best boosted-tree model

Require package

#install.packages("")
library(knitr)
library(rmarkdown)
library(MuMIn)
library(tidyverse)
library(caret)
library(corrplot)
library(readxl)
library(caret)
library(ggiraphExtra)
library(knitr)
library(ggplot2)
library(ggpubr)
library(rpart.plot)
library(rpart)
library(DT)

Read in data

#read in hour data
HourData <- read.csv("hour.csv")
HourData<- HourData %>% select(-casual, -registered)
HourData$yr <- as.factor(HourData$yr)
HourData$holiday <- as.factor(HourData$holiday)
HourData$workingday <- as.factor(HourData$workingday)

#filter data by weekday
HourData <-HourData %>% filter(weekday==params$w)
#showing data
HourData <-HourData %>% select(-weekday, -workingday,-instant)
tbl_df(HourData)

## # A tibble: 2,479 x 12
##    dteday season yr     mnth    hr holiday weathersit  temp atemp   hum
##    <chr>   <int> <fct> <int> <int> <fct>        <int> <dbl> <dbl> <dbl>
##  1 2011-~      1 0         1     0 0                1  0.22 0.197  0.44
##  2 2011-~      1 0         1     1 0                1  0.2  0.167  0.44
##  3 2011-~      1 0         1     4 0                1  0.16 0.136  0.47
##  4 2011-~      1 0         1     5 0                1  0.16 0.136  0.47
##  5 2011-~      1 0         1     6 0                1  0.14 0.106  0.5 
##  6 2011-~      1 0         1     7 0                1  0.14 0.136  0.5 
##  7 2011-~      1 0         1     8 0                1  0.14 0.121  0.5 
##  8 2011-~      1 0         1     9 0                1  0.16 0.136  0.43
##  9 2011-~      1 0         1    10 0                1  0.18 0.167  0.43
## 10 2011-~      1 0         1    11 0                1  0.2  0.182  0.4 
## # ... with 2,469 more rows, and 2 more variables: windspeed <dbl>, cnt <int>

#Separate dataset into train (70%) and test (30%) data set
set.seed(1997)
train <- sample(1:nrow(HourData), size = nrow(HourData)*0.7)
test <- dplyr::setdiff(1:nrow(HourData), train)
HourDataTrain <- HourData[train, ]
HourDataTest <- HourData[test, ]

Summarize the training data

Here I will show you some summary of my training dataset.

1. I conduct a histogram of the rental count, since this is my response variable.
2. I built up a summary table of all the weather measurement.
3. I also showing the weather summary via a boxplot.
4. I plot the rental count distributed by time.
5. I plot the rental count distributed by weather situation.

# plot the histogram of rental count
hist <- ggplot(data=HourDataTrain, aes(x=cnt))+geom_histogram(binwidth = 20, aes(color=yr))
hist <-hist+labs(title="Histogram of the retal count", x="rental count")
hist <-hist+scale_fill_discrete(labels=c(2011,2012))
hist

#prin out summary table for tempature humidity and windspeed
sum <- HourDataTrain%>% select(c(temp, atemp, hum, windspeed))
kable(apply(sum, 2,summary), caption="Numeric Summary for weather measurement")

	temp	atemp	hum	windspeed
Min.	0.0200000	0.0303000	0.2100000	0.0000000
1st Qu.	0.3400000	0.3333000	0.4900000	0.1045000
Median	0.5200000	0.5000000	0.6400000	0.1940000
Mean	0.4961614	0.4763708	0.6329395	0.1894927
3rd Qu.	0.6600000	0.6212000	0.7800000	0.2537000
Max.	0.9200000	0.8485000	1.0000000	0.6866000

#plot the boxplot of tempature humidity and windspeed (not genralized amount)
#plot base
boxplot <- ggplot(data = HourDataTrain, aes(x=season))
#adding 4 variables
tem <-boxplot+geom_boxplot(aes(y=temp*41, group=season))+labs(y="Tempature (c)", title = "boxplot for weather measurement")
fetem <-boxplot+geom_boxplot(aes(y=atemp*50, group=season))+labs(y="Feeling Tempature (c)")
hum <-boxplot+geom_boxplot(aes(y=hum*100, group=season))+labs(y="Humidity")
wind <-boxplot+geom_boxplot(aes(y=windspeed*67, group=season))+labs(y="Wind Speed")
#combine 4 plots into 1
ggarrange(tem, fetem, hum , wind, ncol = 2, nrow = 2)

# plot the count distribution among time and weather
# by time
barplot1<-ggplot(data = HourDataTrain, aes(x=hr))+geom_col(aes(y=cnt, fill=yr))+facet_wrap(~mnth)
barplot1 <- barplot1+labs(x="time", y="Rental Count", title="Retal count distribution by month" )
barplot1+scale_fill_discrete(name="year", labels=c(2011,2012))

# by weather
barplot2 <-ggplot(data = HourDataTrain, aes(x=weathersit))+geom_col(aes(y=cnt, fill=yr))+facet_wrap(~mnth)
barplot2 <- barplot2+labs(x="Weather situation, 1: clear day, 2: misty day, 3:rain or snow", y="Rental Count", title="Retal count distribution by month" )
barplot2+scale_fill_discrete(name="year", labels=c(2011,2012))

Training Model

Here I use two different method to train my model. First method is using a tree-based models with leave one out cross validation. For the second method, I use the boosted tree model with cross validation. Both two training are done using the train function from caret package. The data was cantered and scaled before training.

Tree-based model

Since our respons variable is continuous. I use the regression tree model to training my data. The method= "rpart" was used in train function Moreover, because I want to use the leave-one-out cross validation for this training, therefore,the method= "LOOCV" was used in trainControl. We can adjust the grid parameter by ourselves. Since the default result shows that cp should be very small to have a lowest RMSE. I set a range [0.0001,0.0005] to fit for every weekday. Something to notice, because the cp is too small, when I draw my regression tree, it seems like a mess.

# set up training control, using leave one out cross validation.
set.seed(615)
trctrl <- trainControl(method = "LOOCV", number = 1)

# getModelInfo("rpart")
# training using regression tree models with cp in [0.0001,0.0005]
# since the cp seems have to be really small when I used the default cp to train

model1 <- cnt~season+yr+mnth+hr+holiday+weathersit+temp+atemp+hum+windspeed

RegTree_fit1 <- train(model1, data = HourDataTrain, method = "rpart",
                 trControl=trctrl,
                 preProcess = c("center", "scale"),
                 tuneGrid=expand.grid(cp=seq(0.0001,0.0005,0.00004))
)

# show the training result
RegTree_fit1

## CART 
## 
## 1735 samples
##   10 predictor
## 
## Pre-processing: centered (10), scaled (10) 
## Resampling: Leave-One-Out Cross-Validation 
## Summary of sample sizes: 1734, 1734, 1734, 1734, 1734, 1734, ... 
## Resampling results across tuning parameters:
## 
##   cp       RMSE      Rsquared   MAE     
##   0.00010  74.41616  0.8294464  43.44368
##   0.00014  74.68805  0.8281557  44.10587
##   0.00018  74.66864  0.8282416  44.15918
##   0.00022  74.69121  0.8281267  44.11556
##   0.00026  74.63884  0.8282762  44.10331
##   0.00030  74.79066  0.8275371  44.69143
##   0.00034  74.78446  0.8275581  44.66231
##   0.00038  74.76520  0.8276521  44.63508
##   0.00042  74.50287  0.8287256  44.49473
##   0.00046  74.40986  0.8291284  44.37951
##   0.00050  74.55470  0.8285169  44.60506
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was cp = 0.00046.

# plot the RMSE of selected cp
plot(RegTree_fit1)

# plot my final tree model
rpart.plot(RegTree_fit1$finalModel)

Boosted tree model

Here I want to training my data using boosted tree model. The method= "gbm" was used in train function Because I want to use thecross validation for this training, therefore,the method= "cv" was used in trainControl. We can adjust the grid parameter by ourselves. I set a range of number of tree [100,1250] and interaction 5~11 to fit for every weekday.

# set up training control, using cross validation with 10 folder
set.seed(615)
trctrl <- trainControl(method = "cv", number = 10)

# training using boosted tree models with boosting interation in [700,1250] and try max tree depth 5~9
model2 <- cnt~season+yr+mnth+hr+holiday+weathersit+temp+atemp+hum+windspeed
RegTree_fit2 <- train(model2, data = HourDataTrain, method = "gbm",
                trControl=trctrl,
                preProcess = c("center", "scale"),
                tuneGrid=expand.grid(n.trees=seq(100,1250,25),
                                     interaction.depth=5:11,
                                     shrinkage=0.1, n.minobsinnode=10)
                 )

# show the training result of boosted tree
RegTree_fit2$bestTune

##     n.trees interaction.depth shrinkage n.minobsinnode
## 274    1050                10       0.1             10

# plot the RMSE of different parameters
plot(RegTree_fit2)

Predicting using the best tree-base model

Using the best boosted tree model to testing the data.

# predict use predict function
tree_pred <- predict(RegTree_fit1, newdata = HourDataTest)

#Calculate the Root MSE
RMSE_tree<- sqrt(mean((tree_pred-HourDataTest$cnt)^2))
label <- paste0("RMSE =", RMSE_tree)

# plot the prediction
count <- data.frame(true_count=HourDataTest$cnt,prediction=tree_pred )
predPlot <- ggplot(data=count, aes(x=true_count,y=prediction))
predPlot <- predPlot+labs(title="Prediction V.s. True Count using tree-base model")+geom_point()
predPlot <- predPlot+geom_smooth(color="orange")+geom_abline(aes(intercept=0,slope=1), color="blue")
predPlot <- predPlot+geom_text(x=200, y=600,label=label, color="brown")
predPlot

Predicting using the best boosted-tree model

# predict use predict function
boosted_pred <- predict(RegTree_fit2, newdata = HourDataTest)

#Calculate the Root MSE
RMSE_boosted <- sqrt(mean((boosted_pred-HourDataTest$cnt)^2))
lab <- paste0("RMSE =", RMSE_boosted)
# plot the prediction
count2 <- data.frame(True_count=HourDataTest$cnt,prediction=boosted_pred )
pred_plot <- ggplot(data=count2, aes(x=True_count,y=prediction))
pred_plot <- pred_plot+labs(title="Prediction V.s. True Count using boosted model")+geom_point()
pred_plot <- pred_plot+geom_smooth(color="orange")+geom_abline(aes(intercept=0,slope=1), color="blue")
pred_plot <- pred_plot+geom_text(x=200, y=600,label=lab, color=" brown")
pred_plot

This site is open source. Improve this page.