OliviasHuntAssayAnalysis.Rmd

---
title: "OliviaFeedingAssay"
author: "Konstantinos Lagogiannis"
date: "23/05/2022"
output:
  pdf_document: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
source("config_lib.R")

setEnvFileLocations("LAPTOP") #HOME,OFFICE,#LAPTOP

source("HuntingEventAnalysis_lib.r")
source("TrajectoryAnalysis.r")

library(RColorBrewer);
library(ggplot2)
library(ggpubr) ##install.packages("ggpubr")

mypalette <- brewer.pal(8,"Set1")
mypal_rgb <- col2rgb(mypalette,alpha = 1)
mypal_rgb["alpha",1:8] <-55 ##Opacity

colourH <- list(HET= mypalette[2], #rgb(0.9,0.01,0.01,0.2),
                HOM = mypalette[3], #rgb(0.01,0.01,0.9,0.2))
                WT=mypalette[1],
                ROT="black") #rgb(0.01,0.7,0.01,0.2),

colourHLine <- list(HET= col2hex(mypal_rgb)[2], #rgb(0.9,0.01,0.01,0.2),
                HOM = col2hex(mypal_rgb)[3], #rgb(0.01,0.01,0.9,0.2))
                WT=col2hex(mypal_rgb)[1],
                ROT="black") #rgb(0.01,0.7,0.01,0.2),
 ##Transparent For MCMC Samples (Live) #c(rgb(0.95,0.01,0.01,0.1),rgb(0.01,0.7,0.01,0.1),rgb(0.01,0.01,0.9,0.1),rgb(0.00,0.00,0.0,1.0)) ####Transparent For MCMC Samples (Empty)

pchH <- list(HET= 1,
             HOM = 16, #rgb(0.01,0.01,0.9,0.2))
             WT=22) #rgb(0.01,0.7,0.01,0.2),

##Load Hunt event data
strHuntEventFile <- paste0(strDataStore,"/HB_allHuntEvents.rds")
datHuntEvents <-readRDS(strHuntEventFile)

##load Tracked Frames
strdatTrackedFramesFile = paste0(strDataStore,"/setn1_Dataset_DATA081022.rds")
datAllFrames <- readRDS(strdatTrackedFramesFile)

vGroupIDs <- c("WT","HET","HOM") #(unique(datAllFrames$group)) ##Rearing Group names(groupsrcdatList)
vCondIDs <- as.character(unique(datAllFrames$testCond))

```

# Hunt Events

```{r, re-calcHuntStats, echo = FALSE}

##Calc Hunt Stat
lHuntStat     <- list();
i = 0
for (g in vGroupIDs)
{
  for (c in vCondIDs)
  {
    i = i + 1
    #message(paste("#### Process Rearing Group ",g," Tested in",c, "###############"))
    lHuntStat[[g]] <- calcHuntStat3(datHuntEvents[datHuntEvents$groupID == g,])
    #lHuntStat[[i]] = calcHuntStat3(datHuntEvent)
    stopifnot(length(lHuntStat[[g]]$vHLarvaEventCount) > 0)
  }
}
datHuntStat = do.call(rbind,lHuntStat)#

#message("~~ Finished Summarizing HuntEvent Data For Each Group ~~ ")

```

##  Number of Hunt Events

Box plots comparing number of hunt events per larva across groups. Initial a basic boxplot without the test and p-values, followed by a plot with Wilcox non parametric test and a boxplot with pairwise parametric tests.

 The Frames Are in datAllFrames , the huntEvents are in datAllHuntEvents
## Gaussian Kernel Empirical distribution 

```{r, hunteventsperGroup, echo = FALSE}

  datEventsPerExpID <- getHuntEventsCountsPerExp(datHuntEvents)
  densHuntEventCounts <- list(HOM = density(datHuntStat["HOM",]$vHLarvaEventCount),
                        HET = density(datHuntStat["HET",]$vHLarvaEventCount),
                        WT  =  density(datHuntStat["WT",]$vHLarvaEventCount)               
                        )
  BW = "20"
#pdf(paste0(strPlotExportPath,"/huntcounts_densities_empirical.pdf"),width=7,height=5,title="Gaussian Kernel Estimated Distribution of hunt event counts per group") 
  
    plot(densHuntEventCounts$HOM,col=unlist(colourH["HOM"]), lwd=3,lty=lineTypeL[1],main="Hunt-event counts per group ",xlim=c(0,350),ylim=c(0,0.01))
    lines(densHuntEventCounts$HET,col=(colourH["HET"][[1]] ), lwd=3,lty=lineTypeL[2])
    lines(densHuntEventCounts$WT,col=unlist(colourH["WT"]), lwd=3,lty=lineTypeL[3])
    legend("topright",legend=paste(vGroupIDs),lwd=3,col=c(colourH[["HOM"]], colourH[["HET"]], colourH[["WT" ]]),lty=lineTypeL )
    mtext(side = 1,cex=cex, line = lineAxis, expression(paste("Hunt Events per larva in 20min") )  ) #(",lambda," )

#dev.off()
    # t.test(datEventsPerExpID[datEventsPerExpID$groupID == "HET",]$nHuntEvents,datEventsPerExpID[datEventsPerExpID$groupID == "WT",]$nHuntEvents)
    
    #message("HOM Event Count Stats")
    #print(densHuntEventCounts$HOM)
  
    #message("WT Event Count Stats")
    #print(densHuntEventCounts$WT)


