## Random Forest ##
rm(list=ls()) # clear environment 
Sys.setenv(LANG = "en") # language englisch
library(randomForest)

setwd("~/Documents/Uni/PHD/Statistics/RandomForest/Review")
RFData <- read.csv2("~/Documents/Uni/PHD/Statistics/RandomForest/Review/MelasTraitsMean_Rev.csv", dec = ",")

RFData1=RFData[,c(1:55)] # select columns 

# remove species with only NA 61, 394
RFData1=RFData1[c(1:60,62:393,395:413),]

is.data.frame(RFData1)
str(RFData1)
head(RFData1)

## seperate species with known pollinator from species with unknown pollinator
knownPoll <- subset(RFData1,Syndrome=="Buzz"|Syndrome=="Passerine"|Syndrome=="Nectar foraging vertebrates"|Syndrome=="Generalist"|Syndrome=="Non Buzz")
knownPoll
unknownPoll <- subset(RFData1,Syndrome!="Buzz"&Syndrome!="Passerine"&Syndrome!="Nectar foraging vertebrates"&Syndrome!="Generalist"&Syndrome!="Non Buzz")
unknownPoll


#find species with missing traits
missingKno <- as.data.frame(which(is.na(knownPoll[,c(4:43,45:47,52,53,55)]), arr.ind=TRUE)) # all traits except: 44, 48, 49, 50, 51, 54 
missingunk <- as.data.frame(which(is.na(unknownPoll[,c(4:43,45:47,52,53,55)]), arr.ind=TRUE))

##### Select traits for the model and check for mising data
#remove species with missing traits 
UnkCHA=unknownPoll[c(19,24,90,94),c(4:43,45:47,52:53,55)] 
KnoCHA=knownPoll[c(1:6,8:17,19:22,25:28,30:35,37:44,46:56,59:69,71:76,78:83,86,88,90,91,93:113,115,117:181,183:189,191:215,217,221:236,238,240:253),c(4:43,45:47,52:53,55)] 

## check for NAs in DF
apply(is.na(UnkCHA), 2, which)
apply(is.na(KnoCHA), 2, which)

#
knownPoll1 <- knownPoll[c(1:6,8:17,19:22,25:28,30:35,37:44,46:56,59:69,71:76,78:83,86,88,90,91,93:113,115,117:181,183:189,191:215,217,221:236,238,240:253),c(1,4:43,45:47,52:53,55)] 
unknownPoll1 <- unknownPoll[c(19,24,90,94),c(1,4:43,45:47,52:53,55)] 

spk<-knownPoll1$Taxon
spuk<-unknownPoll1$Taxon
str(knownPoll)

# safe set of selected traits and species with known pollinators
write.csv(KnoCHA, "~/Documents/Uni/PHD/Statistics/RandomForest/Review/knownPoll_NoBuzz.csv")
knownPoll1 <- read.csv("~/Documents/Uni/PHD/Statistics/RandomForest/Review/knownPoll_NoBuzz.csv", row.names=1) ## save file no traits with NA included

# check calsses
sapply(knownPoll1, class)

# transfomre Syndrome from character to factor
knownPoll1Trans <- transform(knownPoll1, Syndrome  = as.factor(as.character(Syndrome)))

####create empty datasets to store data. for all traits with NA removed, all species with known pollinators - knownPoll1##
pred <- data.frame(matrix(nrow=4,ncol=100));colnames(pred)<-btstrap #this is for predicting the unknown
imp_wor <- data.frame(matrix(nrow=227, ncol=2));colnames(imp_wor)<-c("char","Gini") #known character importance - GiniIndex
pred_wor <- data.frame(matrix(nrow=227,ncol=100));colnames(pred_wor)<-btstrap #known, what was predicted, was it wrong?
errrate_wor <-data.frame(matrix(nrow=100,ncol=6));colnames(errrate_wor)<-c("OOB","Buzz","Generalist","Nectar foraging vertebrates","Non Buzz","Passerine") #error rate of known

####make loop to find out most imp characters in differentiating the pollinators - knownPoll1##
for(btstrap in 1:100) 
{
  
  output.forestG <- randomForest(Syndrome~.,data=knownPoll1Trans,ntree=500,proximity=TRUE,importance=TRUE) # syndrome vs. traits
  prered1<-predict(output.forestG,UnkCHA) # use output from syndrome vs traits and predict unknown
  
  pred[,btstrap]<-prered1
  
  pred_wor[,btstrap]<-output.forestG$predicted
  imp_wor <- (importance(output.forestG,type=2)) #Gini index - importance of traits to characterize the syndromes
  write.csv(imp_wor, paste("Imp_wor",btstrap,".csv",sep=""))
  impSyn <- (importance(output.forestG)) # trait importance for all syndromes seperate
  write.csv(impSyn, paste("ImpSyn", btstrap,".csv",sep=""))           
  errrate_wor[btstrap,1]<-median(output.forestG$err.rate[,1]) 
  errrate_wor[btstrap,2]<-median(output.forestG$err.rate[,2]) 
  errrate_wor[btstrap,3]<-median(output.forestG$err.rate[,3]) 
  errrate_wor[btstrap,4]<-median(output.forestG$err.rate[,4])
  errrate_wor[btstrap,5]<-median(output.forestG$err.rate[,5])
  errrate_wor[btstrap,6]<-median(output.forestG$err.rate[,6])
}

write.csv(pred, "bootPredAlltraits_Rev_NoBuzz.csv") ### predicted syndromes for species with unknown pollinator
write.csv(errrate_wor, "Error_predAlltraits_Rev_NoBuzz.csv") ### error rate for predicting known syndromes
write.csv(pred_wor, "KnownPred_Alltraits_Rev_NoBuzz.csv") ### what was predicted for the known pollinators and how often within the 100 btstrap


###### Trait importance
#mean gini overall

imps <- data.frame(matrix(nrow=45,ncol=100)) # 45 = number of traits

code<-1:100
for(j in 1:length(code)){
  a <- read.csv(paste("ImpSyn",code[j],".csv",sep=""),h=T)
  imps[,j]<-a$MeanDecreaseGini
}

code1<-1:45 # number of traits
relimp<-data.frame(matrix(nrow=45,ncol=2));colnames(relimp)<-c("character","meanGini")

relimp[,1]<-a$X

for(i in 1:length(code1)){
  relimp[i,2] <- mean(as.numeric(imps[i,]))
}


OrderNA <-relimp[order(relimp$meanGini, decreasing = TRUE),] ### to write list with traits and importance in numbers
write.csv(OrderNA, "traitImp_MeanGini_Alltraits_Rev.csv")

## plot trait importance in a graph
OrderNA1 <-relimp[order(relimp$meanGini, decreasing = FALSE),]
#colNAMean<-c("gray0","azure3","azure3","azure3","dimgray","gray0","azure3","gray0","azure3","gray0","azure3","gray0","azure3","azure3","azure3","azure3","gray0") # create color vector to distinguish between efficiency, attraction and neutral traits or organ modules

pdf("knownGINI_All_traits_Rev.pdf", width=7, height=8) ### importance of traits to predict syndrome - mean
#quartz()
dotchart(OrderNA1$meanGini, OrderNA1$character, pch=16, cex=0.7, main = "Trait Importance Rev. MeanGini") 
#legend("topleft", legend = c("Efficiency Traits","Attraction Traits","Neutral"),
# pch = 16, cex = 0.5)
dev.size(units = c("cm"))
dev.off()

#mean gini Bee

imps1 <- data.frame(matrix(nrow=45,ncol=100))

code<-1:100
for(j in 1:length(code)){
  a1 <- read.csv(paste("ImpSyn",code[j],".csv",sep=""),h=T)
  imps1[,j]<-a1$Buzz
}

code1<-1:45 # number of traits
relimp1<-data.frame(matrix(nrow=45,ncol=2));colnames(relimp1)<-c("character","meanGini")

relimp1[,1]<-a1$X

for(i in 1:length(code1)){
  relimp1[i,2] <- mean(as.numeric(imps1[i,]))
}

OrderNA2 <-relimp1[order(relimp1$meanGini, decreasing = TRUE),] ### to write
write.csv(OrderNA2, "traitImp_AllTraits_MeanBee_Rev.csv")

OrderNA3 <-relimp1[order(relimp1$meanGini, decreasing = FALSE),]

#colNABee<-c("azure3","gray0","azure3","gray0","azure3","dimgray","azure3","azure3","azure3","gray0","azure3","gray0","gray0","azure3","azure3","azure3","gray0")

