---
title: "Lab 7 Bagging and Random Forest"
output: pdf_document
---
  
#### 1)	`survived` is a numeric value. We need to first transform it to a categorical value and saved it as a new variable survived01. Use `titanic3$survived01 = as.factor(titanic3$survived)` to do so and check that this variable has been included in the dataset. 

```{r import, message=FALSE, warning=FALSE}

library(dplyr)

titanic3 <- read.csv("C:\\Users\\tingting\\Dropbox\\modules\\BigData\\sessions\\labs\\Lab7 RF\\titanic3.csv")
titanic3 <- select(titanic3,-name, -ticket, -boat, -body, -home.dest, -cabin) %>% 
  mutate(embarked = factor(embarked), 
         sex = factor(sex),
         pclass = factor(pclass),
         survived01 = as.factor(titanic3$survived))
#embed survived01=as.factor(titanic3$survived) in the mutate() function
#to simplify the code. Use help to see more about mutate().
#Here, it is same to use survived01 = factor(titanic3$survived)

summary(titanic3)

```

#### 2)	Install the package of randomForest and include this package into your code. In order to call the randomForest() function, all the missing value rows need to be dealt with. The simplest way is to remove those rows. Use titanic3 <- na.omit(titanic3) to do that.
```{r, message=FALSE, warning=FALSE}

library(randomForest)
nrow(titanic3)
titanic3 <- na.omit(titanic3)
nrow(titanic3)
```

#### 3)	Use a seed to set half of the dataset to be training dataset and the other half to be test dataset. 

```{r}
set.seed(8)
train <- sample(1:nrow(titanic3), nrow(titanic3)/2)
test <- titanic3[-train,]  #the test set
x_test <- test[,-c(2,9)]  #the predictors in the test set
#-c(2,9) is to remove survived and survived01
# names(titanic3)
#[1] "pclass"     "survived"   "sex"        "age"        "sibsp"  
#[6] "parch"      "fare"       "embarked"   "survived01"

survived.test <- titanic3$survived[-train]
survived01.test <- titanic3$survived01[-train]
```

#### 4)	Use the training dataset to build a bagged model for
-	`y`: survived 
-	`x`: all the features other than survived and survived01.

Compute the mean error rate on the test dataset.

\textbf{Remark:} You might get a warning message, saying that 

\emph{In randomForest.default(m, y, ...) : The response has five or fewer unique values. Are you sure you want to do regression? }

Ignore the message for now. It's doable and you will get a bagged model anyway. 

```{r, message=FALSE, warning=FALSE}
#set.seed(6)
bag.titanic <- randomForest(survived ~ .-survived01, 
                           data = titanic3, 
                           subset = train,
                           mtry=7, 
                           importance = TRUE)
print(bag.titanic)
bag.pred <- predict(bag.titanic, newdata = test)
bag.pred.class <- ifelse(bag.pred <= 0.5, "0", "1")
print(mean(bag.pred.class!=survived.test))

```


#### 5)	Using the same training and test dataset, build a bagged model for
-	y: survived01 
-	x: all the features other than survived and survived01
a) Find out on how many trees your model is built and the OOB error
b) Compute the mean error rate on the test dataset.


```{r}
bag.titanic01 <- randomForest(survived01 ~ .-survived,
                             data = titanic3, 
                             subset = train,
                             mtry = 7, 
                             importance = TRUE)
print(bag.titanic01)
bag.pred01 <- predict(bag.titanic01, newdata = test, type="class")
print(mean(bag.pred01!=survived.test))
bag.titanic01
```

The OOB error is 23.18% and the number of trees grown is 500. Note that OOB error rate is only available for classification trees. 

#### 6) Plot the variable importance plot for the two bagged models you built in 4) and 5) and comment whether the importance coincides.
```{r}
importance(bag.titanic)
varImpPlot(bag.titanic)
importance(bag.titanic01)
varImpPlot(bag.titanic01)

#The barplot can be:

barplot(sort(importance(bag.titanic)[,1], decreasing = TRUE),
        xlab = "Relative Importance",
        horiz = TRUE, 
        col = "red",
        las=1 #The las argument will allow rotation of 90 degrees for labels
        )
```

Yes, the importance of both models coincide if we look at the %IncMSE. 

#### 7)	Plot a graph that shows the test error rate of a single tree (red dashed line), the mean test error rates for majority vote (black curve) and the test error rates for averaging the probabilities (blue curve), both in relation to the number of trees. Add a legend if you can.

<!--
```{r, message=FALSE, warning=FALSE}
#The following is the inefficient approach:
# calculate the black line:
yhat.ter.mv=rep(0,200)      # a vector for Test Error Rate using majority vote (by default)

for(i in 1:200){
  set.seed(6)
  bag.titanic01 = randomForest(survived01 ~ .-survived, data=titanic3, subset = train, mtry=7, importance = TRUE, ntree=i)
  bag.pred01 = predict(bag.titanic01, newdata = test, type="class")
  yhat.ter.mv[i] = mean(bag.pred01!=survived.test)
}
#plot the black line
plot(yhat.ter.mv,xlab="Number of Bootstrap Data Sets",ylab="Test Error Rate", type="l",ylim=c(0.20,0.30), xlim=c(0,205))

#plot the red dashed line
abline(h=yhat.ter.mv[1],lty=2,col="red")


yhat.ter.ave=rep(0,200)     # a vector for Test Error Rate using averaging
for(j in 1:200){
  set.seed(6)
  bag.titanic = randomForest(survived ~ .-survived01, data=titanic3, subset = train, mtry=7, importance = TRUE, ntree=j)
  bag.pred = predict(bag.titanic, newdata = test)
  bag.pred.class = ifelse(bag.pred<=0.5, "0", "1")
  yhat.ter.ave[j] = mean(bag.pred.class!=survived.test)
}


lines(yhat.ter.ave,col="blue")

legend("topright",c("single tree","majority vote","averaging prob"),lty=c(2,1,1),col=c("red","black","blue"))
```