```
  
### Statistical Model comparison

```{r huntrate statistical model inference-setup}
library(rjags)


## Run Baysian Inference on Model for Hunt Event Counts IN A Group/ Test Condition
## Return Samples Drawn structure
mcmc_drawEventCountModels <- function(nGroupEvents,strModelFilename)
{
  varnames1=c("n","q","r")
  burn_in=1000;
  steps=15000;
  plotsamples = 1000
  thin=2;
  chains = 3
  
  
  ##Select the event Count Column 2
  nEvents=as.numeric(nGroupEvents)
  nDat = length(nEvents);
  dataG=list(n=nEvents,NTOT=nDat);
  
  model=jags.model(file=strModelFilename,data=dataG,n.chains=chains);
  
  update(model,burn_in)
  
  drawSamples=jags.samples(model,steps,thin=thin,variable.names=varnames1)
  
  return(drawSamples) 
}


## Compare Model TO Data Using CDF ##
plotEventCountDistribution_cdf <- function(datHEventCount,drawHEvent,lcolour="black",lpch=16,lty=1,Plim,nplotSamples=100,newPlot = FALSE)
{
  XLim <- 400
  x <- seq(0,XLim,1)

  cdfD_N <- ecdf(datHEventCount$nHuntEvents)

  plot(cdfD_N,col=lcolour,pch=lpch,xlab=NA,ylab=NA,main="",xlim=c(0,XLim),ylim=c(0,1),cex=cex+0.1,cex.axis=cex,cex.lab=cex,add=!newPlot)
  ##Construct CDF of Model by Sampling randomly from Model distribution for exp rate parameter
  for (c in 1:NROW(drawHEvent$q[1,1,])) {
    for (j in (NROW(drawHEvent$q[,,c])-nplotSamples):NROW(drawHEvent$q[,,c]) )
    {
      cdfM <- dnbinom(x,size=drawHEvent$r[,j,c],prob=  drawHEvent$q[,j,c]  )##1-exp(-q*x) ##ecdf(  dexp( x, q  ) )
      lines(x,cumsum(cdfM),col=lcolour,lty=lty) #add=TRUE,
    }
  }
  plot(cdfD_N,col=colourP[4],pch=lpch,xlab=NA,ylab=NA,main="",xlim=c(0,XLim),ylim=c(0,1),cex=cex+0.1,cex.axis=cex,cex.lab=cex,add=TRUE)
  #axis(side = 4)
  #mtext(side = 4, line = 2.1, 'Counts')
  ### Draw Distribution oF Hunt Rates - 
  ## For Poisson- gamma mixture we recover Gamma Params from nbinomial as shape r>shape r and scale theta:(1-p)/p
  
  ##For the EXP Mixture with Poisson The Exp Rate was recovered as :
  ##(z= p/(1-p))
#  hist( (1-tail(drawHEvent$q[,,c],nplotSamples))/tail(drawHEvent$q[,,c],nplotSamples)  )
  
  
}

## Discrete - Geometric Cause Mixture of rates - assuming rates drawn from most informative Prior distribution (EXP)
## Assuming nbinom(r,p) Poisson(L|a,b) Gamma(a,b) then r=a, p=1/(b+1) -> b=(1-p)/p
## Give geometric
modelGEventRateGeom="model { 
q ~ dunif(0.0,1)
r ~ dgamma(1,1)

for(j in 1:NTOT){
  n[j] ~  dnegbin(q,r) ##Model Number Of Hunt Events Per Larvae
  }
}"



strModelName = "modelGroupEventRate.tmp"
fileConn=file(strModelName)
writeLines(modelGEventRateGeom,fileConn);
close(fileConn)

### Cut And Examine The data Where There Are Between L and M rotifers Initially
preyCntRange <- c(0,1100)
plotsamples <- 20
drawSamples_HOM <- mcmc_drawEventCountModels(datHuntStat["HOM",]$vHLarvaEventCount,"modelGroupEventRate.tmp")
drawSamples_HET <- mcmc_drawEventCountModels(datHuntStat["HET",]$vHLarvaEventCount,"modelGroupEventRate.tmp")
drawSamples_WT <- mcmc_drawEventCountModels(datHuntStat["WT",]$vHLarvaEventCount,"modelGroupEventRate.tmp")
```

Plot CDS of empirical and Statistical model to validate fit  

```{r huntrate data vs stat model CDF }

#pdf(paste0(strPlotExportPath,"/huntcounts_modelVsData_CDF_HOM.pdf"),width=7,height=5,title="CDF Model hunt event counts per group") 

  plotEventCountDistribution_cdf(datEventsPerExpID[datEventsPerExpID$groupID == "HOM",],
                                 drawSamples_HOM,colourHLine[["HOM" ]],pchL[1],lineTypeL[1],Plim,plotsamples,newPlot=TRUE )
   mtext(side = 1,cex=cex, line = lineAxis, expression(paste("Number of hunt events in 20min") )  ) #(",lambda," )
  title("HOM Huntrate CDF model")  
#dev.off()

