library(nnet)
library(class)
x <- read.table("https://www.uio.no/studier/emner/matnat/math/STK2100/data/phoneme.data",header=T,sep=",")
#remove row.names
head(x[,c(1:3,256:258)])
x = x[,-1]

Err = matrix(nrow=11,ncol=1)

#Only considering aa and a0
x2 = x[x$g=="aa" | x$g=="ao",]
x2$g = as.factor(x2$g)
x2$g = droplevels(x2$g)
n <- nrow(x2)
set.seed(123)
train <- sample(1:n,as.integer(0.5*n),replace=FALSE)
test <- c(1:n)[-train]
n.train = length(train)
n.test = length(test)
#Remove speaker
X2 = x2[,1:257]
X2$g = as.factor(X2$g)

##TREE
library(tree)
 par(mfrow=c(1,2))

large.tree <- tree(g ~ ., data = X2[train,],
                   control = tree.control(nobs = n.train, mindev = 0.001))
cv.err <- cv.tree(large.tree, K = 10, method = 'deviance')
plot(cv.err)
mod.tree <- prune.tree(large.tree, best = 6)
plot(mod.tree)
text(mod.tree)

train2 <- sample(1:n,as.integer(0.5*n),replace=FALSE)
test2 <- c(1:n)[-train2]
n.train2 = length(train2)
n.test2 = length(test2)
large.tree2 <- tree(g ~ ., data = X2[train2,],
                   control = tree.control(nobs = n.train2, mindev = 0.001))
cv.err2 <- cv.tree(large.tree2, K = 10, method = 'deviance')
plot(cv.err2)
mod.tree2 <- prune.tree(large.tree2, best = 6)
plot(mod.tree2)
text(mod.tree2)

#Bagging
B = 50
pred.bag = matrix(nrow=B,ncol=n.test)
for(b in 1:B)
{
 ind = sample(train,n.train,replace=TRUE)
 large.tree <- tree(g ~ ., data = X2[ind,],
                     control = tree.control(nobs = n.train2, mindev = 0.001))
 mod.tree <- prune.tree(large.tree, best = 6)
 pred.bag[b,] = predict(mod.tree,newdata=x2[test,], type = 'class')
}
pred.bag = pred.bag-1
pred.comb = colMeans(pred.bag)>0.5
cl.comb = rep("aa",n.test)
cl.comb[pred.comb==1] = "ao"
show(table(x2$g[test],cl.comb))
cat("Error rate bagging=",sum(x2$g[test]!=cl.comb)/n.test,"\n")

#Bagging using Random Forest
library(randomForest)
mod.bag = randomForest(g ~ ., data = X2,mtry=256,
                       importance=TRUE,subset=train)
pred.bag = predict(mod.bag,newdata=x2[test,])
show(table(x2$g[test],pred.bag))
cat("Error rate bagging=",sum(x2$g[test]!=pred.bag)/n.test,"\n")
varImpPlot(mod.bag)

#Boosting
# boosting
library(mboost)
##NB: The following two commands are VERY time-consuming
mod.boost <- mboost(g~., data = X2[train,], 
                       family = AdaExp(), baselearner = 'btree')
cv.mstop <- cvrisk(mod.boost)
plot(cv.mstop)
show(mstop(cv.mstop))
pred.boost <- predict(mod.boost, newdata = X2[test,], type = 'class')
show(table(x2$g[test],pred.boost))
cat("Error rate boosting=",sum(X2$g[test]!=pred.boost)/n.test,"\n")
pred.boost2 <- predict(mod.boost[mstop(cv.mstop)], 
                       newdata = X2[test,], type = 'class')
show(table(x2$g[test],pred.boost2))
cat("Error rate boosting2=",sum(X2$g[test]!=pred.boost2)/n.test,"\n")


#Random forest
# random forests
library(randomForest)
library(ranger)
mod.randomForest <- ranger(g~., data = X2[train,],mtry = sqrt(ncol(X2)))
pred.randomForest <- predict(mod.randomForest, data = X2[test,],type='response')$predictions
show(table(x2$g[test],pred.randomForest))
cat("Error rate random Forest=",sum(X2$g[test]!=pred.randomForest)/n.test,"\n")
mod.randomForest <- randomForest(g~., data = X2[train,],mtry = sqrt(ncol(X2)))
pred.randomForest <- predict(mod.randomForest, newdata = X2[test,],type='class')
show(table(x2$g[test],pred.randomForest))
cat("Error rate random Forest=",sum(X2$g[test]!=pred.randomForest)/n.test,"\n")
varImpPlot(mod.randomForest)