library(rpart)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.2.3
library(caret)
## Warning: package 'caret' was built under R version 3.2.2
## Loading required package: lattice
library(randomForest)
## randomForest 4.6-12
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## 
## The following object is masked from 'package:ggplot2':
## 
##     margin
library(e1071)

Questions 1 & 2

data = read.csv('seaflow_21min.csv')
nrow(data)
## [1] 72343
summary(data)
##     file_id           time          cell_id            d1       
##  Min.   :203.0   Min.   : 12.0   Min.   :    0   Min.   : 1328  
##  1st Qu.:204.0   1st Qu.:174.0   1st Qu.: 7486   1st Qu.: 7296  
##  Median :206.0   Median :362.0   Median :14995   Median :17728  
##  Mean   :206.2   Mean   :341.5   Mean   :15008   Mean   :17039  
##  3rd Qu.:208.0   3rd Qu.:503.0   3rd Qu.:22401   3rd Qu.:24512  
##  Max.   :209.0   Max.   :643.0   Max.   :32081   Max.   :54048  
##        d2          fsc_small        fsc_perp        fsc_big     
##  Min.   :   32   Min.   :10005   Min.   :    0   Min.   :32384  
##  1st Qu.: 9584   1st Qu.:31341   1st Qu.:13496   1st Qu.:32400  
##  Median :18512   Median :35483   Median :18069   Median :32400  
##  Mean   :17437   Mean   :34919   Mean   :17646   Mean   :32405  
##  3rd Qu.:24656   3rd Qu.:39184   3rd Qu.:22243   3rd Qu.:32416  
##  Max.   :54688   Max.   :65424   Max.   :63456   Max.   :32464  
##        pe          chl_small        chl_big           pop       
##  Min.   :    0   Min.   : 3485   Min.   :    0   crypto :  102  
##  1st Qu.: 1635   1st Qu.:22525   1st Qu.: 2800   nano   :12698  
##  Median : 2421   Median :30512   Median : 7744   pico   :20860  
##  Mean   : 5325   Mean   :30164   Mean   : 8328   synecho:18146  
##  3rd Qu.: 5854   3rd Qu.:38299   3rd Qu.:12880   ultra  :20537  
##  Max.   :58675   Max.   :64832   Max.   :57184

Question 3

set.seed(000)
train_ind = createDataPartition(data$pop,p=0.5)[[1]]#must access the first (and only) element
#of the list to be used as an index array
train = data[train_ind,]
test = data[-train_ind,]

ans = mean(train[,'time'])
ans
## [1] 340.7668
#340.7668

Question 4

#png("plot1.png",width = 480, height = 480)
ggplot(train, aes(x=chl_small,y=pe,color=pop)) + geom_point(size=1)

#dev.off()
#pico & nano

Questions 5,6,7

fol <- pop ~ fsc_small + fsc_perp + fsc_big + pe + chl_big + chl_small
model <- rpart(fol, method="class", data=train)
print(model)
## n= 36172 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
##  1) root 36172 25742 pico (0.0014 0.18 0.29 0.25 0.28)  
##    2) pe< 5001.5 26341 15958 pico (0 0.22 0.39 0 0.38)  
##      4) chl_small< 32542.5 11789  2173 pico (0 0.00025 0.82 0 0.18) *
##      5) chl_small>=32542.5 14552  6608 ultra (0 0.4 0.053 0 0.55)  
##       10) chl_small>=41057.5 5485   807 nano (0 0.85 0.00018 0 0.15) *
##       11) chl_small< 41057.5 9067  1929 ultra (0 0.13 0.084 0 0.79) *
##    3) pe>=5001.5 9831   758 synecho (0.0052 0.051 0.0048 0.92 0.016)  
##      6) chl_small>=38125.5 642   143 nano (0.079 0.78 0 0.062 0.081) *
##      7) chl_small< 38125.5 9189   156 synecho (0 0.00065 0.0051 0.98 0.011) *
#crypto
#5001.5
#pe, chl_small

Question 8

pred = predict(model, newdata = test)
prediction = (colnames(pred)[max.col(pred,ties.method="first")])
acc = sum(prediction== as.character(test$pop))/length(prediction)
acc
## [1] 0.8544967
#0.8544967

Questions 9

model2 <- randomForest(fol, data=train)
print(model2)
## 
## Call:
##  randomForest(formula = fol, data = train) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 2
## 
##         OOB estimate of  error rate: 7.93%
## Confusion matrix:
##         crypto nano  pico synecho ultra  class.error
## crypto      50    1     0       0     0 0.0196078431
## nano         1 5546     0       4   798 0.1264766105
## pico         0    0 10074      13   343 0.0341323106
## synecho      1    2     0    9068     2 0.0005510856
## ultra        0  357  1340       7  8565 0.1659363132
pred2 = predict(model2, newdata = test)
acc2 = sum(pred2== test$pop)/length(pred2)
acc2
## [1] 0.9199635
#0.9198806