#pdf(paste0(strPlotExportPath,"/huntcounts_modelVsData_CDF_HET.pdf"),width=7,height=5,title="CDF Model hunt event counts per group") 

  plotEventCountDistribution_cdf(datEventsPerExpID[datEventsPerExpID$groupID == "HET",],drawSamples_HET, colourHLine[["HET"]],pchL[2],lineTypeL[2],Plim,plotsamples,newPlot=TRUE )
 mtext(side = 1,cex=cex, line = lineAxis, expression(paste("Number of hunt events in 20min") )  ) #(",lambda," )

   title("HET Huntrate CDF model")  
#dev.off()
 
 
#pdf(paste0(strPlotExportPath,"/huntcounts_modelVsData_CDF_WT.pdf"),width=7,height=5,title="CDF Model hunt event counts per group") 

  plotEventCountDistribution_cdf(datEventsPerExpID[datEventsPerExpID$groupID == "WT",],drawSamples_WT, colourHLine[["WT"]],pchL[3],lineTypeL[3],Plim,plotsamples,newPlot=TRUE )
   mtext(side = 1,cex=cex, line = lineAxis, expression(paste("Number of hunt events in 20min") )  ) #(",lambda," )

   title("WT Huntrate CDF model")  
 
#dev.off()


```
  
After checking model fit with CDF, then compare (Poisson) inferred huntrates per group using model 

```{r huntrate model compare}


### Draw Distribution oF Hunt Rates - 
## for the exp draw (z= p/(1-p)) ## But it is the same for Rate Of Gamma Too / Or inverse for scale
plotsamples <- 5000
schain <-1:3

### The Prob Of Success p from NegBinom translates to Gamma Rate p/(1-p), or scale: (1-p)/p
HEventHuntGammaRate_HOM <-tail(drawSamples_HOM$q[,,schain],plotsamples)/(1-tail(drawSamples_HOM$q[,,schain],plotsamples));
HEventHuntGammaRate_HET <-tail(drawSamples_HET$q[,,schain],plotsamples)/(1-tail(drawSamples_HET$q[,,schain],plotsamples));
HEventHuntGammaRate_WT <- tail(drawSamples_WT$q[,,schain],plotsamples)/(1-tail(drawSamples_WT$q[,,schain],plotsamples));
HEventHuntGammaShape_HOM <- tail(drawSamples_HOM$r[,,schain],plotsamples);
HEventHuntGammaShape_HET <- tail(drawSamples_HET$r[,,schain],plotsamples)
HEventHuntGammaShape_WT <- tail(drawSamples_WT$r[,,schain],plotsamples);

## Calc and Save POisson Rate Samples
lPoissonHuntRates <- list(HOM=HEventHuntGammaShape_HOM*1/HEventHuntGammaRate_HOM,
                          HET=HEventHuntGammaShape_HET*1/HEventHuntGammaRate_HET,
                          WT=HEventHuntGammaShape_WT*1/HEventHuntGammaRate_WT)
## Calculate Hunt Rates ###
pBW = 2
densHPoissonRate_HOM <- density( lPoissonHuntRates$HOM,bw=pBW)
densHPoissonRate_HET <- density( lPoissonHuntRates$HET,bw=pBW)
densHPoissonRate_WT <- density( lPoissonHuntRates$WT,bw=pBW)

##  Estimate  Hunt Rates Distributions ###
densHPoissonRate <- list(HOM= densHPoissonRate_HOM,
                         HET= densHPoissonRate_HET,
                         WT= densHPoissonRate_WT
                         )


Ylim <- 0.04
#pdf(paste0(strPlotExportPath,"/huntPoissonRate_modelComparison.pdf"),width=7,height=5,title=" Model Inferred hunt rates per group") 

    plot(densHPoissonRate$HOM$x, densHPoissonRate$HOM$y,type='l',lty=lineTypeL[1],col=colourH["HOM"][[1]] ,lwd=4,ylab=NA,xlab=NA,xlim=c(0,230),ylim=c(0,Ylim),cex=cex,cex.axis=cex )
  lines(densHPoissonRate$HET$x, densHPoissonRate$HET$y,type='l',lty=lineTypeL[2],col=colourH["HET"][[1]],lwd=4,ylab=NA,xlab=NA)
  lines(densHPoissonRate$WT$x, densHPoissonRate$WT$y,type='l',lty=lineTypeL[3],col=colourH["WT"][[1]],lwd=4,ylab=NA,xlab=NA)
  
  legend("topright",legend=paste(vGroupIDs),lwd=3,
         col=c(colourH[["HOM"]], colourH[["HET" ]], colourH[["WT" ]] ),
         lty=lineTypeL)

    
  mtext(side = 1,cex=cex, line = lineAxis, expression(paste("Inferred hunt-rate from statistical model  (N/20 min.)") )  ) #(",lambda," )
  mtext(side = 2,cex=cex, line = lineAxis, " Density function ")
 title("Distribution of model-inferred hunt rates ")

 message("test HOM vs WT PoissonRate model param")
 print(t.test(lPoissonHuntRates$HOM,lPoissonHuntRates$WT))
#dev.off()
 
 message("test HOM vs WT PoissonRate model param")
 print(t.test(lPoissonHuntRates$HET,lPoissonHuntRates$WT))