%or -->

```{r, message=FALSE, warning=FALSE}
# calculate the black line:
#Here we insert the test set (xtest and ytest) when building the model. 
#Please read the last three slides of Session 7 for more information.
bag.titanic01=randomForest(survived01 ~ .-survived, 
                           data=titanic3, 
                           subset = train,
                           importance = TRUE, 
                           xtest=x_test, 
                           ytest=survived01.test,
                           mtry=7,
                           ntree=200)


#plot the black line
plot(1:200, bag.titanic01$test$err.rate[,1],
     type="l",
     xlab="Number of Bootstrap Data Sets",
     ylab="Test Error Rate", 
     ylim=c(0.17,0.30), xlim=c(0,205))

#plot the red dashed line
abline(h = bag.titanic01$test$err.rate[1,1],
       lty=2,col="red")


yhat.ter.ave <- rep(0,200)     # a vector for Test Error Rate using averaging
for(j in 1:200){
  #set.seed(6)
  bag.titanic <- randomForest(survived ~ .-survived01,
                             data=titanic3, 
                             subset = train, 
                             mtry=7,
                             importance = TRUE, 
                             ntree=j)
  bag.pred <- predict(bag.titanic, newdata = test)
  bag.pred.class <- ifelse(bag.pred<=0.5, "0", "1")
  yhat.ter.ave[j] <- mean(bag.pred.class!=survived.test)
}

lines(yhat.ter.ave,col="blue")

legend("topright",
       c("single tree","majority vote","averaging prob"),
       lty=c(2,1,1),
       col=c("red","black","blue"))
```

#### 8)	Plot a graph that shows the best value of mtry for the random forest model
-	y: survived01 
-	x: all the features other than survived and survived01
-	mtry: range from 1 to 7

<!--
```{r}
testErrorRate=rep(0,7)
for(i in 1:7){
  bag.titanic01 = randomForest(survived01 ~ .-survived,
                               data=titanic3, subset = train,
                               mtry=i, importance = TRUE, 
                               ntree= 500)
  bag.pred01 = predict(bag.titanic01, newdata = test, type="class")
  testErrorRate[i] = mean(bag.pred01!=survived.test)
}
plot(testErrorRate,type="b",xlab="mtry",ylab="Test Error Rate")
```
%-->

```{r}
testErrorRate <- rep(0,7)
for(i in 1:7){
  set.seed(6)
  bag.titanic01 <- randomForest(survived01 ~ .-survived,
                               data=titanic3, 
                               subset=train,
                               mtry=i, 
                               importance=TRUE,
                               xtest=x_test,
                               ytest=survived01.test,
                               ntree=500)
  testErrorRate[i] <- bag.titanic01$test$err.rate[500,1]
}
plot(testErrorRate,type="b",xlab="mtry",ylab="Test Error Rate")

```

#### 9)	Play with mtry and ntree, plot a graph that shows test error rate vs ntree for different mtry, and find the best/reasonably good combination of mtry and ntree from the plot. Add a legend if you can.

<!--```{r}
plot(0,xlab="Number of Trees",ylab="Test Error Rate", xlim=c(1,540), ylim =c(0.17,0.25))
testER = matrix(nrow=7,ncol=500)

for(i in 1:7){
  for(j in 1:500){
    set.seed(1)
    bag.titanic01 = randomForest(survived01 ~ .-survived,
                                data=titanic3, 
                                subset = train,
                                mtry=i, 
                                ntree=j, 
                                importance = TRUE)
    bag.pred01 = predict(bag.titanic01, 
                         newdata = test, 
                         type="class")
    testER[i,j] = mean(bag.pred01!=survived.test)
  }
  lines(testER[i,],col=i,type="l")
}

legend(title = "mtry", "bottomright",c("1","2","3","4","5","6","7"),lty=rep(1,7),col=1:7)
```

or --->

```{r}
plot(0,
     xlab="Number of Trees",ylab="Test Error Rate", 
     xlim=c(1,540), ylim =c(0.18,0.28))


for(i in 1:7){
  # It's also possible to call randomForest using x and y as the training set.
  # bag.titanic01=randomForest(x=x_train, y=y_train01,
  #                            importance = TRUE, 
  #                            xtest=x_test, ytest=survived01.test,
  #                            mtry=7,ntree=500)
  bag.titanic01=randomForest(survived01 ~ .-survived,
                             data=titanic3, 
                             subset=train,
                             importance = TRUE, 
                             xtest=x_test, 
                             ytest=survived01.test,
                             mtry=i,
                             ntree=500)
  lines(bag.titanic01$test$err.rate[,1],col=i,type="l")
}

legend(title = "mtry",
       "topright",
       c("1","2","3","4","5","6","7"),
       lty=rep(1,7),col=1:7)
```


From the plot, we can see that mtry = 2 and 3 are better than the others in the long run. This result coincides with the empirical result: pick `mtry = sqrt(p)` when it is a classification tree. Here `p = 7` and `sqrt(p) = sqrt(7) = 2.65` $\in [2,3]$.