Question10

importance(model2)
##           MeanDecreaseGini
## fsc_small        2775.7407
## fsc_perp         2086.7433
## fsc_big           193.9758
## pe               8981.5464
## chl_big          4762.6839
## chl_small        8054.9010
#pe,chl_small

Question 11

model3 = svm(fol, data=train)
print(model3)
## 
## Call:
## svm(formula = fol, data = train)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  1 
##       gamma:  0.1666667 
## 
## Number of Support Vectors:  9196
pred3 = predict(model3, newdata = test)
acc3 = sum(pred3== test$pop)/length(pred3)
acc3
## [1] 0.9198253
#0.9198253

Question 12

confmat1 = table(pred = prediction, true = test$pop)
confmat2 = table(pred = pred2, true = test$pop)
confmat3 = table(pred = pred3, true = test$pop)

confmat1
##          true
## pred      crypto nano pico synecho ultra
##   nano        51 5153    0      42   818
##   pico         0    4 9614       0  2239
##   synecho      0   19   45    9031   101
##   ultra        0 1173  771       0  7110
confmat2
##          true
## pred      crypto  nano  pico synecho ultra
##   crypto      47     1     0       0     0
##   nano         0  5575     0       0   346
##   pico         0     0 10084       0  1418
##   synecho      4     5     9    9073     7
##   ultra        0   768   337       0  8497
confmat3
##          true
## pred      crypto  nano  pico synecho ultra
##   crypto      44     1     0       0     0
##   nano         2  5628     0       0   370
##   pico         0     0 10065      36  1393
##   synecho      5     8    49    9033     4
##   ultra        0   712   316       4  8501
#ultra is mistaken for pico

Question 13

From

summary(data)
##     file_id           time          cell_id            d1       
##  Min.   :203.0   Min.   : 12.0   Min.   :    0   Min.   : 1328  
##  1st Qu.:204.0   1st Qu.:174.0   1st Qu.: 7486   1st Qu.: 7296  
##  Median :206.0   Median :362.0   Median :14995   Median :17728  
##  Mean   :206.2   Mean   :341.5   Mean   :15008   Mean   :17039  
##  3rd Qu.:208.0   3rd Qu.:503.0   3rd Qu.:22401   3rd Qu.:24512  
##  Max.   :209.0   Max.   :643.0   Max.   :32081   Max.   :54048  
##        d2          fsc_small        fsc_perp        fsc_big     
##  Min.   :   32   Min.   :10005   Min.   :    0   Min.   :32384  
##  1st Qu.: 9584   1st Qu.:31341   1st Qu.:13496   1st Qu.:32400  
##  Median :18512   Median :35483   Median :18069   Median :32400  
##  Mean   :17437   Mean   :34919   Mean   :17646   Mean   :32405  
##  3rd Qu.:24656   3rd Qu.:39184   3rd Qu.:22243   3rd Qu.:32416  
##  Max.   :54688   Max.   :65424   Max.   :63456   Max.   :32464  
##        pe          chl_small        chl_big           pop       
##  Min.   :    0   Min.   : 3485   Min.   :    0   crypto :  102  
##  1st Qu.: 1635   1st Qu.:22525   1st Qu.: 2800   nano   :12698  
##  Median : 2421   Median :30512   Median : 7744   pico   :20860  
##  Mean   : 5325   Mean   :30164   Mean   : 8328   synecho:18146  
##  3rd Qu.: 5854   3rd Qu.:38299   3rd Qu.:12880   ultra  :20537  
##  Max.   :58675   Max.   :64832   Max.   :57184

we can see that the most likely candidate is fsc_big as it spans quite a small range of values. In order to confirm or disprove the assumption we check how many values this variable assumes

levels(as.factor(data$fsc_big))
## [1] "32384" "32400" "32416" "32432" "32448" "32464"
#answer = fsc_big

Question 14

#png("plot2.png",width = 480, height = 480)
ggplot(data, aes(x=time,y=chl_big,color=file_id)) + geom_point(size=1)

#dev.off()
new_data = data[-which(data$file_id == 208),]
set.seed(000)
new_train_ind = createDataPartition(new_data$pop,p=0.5)[[1]]
new_train = new_data[new_train_ind,]
new_test = new_data[-new_train_ind,]
#new_fol <- pop ~ fsc_small + fsc_perp + fsc_big + pe  + chl_small
model3bis = svm(fol, data=new_train)
print(model3bis)
## 
## Call:
## svm(formula = fol, data = new_train)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  1 
##       gamma:  0.1666667 
## 
## Number of Support Vectors:  3765
pred3bis = predict(model3bis, newdata = new_test)
acc3bis = sum(pred3bis== new_test$pop)/length(pred3bis)
acc3bis
## [1] 0.9721605
#0.9721605
change = -(acc3-acc3bis)
change
## [1] 0.05233527
#0.05233527