#dev.off()
```
 It is clear that HETs have higher huntrates to Control (WT) and HOM. From this analysis we may further infer probabilities of a HET larvae having higher huntrates 
 
 
The empirical mean number vs the model based mean of of hunt-events estimated per group is :

```{r mean-huntevent_summary,echo=FALSE}
lmuModelRates <- list(HET=mean(lPoissonHuntRates$HET),
                      HOM=mean(lPoissonHuntRates$HOM),
                      WT=mean(lPoissonHuntRates$WT)
                      )
print(tapply(datHuntEventCountsPerGroup$LarvalEventCount,datHuntEventCountsPerGroup$Group,mean))
print(data.frame(lmuModelRates) )
message("with SD :")

print(tapply(datHuntEventCountsPerGroup$LarvalEventCount,datHuntEventCountsPerGroup$Group,sd))
```
 
```{r, hunteventsperGroup_huntStat, echo = FALSE}

#print(table(datAllHuntEvent$groupID))


boxplot(as.numeric(datHuntStat[vGroupIDs[1],]$vHLarvaEventCount),
        as.numeric(datHuntStat[vGroupIDs[2],]$vHLarvaEventCount),
        as.numeric(datHuntStat[vGroupIDs[3],]$vHLarvaEventCount),
        names = c(vGroupIDs[1],vGroupIDs[2],vGroupIDs[3]),
        col=c(colourH[[vGroupIDs[1] ]], colourH[[vGroupIDs[2] ]], colourH[[vGroupIDs[3] ]] ),
        ylab="Number of Hunt events per larva ")

testHEvents =  list(
            WTvsHET = t.test(as.numeric(datHuntStat["WT",]$vHLarvaEventCount),
            as.numeric(datHuntStat["HET",]$vHLarvaEventCount) ),
            WTvsHOM = t.test(as.numeric(datHuntStat["WT",]$vHLarvaEventCount),
            as.numeric(datHuntStat["HOM",]$vHLarvaEventCount) ),
            HOMvsHET = t.test(as.numeric(datHuntStat["HOM",]$vHLarvaEventCount),
            as.numeric(datHuntStat["HET",]$vHLarvaEventCount) )
          )


message("Significance WT vs HET: ", testHEvents$WTvsHET$p.value,
      " mean :", as.numeric(testHEvents$WTvsHET$estimate[1]), 
      " vs ", as.numeric(testHEvents$WTvsHET$estimate[2]) )

message("Significance WT vs HOM: ", testHEvents$WTvsHOM$p.value, " mean :", 
      testHEvents$WTvsHOM$estimate[1]," vs ", testHEvents$WTvsHOM$estimate[2] )
message( "Significance HOM vs HET: ", testHEvents$HOMvsHET$p.value, 
      " mean :", testHEvents$HOMvsHET$estimate[1]," vs ",testHEvents$HOMvsHET$estimate[2])

## GGPlot Way ##
#Prep dat
datHuntEventCountsPerGroup <- rbind.data.frame(cbind(as.numeric(datHuntStat["WT",]$vHLarvaEventCount),"WT"),
      cbind(as.numeric(datHuntStat["HET",]$vHLarvaEventCount),"HET"),
      cbind(as.numeric(datHuntStat["HOM",]$vHLarvaEventCount),"HOM"),stringsAsFactors=FALSE
      )
datHuntEventCountsPerGroup$V1 <- as.numeric(datHuntEventCountsPerGroup$V1)
names(datHuntEventCountsPerGroup)=c("LarvalEventCount","Group")

p <- ggboxplot(datHuntEventCountsPerGroup, x = "Group", y = "LarvalEventCount",
          color = "Group", palette = "jco",
          add = "jitter",ylab="Number of hunt events per larva")

comp_pairs <- list( c("WT", "HET"), c("WT", "HOM"), c("HOM", "HET") )
#  Add p-value
p + stat_compare_means(
                       method = "wilcox.test",comparisons = comp_pairs)+ggtitle("Non-parametric comparison on number of hunt events")
# Change method
p + stat_compare_means(method = "t.test",comparisons = comp_pairs, label="p.format",label.y=c(15,20,30)) + ggtitle("Parametric comparison on number of hunt events")

#p + stat_compare_means( aes(label = paste0(..method.., "\n", "p =", ..p.format..)), 
#                        label.x = 1.5, label.y = 40,comparisons = comp_pairs) 
```





## Timing of Hunt events

How are hunt events distributed during the 20min of prey encounters?
We pool hunt event occurrences withing each group and plot the cumulative distribution function across the 20min

```{r hunteventsTiming, echo = FALSE}

plot(ecdf(datHuntEvents[datHuntEvents$groupID == "HET",]$startFrame/(50*60)),col=colourH$HET,pch=pchH$HET,
     main="Cumulative distribution of hunt events through time",xlab="Time (min)")
lines(ecdf(datHuntEvents[datHuntEvents$groupID == "HOM",]$startFrame/(50*60)),col=colourH$HOM,pch=pchH$HOM)
lines(ecdf(datHuntEvents[datHuntEvents$groupID == "WT",]$startFrame/(50*60)),col=colourH$WT,pch=pchH$WT,)
legend("topleft",legend=names(colourH), col=unlist(colourH),pch=unlist(pchH) )
```

##   Hunt-event Duration 

### Total time spent hunting 

```{r hunteventsDuration,echo = FALSE}