pdf("knownGINI_All_traits_Bee_Rev.pdf", width=7, height=7) ### importance of traits to predict syndrome
#quartz()
dotchart(OrderNA3$meanGini, OrderNA3$character, pch=16, cex=0.7, main = "Trait Importance Rev. Bee") #geil! :D
#legend("topleft", legend = c("Efficiency Traits","Attraction Traits","Neutral"),
# col =  c("azure3","gray0","dimgray"),
# pch = 16, cex = 0.5)

dev.size(units = c("cm"))
dev.off()

######################
#mean gini Generalist

imps2 <- data.frame(matrix(nrow=45,ncol=100))

code<-1:100
for(j in 1:length(code)){
  a2 <- read.csv(paste("ImpSyn",code[j],".csv",sep=""),h=T)
  imps2[,j]<-a2$Generalist
}

code1<-1:45 # number of traits
relimp2<-data.frame(matrix(nrow=45,ncol=2));colnames(relimp2)<-c("character","meanGini")

relimp2[,1]<-a2$X

for(i in 1:length(code1)){
  relimp2[i,2] <- mean(as.numeric(imps2[i,]))
}

OrderNA4 <-relimp2[order(relimp2$meanGini, decreasing = TRUE),] ### to write
write.csv(OrderNA4, "traitImp_AllTraits_Generalist_Rev.csv")

OrderNA5 <-relimp2[order(relimp2$meanGini, decreasing = FALSE),]

#colNAPas<-c("gray0","azure3","azure3","azure3","gray0","gray0","gray0","azure3","dimgray","gray0","azure3","azure3","azure3","azure3","azure3","azure3","gray0")

pdf("knownGINI_All_traits_Generalist_Rev.pdf", width=7, height=7) ### importance of traits to predict syndrome
#quartz()
dotchart(OrderNA5$meanGini, OrderNA5$character, pch=16, cex=0.7, main = "Trait Importance Rev. Generalist") 
#legend("topleft", legend = c("Efficiency Traits","Attraction Traits","Neutral"),
#   col =  c("azure3","gray0","dimgray"),
#   pch = 16, cex = 0.5)
dev.size(units = c("cm"))
dev.off()

#mean gini Nectar

imps3 <- data.frame(matrix(nrow=45,ncol=100))

code<-1:100
for(j in 1:length(code)){
  a3 <- read.csv(paste("ImpSyn",code[j],".csv",sep=""),h=T)
  imps3[,j]<-a3$Nectar.foraging.vertebrates
}

code1<-1:45 # number of traits
relimp3<-data.frame(matrix(nrow=45,ncol=2));colnames(relimp3)<-c("character","meanGini")

relimp3[,1]<-a3$X

for(i in 1:length(code1)){
  relimp3[i,2] <- mean(as.numeric(imps3[i,]))
}

OrderNA6 <-relimp3[order(relimp3$meanGini, decreasing = TRUE),] ### to write
write.csv(OrderNA6, "traitImp_AllTraits_Nectar_Rev.csv")

OrderNA7 <-relimp3[order(relimp3$meanGini, decreasing = FALSE),]

#colNAMix<-c("azure3","gray0","azure3","gray0","gray0","azure3","azure3","dimgray","gray0","gray0","azure3","azure3","azure3","azure3","azure3","azure3","gray0")

pdf("knownGINI_All_traits_Nectar_Rev.pdf", width=7, height=7) ### importance of traits to predict syndrome
#quartz()
dotchart(OrderNA7$meanGini, OrderNA7$character,  pch=16, cex=0.7, main = "Trait Importance Rev. Nectar") 
#legend("topleft", legend = c("Efficiency Traits","Attraction Traits","Neutral"),
#  col =  c("azure3","gray0","dimgray"),
# pch = 16, cex = 0.5)
dev.size(units = c("cm"))
dev.off()


#mean gini Passerine

imps4 <- data.frame(matrix(nrow=45,ncol=100))

code<-1:100
for(j in 1:length(code)){
  a4 <- read.csv(paste("ImpSyn",code[j],".csv",sep=""),h=T)
  imps4[,j]<-a4$Passerine
}

code1<-1:45 # number of traits
relimp4<-data.frame(matrix(nrow=45,ncol=2));colnames(relimp4)<-c("character","meanGini")

relimp4[,1]<-a4$X

for(i in 1:length(code1)){
  relimp4[i,2] <- mean(as.numeric(imps4[i,]))
}

OrderNA7 <-relimp4[order(relimp4$meanGini, decreasing = TRUE),] ### to write
write.csv(OrderNA7, "traitImp_AllTraits_Passerine_Rev.csv")

OrderNA8 <-relimp4[order(relimp4$meanGini, decreasing = FALSE),]

#colNAMix<-c("azure3","gray0","azure3","gray0","gray0","azure3","azure3","dimgray","gray0","gray0","azure3","azure3","azure3","azure3","azure3","azure3","gray0")

pdf("knownGINI_All_traits_Passerine_Rev.pdf", width=7, height=7) ### importance of traits to predict syndrome
#quartz()
dotchart(OrderNA8$meanGini, OrderNA8$character,  pch=16, cex=0.7, main = "Trait Importance Rev. Passerine") 
#legend("topleft", legend = c("Efficiency Traits","Attraction Traits","Neutral"),
#  col =  c("azure3","gray0","dimgray"),
# pch = 16, cex = 0.5)
dev.size(units = c("cm"))
dev.off()

######### SummarizePrediction

### summary for unknown syndrome predictions

#import prediction files for unknown pollinators and check how often one value is included
bootPred <- read.csv("~/Documents/Uni/PHD/Statistics/RandomForest/Review/bootPredAlltraits_Rev.csv", row.names=1) 

code<-c("Buzz","Generalist","Nectar foraging vertebrates","Passerine")
pollpreds <- data.frame(matrix(nrow=5, ncol=4));colnames(pollpreds)<-c("Buzz","Generalist","Nectar foraging vertebrates","Passerine");rownames(pollpreds)<-spuk
zeile<-1:5

for(j in 1:length(zeile)){
  for(i in 1:length(code)){
    pollpreds[j,1]<-length(which(bootPred[j,]=="Buzz"))/100
    pollpreds[j,2]<-length(which(bootPred[j,]=="Generalist"))/100
    pollpreds[j,3]<-length(which(bootPred[j,]=="Nectar foraging vertebrates"))/100
    pollpreds[j,4]<-length(which(bootPred[j,]=="Passerine"))/100
    
  }}
pollpreds

write.csv(pollpreds, "sumPredAlltraits_Rev.csv") 

#### summary for known pollinator predictions

KnownPred <- read.csv("~/Documents/Uni/PHD/Statistics/RandomForest/Review/KnownPred_Alltraits_Rev.csv", row.names=1)
polls<-knownPoll1$Syndrome

pollpreds2 <- data.frame(matrix(nrow=225, ncol=2));colnames(pollpreds2)<-c("poll","p");rownames(pollpreds2)<-spk ## create empty data frame
pollpreds2[1]<-polls 

