library(readr)
library(grid)
library(gridExtra)
library(ggpubr)
library(EnvStats)
library(boot)
library(FSA)

# Partition data
UBAT_FlashKinetics <- read_csv("UBAT_FlashKinetics.csv")
fusi<-subset(UBAT_FlashKinetics, type=="fusi")
nocti<-subset(UBAT_FlashKinetics, type=="nocti")
baha<-subset(UBAT_FlashKinetics, type=="baha")
poly<-subset(UBAT_FlashKinetics, type=="poly")
moni<-subset(UBAT_FlashKinetics, type=="moni")
leidyi<-subset(UBAT_FlashKinetics, type=="leidyi")
oiko<-subset(UBAT_FlashKinetics, type=="oiko")


#correlation of variables
cor(UBAT_FlashKinetics$`peak intensity`,UBAT_FlashKinetics$`integrated flash`,use="complete.obs")

# Kinetics of all
UBAT_FlashKinetics_new<-UBAT_FlashKinetics
UBAT_FlashKinetics_new$type<-factor(UBAT_FlashKinetics_new$type,c("baha","poly","moni","nocti","fusi","leidyi","oiko"))

boxplotFD<-ggplot(UBAT_FlashKinetics_new, aes(x=UBAT_FlashKinetics_new$type,y=UBAT_FlashKinetics$`flash duration`))+geom_boxplot(fill="9FF",alpha=0.2)+xlab("Species")+ylab("Flash duration (ms)")+ stat_summary(fun.y = mean,geom="point",shape=20,size=3,color="red",fill="red") 
boxplotFD


boxplotRT<-ggplot(UBAT_FlashKinetics_new, aes(x=UBAT_FlashKinetics_new$type,y=UBAT_FlashKinetics_new$`rise time`))+geom_boxplot(fill="9FF",alpha=0.2)+xlab("Species")+ylab("Rise time (ms)")+stat_summary(fun.y = mean,geom="point",shape=20,size=3,color="red",fill="red")+annotate("text",x=1:length(table(UBAT_FlashKinetics_new$type)),y=aggregate(UBAT_FlashKinetics_new$`flash duration`~UBAT_FlashKinetics_new$type, UBAT_FlashKinetics_new,median)[,2],label=table(UBAT_FlashKinetics_new$type),col="red")+ylim(0,500)
boxplotRT


boxplotDT<-ggplot(UBAT_FlashKinetics_new, aes(x=UBAT_FlashKinetics_new$type,y=UBAT_FlashKinetics_new$`decay time`))+geom_boxplot(fill="9FF",alpha=0.2)+xlab("Species")+ylab("Decay time (ms)")+ylim(0,1600)+stat_summary(fun.y = mean,geom="point",shape=20,size=3,color="red",fill="red") 
boxplotDT

boxplotEF<-ggplot(UBAT_FlashKinetics_new, aes(x=UBAT_FlashKinetics_new$type,y=UBAT_FlashKinetics_new$`efold`))+geom_boxplot(fill="9FF",alpha=0.2)+xlab("Species")+ylab("e-folding time (ms)")+ylim(0,500)+stat_summary(fun.y = mean,geom="point",shape=20,size=3,color="red",fill="red")
boxplotEF
library(scales)
library(ggbreak)
boxplotIE<-ggplot(UBAT_FlashKinetics_new, aes(x=UBAT_FlashKinetics_new$type,y=UBAT_FlashKinetics$`integrated flash`))+geom_boxplot(fill="9FF",alpha=0.2)+xlab("Species")+ylab("Integrated emission (photons)")+ stat_summary(fun.y = mean,geom="point",shape=20,size=3,color="red",fill="red")+scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x))
boxplotIE

boxplotIE<-ggplot(UBAT_FlashKinetics_new, aes(x=UBAT_FlashKinetics_new$type,y=UBAT_FlashKinetics$`integrated flash`))+geom_boxplot(fill="9FF",alpha=0.2)+xlab("Species")+ylab("Integrated emission (photons)")+ stat_summary(fun.y = mean,geom="point",shape=20,size=3,color="red",fill="red")+scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x))
boxplotIE
grid.arrange(boxplotRT,boxplotDT,boxplotEF,boxplotFD,boxplotPI,boxplotIE,ncol=2)