t.test(datHuntStat["WT",]$vHDurationPerLarva,datHuntStat["HET",]$vHDurationPerLarva)

c_FrmToMin <- 1/(60*mean(datAllFrames$fps))
boxplot(datHuntStat[vGroupIDs[1],]$vHDurationPerLarva*c_FrmToMin,
        datHuntStat[vGroupIDs[2],]$vHDurationPerLarva*c_FrmToMin,
        datHuntStat[vGroupIDs[3],]$vHDurationPerLarva*c_FrmToMin,
        names = c(vGroupIDs[1],vGroupIDs[2],vGroupIDs[3]),
        col=c(colourH[[vGroupIDs[1]]],colourH[[vGroupIDs[2] ]],colourH[[vGroupIDs[3] ]] ),
        ylab="Total time spent hunting per larva (min)")


testHEvents =  list(
            WTvsHET = t.test(as.numeric(datHuntStat["WT",]$vHDurationPerLarva*c_FrmToMin),
            as.numeric(datHuntStat["HET",]$vHDurationPerLarva*c_FrmToMin) ),
            WTvsHOM = t.test(as.numeric(datHuntStat["WT",]$vHDurationPerLarva*c_FrmToMin),
            as.numeric(datHuntStat["HOM",]$vHDurationPerLarva*c_FrmToMin) ),
            HOMvsHET = t.test(as.numeric(datHuntStat["HOM",]$vHDurationPerLarva*c_FrmToMin),
            as.numeric(datHuntStat["HET",]$vHDurationPerLarva*c_FrmToMin) )
          )


message("Significance WT vs HET: ", testHEvents$WTvsHET$p.value, " mean :",testHEvents$WTvsHET$estimate[1]," vs ",testHEvents$WTvsHET$estimate[2] )
message("Significance WT vs HOM: ", testHEvents$WTvsHOM$p.value, " mean:",testHEvents$WTvsHOM$estimate[1]," vs ",testHEvents$WTvsHOM$estimate[2] )
message("Significance HOM vs HET: ", testHEvents$HOMvsHET$p.value, " mean:",testHEvents$HOMvsHET$estimate[1]," vs ",testHEvents$HOMvsHET$estimate[2] )


## GGPlot Way ##
#Prep dat
datHuntEventDurationPerGroup <- rbind.data.frame(cbind(as.numeric(datHuntStat["WT",]$vHDurationPerLarva),"WT"),
      cbind(as.numeric(datHuntStat["HET",]$vHDurationPerLarva),"HET"),
      cbind(as.numeric(datHuntStat["HOM",]$vHDurationPerLarva),"HOM"),stringsAsFactors=FALSE
      )
datHuntEventDurationPerGroup$V1 <- as.numeric(datHuntEventDurationPerGroup$V1)/(50) #Convert to seconds 50 fps
names(datHuntEventDurationPerGroup)=c("HuntDurationPerLarva","Group")

p <- ggboxplot(datHuntEventDurationPerGroup, x = "Group", y = "HuntDurationPerLarva",
          color = "Group", palette = "jco",
          add = "jitter",ylab="time spent hunting (sec)")

comp_pairs <- list( c("WT", "HET"), c("WT", "HOM"), c("HOM", "HET") )
#  Add p-value
p + stat_compare_means(
                       method = "wilcox.test",comparisons = comp_pairs) + ggtitle("Non-parametric comparison of total hunt duration")
# Change method
p + stat_compare_means(method = "t.test",comparisons = comp_pairs, label="p.format") + ggtitle("Parametric comparison of hunt duration")

#p + stat_compare_means( aes(label = paste0(..method.., "\n", "p =", ..p.format..)), 
#                        label.x = 1.5, label.y = 40,comparisons = comp_pairs) 

```


The mean hunt duration (sec) estimated for each group is :
```{r mean-huntduration,echo=FALSE}

print(tapply(datHuntEventDurationPerGroup$HuntDurationPerLarva,datHuntEventDurationPerGroup$Group,mean))
message("with SD:")
print(tapply(datHuntEventDurationPerGroup$HuntDurationPerLarva,datHuntEventDurationPerGroup$Group,sd))

```


## Mean hunt-event episode duration

The duration of each hunt event may be different between genotypes. In previous work I have found that the duration of a hunt episode appears to be intrinsically controlled - like an internal timer- clock that is not strongly influenced by external stimuli.
We compare the mean duration of hunt-events between groups, by comparing the mean hunt-episode duration between larva of each group

```{r durationperhuntevent,echo = FALSE}

datHuntStat["HET",]$meanEpisodeDuration
datHuntStat["HOM",]$meanEpisodeDuration
datHuntStat["WT",]$meanEpisodeDuration

datHuntEvents$duration_sec <- (datHuntEvents$endFrame - datHuntEvents$startFrame)/mean(datAllFrames$fps)