pollpreds2[1,2]<-length(which(KnownPred[1,]=="Buzz"))/100
pollpreds2[2,2]<-length(which(KnownPred[2,]=="Buzz"))/100
pollpreds2[3,2]<-length(which(KnownPred[3,]=="Buzz"))/100
pollpreds2[4,2]<-length(which(KnownPred[4,]=="Buzz"))/100
pollpreds2[5,2]<-length(which(KnownPred[5,]=="Buzz"))/100
pollpreds2[6,2]<-length(which(KnownPred[6,]=="Buzz"))/100
pollpreds2[7,2]<-length(which(KnownPred[7,]=="Passerine"))/100
pollpreds2[8,2]<-length(which(KnownPred[8,]=="Passerine"))/100
pollpreds2[9,2]<-length(which(KnownPred[9,]=="Passerine"))/100
pollpreds2[10,2]<-length(which(KnownPred[10,]=="Passerine"))/100
pollpreds2[11,2]<-length(which(KnownPred[11,]=="Buzz"))/100
pollpreds2[12,2]<-length(which(KnownPred[12,]=="Buzz"))/100
pollpreds2[13,2]<-length(which(KnownPred[13,]=="Buzz"))/100
pollpreds2[14,2]<-length(which(KnownPred[14,]=="Buzz"))/100
pollpreds2[15,2]<-length(which(KnownPred[15,]=="Buzz"))/100
pollpreds2[16,2]<-length(which(KnownPred[16,]=="Buzz"))/100
pollpreds2[17,2]<-length(which(KnownPred[17,]=="Buzz"))/100
pollpreds2[18,2]<-length(which(KnownPred[18,]=="Buzz"))/100
pollpreds2[19,2]<-length(which(KnownPred[19,]=="Buzz"))/100
pollpreds2[20,2]<-length(which(KnownPred[20,]=="Buzz"))/100
pollpreds2[21,2]<-length(which(KnownPred[21,]=="Nectar foraging vertebrates"))/100
pollpreds2[22,2]<-length(which(KnownPred[22,]=="Buzz"))/100
pollpreds2[23,2]<-length(which(KnownPred[23,]=="Buzz"))/100
pollpreds2[24,2]<-length(which(KnownPred[24,]=="Buzz"))/100
pollpreds2[25,2]<-length(which(KnownPred[25,]=="Buzz"))/100
pollpreds2[26,2]<-length(which(KnownPred[26,]=="Buzz"))/100
pollpreds2[27,2]<-length(which(KnownPred[27,]=="Buzz"))/100
pollpreds2[28,2]<-length(which(KnownPred[28,]=="Nectar foraging vertebrates"))/100
pollpreds2[29,2]<-length(which(KnownPred[29,]=="Buzz"))/100
pollpreds2[30,2]<-length(which(KnownPred[30,]=="Buzz"))/100
pollpreds2[31,2]<-length(which(KnownPred[31,]=="Nectar foraging vertebrates"))/100
pollpreds2[32,2]<-length(which(KnownPred[32,]=="Nectar foraging vertebrates"))/100
pollpreds2[33,2]<-length(which(KnownPred[33,]=="Nectar foraging vertebrates"))/100
pollpreds2[34,2]<-length(which(KnownPred[34,]=="Nectar foraging vertebrates"))/100
pollpreds2[35,2]<-length(which(KnownPred[35,]=="Buzz"))/100
pollpreds2[36,2]<-length(which(KnownPred[36,]=="Buzz"))/100
pollpreds2[37,2]<-length(which(KnownPred[37,]=="Buzz"))/100
pollpreds2[38,2]<-length(which(KnownPred[38,]=="Buzz"))/100
pollpreds2[39,2]<-length(which(KnownPred[39,]=="Buzz"))/100
pollpreds2[40,2]<-length(which(KnownPred[40,]=="Buzz"))/100
pollpreds2[41,2]<-length(which(KnownPred[41,]=="Buzz"))/100
pollpreds2[42,2]<-length(which(KnownPred[42,]=="Nectar foraging vertebrates"))/100
pollpreds2[43,2]<-length(which(KnownPred[43,]=="Buzz"))/100
pollpreds2[44,2]<-length(which(KnownPred[44,]=="Buzz"))/100
pollpreds2[45,2]<-length(which(KnownPred[45,]=="Buzz"))/100
pollpreds2[46,2]<-length(which(KnownPred[46,]=="Buzz"))/100
pollpreds2[47,2]<-length(which(KnownPred[47,]=="Buzz"))/100
pollpreds2[48,2]<-length(which(KnownPred[48,]=="Buzz"))/100
pollpreds2[49,2]<-length(which(KnownPred[49,]=="Nectar foraging vertebrates"))/100
pollpreds2[50,2]<-length(which(KnownPred[50,]=="Buzz"))/100
pollpreds2[51,2]<-length(which(KnownPred[51,]=="Buzz"))/100
pollpreds2[52,2]<-length(which(KnownPred[52,]=="Buzz"))/100
pollpreds2[53,2]<-length(which(KnownPred[53,]=="Buzz"))/100
pollpreds2[54,2]<-length(which(KnownPred[54,]=="Buzz"))/100
pollpreds2[55,2]<-length(which(KnownPred[55,]=="Buzz"))/100
pollpreds2[56,2]<-length(which(KnownPred[56,]=="Buzz"))/100
pollpreds2[57,2]<-length(which(KnownPred[57,]=="Buzz"))/100
pollpreds2[58,2]<-length(which(KnownPred[58,]=="Buzz"))/100
pollpreds2[59,2]<-length(which(KnownPred[59,]=="Buzz"))/100
pollpreds2[60,2]<-length(which(KnownPred[60,]=="Buzz"))/100
pollpreds2[61,2]<-length(which(KnownPred[61,]=="Buzz"))/100
pollpreds2[62,2]<-length(which(KnownPred[62,]=="Buzz"))/100
pollpreds2[63,2]<-length(which(KnownPred[63,]=="Buzz"))/100
pollpreds2[64,2]<-length(which(KnownPred[64,]=="Buzz"))/100
pollpreds2[65,2]<-length(which(KnownPred[65,]=="Buzz"))/100
pollpreds2[66,2]<-length(which(KnownPred[66,]=="Buzz"))/100
pollpreds2[67,2]<-length(which(KnownPred[67,]=="Buzz"))/100
pollpreds2[68,2]<-length(which(KnownPred[68,]=="Buzz"))/100
pollpreds2[69,2]<-length(which(KnownPred[69,]=="Buzz"))/100
pollpreds2[70,2]<-length(which(KnownPred[70,]=="Buzz"))/100
pollpreds2[71,2]<-length(which(KnownPred[71,]=="Nectar foraging vertebrates"))/100
pollpreds2[72,2]<-length(which(KnownPred[72,]=="Buzz"))/100
pollpreds2[73,2]<-length(which(KnownPred[73,]=="Buzz"))/100
pollpreds2[74,2]<-length(which(KnownPred[74,]=="Buzz"))/100
pollpreds2[75,2]<-length(which(KnownPred[75,]=="Nectar foraging vertebrates"))/100
pollpreds2[76,2]<-length(which(KnownPred[76,]=="Buzz"))/100
pollpreds2[77,2]<-length(which(KnownPred[77,]=="Nectar foraging vertebrates"))/100
pollpreds2[78,2]<-length(which(KnownPred[78,]=="Buzz"))/100
pollpreds2[79,2]<-length(which(KnownPred[79,]=="Buzz"))/100
pollpreds2[80,2]<-length(which(KnownPred[80,]=="Passerine"))/100
pollpreds2[81,2]<-length(which(KnownPred[81,]=="Buzz"))/100
pollpreds2[82,2]<-length(which(KnownPred[82,]=="Buzz"))/100
pollpreds2[83,2]<-length(which(KnownPred[83,]=="Buzz"))/100
pollpreds2[84,2]<-length(which(KnownPred[84,]=="Passerine"))/100
pollpreds2[85,2]<-length(which(KnownPred[85,]=="Nectar foraging vertebrates"))/100
pollpreds2[86,2]<-length(which(KnownPred[86,]=="Nectar foraging vertebrates"))/100
pollpreds2[87,2]<-length(which(KnownPred[87,]=="Nectar foraging vertebrates"))/100
pollpreds2[88,2]<-length(which(KnownPred[88,]=="Nectar foraging vertebrates"))/100
pollpreds2[89,2]<-length(which(KnownPred[89,]=="Buzz"))/100
pollpreds2[90,2]<-length(which(KnownPred[90,]=="Nectar foraging vertebrates"))/100
pollpreds2[91,2]<-length(which(KnownPred[91,]=="Buzz"))/100
pollpreds2[92,2]<-length(which(KnownPred[92,]=="Buzz"))/100
pollpreds2[93,2]<-length(which(KnownPred[93,]=="Buzz"))/100
pollpreds2[94,2]<-length(which(KnownPred[94,]=="Buzz"))/100
pollpreds2[95,2]<-length(which(KnownPred[95,]=="Buzz"))/100
pollpreds2[96,2]<-length(which(KnownPred[96,]=="Buzz"))/100
pollpreds2[97,2]<-length(which(KnownPred[97,]=="Buzz"))/100
pollpreds2[98,2]<-length(which(KnownPred[98,]=="Buzz"))/100
pollpreds2[99,2]<-length(which(KnownPred[99,]=="Generalist"))/100
pollpreds2[100,2]<-length(which(KnownPred[100,]=="Buzz"))/100
pollpreds2[101,2]<-length(which(KnownPred[101,]=="Buzz"))/100
pollpreds2[102,2]<-length(which(KnownPred[102,]=="Buzz"))/100
pollpreds2[103,2]<-length(which(KnownPred[103,]=="Buzz"))/100
pollpreds2[104,2]<-length(which(KnownPred[104,]=="Buzz"))/100
pollpreds2[105,2]<-length(which(KnownPred[105,]=="Nectar foraging vertebrates"))/100
pollpreds2[106,2]<-length(which(KnownPred[106,]=="Buzz"))/100
pollpreds2[107,2]<-length(which(KnownPred[107,]=="Buzz"))/100
pollpreds2[108,2]<-length(which(KnownPred[108,]=="Generalist"))/100
pollpreds2[109,2]<-length(which(KnownPred[109,]=="Nectar foraging vertebrates"))/100
pollpreds2[110,2]<-length(which(KnownPred[110,]=="Buzz"))/100
pollpreds2[111,2]<-length(which(KnownPred[111,]=="Buzz"))/100
pollpreds2[112,2]<-length(which(KnownPred[112,]=="Buzz"))/100
pollpreds2[113,2]<-length(which(KnownPred[113,]=="Buzz"))/100
pollpreds2[114,2]<-length(which(KnownPred[114,]=="Buzz"))/100
pollpreds2[115,2]<-length(which(KnownPred[115,]=="Nectar foraging vertebrates"))/100
pollpreds2[116,2]<-length(which(KnownPred[116,]=="Buzz"))/100
pollpreds2[117,2]<-length(which(KnownPred[117,]=="Buzz"))/100
pollpreds2[118,2]<-length(which(KnownPred[118,]=="Buzz"))/100
pollpreds2[119,2]<-length(which(KnownPred[119,]=="Nectar foraging vertebrates"))/100
pollpreds2[120,2]<-length(which(KnownPred[120,]=="Buzz"))/100
pollpreds2[121,2]<-length(which(KnownPred[121,]=="Buzz"))/100
pollpreds2[122,2]<-length(which(KnownPred[122,]=="Buzz"))/100
pollpreds2[123,2]<-length(which(KnownPred[123,]=="Buzz"))/100
pollpreds2[124,2]<-length(which(KnownPred[124,]=="Nectar foraging vertebrates"))/100
pollpreds2[125,2]<-length(which(KnownPred[125,]=="Buzz"))/100
pollpreds2[126,2]<-length(which(KnownPred[126,]=="Generalist"))/100
pollpreds2[127,2]<-length(which(KnownPred[127,]=="Buzz"))/100
pollpreds2[128,2]<-length(which(KnownPred[128,]=="Buzz"))/100
pollpreds2[129,2]<-length(which(KnownPred[129,]=="Buzz"))/100
pollpreds2[130,2]<-length(which(KnownPred[130,]=="Buzz"))/100
pollpreds2[131,2]<-length(which(KnownPred[131,]=="Buzz"))/100
pollpreds2[132,2]<-length(which(KnownPred[132,]=="Generalist"))/100
pollpreds2[133,2]<-length(which(KnownPred[133,]=="Buzz"))/100
pollpreds2[134,2]<-length(which(KnownPred[134,]=="Nectar foraging vertebrates"))/100
pollpreds2[135,2]<-length(which(KnownPred[135,]=="Buzz"))/100
pollpreds2[136,2]<-length(which(KnownPred[136,]=="Generalist"))/100
pollpreds2[137,2]<-length(which(KnownPred[137,]=="Nectar foraging vertebrates"))/100
pollpreds2[138,2]<-length(which(KnownPred[138,]=="Buzz"))/100
pollpreds2[139,2]<-length(which(KnownPred[139,]=="Buzz"))/100
pollpreds2[140,2]<-length(which(KnownPred[140,]=="Buzz"))/100
pollpreds2[141,2]<-length(which(KnownPred[141,]=="Buzz"))/100
pollpreds2[142,2]<-length(which(KnownPred[142,]=="Buzz"))/100
pollpreds2[143,2]<-length(which(KnownPred[143,]=="Buzz"))/100
pollpreds2[144,2]<-length(which(KnownPred[144,]=="Buzz"))/100
pollpreds2[145,2]<-length(which(KnownPred[145,]=="Buzz"))/100
pollpreds2[146,2]<-length(which(KnownPred[146,]=="Buzz"))/100
pollpreds2[147,2]<-length(which(KnownPred[147,]=="Buzz"))/100
pollpreds2[148,2]<-length(which(KnownPred[148,]=="Buzz"))/100
pollpreds2[149,2]<-length(which(KnownPred[149,]=="Buzz"))/100
pollpreds2[150,2]<-length(which(KnownPred[150,]=="Generalist"))/100
pollpreds2[151,2]<-length(which(KnownPred[151,]=="Buzz"))/100
pollpreds2[152,2]<-length(which(KnownPred[152,]=="Buzz"))/100
pollpreds2[153,2]<-length(which(KnownPred[153,]=="Buzz"))/100
pollpreds2[154,2]<-length(which(KnownPred[154,]=="Buzz"))/100
pollpreds2[155,2]<-length(which(KnownPred[155,]=="Buzz"))/100
pollpreds2[156,2]<-length(which(KnownPred[156,]=="Buzz"))/100
pollpreds2[157,2]<-length(which(KnownPred[157,]=="Buzz"))/100
pollpreds2[158,2]<-length(which(KnownPred[158,]=="Buzz"))/100
pollpreds2[159,2]<-length(which(KnownPred[159,]=="Buzz"))/100
pollpreds2[160,2]<-length(which(KnownPred[160,]=="Buzz"))/100
pollpreds2[161,2]<-length(which(KnownPred[161,]=="Buzz"))/100
pollpreds2[162,2]<-length(which(KnownPred[162,]=="Buzz"))/100
pollpreds2[163,2]<-length(which(KnownPred[163,]=="Buzz"))/100
pollpreds2[164,2]<-length(which(KnownPred[164,]=="Buzz"))/100
pollpreds2[165,2]<-length(which(KnownPred[165,]=="Generalist"))/100
pollpreds2[166,2]<-length(which(KnownPred[166,]=="Buzz"))/100
pollpreds2[167,2]<-length(which(KnownPred[167,]=="Buzz"))/100
pollpreds2[168,2]<-length(which(KnownPred[168,]=="Buzz"))/100
pollpreds2[169,2]<-length(which(KnownPred[169,]=="Generalist"))/100
pollpreds2[170,2]<-length(which(KnownPred[170,]=="Buzz"))/100
pollpreds2[171,2]<-length(which(KnownPred[171,]=="Buzz"))/100
pollpreds2[172,2]<-length(which(KnownPred[172,]=="Buzz"))/100
pollpreds2[173,2]<-length(which(KnownPred[173,]=="Buzz"))/100
pollpreds2[174,2]<-length(which(KnownPred[174,]=="Buzz"))/100
pollpreds2[175,2]<-length(which(KnownPred[175,]=="Buzz"))/100
pollpreds2[176,2]<-length(which(KnownPred[176,]=="Buzz"))/100
pollpreds2[177,2]<-length(which(KnownPred[177,]=="Buzz"))/100
pollpreds2[178,2]<-length(which(KnownPred[178,]=="Buzz"))/100
pollpreds2[179,2]<-length(which(KnownPred[179,]=="Buzz"))/100
pollpreds2[180,2]<-length(which(KnownPred[180,]=="Buzz"))/100
pollpreds2[181,2]<-length(which(KnownPred[181,]=="Buzz"))/100
pollpreds2[182,2]<-length(which(KnownPred[182,]=="Buzz"))/100
pollpreds2[183,2]<-length(which(KnownPred[183,]=="Buzz"))/100
pollpreds2[184,2]<-length(which(KnownPred[184,]=="Buzz"))/100
pollpreds2[185,2]<-length(which(KnownPred[185,]=="Buzz"))/100
pollpreds2[186,2]<-length(which(KnownPred[186,]=="Buzz"))/100
pollpreds2[187,2]<-length(which(KnownPred[187,]=="Buzz"))/100
pollpreds2[188,2]<-length(which(KnownPred[188,]=="Buzz"))/100
pollpreds2[189,2]<-length(which(KnownPred[189,]=="Buzz"))/100
pollpreds2[190,2]<-length(which(KnownPred[190,]=="Buzz"))/100
pollpreds2[191,2]<-length(which(KnownPred[191,]=="Buzz"))/100
pollpreds2[192,2]<-length(which(KnownPred[192,]=="Buzz"))/100
pollpreds2[193,2]<-length(which(KnownPred[193,]=="Buzz"))/100
pollpreds2[194,2]<-length(which(KnownPred[194,]=="Buzz"))/100
pollpreds2[195,2]<-length(which(KnownPred[195,]=="Buzz"))/100
pollpreds2[196,2]<-length(which(KnownPred[196,]=="Buzz"))/100
pollpreds2[197,2]<-length(which(KnownPred[197,]=="Buzz"))/100
pollpreds2[198,2]<-length(which(KnownPred[198,]=="Buzz"))/100
pollpreds2[199,2]<-length(which(KnownPred[199,]=="Buzz"))/100
pollpreds2[200,2]<-length(which(KnownPred[200,]=="Buzz"))/100
pollpreds2[201,2]<-length(which(KnownPred[201,]=="Buzz"))/100
pollpreds2[202,2]<-length(which(KnownPred[202,]=="Buzz"))/100
pollpreds2[203,2]<-length(which(KnownPred[203,]=="Buzz"))/100
pollpreds2[204,2]<-length(which(KnownPred[204,]=="Buzz"))/100
pollpreds2[205,2]<-length(which(KnownPred[205,]=="Buzz"))/100
pollpreds2[206,2]<-length(which(KnownPred[206,]=="Buzz"))/100
pollpreds2[207,2]<-length(which(KnownPred[207,]=="Buzz"))/100
pollpreds2[208,2]<-length(which(KnownPred[208,]=="Buzz"))/100
pollpreds2[209,2]<-length(which(KnownPred[209,]=="Buzz"))/100
pollpreds2[210,2]<-length(which(KnownPred[210,]=="Buzz"))/100
pollpreds2[211,2]<-length(which(KnownPred[211,]=="Buzz"))/100
pollpreds2[212,2]<-length(which(KnownPred[212,]=="Buzz"))/100
pollpreds2[213,2]<-length(which(KnownPred[213,]=="Buzz"))/100
pollpreds2[214,2]<-length(which(KnownPred[214,]=="Buzz"))/100
pollpreds2[215,2]<-length(which(KnownPred[215,]=="Buzz"))/100
pollpreds2[216,2]<-length(which(KnownPred[216,]=="Buzz"))/100
pollpreds2[217,2]<-length(which(KnownPred[217,]=="Buzz"))/100
pollpreds2[218,2]<-length(which(KnownPred[218,]=="Buzz"))/100
pollpreds2[219,2]<-length(which(KnownPred[219,]=="Buzz"))/100
pollpreds2[220,2]<-length(which(KnownPred[220,]=="Buzz"))/100
pollpreds2[221,2]<-length(which(KnownPred[221,]=="Buzz"))/100
pollpreds2[222,2]<-length(which(KnownPred[222,]=="Buzz"))/100
pollpreds2[223,2]<-length(which(KnownPred[223,]=="Buzz"))/100
pollpreds2[224,2]<-length(which(KnownPred[224,]=="Buzz"))/100
pollpreds2[225,2]<-length(which(KnownPred[225,]=="Buzz"))/100