# P. fusiformis
fusiPI<-ggplot(fusi, aes(x=fusi$`peak intensity`))+geom_histogram(binwidth = 5000000000,color="grey",alpha=0.7)+labs(x="Peak intensity (photons/s)")+labs(y="Frequency")+theme_bw()+expand_limits(x=c(1e9,9e10))+scale_x_continuous(breaks=c(1e9,1e10,2e10,3e10,4e10,5e10,6e10,7e10,8e10,9e10))+scale_y_continuous(breaks = c(0,1,2,3,4,5,6,7,8,9,10,11,12))+ggtitle("P. fusiformis Peak intensity")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) 
fusiPI

fusiFD<-ggplot(fusi, aes(x=fusi$`flash duration`))+geom_histogram(binwidth = 50,color="grey",alpha=0.7)+labs(x="Flash duration (msec)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(200,300, 400,500, 600,700, 800,900, 1000,1100,1200,1300,1400,1500,1600,1700))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16))+ggtitle("P. fusiformis Flash duration")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))  
fusiFD

fusiRT<-ggplot(fusi, aes(x=fusi$`rise time`))+geom_histogram(binwidth = 20,color="grey",alpha=0.7)+labs(x="Rise time (msec)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(20,40,60,80,100,120,140,160,180,200))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18,20))+ggtitle("P. fusiformis Rise time") 
fusiRT

fusiIE<-ggplot(fusi, aes(x=fusi$`integrated flash`))+geom_histogram(binwidth = 2500000000,color="grey",alpha=0.7)+labs(x="Integrated emission (photons)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(1e9,5e9,1e10,1.5e10,2e10,2.5e10,3e10,3.5e10,4e10))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18,20))+ggtitle("P. fusiformis Integrated emission")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
fusiIE

fusiDT<-ggplot(fusi, aes(x=fusi$`decay time`))+geom_histogram(binwidth = 100,color="grey",alpha=0.7)+labs(x="Decay time (msec)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(200,300, 400,500, 600,700, 800,900, 1000,1100,1200,1300,1400,1500,1600,1700))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18,20))+ggtitle("P. fusiformis Decay time")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) 
fusiDT
grid.arrange(fusiPI,fusiFD,fusiRT,fusiDT,fusiIE,ncol=2)

# P. noctiluca

noctiPI<-ggplot(nocti, aes(x=nocti$`peak intensity`))+geom_histogram(binwidth = 2500000000,color="grey",alpha=0.7)+labs(x="Peak intensity (photons/s)")+labs(y="Frequency")+theme_bw()+expand_limits(x=c(1e9,3e10))+scale_x_continuous(breaks=c(1e9,5e9,1e10,1.5e10,2e10,2.5e10,3e10))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18,20))+ggtitle("P. noctiluca Peak intensity")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
noctiPI

noctiFD<-ggplot(nocti, aes(x=nocti$`flash duration`))+geom_histogram(binwidth = 50,color="grey",alpha=0.7)+labs(x="Flash duration (msec)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(200,300, 400,500, 600,700, 800,900, 1000,1100,1200,1300,1400,1500,1600,1700))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16))+ggtitle("P. noctiluca Flash duration") 
noctiFD

noctiRT<-ggplot(nocti, aes(x=nocti$`rise time`))+geom_histogram(binwidth = 10,color="grey",alpha=0.7)+labs(x="Rise time (msec)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(20,30,40,50,60,70,80,90,100,110,120,140,160,180,200))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18,20,22,24,26,28,30,32))+ggtitle("P. noctiluca Rise time") 
noctiRT

noctiIE<-ggplot(nocti, aes(x=nocti$`integrated flash`))+geom_histogram(binwidth = 1250000000,color="grey",alpha=0.7)+labs(x="Integrated emission (photons)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(1e8,2.5e9,5e9,7.5e9,1e10,1.25e10,1.5e10,1.75e10,2e10,2.5e10,3e10,3.5e10,4e10))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18,20))+ggtitle("P. noctiluca Integrated emission")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
noctiIE

noctiDT<-ggplot(nocti, aes(x=nocti$`decay time`))+geom_histogram(binwidth = 50,color="grey",alpha=0.7)+labs(x="Decay time (msec)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(200,300, 400,500, 600,700, 800,900, 1000,1100,1200,1300,1400,1500,1600,1700))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18,20))+ggtitle("P. noctiluca Decay time") 
noctiDT
grid.arrange(noctiPI,noctiFD,noctiRT,noctiDT,noctiIE,ncol=2)


# P. bahamense

bahaPI<-ggplot(baha, aes(x=baha$`peak intensity`))+geom_histogram(binwidth = 50000000,color="grey",alpha=0.7)+labs(x="Peak intensity (photons/s)")+labs(y="Frequency")+theme_bw()+expand_limits(x=c(1e8,8e8))+scale_x_continuous(breaks=c(1e8,2e8,3e8,4e8,5e8,6e8,7e8,8e8))+scale_y_continuous(breaks = c(0,1,2,3,4,5,6,7,8,9,10,11,12))+ggtitle("P. bahamense Peak intensity")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
bahaPI

bahaFD<-ggplot(baha, aes(x=baha$`flash duration`))+geom_histogram(binwidth = 25,color="grey",alpha=0.7)+labs(x="Flash duration (msec)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(0,50,100,150,200,250,300,350,400,450,500,550,600))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18,20,22,24,26))+ggtitle("P. bahamense Flash duration") 
bahaFD

bahaRT<-ggplot(baha, aes(x=baha$`rise time`))+geom_histogram(binwidth = 10,color="grey",alpha=0.7)+labs(x="Rise time (msec)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(20,30,40,50,60,70,80,90,100,110,120,140,160,180,200))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18,20,22,24,26,28,30,32))+ggtitle("P. bahamense Rise time") 
bahaRT

bahaIE<-ggplot(baha, aes(x=baha$`integrated flash`))+geom_histogram(binwidth = 50000000,color="grey",alpha=0.7)+labs(x="Integrated emission (photons)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(1e7,1e8,2e8,3e8,4e8,5e8,6e8))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18,20))+ggtitle("P. bahamense Integrated emission")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
bahaIE

bahaDT<-ggplot(baha, aes(x=baha$`decay time`))+geom_histogram(binwidth = 25,color="grey",alpha=0.7)+labs(x="Decay time (msec)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(0,25,50,75,100,125,150,175,200,225,250))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18,20,22,24,26,28,30))+ggtitle("P. bahamense Decay time") 
bahaDT
grid.arrange(bahaPI,bahaFD,bahaRT,bahaDT,bahaIE,ncol=2)

# L. polyedrum

polyPI<-ggplot(poly, aes(x=poly$`peak intensity`))+geom_histogram(binwidth = 25000000,color="grey",alpha=0.7)+labs(x="Peak intensity (photons/s)")+labs(y="Frequency")+theme_bw()+expand_limits(x=c(1e7,3e8))+scale_x_continuous(breaks=c(5e7,1e8,1.5e8,2e8,2.5e8))+scale_y_continuous(breaks = c(0,1,2,3,4,5,6,7,8,9,10,11,12))+ggtitle("L. polyedrum Peak intensity")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
polyPI

polyFD<-ggplot(poly, aes(x=poly$`flash duration`))+geom_histogram(binwidth = 25,color="grey",alpha=0.7)+labs(x="Flash duration (msec)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(0,50,100,150,200,250,300,350,400,450,500,550,600))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16))+ggtitle("L. polyedrum Flash duration") 
polyFD

polyRT<-ggplot(poly, aes(x=poly$`rise time`))+geom_histogram(binwidth = 10,color="grey",alpha=0.7)+labs(x="Rise time (msec)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(20,30,40,50,60,70,80,90,100,110,120,140,160,180,200))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18,20,22,24,26,28,30,32))+ggtitle("L. polyedrum Rise time") 
polyRT

polyIE<-ggplot(poly, aes(x=poly$`integrated flash`))+geom_histogram(binwidth = 10000000,color="grey",alpha=0.7)+labs(x="Integrated emission (photons)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(1e7,2e7,3e7,4e7,5e7,6e7,7e7,8e7,9e7,1e8,1.1e8))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18,20))+ggtitle("L. polyedrum Integrated emission")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
polyIE

polyDT<-ggplot(poly, aes(x=poly$`decay time`))+geom_histogram(binwidth = 25,color="grey",alpha=0.7)+labs(x="Decay time (msec)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(0,25,50,75,100,125,150,175,200,225,250))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18,20,22,24,26,28,30))+ggtitle("L. polyedrum Decay time") 
polyDT
grid.arrange(polyPI,polyFD,polyRT,polyDT,polyIE,ncol=2)

# A. monilatum
moniPI<-ggplot(moni, aes(x=moni$`peak intensity`))+geom_histogram(binwidth = 100000000,color="grey",alpha=0.7)+labs(x="Peak intensity (photons/s)")+labs(y="Frequency")+theme_bw()+expand_limits(x=c(1e8,1e9))+scale_x_continuous(breaks=c(1e8,3e8,5e8,7e8,9e8,1.1e9,1.3e9,1.5e9,1.7e9,1.9e9))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16))+ggtitle("A. monilatum Peak intensity")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
moniPI

moniFD<-ggplot(moni, aes(x=moni$`flash duration`))+geom_histogram(binwidth = 25,color="grey",alpha=0.7)+labs(x="Flash duration (msec)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(0,50,100,150,200,250,300,350,400,450,500,550,600))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16))+ggtitle("A. monilatum Flash duration") 
moniFD

moniRT<-ggplot(moni, aes(x=moni$`rise time`))+geom_histogram(binwidth = 10,color="grey",alpha=0.7)+labs(x="Rise time (msec)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(20,30,40,50,60,70,80,90,100,110,120,140,160,180,200))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18,20,22,24,26,28,30,32))+ggtitle("A monilatum Rise time") 
moniRT

moniIE<-ggplot(moni, aes(x=moni$`integrated flash`))+geom_histogram(binwidth = 50000000,color="grey",alpha=0.7)+labs(x="Integrated emission (photons)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(1e7,1e8,2e8,3e8,4e8,5e8,6e8,7e8,8e8,9e8))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18,20))+ggtitle("A monilatum Integrated emission")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
moniIE

moniDT<-ggplot(moni, aes(x=moni$`decay time`))+geom_histogram(binwidth = 25,color="grey",alpha=0.7)+labs(x="Decay time (msec)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(0,25,50,75,100,125,150,175,200,225,250,275,300,325))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18,20,22,24,26,28,30))+ggtitle("A monilatum Decay time") 
moniDT
grid.arrange(moniPI,moniFD,moniRT,moniDT,moniIE,ncol=2)


# M. leidyi
leidyiPI<-ggplot(leidyi, aes(x=leidyi$`peak intensity`))+geom_histogram(binwidth = 100000000000,color="grey",alpha=0.7)+labs(x="Peak intensity (photons/s)")+labs(y="Frequency")+theme_bw()+expand_limits(x=c(1e9,9e10))+scale_x_continuous(breaks=c(1e10,1e11,2e11,3e11,4e11,5e11,6e11,7e11,8e11,9e11,1e12))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18,20,22,24,26))+ggtitle("M. leidyi Peak intensity")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) 
leidyiPI

leidyiFD<-ggplot(leidyi, aes(x=leidyi$`flash duration`))+geom_histogram(binwidth = 200,color="grey",alpha=0.7)+labs(x="Flash duration (msec)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(200,400,600,800,1000,1200,1400,1600,1800,2000,2200,2400))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16))+ggtitle("M. leidyi Flash duration")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))  
leidyiFD

leidyiRT<-ggplot(leidyi, aes(x=leidyi$`rise time`))+geom_histogram(binwidth = 50,color="grey",alpha=0.7)+labs(x="Rise time (msec)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(100,150,200,250,300,350,400,450,500,600))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18,20))+ggtitle("M. leidyi Rise time") 
leidyiRT

leidyiIE<-ggplot(leidyi, aes(x=leidyi$`integrated flash`))+geom_histogram(binwidth = 50000000000,color="grey",alpha=0.7)+labs(x="Integrated emission (photons)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(1e10,1e11,2e11,3e11,4e11,5e11))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18,20,22,24,26))+ggtitle("M. leidyi Integrated emission")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
leidyiIE

leidyiDT<-ggplot(leidyi, aes(x=leidyi$`decay time`))+geom_histogram(binwidth = 200,color="grey",alpha=0.7)+labs(x="Decay time (msec)")+labs(y="Frequency")+theme_bw()+scale_x_continuous(breaks=c(200,400,600,800,1000,1200,1400,1600,1800,2000,2200))+scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18,20))+ggtitle("M. leidyi Decay time")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) 
leidyiDT
grid.arrange(leidyiPI,leidyiFD,leidyiRT,leidyiDT,leidyiIE,ncol=2)

# Peak intensity testing
shapiro.test(nocti$`peak intensity`)
library(FSA)
kruskal.test(UBAT_FlashKinetics$`peak intensity`~ UBAT_FlashKinetics$type, data=UBAT_FlashKinetics)
dunnTest(UBAT_FlashKinetics$`peak intensity`~UBAT_FlashKinetics$type,data = UBAT_FlashKinetics,method = "holm")

# Flash duration testing
kruskal.test(UBAT_FlashKinetics$`flash duration`~ UBAT_FlashKinetics$type, data=UBAT_FlashKinetics)
dunnTest(UBAT_FlashKinetics$`flash duration`~UBAT_FlashKinetics$type,data = UBAT_FlashKinetics,method = "holm")

# Decay time testing
kruskal.test(UBAT_FlashKinetics$`decay time`~ UBAT_FlashKinetics$type, data=UBAT_FlashKinetics)
dunnTest(UBAT_FlashKinetics$`decay time`~UBAT_FlashKinetics$type,data = UBAT_FlashKinetics,method = "holm")

# Rise time testing
kruskal.test(UBAT_FlashKinetics$`rise time`~ UBAT_FlashKinetics$type, data=UBAT_FlashKinetics)
dunnTest(UBAT_FlashKinetics$`rise time`~UBAT_FlashKinetics$type,data = UBAT_FlashKinetics,method = "holm")

# Integrated emission testing
kruskal.test(UBAT_FlashKinetics$`integrated flash`~ UBAT_FlashKinetics$type, data=UBAT_FlashKinetics)
dunnTest(UBAT_FlashKinetics$`integrated flash`~UBAT_FlashKinetics$type,data = UBAT_FlashKinetics,method = "holm")

#E-folding testing
kruskal.test(UBAT_FlashKinetics$`efold`~ UBAT_FlashKinetics$type, data=UBAT_FlashKinetics)
dunnTest(UBAT_FlashKinetics$`efold`~UBAT_FlashKinetics$type,data = UBAT_FlashKinetics,method = "holm")

# bootstrapping for CI
samplemean<-function(fusi, i) 
{    m <- mean(fusi$`peak intensity`[i],na.rm=TRUE)
n <- length(i)
v <- (n-1)*var(fusi$`peak intensity`[i])/n^2
c(m, v)}
fusibootPI<-boot(fusi,samplemean,R=5000)
fusiPICI<-boot.ci(fusibootPI,type=c("norm", "basic", "perc", "stud"))
fusiPICI

# Mean of FK
meanfusi<-c(mean(fusi$`peak intensity`,na.rm=TRUE),mean(fusi$`rise time`,na.rm=TRUE),mean(fusi$`decay time`,na.rm=TRUE),mean(fusi$`flash duration`,na.rm=TRUE),mean(fusi$`integrated flash`,na.rm=TRUE))
meanfusi

#  SD of FK
sd(leidyi$`peak intensity`,na.rm=TRUE)


#Correlation 
ggplot(UBAT_FlashKinetics, aes(x=UBAT_FlashKinetics$`peak intensity`, y=UBAT_FlashKinetics$`integrated flash`)) + geom_point()+geom_smooth(method="lm")+geom_point() +stat_regline_equation()+stat_cor(aes(label=..rr.label..))

gplot(UBAT_FlashKinetics, aes(x=log(UBAT_FlashKinetics$`flash duration`), y=log(UBAT_FlashKinetics$`integrated flash`))) + geom_point()+geom_smooth(method="lm") + geom_point() +stat_cor(aes(label=..rr.label..))+stat_regline_equation(label.x=4,label.y=25)

ggplot(UBAT_FlashKinetics, aes(x=UBAT_FlashKinetics$`rise time`, y=log(UBAT_FlashKinetics$`peak intensity`))) + geom_point()+geom_smooth(method="lm") + geom_point() +stat_cor(aes(label=..rr.label..))+stat_regline_equation(label.x=4,label.y=25)

#Pearson test
pearson<-cor.test(UBAT_FlashKinetics$`peak intensity`, UBAT_FlashKinetics$`integrated flash`,method = "pearson")
pearson

#scatter plots
library(ggplot2)
ggplot(UBAT_FlashKinetics, aes(UBAT_FlashKinetics$`rise time`,UBAT_FlashKinetics$`efold`,col=UBAT_FlashKinetics$type))+geom_point()

#Cluster analysis
library(cluster)
  #for all observations
z <- UBAT_FlashKinetics[,-c(1,1)]
means <- apply(z,2,mean)
sds <- apply(z,2,sd)
nor <- scale(z,center=means,scale=sds)
distance=dist(nor)
cluster=hclust(distance,method='average')
plot(cluster)
plot(cluster, labels=UBAT_FlashKinetics$type,main = 'Cluster')
rect.hclust(cluster,k=7)

  #using sample means
UBAT_means <- read_csv("UBAT_means.csv")
z_means <- UBAT_means[,-c(1,1)]
means_means <- apply(z_means,2,mean)
sds_means <- apply(z_means,2,sd)
nor_means <- scale(z_means,center=means_means,scale=sds_means)
distance=dist(nor_means)
cluster=hclust(distance,method='average')
plot(cluster)
plot(cluster, labels=UBAT_means$type,main = 'Cluster')


#Linear Discriminant Analysis and multivariate model
library(ggplot2)
library(MASS)
UBAT_LDA<-read_csv("UBAT_LDA.csv")
UBAT_LDA[2:7]<-scale(UBAT_LDA[2:7])
UBAT_LDA<-UBAT_LDA[,-5]
UBAT_LDA<-UBAT_LDA[,-5]
apply(UBAT_LDA[2:5],2,mean,na.rm=TRUE)
apply(UBAT_LDA[2:5],2,sd,na.rm=TRUE)
set.seed(5)
sample<-sample(c(TRUE,FALSE),nrow(UBAT_LDA),replace=TRUE,prob = c(0.7,0.3))
train<-UBAT_LDA[sample,]
test<-UBAT_LDA[!sample,]
model<-lda(type~.,data=train)
model
predicted<-predict(model,test)
head(predicted$class)
head(predicted$posterior)
mean(predicted$class==test$type,na.rm=TRUE)

lda_plot<-cbind(train,predict(model)$x)
ggplot(lda_plot,aes(LD1,LD2))+geom_point(aes(color=type))

lda.train<-predict(model,train)
train$lda<-lda.train$class
table(train$lda,train$type)

lda.test<-predict(model,test)
test$lda<-lda.test$class
table(test$lda,test$type)

#factor analysis
library(psych)
fa <- fa(r = UBAT_LDA[2:7], nfactors = 6, rotate = "varimax")
summary(fa)
fa$loadings

#without baha
UBAT_LDA<-read_csv("UBAT_LDA_nobaha.csv")
UBAT_LDA<-UBAT_LDA[,-5]
UBAT_LDA<-UBAT_LDA[,-5]
UBAT_LDA_scale<-scale(UBAT_LDA[2:5])
apply(UBAT_LDA[2:5],2,mean,na.rm=TRUE)
apply(UBAT_LDA[2:5],2,sd,na.rm=TRUE)
set.seed(1)
sample<-sample(c(TRUE,FALSE),nrow(UBAT_LDA),replace=TRUE,prob = c(0.3,0.7))
train<-UBAT_LDA[sample,]
test<-UBAT_LDA[!sample,]
model<-lda(type~.,data=train)
model
predicted<-predict(model,test)
head(predicted$class)
head(predicted$posterior)
mean(predicted$class==test$type,na.rm=TRUE)

lda_plot<-cbind(train,predict(model)$x)
ggplot(lda_plot,aes(LD1,LD2))+geom_point(aes(color=type))