vHuntEventDurationPerLarva <-tapply(datHuntEvents$duration_sec,datHuntEvents$expID,mean)
vGroupForLarva <- tapply(as.character(datHuntEvents$groupID),datHuntEvents$expID,head,1)
datEventDurationPerLarva = cbind.data.frame(meanHuntEpisodeDuration_sec=as.numeric(vHuntEventDurationPerLarva),
                                            Group=as.character(vGroupForLarva),stringsAsFactors=FALSE )

## GGPLOT Episode Duration Stats ##

p <- ggboxplot(datEventDurationPerLarva, x = "Group", y = "meanHuntEpisodeDuration_sec",
          color = "Group", palette = "jco",
          add = "jitter",ylab="mean hunt-even duration per larva (sec)")

comp_pairs <- list( c("WT", "HET"), c("WT", "HOM"), c("HOM", "HET") )
#  Add p-value
p + stat_compare_means(
                       method = "wilcox.test",comparisons = comp_pairs) + ggtitle("Non-parametric comparison of total hunt-event duration")
# Change method
p + stat_compare_means(method = "t.test",comparisons = comp_pairs, label="p.format") + ggtitle("Parametric comparison of hunt-event duration")


```
The mean hunt-episode duration (sec) is estimated for each group as :
```{r mean-huntEventduration,echo=FALSE}

print(tapply(datEventDurationPerLarva$meanHuntEpisodeDuration_sec,datEventDurationPerLarva$Group,mean))
message("with SD :")
print(tapply(datEventDurationPerLarva$meanHuntEpisodeDuration_sec,datEventDurationPerLarva$Group,sd))

```

## Prey Density Change

Plot prey Trajectories per group

```{r preyDynamics, echo = FALSE, eval=TRUE}

  ## Plot Prey Dynamics ##
  mat_y <- list()
  PreyCountSampleFrames <- list()
  muPreyCount <- list()
  sdPreyCount <- list()
  errPreyCount <- list()
  for (grpID in unique(datAllFrames$groupID) )
  {
    datGroupFrames <- datAllFrames[datAllFrames$groupID== grpID,]
    bFirst <- TRUE
    nFrames <- max(datAllFrames$frameN)
    mat_y[[grpID]] <- matrix(NA,nrow=NROW(unique(datGroupFrames$expID)),
                               ncol=nFrames)
    i = 1
    ## This Is a very Slow Way to ensure sampling of Prey Counts across exp on same frameN 
    for (fID in unique(datGroupFrames$expID) )
    {
      ## Cautionary way because dublicate frameN records may exist
      x <- datGroupFrames[datGroupFrames$expID == fID & datGroupFrames$frameN <= nFrames,"frameN"]
      y <- datGroupFrames[datGroupFrames$expID == fID &
                                         datGroupFrames$frameN  <= nFrames ,]$PreyCount
      ##fill specific frameIdx Colums for Vector - in case Frames are missing
      mat_y[[grpID]][i,x] <- y
      # 
      # if (bFirst){
      #   
      #   bFirst    <- FALSE
      #   plot(x,y, type="l",ylim=c(0,100), #xlim=c(1,20)
      #        main=grpID,col=colourH[[grpID]],ylab="Prey count",xlab="Time (min)" )  
      # }else{
      #   
      #   lines(x,y,col=colourH[[grpID]])
      # }
    i = i + 1
    } ##For Each ExpID
    
    
    muPreyCount[[grpID]] = apply(mat_y[[grpID]][],2,mean,na.rm=T)
    PreyCountSampleFrames[[grpID]] <- which(!is.na(muPreyCount[[grpID]]))##Because Rot tracks have missing frames /recorded at intervals
    muPreyCount[[grpID]] <- meanf(muPreyCount[[grpID]][PreyCountSampleFrames[[grpID]]],3500)
    sdPreyCount[[grpID]] = meanf(apply(mat_y[[grpID]],2,sd,na.rm=T),3500)
    errPreyCount[[grpID]] = sdPreyCount[[grpID]]/sqrt(NROW(mat_y[[grpID]]) )

    
  }##For Each GRoup

x <- PreyCountSampleFrames$WT/(50*60)
plot(x,muPreyCount$WT,col=colourH$WT,type="l",xlim=c(1,18),xlab="Time (min)",ylab="Mean prey count",lwd=1,ylim=c(40,70))
x <- PreyCountSampleFrames$HET/(50*60)
lines(x,muPreyCount$HET,col=colourH$HET,type="l",lwd=1)
x <- PreyCountSampleFrames$HOM/(50*60)
lines(x,muPreyCount$HOM-7,col=colourH$HOM,type="l",lwd=3)
x <- PreyCountSampleFrames$ROT/(50*60)
lines(x,muPreyCount$ROT-40,col=colourH$ROT,type="l",lwd=3)
legend("topright",legend=c("WT","HET","HOM","ROT"), col=c(colourH$WT,colourH$HET,colourH$HOM,colourH$ROT),lty=1,lwd=2 )

  
  
  
x <- PreyCountSampleFrames$WT/(50*60)
plot(x,muPreyCount$WT,col=colourH$WT,type="l",xlim=c(1,18),xlab="Time (min)",ylab="Mean prey count",lwd=1,ylim=c(40,70))
#lines(x,muPreyCount$WT+errPreyCount$WT,col=colourH$WT,lwd=2)
polygon(c(rev(x),x), c(rev(muPreyCount$WT-errPreyCount$WT), (muPreyCount$WT+errPreyCount$WT) ),
        density = 40,angle = -45, col=colourH$WT,lwd=1)