write.csv(pollpreds2, "PollObs_pred_Alltraites_rev.csv") 

# get alternative predictions for species with lower than 100% (Memecylon umbellatum, Miconia brevitheca, Miconia corymbiformis, Miconia ligustroides, Miconia pepericarpa, Miconia pusilliflora)
taxon_name <- c("Miconia brevitheca",
  "Miconia corymbiformis",
  "Miconia ligustroides",
  "Miconia pusilliflora"
)

# Get row number(s) where the taxa matches 'taxon_name'
matching_rows <- which(rownames(pollpreds2) %in% taxon_name)

# extract species with predictions from KnownPred
altPred <- KnownPred[matching_rows, ]
  
# Define the terms you want to count
terms <- c("Buzz", "Generalist", "Nectar foraging vertebrates", "Passerine")

# Create an empty dataframe with column names
percent_df <- data.frame(matrix(NA, nrow = 4, ncol = length(terms)))
colnames(percent_df) <- terms

# Miconia brevitheca
percent_df[1,1]<-length(which(altPred[1,]=="Buzz"))/100
percent_df[1,2]<- ""
percent_df[1,3]<-length(which(altPred[1,]=="Nectar foraging vertebrates"))/100
percent_df[1,4]<-length(which(altPred[1,]=="Passerine"))/100