x <- PreyCountSampleFrames$HET/(50*60)
lines(x,muPreyCount$HET,col=colourH$HET,type="l",lwd=1)
polygon(c(rev(x),x), c(rev(muPreyCount$HET-errPreyCount$HET), (muPreyCount$HET+errPreyCount$HET) ),
        density = 30,angle = 15, col=colourH$HET,lwd=1)


x <- PreyCountSampleFrames$HOM/(50*60)
lines(x,muPreyCount$HOM,col=colourH$HOM,type="l",lwd=3)
polygon(c(rev(x),x), c(rev(muPreyCount$HOM-errPreyCount$HOM), (muPreyCount$HOM+errPreyCount$HOM) ),
        density = 22,angle = 45, col=colourH$HOM,lwd=1)



##Add The Mean Rotifer Trend
x <- PreyCountSampleFrames$ROT/(50*60) #which(!is.na(muPreyCount$ROT))##Because Rot tracks have missing frames /recorded at intervals
lines(x,muPreyCount$ROT,col=colourH$ROT,type="l",lwd=3,cex=0.5)
#polygon(c(rev(x),x), c(rev(muPreyCount$ROT[x]-errPreyCount$ROT[x]), (muPreyCount$ROT[x]+errPreyCount$ROT[x]) ),
        #density = 22,angle = 45, col=colourH$ROT,lwd=1)
##Linear Regression 
lm.ROT <- lm(x ~ muPreyCount$ROT)
abline(lm.ROT,lwd=3,lty=2,col=colourH$ROT)

legend("topright",legend=c("WT","HET","HOM","ROT"), col=c(colourH$WT,colourH$HET,colourH$HOM,colourH$ROT),lty=1,lwd=2 )

## Prey Reduction 
muPreyReduction<- list(
          HET=mean(muPreyCount$HET[120:2000])-mean(muPreyCount$HET[(nFrames-6000):(nFrames-5000)]),
          WT=mean(muPreyCount$WT[120:2000])- mean(muPreyCount$WT[(nFrames-6000):(nFrames-5000)]),
          HOM=mean(muPreyCount$HOM[120:2000])- mean(muPreyCount$HOM[(nFrames-6000):(nFrames-5000)]),
          ROT=mean(muPreyCount$ROT[120:2000],na.rm=T)-mean(muPreyCount$ROT[(nFrames-6000):(nFrames-5000)],na.rm=T)
          )

message("Mean prey reduction for WT:",format(muPreyReduction$WT,digits =3),
        " HET:",format(muPreyReduction$HET,digits =3)," HOM:",format(muPreyReduction$HOM,digits =3))

##Prep Prey dataInto Single Data Frame
lfit <- list()
ldatFrm <- list()
for (grpID in names(muPreyCount))
{
 ##nEed to concatenate all timeseries into one long vector, and make a corresponding time vector 
  x <- c(matrix(matrix( rep(seq(1,nFrames),nrow(mat_y[[grpID]]))/(50*60),nrow=nrow(mat_y[[grpID]]),byrow=F),nrow=1,byrow=T))
  y <-   unlist(as.data.frame(t(mat_y[[grpID]]) ))
  
  ldatFrm[[grpID]] <- cbind.data.frame(x,y,group=grpID)
  # ldatRot[[grpID]] <- cbind.data.frame(frameN = x,
  #                                      PreyCount = y)  #matrix(mat_y[[grpID]],nrow=1,byrow=T)
  ## Do Fit On All Timeseries
  lfit[[grpID]] <- lm(y ~ x)
}
datPreyReduction <-  do.call(rbind,ldatFrm)




## GGPlot Way ##
#Prep dat of Reduction Per experiment

##Special Treamtmet for ROT as data does not come from Hunt Event analysis
## Calc Prey Reduction -
lvInitPreyPerExp  <- list()
lvFinalPreyPerExp <- list()
lvNaturalPreyReduction <- list()
for (grpID in names(muPreyCount) )
{
  lvInitPreyPerExp[[grpID]] <- apply(mat_y[[grpID]] [,1:2000],1,mean,na.rm=T)
  lvFinalPreyPerExp[[grpID]]  <- apply(mat_y[[grpID]] [,(nFrames-6000):(nFrames-5000)],1,mean,na.rm=T)
  lvNaturalPreyReduction[[grpID]]  <- lvInitPreyPerExp[[grpID]] - lvFinalPreyPerExp[[grpID]] 
  lvNaturalPreyReduction[[grpID]]  <- lvNaturalPreyReduction[[grpID]][!is.nan(lvNaturalPreyReduction[[grpID]] )]
}


datPreyReductionPerGroup <- rbind.data.frame(
      cbind(as.numeric(lvNaturalPreyReduction[["WT"]]),"WT"),
      cbind(as.numeric(lvNaturalPreyReduction[["HET"]]),"HET"),
      cbind(as.numeric(lvNaturalPreyReduction[["HOM"]]),"HOM"),
      cbind(as.numeric(lvNaturalPreyReduction[["ROT"]]),"ROT"),
      stringsAsFactors=FALSE
      )

datPreyReductionPerGroup$V1 <- as.numeric(datPreyReductionPerGroup$V1) #Convert to seconds 50 fps
names(datPreyReductionPerGroup)=c("PreyReductionPerLarva","Group")

p <- ggboxplot(datPreyReductionPerGroup, x = "Group", y = "PreyReductionPerLarva",
          color = "Group", palette = "jco",
          add = "jitter",ylab="Number of prey ")

comp_pairs <- list( c("WT", "HET"), c("WT", "HOM"), c("HOM", "HET"),c("WT", "ROT") )
#  Add p-value
p + stat_compare_means(
                       method = "wilcox.test",comparisons = comp_pairs) + ggtitle("Non-parametric comparison of prey reduction")
# Change method
p + stat_compare_means(method = "t.test",comparisons = comp_pairs, label="p.format") + ggtitle("Parametric comparison of prey reduction")


```



The mean prey reduction per larva of each group is :
```{r mean-preyreduction, echo=FALSE}

print(tapply(datPreyReductionPerGroup$PreyReductionPerLarva, datPreyReductionPerGroup$Group,mean))
message("with SD:")
print(tapply(datPreyReductionPerGroup$PreyReductionPerLarva, datPreyReductionPerGroup$Group,sd))
```

### Compare consumption using linear regression 

```{r prey-linear-regression}
lfit_group <- lm(x~y*group,data=datPreyReduction)
anova(lfit_group)

# Compare slopes #
library(lsmeans)
lfit_group$coefficients
lfit_group.lst <- lstrends(lfit_group,"group",var="y")
pairs(lfit_group.lst)

##Show Data With Linear Fit ##Fix this To Overall
for (grpID in names(muPreyCount))
{
 
  message("Linear fit summary for group ",grpID)
  print(summary(lfit[[grpID]]))
    
  plot(lfit[[grpID]]$model$x,lfit[[grpID]]$model$y, type="p",ylim=c(0,100), #xlim=c(1,20)
             main=grpID,col=colourH[[grpID]],ylab="Prey count",xlab="Time (min)",cex=0.05 )  
  
  abline(lfit[[grpID]],lwd=4,col="black",lty=2)
  legend("topright",
         legend= paste(grpID,"r2=",format(summary(lfit[[grpID]])$r.squared,digits=4), "\n c=",format(lfit[[grpID]]$coefficients[1],digits=4),", w=",format(lfit[[grpID]]$coefficients[2],digits=4)))
}


```



## Eye Motion

We summarize the left-right eye motion from each group plot 2D densities. One method produces a 2D distribution by counting the number of larvae whose eyes have occupied a particular point in L-R eye angle phase. The second method uses kernel density estimation methods to estimate the distirbution of eye movements.
These plots usually serve to validate that eye tracking works as expected. A common resting point for the eyes is around -10 and 10 ie Vergence of 20 degrees, while Vergence above 45 degrees indicates larva has initiated hunting.


```{r sampleHuntEventDetect-validation,echo=FALSE}
 ### Make Eye Phase Space Density Plots ##########
  s_expID <- 10
  datSampleFish <- datAllFrames[datAllFrames$expID == s_expID,]
  datSampleFish_EventFrames <- datHuntEvents[datHuntEvents$expID ==s_expID, ]
  plot(datSampleFish$frameN/(datSampleFish$fps*60), #
    (datSampleFish$LEyeAngle-datSampleFish$REyeAngle),type="l",#xlim=c(0.1,10),
    xlab="Time (min)",ylab="Eye vergence angle (deg)",
    main=paste(" fID:",s_expID, "hunt events:",NROW(datHuntEvents[datHuntEvents$expID == s_expID,] )),
    ylim=c(-10,60))
  abline(h=G_THRESHUNTVERGENCEANGLE,lty=2)
  points(datSampleFish_EventFrames$startFrame/(mean(datSampleFish$fps)*60), datSampleFish_EventFrames$eyeVergence ,col="red" )

```

The above figure shows a sample where hunt events where detected for fish ID `r print(s_expID)` based on eye-vergence data.



```{r eyemovements,echo = FALSE, eval=TRUE,cache=TRUE}
  
#  datSampleFish[datSampleFish$frameN %in% datSampleFish_EventFrames$startFrame,"LEyeAngle"]
   
  for (i in vGroupIDs)
  {
    message(paste("#### Eye ProcessGroup ",i," ###############"))
    
    
    ##Take All larva IDs recorded - Regardless of Data Produced - No Tracks Is Also Data
    #vexpID = unique(filtereddatAllFrames$expID)
    ##Select Larvaof this Group
    
    datAllGroupFrames <- datAllFrames[which(datAllFrames$groupID == i),]
    #Note:A Larva ID Corresponds to A specific Condition ex. NF1E (Same Fish Is tested in 2 conditions tho ex. NF1E, NF1L)
    vexpID = unique(datAllGroupFrames$expID)
    #plotGroupMotion(datAllGroupFrames,lHuntStat[[i]],vexpID)
    #######################################################################
    ###  EYE - PLOT Scatter and Eye Densities #####
    strCond = i;
    source("EyeScatterAndDensities.r")
    #####
  }
```

The distribution of WT eye-motion shows less skew than HET and HOM, which suggests that WT eyes are more free to move independently of each other.