# Miconia corymbiformis
percent_df[2,1]<-length(which(altPred[2,]=="Buzz"))/100
percent_df[2,2]<- ""
percent_df[2,3]<-length(which(altPred[2,]=="Nectar foraging vertebrates"))/100
percent_df[2,4]<-length(which(altPred[2,]=="Passerine"))/100

# Miconia ligustroides
percent_df[3,1]<-length(which(altPred[3,]=="Buzz"))/100
percent_df[3,2]<- ""
percent_df[3,3]<-length(which(altPred[3,]=="Nectar foraging vertebrates"))/100
percent_df[3,4]<-length(which(altPred[3,]=="Passerine"))/100

# Miconia pusilliflora
percent_df[4,1]<-length(which(altPred[4,]=="Buzz"))/100
percent_df[4,2]<- ""
percent_df[4,3]<-length(which(altPred[4,]=="Nectar foraging vertebrates"))/100
percent_df[4,4]<-length(which(altPred[4,]=="Passerine"))/100

rownames(percent_df) <- taxon_name

write.csv(percent_df, "~/Documents/Uni/PHD/Statistics/RandomForest/Review/AlternativePredictions_Rev.csv")

#Error file

errrate_wor<-read.csv("~/Documents/Uni/PHD/Statistics/RandomForest/Review/Error_predAlltraits_Rev_NoBuzz.csv")

errorRate <- data.frame(matrix(nrow=1,ncol=6));colnames(errorRate)<-c("OOB","Buzz","Generalist","Nectar foraging vertebrates","Non Buzz","Passerine")

errorRate[1,1] <- median(errrate_wor$OOB) 
errorRate[1,2] <- median(errrate_wor$Buzz,na.rm=TRUE) # median error rate for predicting the bee syndrome right for species with known bee syndrome
errorRate[1,3] <- median(errrate_wor$Generalist,na.rm=TRUE) # median error rate for predicting the generalist syndrome right for species with known bee syndrome 
errorRate[1,4] <- median(errrate_wor$Nectar.foraging.vertebrates,na.rm=TRUE) # median error rate for predicting the nectar foraging vertebrates syndrome right for species with known bee syndrome 
errorRate[1,5] <- median(errrate_wor$Non.Buzz,na.rm=TRUE)
errorRate[1,6] <- median(errrate_wor$Passerine,na.rm=TRUE) # median error rate for predicting the Passerine syndrome right for species with known bee syndrome 


write.csv(errorRate, "~/Documents/Uni/PHD/Statistics/RandomForest/Review/Error_perSyn_Alltraits_Rev.csv")

### test if generalists are a syndrome
# move all generalists in the buzz and compare pred acc. 

# rename generalist to buzz
library(dplyr)

RFData2 <- RFData1 %>%
  mutate(Syndrome = recode(Syndrome, "Generalist" = "Buzz"))

## seperate species with known pollinator from species with unknown pollinator

knownPoll <- subset(RFData2,Syndrome=="Buzz"|Syndrome=="Passerine"|Syndrome=="Nectar foraging vertebrates")
knownPoll
unknownPoll <- subset(RFData2,Syndrome!="Buzz"&Syndrome!="Passerine"&Syndrome!="Nectar foraging vertebrates")
unknownPoll

##### Select traits for the model and check for mising data

#find species with missing traits
missingKno <- as.data.frame(which(is.na(knownPoll[,c(4:43,45:47,52,53,55)]), arr.ind=TRUE)) # all traits except: 44, 48, 49, 50, 51, 54 
missingunk <- as.data.frame(which(is.na(unknownPoll[,c(4:43,45:47,52,53,55)]), arr.ind=TRUE))

#remove species with missing traits 
UnkCHA=unknownPoll[c(19,24,88,91,95),c(4:43,45:47,52:53,55)] 
KnoCHA=knownPoll[c(1:6,8:17,19:22,25:28,30:35,37:44,46:56,59:69,71:76,78:83,86,88,90,92:111,114,116:180,182:188,190:214,216,220:235,237,239:252),c(4:43,45:47,52:53,55)] 

knownPoll1 <- knownPoll[c(1:6,8:17,19:22,25:28,30:35,37:44,46:56,59:69,71:76,78:83,86,88,90,92:111,114,116:180,182:188,190:214,216,220:235,237,239:252),c(1,4:43,45:47,52:53,55)] 
unknownPoll1 <- unknownPoll[c(19,24,88,91,95),c(1,4:43,45:47,52:53,55)] 

spk<-knownPoll1$Taxon
spuk<-unknownPoll1$Taxon
str(knownPoll)

## check for NAs in DF
apply(is.na(UnkCHA), 2, which)
apply(is.na(KnoCHA), 2, which)

# safe set of selected traits and species with known pollinators
write.csv(KnoCHA, "~/Documents/Uni/PHD/Statistics/RandomForest//knownPoll_GentoBuzz.csv")
knownPoll1 <- read.csv("~/Documents/Uni/PHD/Statistics/RandomForest//knownPoll_GentoBuzz.csv", row.names=1) ## save file no traits with NA included

# check calsses
sapply(knownPoll1, class)

# transfomre Syndrome from character to factor
knownPoll1Trans <- transform(knownPoll1, Syndrome  = as.factor(as.character(Syndrome)))

####create empty datasets to store data. for all traits with NA removed, all species with known pollinators - knownPoll1##
pred <- data.frame(matrix(nrow=5,ncol=100));colnames(pred)<-btstrap #this is for predicting the unknown
imp_wor <- data.frame(matrix(nrow=225, ncol=2));colnames(imp_wor)<-c("char","Gini") #known character importance - GiniIndex
pred_wor <- data.frame(matrix(nrow=225,ncol=100));colnames(pred_wor)<-btstrap #known, what was predicted, was it wrong?
errrate_wor <-data.frame(matrix(nrow=100,ncol=4));colnames(errrate_wor)<-c("OOB","Buzz","Nectar foraging vertebrates","Passerine") #error rate of known

####make loop to find out most imp characters in differentiating the pollinators - knownPoll1##
for(btstrap in 1:100) 
{
  
  output.forestG <- randomForest(Syndrome~.,data=knownPoll1Trans,ntree=500,proximity=TRUE,importance=TRUE) # syndrome vs. traits
  prered1<-predict(output.forestG,UnkCHA) # use output from syndrome vs traits and predict unknown
  
  pred[,btstrap]<-prered1
  
  pred_wor[,btstrap]<-output.forestG$predicted
  imp_wor <- (importance(output.forestG,type=2)) #Gini index - importance of traits to characterize the syndromes
  write.csv(imp_wor, paste("Imp_wor",btstrap,".csv",sep=""))
  impSyn <- (importance(output.forestG)) # trait importance for all syndromes seperate
  write.csv(impSyn, paste("ImpSyn", btstrap,".csv",sep=""))           
  errrate_wor[btstrap,1]<-median(output.forestG$err.rate[,1]) 
  errrate_wor[btstrap,2]<-median(output.forestG$err.rate[,2]) 
  errrate_wor[btstrap,3]<-median(output.forestG$err.rate[,3]) 
  errrate_wor[btstrap,4]<-median(output.forestG$err.rate[,4])
}

write.csv(pred, "bootPredAlltraits_GentoBuzz.csv") ### predicted syndromes for species with unknown pollinator
write.csv(errrate_wor, "Error_predAlltraits_GentoBuzz.csv") ### error rate for predicting known syndromes
write.csv(pred_wor, "KnownPred_Alltraits_GentoBuzz.csv") ### what was predicted for the known pollinators and how often within the 100 btstrap

#Error file

errrate_wor<-read.csv("~/Documents/Uni/PHD/Statistics/RandomForest/Error_predAlltraits_GentoBuzz.csv")

errorRate <- data.frame(matrix(nrow=1,ncol=4));colnames(errorRate)<-c("OOB","Buzz","Nectar foraging vertebrates","Passerine")

errorRate[1,1] <- median(errrate_wor$OOB) 
errorRate[1,2] <- median(errrate_wor$Buzz,na.rm=TRUE) # median error rate for predicting the bee syndrome right for species with known bee syndrome
errorRate[1,3] <- median(errrate_wor$Nectar.foraging.vertebrates,na.rm=TRUE) # median error rate for predicting the nectar foraging vertebrates syndrome right for species with known bee syndrome 
errorRate[1,4] <- median(errrate_wor$Passerine,na.rm=TRUE) # median error rate for predicting the Passerine syndrome right for species with known bee syndrome 

write.csv(errorRate, "~/Documents/Uni/PHD/Statistics/RandomForest/Output/Error_perSyn_Alltraits_GentoBuzz.csv")
 
# rename generalist to NF
library(dplyr)

RFData2 <- RFData1 %>%
  mutate(Syndrome = recode(Syndrome, "Generalist" = "Nectar foraging vertebrates"))

## seperate species with known pollinator from species with unknown pollinator

knownPoll <- subset(RFData2,Syndrome=="Buzz"|Syndrome=="Passerine"|Syndrome=="Nectar foraging vertebrates")
knownPoll
unknownPoll <- subset(RFData2,Syndrome!="Buzz"&Syndrome!="Passerine"&Syndrome!="Nectar foraging vertebrates")
unknownPoll

#find species with missing traits
missingKno <- as.data.frame(which(is.na(knownPoll[,c(4:43,45:47,52,53,55)]), arr.ind=TRUE)) # all traits except: 44, 48, 49, 50, 51, 54 
missingunk <- as.data.frame(which(is.na(unknownPoll[,c(4:43,45:47,52,53,55)]), arr.ind=TRUE))

#remove species with missing traits 
UnkCHA=unknownPoll[c(19,24,88,91,95),c(4:43,45:47,52:53,55)] 
KnoCHA=knownPoll[c(1:6,8:17,19:22,25:28,30:35,37:44,46:56,59:69,71:76,78:83,86,88,90,92:111,114,116:180,182:188,190:214,216,220:235,237,239:252),c(4:43,45:47,52:53,55)] 

knownPoll1 <- knownPoll[c(1:6,8:17,19:22,25:28,30:35,37:44,46:56,59:69,71:76,78:83,86,88,90,92:111,114,116:180,182:188,190:214,216,220:235,237,239:252),c(1,4:43,45:47,52:53,55)] 
unknownPoll1 <- unknownPoll[c(19,24,88,91,95),c(1,4:43,45:47,52:53,55)] 

spk<-knownPoll1$Taxon
spuk<-unknownPoll1$Taxon
str(knownPoll)

##### Select traits for the model and check for mising data

## check for NAs in DF
apply(is.na(UnkCHA), 2, which)
apply(is.na(KnoCHA), 2, which)

# safe set of selected traits and species with known pollinators
write.csv(KnoCHA, "~/Documents/Uni/PHD/Statistics/RandomForest//knownPoll_GentoNF.csv")
knownPoll1 <- read.csv("~/Documents/Uni/PHD/Statistics/RandomForest//knownPoll_GentoNF.csv", row.names=1) ## save file no traits with NA included

# check calsses
sapply(knownPoll1, class)

# transfomre Syndrome from character to factor
knownPoll1Trans <- transform(knownPoll1, Syndrome  = as.factor(as.character(Syndrome)))

####create empty datasets to store data. for all traits with NA removed, all species with known pollinators - knownPoll1##
pred <- data.frame(matrix(nrow=5,ncol=100));colnames(pred)<-btstrap #this is for predicting the unknown
imp_wor <- data.frame(matrix(nrow=225, ncol=2));colnames(imp_wor)<-c("char","Gini") #known character importance - GiniIndex
pred_wor <- data.frame(matrix(nrow=225,ncol=100));colnames(pred_wor)<-btstrap #known, what was predicted, was it wrong?
errrate_wor <-data.frame(matrix(nrow=100,ncol=4));colnames(errrate_wor)<-c("OOB","Buzz","Nectar foraging vertebrates","Passerine") #error rate of known

####make loop to find out most imp characters in differentiating the pollinators - knownPoll1##
for(btstrap in 1:100) 
{
  
  output.forestG <- randomForest(Syndrome~.,data=knownPoll1Trans,ntree=500,proximity=TRUE,importance=TRUE) # syndrome vs. traits
  prered1<-predict(output.forestG,UnkCHA) # use output from syndrome vs traits and predict unknown
  
  pred[,btstrap]<-prered1
  
  pred_wor[,btstrap]<-output.forestG$predicted
  imp_wor <- (importance(output.forestG,type=2)) #Gini index - importance of traits to characterize the syndromes
  write.csv(imp_wor, paste("Imp_wor",btstrap,".csv",sep=""))
  impSyn <- (importance(output.forestG)) # trait importance for all syndromes seperate
  write.csv(impSyn, paste("ImpSyn", btstrap,".csv",sep=""))           
  errrate_wor[btstrap,1]<-median(output.forestG$err.rate[,1]) 
  errrate_wor[btstrap,2]<-median(output.forestG$err.rate[,2]) 
  errrate_wor[btstrap,3]<-median(output.forestG$err.rate[,3]) 
  errrate_wor[btstrap,4]<-median(output.forestG$err.rate[,4])
}

write.csv(pred, "bootPredAlltraits_GentoNF.csv") ### predicted syndromes for species with unknown pollinator
write.csv(errrate_wor, "Error_predAlltraits_GentoNF.csv") ### error rate for predicting known syndromes
write.csv(pred_wor, "KnownPred_Alltraits_GentoNF.csv") ### what was predicted for the known pollinators and how often within the 100 btstrap

#Error file

errrate_wor<-read.csv("~/Documents/Uni/PHD/Statistics/RandomForest/Error_predAlltraits_GentoNF.csv")

errorRate <- data.frame(matrix(nrow=1,ncol=4));colnames(errorRate)<-c("OOB","Buzz","Nectar foraging vertebrates","Passerine")

errorRate[1,1] <- median(errrate_wor$OOB) 
errorRate[1,2] <- median(errrate_wor$Buzz,na.rm=TRUE) # median error rate for predicting the bee syndrome right for species with known bee syndrome
errorRate[1,3] <- median(errrate_wor$Nectar.foraging.vertebrates,na.rm=TRUE) # median error rate for predicting the nectar foraging vertebrates syndrome right for species with known bee syndrome 
errorRate[1,4] <- median(errrate_wor$Passerine,na.rm=TRUE) # median error rate for predicting the Passerine syndrome right for species with known bee syndrome 

write.csv(errorRate, "~/Documents/Uni/PHD/Statistics/RandomForest/Output/Error_perSyn_Alltraits_GentoNF.csv")

# rename generalist to Pass
library(dplyr)

RFData2 <- RFData1 %>%
  mutate(Syndrome = recode(Syndrome, "Generalist" = "Passerine"))

## seperate species with known pollinator from species with unknown pollinator
knownPoll <- subset(RFData2,Syndrome=="Buzz"|Syndrome=="Passerine"|Syndrome=="Nectar foraging vertebrates")
knownPoll
unknownPoll <- subset(RFData2,Syndrome!="Buzz"&Syndrome!="Passerine"&Syndrome!="Nectar foraging vertebrates")
unknownPoll

#find species with missing traits
missingKno <- as.data.frame(which(is.na(knownPoll[,c(4:43,45:47,52,53,55)]), arr.ind=TRUE)) # all traits except: 44, 48, 49, 50, 51, 54 
missingunk <- as.data.frame(which(is.na(unknownPoll[,c(4:43,45:47,52,53,55)]), arr.ind=TRUE))

#remove species with missing traits 
UnkCHA=unknownPoll[c(19,24,88,91,95),c(4:43,45:47,52:53,55)] 
KnoCHA=knownPoll[c(1:6,8:17,19:22,25:28,30:35,37:44,46:56,59:69,71:76,78:83,86,88,90,92:111,114,116:180,182:188,190:214,216,220:235,237,239:252),c(4:43,45:47,52:53,55)] 

knownPoll1 <- knownPoll[c(1:6,8:17,19:22,25:28,30:35,37:44,46:56,59:69,71:76,78:83,86,88,90,92:111,114,116:180,182:188,190:214,216,220:235,237,239:252),c(1,4:43,45:47,52:53,55)] 
unknownPoll1 <- unknownPoll[c(19,24,88,91,95),c(1,4:43,45:47,52:53,55)] 

spk<-knownPoll1$Taxon
spuk<-unknownPoll1$Taxon
str(knownPoll)

##### Select traits for the model and check for mising data

## check for NAs in DF
apply(is.na(UnkCHA), 2, which)
apply(is.na(KnoCHA), 2, which)

# safe set of selected traits and species with known pollinators
write.csv(KnoCHA, "~/Documents/Uni/PHD/Statistics/RandomForest//knownPoll_GentoPass.csv")
knownPoll1 <- read.csv("~/Documents/Uni/PHD/Statistics/RandomForest//knownPoll_GentoPass.csv", row.names=1) ## save file no traits with NA included

# check calsses
sapply(knownPoll1, class)

# transfomre Syndrome from character to factor
knownPoll1Trans <- transform(knownPoll1, Syndrome  = as.factor(as.character(Syndrome)))

####create empty datasets to store data. for all traits with NA removed, all species with known pollinators - knownPoll1##
pred <- data.frame(matrix(nrow=5,ncol=100));colnames(pred)<-btstrap #this is for predicting the unknown
imp_wor <- data.frame(matrix(nrow=225, ncol=2));colnames(imp_wor)<-c("char","Gini") #known character importance - GiniIndex
pred_wor <- data.frame(matrix(nrow=225,ncol=100));colnames(pred_wor)<-btstrap #known, what was predicted, was it wrong?
errrate_wor <-data.frame(matrix(nrow=100,ncol=4));colnames(errrate_wor)<-c("OOB","Buzz","Nectar foraging vertebrates","Passerine") #error rate of known

####make loop to find out most imp characters in differentiating the pollinators - knownPoll1##
for(btstrap in 1:100) 
{
  
  output.forestG <- randomForest(Syndrome~.,data=knownPoll1Trans,ntree=500,proximity=TRUE,importance=TRUE) # syndrome vs. traits
  prered1<-predict(output.forestG,UnkCHA) # use output from syndrome vs traits and predict unknown
  
  pred[,btstrap]<-prered1
  
  pred_wor[,btstrap]<-output.forestG$predicted
  imp_wor <- (importance(output.forestG,type=2)) #Gini index - importance of traits to characterize the syndromes
  write.csv(imp_wor, paste("Imp_wor",btstrap,".csv",sep=""))
  impSyn <- (importance(output.forestG)) # trait importance for all syndromes seperate
  write.csv(impSyn, paste("ImpSyn", btstrap,".csv",sep=""))           
  errrate_wor[btstrap,1]<-median(output.forestG$err.rate[,1]) 
  errrate_wor[btstrap,2]<-median(output.forestG$err.rate[,2]) 
  errrate_wor[btstrap,3]<-median(output.forestG$err.rate[,3]) 
  errrate_wor[btstrap,4]<-median(output.forestG$err.rate[,4])
}

write.csv(pred, "bootPredAlltraits_GentoPass.csv") ### predicted syndromes for species with unknown pollinator
write.csv(errrate_wor, "Error_predAlltraits_GentoPass.csv") ### error rate for predicting known syndromes
write.csv(pred_wor, "KnownPred_Alltraits_GentoPass.csv") ### what was predicted for the known pollinators and how often within the 100 btstrap

#Error file

errrate_wor<-read.csv("~/Documents/Uni/PHD/Statistics/RandomForest/Error_predAlltraits_GentoPass.csv")

errorRate <- data.frame(matrix(nrow=1,ncol=4));colnames(errorRate)<-c("OOB","Buzz","Nectar foraging vertebrates","Passerine")

errorRate[1,1] <- median(errrate_wor$OOB) 
errorRate[1,2] <- median(errrate_wor$Buzz,na.rm=TRUE) # median error rate for predicting the bee syndrome right for species with known bee syndrome
errorRate[1,3] <- median(errrate_wor$Nectar.foraging.vertebrates,na.rm=TRUE) # median error rate for predicting the nectar foraging vertebrates syndrome right for species with known bee syndrome 
errorRate[1,4] <- median(errrate_wor$Passerine,na.rm=TRUE) # median error rate for predicting the Passerine syndrome right for species with known bee syndrome 

write.csv(errorRate, "~/Documents/Uni/PHD/Statistics/RandomForest/Output/Error_perSyn_Alltraits_GentoPass.csv")

##############
# count missing

missing_counts <- colSums(is.na(RFData1))
empty_df <- data.frame(matrix(NA, nrow = 1, ncol = ncol(RFData1)))
empty_df[1, ] <- missing_counts

colnames(empty_df) <- colnames(RFData1) 

empty_df1 <- empty_df[c(1),c(4:43,45:47,52:53,55)] 

write.csv(empty_df, "~/Documents/Uni/PHD/Statistics/RandomForest/Output/SpecieswithMissingTraits.csv")

row_means <- rowMeans(empty_df[, 2:46])

#### run RF without most importatnt traits - due to lots of NA in reward and pollen exp.
#find missing traits
missingKno1 <- as.data.frame(which(is.na(knownPoll[,c(4,8:43,45:47,54,55)]), arr.ind=TRUE)) 
missingunk1 <- as.data.frame(which(is.na(unknownPoll[,c(4,8:43,45:47,54,55)]), arr.ind=TRUE))

a <- levels(factor(missingKno1$row))
b <- levels(factor(missingunk1$row))

missing_numbersKno <- setdiff(1:252, a)
missing_numbersUnk <- setdiff(1:159, b)

#remove species with missing traits 
UnkCHA=unknownPoll[missing_numbersUnk,c(4,8:43,45:47,54,55)] 
KnoCHA=knownPoll[missing_numbersKno,c(4,8:43,45:47,54,55)] 

## check for NAs in DF
apply(is.na(UnkCHA), 2, which)
apply(is.na(KnoCHA), 2, which)

#
knownPoll1 <- knownPoll[missing_numbersKno,c(1,4:43,45:47,52:53,55)] 
unknownPoll1 <- unknownPoll[missing_numbersUnk,c(1,4:43,45:47,52:53,55)] 

spk<-knownPoll1$Taxon
spuk<-unknownPoll1$Taxon
str(knownPoll)

write.csv(KnoCHA, "~/Documents/Uni/PHD/Statistics/RandomForest//knownPoll1.csv")
knownPoll2 <- read.csv("~/Documents/Uni/PHD/Statistics/RandomForest//knownPoll1.csv", row.names=1) ## save file no traits with NA included

knownPoll2Trans <- transform(knownPoll2, Syndrome  = as.factor(as.character(Syndrome))) # transfomre Syndrome from character to factor

####create empty datasets to store data. for all traits with NA removed, all species with known pollinators - knownPoll1##
pred <- data.frame(matrix(nrow=62,ncol=100));colnames(pred)<-btstrap #this is for prediciting the unknown, 3 weil 3 unknown 100 weil 100 bootstraps
imp_wor <- data.frame(matrix(nrow=228, ncol=2));colnames(imp_wor)<-c("char","Gini") #known character importance, 234 weil 234 known pollinators
pred_wor <- data.frame(matrix(nrow=228,ncol=100));colnames(pred_wor)<-btstrap #known, what was predicted, was it wrong?
errrate_wor <-data.frame(matrix(nrow=100,ncol=5));colnames(errrate_wor)<-c("OOB","Buzz","Generalist","Nectar foraging vertebrates","Passerine") #error rate of known

####make loop to find out most imp characters in differentiating the pollinators - knownPoll1##
for(btstrap in 1:100) 
{
  
  output.forestG <- randomForest(Syndrome~.,data=knownPoll2Trans,ntree=500,proximity=TRUE,importance=TRUE) # syndrome vs. traits
  prered1<-predict(output.forestG,UnkCHA) # use output from syndrome vs traits and predict unknown
  
  pred[,btstrap]<-prered1
  
  pred_wor[,btstrap]<-output.forestG$predicted
  imp_wor <- (importance(output.forestG,type=2)) #importance of Gini index - importance of traits to characterize the syndromes
  write.csv(imp_wor, paste("Imp_wor",btstrap,".csv",sep=""))
  impSyn <- (importance(output.forestG)) # trait importance for all syndromes seperate
  write.csv(impSyn, paste("ImpSyn", btstrap,".csv",sep=""))           
  errrate_wor[btstrap,1]<-median(output.forestG$err.rate[,1]) #saving out error rate for prediction of known species
  errrate_wor[btstrap,2]<-median(output.forestG$err.rate[,2]) #### WIESO 5 #### OOB + 4 Syndrome?
  errrate_wor[btstrap,3]<-median(output.forestG$err.rate[,3]) #################
  errrate_wor[btstrap,4]<-median(output.forestG$err.rate[,4])
  errrate_wor[btstrap,5]<-median(output.forestG$err.rate[,5])
}

write.csv(pred, "bootPredNoPolltraitsM2.csv") ### predicted syndromes for species with unknown pollinator
write.csv(errrate_wor, "Error_predNoPolltraitsM2.csv") ### error rate for predicting known syndromes
write.csv(pred_wor, "KnownPred_NoPolltraitsM2.csv") ### what was predicted for the known pollinators and how often within the 100 btstrap

######### SummarizePrediction - knownPoll1##

### summary for unknown syndrome predictions

#import prediction files for unknown pollinators and check how often one value is included
bootPred <- read.csv("~/Documents/Uni/PHD/Statistics/RandomForest/bootPredNoPolltraitsM2.csv", row.names=1) ### ordner ändern

code<-c("Buzz","Generalist","Nectar foraging vertebrates","Passerine")
pollpreds <- data.frame(matrix(nrow=62, ncol=4));colnames(pollpreds)<-c("Buzz","Generalist","Nectar foraging vertebrates","Passerine");rownames(pollpreds)<-spuk
zeile<-1:62

for(j in 1:length(zeile)){
  for(i in 1:length(code)){
    pollpreds[j,1]<-length(which(bootPred[j,]=="Buzz"))/100
    pollpreds[j,2]<-length(which(bootPred[j,]=="Generalist"))/100
    pollpreds[j,3]<-length(which(bootPred[j,]=="Nectar foraging vertebrates"))/100
    pollpreds[j,4]<-length(which(bootPred[j,]=="Passerine"))/100
    
  }}
pollpreds

write.csv(pollpreds, "Unk_PredNoPolltraitsM2.csv") ### wie oft wurden die unbekannten syndrome jeder art in den 100 btstaps predicted (1-immer der selbe, <1 nicht immer der selbe)

###################
####### Check which species are now new in KnoCha (KnoCha vs PollObs_pred_Alltraites_rev)
M1Pred <-read.csv("~/Documents/Uni/PHD/Statistics/RandomForest/Review/PollObs_pred_Alltraites_rev.csv")
M2Pred <- knownPoll1

NewSpeciesKnownPoll <- setdiff(M2Pred$Taxon,M1Pred$X) # get missing species
Syndromes <- M2Pred %>%
  filter(Taxon %in% c("Blastus cochinchinensis", "Macairea pachyphylla", "Pachycentria constricta"))

sx <- Syndromes$Syndrome
nam <- Syndromes$Taxon

## predict for new species in the df "Blastus cochinchinensis" "Macairea pachyphylla"    "Pachycentria constricta"
KnownPred <- read.csv("~/Documents/Uni/PHD/Statistics/RandomForest/Review/KnownPred_NoPolltraitsM2.csv", row.names=1)

# select new species
KnownPred <- cbind(M2Pred$Taxon,KnownPred)
selected_species <- KnownPred %>%
  filter(`M2Pred$Taxon` %in% c("Blastus cochinchinensis", "Macairea pachyphylla", "Pachycentria constricta"))

KnownPred <- selected_species[,2:101]

# predict species
pollpreds2 <- data.frame(matrix(nrow=3, ncol=2));colnames(pollpreds2)<-c("poll","p");rownames(pollpreds2)<-nam ## create empty data frame
pollpreds2[1]<-sx## in zeile eins von pollpreds2 wird polls geladen

pollpreds2[1,2]<-length(which(KnownPred[1,]=="Buzz"))/100
pollpreds2[2,2]<-length(which(KnownPred[2,]=="Buzz"))/100
pollpreds2[3,2]<-length(which(KnownPred[3,]=="Buzz"))/100

write.csv(pollpreds2, "PollObs_pred_NoPolltraitesM2.csv")

#file Error3
#sometimes no value, I guess if no bee or hb or pass was collected at that time!
errrate_wor<-read.csv("~/Documents/Uni/PHD/Statistics/RandomForest/Review/Error_predNoPolltraitsM2.csv")

errorRate1 <- data.frame(matrix(nrow=1,ncol=5));colnames(errorRate1)<-c("OOB","Buzz","Generalist","Nectar foraging vertebrates","Passerine")

errorRate1[1,1] <- median(errrate_wor$OOB) # 0.02
errorRate1[1,2] <- median(errrate_wor$Buzz,na.rm=TRUE) # median error rate for predicting the bee syndrome right for species with known bee syndrome = 0
errorRate1[1,3] <- median(errrate_wor$Generalist,na.rm=TRUE) # median error rate for predicting the generalist syndrome right for species with known bee syndrome = 0.3
errorRate1[1,4] <- median(errrate_wor$Nectar.foraging.vertebrates,na.rm=TRUE) # median error rate for predicting the nectar foraging vertebrates syndrome right for species with known bee syndrome = 0.04
errorRate1[1,5] <- median(errrate_wor$Passerine,na.rm=TRUE) # median error rate for predicting the Passerine syndrome right for species with known bee syndrome = 0

write.csv(errorRate1, "~/Documents/Uni/PHD/Statistics/RandomForest/Output/Error_perSyn_NoPolltraitsM2.csv")
####################

check<-read.csv("~/Documents/Uni/PHD/Statistics/RandomForest/Review/Output/Pred/Sum_PredObs_Rev_Fin.csv", sep = ";")
levels(factor(check$Spalte1))
duplicates <- check$Spalte1[duplicated(check$Spalte1)]

missing <- setdiff(RFData$Taxon,check$Spalte1)
