---
title: "Darting_Analysis"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## R Markdown

It is a good idea to run through this workbook one chunk at a time to ensure that things are making sense as it goes along. Technically, on day 1, the only things you need to add before being able to run all chunks is the name of the file at the beginning and the name of your output in the last line of the last chunk.

First, name and save this file in the working directory your data is in.

```{r}
getwd() #should return the working directory where your data and R Markdown file are
funcSE<-function(x) sqrt(var(x,na.rm=TRUE)/length(na.omit(x))) #useful function to get SEM
## Gives count, mean, standard deviation, standard error of the mean, and confidence interval (default 95%).
##   data: a data frame.
##   measurevar: the name of a column that contains the variable to be summariezed
##   groupvars: a vector containing names of columns that contain grouping variables
##   na.rm: a boolean that indicates whether to ignore NA's
##   conf.interval: the percent range of the confidence interval (default is 95%)
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {
    library(plyr)

    # New version of length which can handle NA's: if na.rm==T, don't count them
    length2 <- function (x, na.rm=FALSE) {
        if (na.rm) sum(!is.na(x))
        else       length(x)
    }

    # This does the summary. For each group's data frame, return a vector with
    # N, mean, and sd
    datac <- ddply(data, groupvars, .drop=.drop,
      .fun = function(xx, col) {
        c(N    = length2(xx[[col]], na.rm=na.rm),
          mean = mean   (xx[[col]], na.rm=na.rm),
          sd   = sd     (xx[[col]], na.rm=na.rm)
        )
      },
      measurevar
    )

    # Rename the "mean" column    
    datac <- plyr::rename(datac, c("mean" = measurevar))

    datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean

    # Confidence interval multiplier for standard error
    # Calculate t-statistic for confidence interval: 
    # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
    ciMult <- qt(conf.interval/2 + .5, datac$N-1)
    datac$ci <- datac$se * ciMult

    return(datac)
}
## Norms the data within specified groups in a data frame; it normalizes each
## subject (identified by idvar) so that they have the same mean, within each group
## specified by betweenvars.
##   data: a data frame.
##   idvar: the name of a column that identifies each subject (or matched subjects)
##   measurevar: the name of a column that contains the variable to be summariezed
##   betweenvars: a vector containing names of columns that are between-subjects variables
##   na.rm: a boolean that indicates whether to ignore NA's
normDataWithin <- function(data=NULL, idvar, measurevar, betweenvars=NULL,
                           na.rm=FALSE, .drop=TRUE) {
    library(plyr)

    # Measure var on left, idvar + between vars on right of formula.
    data.subjMean <- ddply(data, c(idvar, betweenvars), .drop=.drop,
     .fun = function(xx, col, na.rm) {
        c(subjMean = mean(xx[,col], na.rm=na.rm))
      },
      measurevar,
      na.rm
    )

    # Put the subject means with original data
    data <- merge(data, data.subjMean)

    # Get the normalized data in a new column
    measureNormedVar <- paste(measurevar, "_norm", sep="")
    data[,measureNormedVar] <- data[,measurevar] - data[,"subjMean"] +
                               mean(data[,measurevar], na.rm=na.rm)

    # Remove this subject mean column
    data$subjMean <- NULL

    return(data)
}
## Summarizes data, handling within-subjects variables by removing inter-subject variability.
## It will still work if there are no within-S variables.
## Gives count, un-normed mean, normed mean (with same between-group mean),
##   standard deviation, standard error of the mean, and confidence interval.
## If there are within-subject variables, calculate adjusted values using method from Morey (2008).
##   data: a data frame.
##   measurevar: the name of a column that contains the variable to be summariezed
##   betweenvars: a vector containing names of columns that are between-subjects variables
##   withinvars: a vector containing names of columns that are within-subjects variables
##   idvar: the name of a column that identifies each subject (or matched subjects)
##   na.rm: a boolean that indicates whether to ignore NA's
##   conf.interval: the percent range of the confidence interval (default is 95%)
summarySEwithin <- function(data=NULL, measurevar, betweenvars=NULL, withinvars=NULL,
                            idvar=NULL, na.rm=FALSE, conf.interval=.95, .drop=TRUE) {

  # Ensure that the betweenvars and withinvars are factors
  factorvars <- vapply(data[, c(betweenvars, withinvars), drop=FALSE],
    FUN=is.factor, FUN.VALUE=logical(1))

  if (!all(factorvars)) {
    nonfactorvars <- names(factorvars)[!factorvars]
    message("Automatically converting the following non-factors to factors: ",
            paste(nonfactorvars, collapse = ", "))
    data[nonfactorvars] <- lapply(data[nonfactorvars], factor)
  }

  # Get the means from the un-normed data
  datac <- summarySE(data, measurevar, groupvars=c(betweenvars, withinvars),
                     na.rm=na.rm, conf.interval=conf.interval, .drop=.drop)

  # Drop all the unused columns (these will be calculated with normed data)
  datac$sd <- NULL
  datac$se <- NULL
  datac$ci <- NULL

  # Norm each subject's data
  ndata <- normDataWithin(data, idvar, measurevar, betweenvars, na.rm, .drop=.drop)

  # This is the name of the new column
  measurevar_n <- paste(measurevar, "_norm", sep="")

  # Collapse the normed data - now we can treat between and within vars the same
  ndatac <- summarySE(ndata, measurevar_n, groupvars=c(betweenvars, withinvars),
                      na.rm=na.rm, conf.interval=conf.interval, .drop=.drop)

  # Apply correction from Morey (2008) to the standard error and confidence interval
  #  Get the product of the number of conditions of within-S variables
  nWithinGroups    <- prod(vapply(ndatac[,withinvars, drop=FALSE], FUN=nlevels,
                           FUN.VALUE=numeric(1)))
  correctionFactor <- sqrt( nWithinGroups / (nWithinGroups-1) )

  # Apply the correction factor
  ndatac$sd <- ndatac$sd * correctionFactor
  ndatac$se <- ndatac$se * correctionFactor
  ndatac$ci <- ndatac$ci * correctionFactor

  # Combine the un-normed means with the normed results
  merge(datac, ndatac)
}
##below are useful packages for plotting and analysis
library(ggplot2)
library(plyr)
library(dplyr)
library(tidyr)
library(data.table)
library(tidyverse)
library(ggpubr)
library(rstatix)
library(WRS2)
pd<-position_dodge(0.3) #common variable for plotting
```


Load in your data from csv's (can save an Excel file as a CSV)
[For now] with: 
titles in row 1
Trial in columnn 1(A)
Animal (Mouse) in column 2(B)
Group in column 3(C)
Sex in column 4(D)
TimeStamp in column 5(E)
StimPeriod in column 6(F)   (Every single row of data will need a defined stimulus period)
Distance_cm in column 7(G)
Samples in column 8(H)
Velocity in colunn 9(I)

Missing data is OK at this point (but you want >99% samples if possible)

**The above can be adjusted, you just may have to adjust the code in some places accordingly.
The titles below are examples. "mydata" data will be used for all following examples.
```{r}
#mydata_Acq1<-read.csv("Statistics-Acq1_all trials.csv")
#mydata_Acq2<-read.csv("Statistics-Acq2_all trials.csv")
#mydata_Ext<-read.csv("Statistics-Ext1_all trials.csv")
#mydata_Corr<-read.csv("CorrelationsExpt1ReplicationGroup.csv")

mydata<-read.csv("Statistics-07-18_Experiment1_Day1Acq.csv") # Include the name of your desired data here. Above are examples of more specific names one can use as well, and would be necessary if you wanted to analyze multiple data files in one R Markdown file. The title of your CSV should be directly copied into this line of code, with a .csv at the end. 

load("mouse_overallThresh.Rda") #load in the overall Threshold data for the relevant species here. It will load in as variable -- threshDF
#load("rat_overallThresh.Rda") #load in the overall Threshold data for the relevant species here. It will load in as variable -- threshDF
```

Run this chunk to get a continuous variable for time @ sampling rate (3.75 Hz)
Will need to do for each session of different lengths.
This means that all video files need to be exactly the same length. This may require trimming in Excel after Ethovision.
```{r}
nObs<-as.numeric(data.table(mydata)[,.N,by=Animal][1,2])
samplen<-length(mydata$Trial)/nObs #This is a roundabout way to determine sample size for any given data
nObsString<-rep(c(1:nObs),samplen)
nObs
head(nObsString);tail(nObsString)
mydata<-cbind(mydata,nObsString)
samplen
```

Identify factors and numeric/dbl values.
```{r}
mydata$Distance_cm<-as.numeric(as.character(mydata$Distance_cm))
mydata$Velocity<-as.numeric(as.character(mydata$Velocity))
mydata$Animal<-as.factor(mydata$Animal)
mydata[is.na(mydata)]<-0
```

View your data
```{r}
head(mydata)
```

Here's a semi-useful easy to make visualization of EVERYTHING
```{r}
ggwholesession<-ggplot(mydata, aes(nObsString,Velocity,color=StimPeriod)) + geom_line()+facet_wrap(Animal~Group,scales='free',labeller = labeller(.cols = label_parsed, .multi_line = FALSE))
ggwholesession
```

This is a semi-tedious way to get markers for the pre and post timestamps for each trial. If there are fewer trials, like in acquisition, remove teh unnecessary markers by putting a comment '#' before the line.
```{r}
PRE1<-mydata$nObsString[min(which(mydata$StimPeriod == "CS1a"))-38]
POST1<-mydata$nObsString[min(which(mydata$StimPeriod == "CS1b"))+76]
PRE2<-mydata$nObsString[min(which(mydata$StimPeriod == "CS2a"))-38]
POST2<-mydata$nObsString[min(which(mydata$StimPeriod == "CS2b"))+76]
PRE3<-mydata$nObsString[min(which(mydata$StimPeriod == "CS3a"))-38]
POST3<-mydata$nObsString[min(which(mydata$StimPeriod == "CS3b"))+76]
PRE4<-mydata$nObsString[min(which(mydata$StimPeriod == "CS4a"))-38]
POST4<-mydata$nObsString[min(which(mydata$StimPeriod == "CS4b"))+76]
PRE5<-mydata$nObsString[min(which(mydata$StimPeriod == "CS5a"))-38]
POST5<-mydata$nObsString[min(which(mydata$StimPeriod == "CS5b"))+76]
PRE6<-mydata$nObsString[min(which(mydata$StimPeriod == "CS6a"))-38]
POST6<-mydata$nObsString[min(which(mydata$StimPeriod == "CS6b"))+76]
PRE7<-mydata$nObsString[min(which(mydata$StimPeriod == "CS7a"))-38]
POST7<-mydata$nObsString[min(which(mydata$StimPeriod == "CS7b"))+76]
PRE8<-mydata$nObsString[min(which(mydata$StimPeriod == "CS8a"))-38]
POST8<-mydata$nObsString[min(which(mydata$StimPeriod == "CS8b"))+76]
PRE9<-mydata$nObsString[min(which(mydata$StimPeriod == "CS9a"))-38]
POST9<-mydata$nObsString[min(which(mydata$StimPeriod == "CS9b"))+76]
PRE10<-mydata$nObsString[min(which(mydata$StimPeriod == "CS10a"))-38]
POST10<-mydata$nObsString[min(which(mydata$StimPeriod == "CS10b"))+76]
PRE11<-mydata$nObsString[min(which(mydata$StimPeriod == "CS11a"))-38]
POST11<-mydata$nObsString[min(which(mydata$StimPeriod == "CS11b"))+76]
PRE12<-mydata$nObsString[min(which(mydata$StimPeriod == "CS12a"))-38]
POST12<-mydata$nObsString[min(which(mydata$StimPeriod == "CS12b"))+76]
PRE13<-mydata$nObsString[min(which(mydata$StimPeriod == "CS13a"))-38]
POST13<-mydata$nObsString[min(which(mydata$StimPeriod == "CS13b"))+76]
PRE14<-mydata$nObsString[min(which(mydata$StimPeriod == "CS14a"))-38]
POST14<-mydata$nObsString[min(which(mydata$StimPeriod == "CS14b"))+76]
PRE15<-mydata$nObsString[min(which(mydata$StimPeriod == "CS15a"))-38]
POST15<-mydata$nObsString[min(which(mydata$StimPeriod == "CS15b"))+76]
PRE16<-mydata$nObsString[min(which(mydata$StimPeriod == "CS16a"))-38]
POST16<-mydata$nObsString[min(which(mydata$StimPeriod == "CS16b"))+76]
PRE1;POST1
PRE2;POST2
PRE3;POST3
PRE4;POST4
PRE5;POST5
PRE6;POST6
PRE7;POST7
PRE8;POST8
PRE9;POST9
PRE10;POST10
PRE11;POST11
PRE12;POST12
PRE13;POST13
PRE14;POST14
PRE15;POST15
PRE16;POST16
```


These are the 'basics' of what we'd like to visualize, with MANY options to customize.
You need to adjust the Trial numbers and PRE and POST numbers for each trial (already done here). Of course, if there are more trials (like in extinciton, you'll have to add those graphs. Copying and pasting and adjusting trial numbers will do the trick; but you will have to have the PRE and POST labels calculated above as well.)
```{r}
gg_T1<-ggplot(subset(mydata,nObsString %in% c(PRE1:POST1)), aes(nObsString,Velocity,color=StimPeriod)) + geom_line()+facet_wrap(Animal~Group,scales='free',labeller = labeller(.cols = label_parsed, .multi_line = FALSE))
gg_T1
```
```{r}
gg_T2<-ggplot(subset(mydata,nObsString %in% c(PRE2:POST2)), aes(nObsString,Velocity,color=StimPeriod)) + geom_line()+facet_wrap(Animal~Group,scales='free',labeller = labeller(.cols = label_parsed, .multi_line = FALSE))
gg_T2
```

```{r}
gg_T3<-ggplot(subset(mydata,nObsString %in% c(PRE3:POST3)), aes(nObsString,Velocity,color=StimPeriod)) + geom_line()+facet_wrap(Animal~Group,scales='free',labeller = labeller(.cols = label_parsed, .multi_line = FALSE))
gg_T3
```
```{r}
gg_T4<-ggplot(subset(mydata,nObsString %in% c(PRE4:POST4)), aes(nObsString,Velocity,color=StimPeriod)) + geom_line()+facet_wrap(Animal~Group,scales='free',labeller = labeller(.cols = label_parsed, .multi_line = FALSE))
gg_T4
```
```{r}
gg_T5<-ggplot(subset(mydata,nObsString %in% c(PRE5:POST5)), aes(nObsString,Velocity,color=StimPeriod)) + geom_line()+facet_wrap(Animal~Group,scales='free',labeller = labeller(.cols = label_parsed, .multi_line = FALSE))
gg_T5
```

So far, in the above plots, each trial (for Acquisition Day 1 from a Luthi mouse replication group as an example) is visualized. They need some aesthetic work if ever to be published, but are a nice way to visualize what's going on and check that everything looks right.

Now for quantifying the darts. Instead of using a findpeaks package and call, I'll try to do it in a more 'simple' logical way using nested for and if loops. FindPeaks is potentially a faster method computationally, but the way the code is currently written is still very efficient.

First, to calculate some potential darting thresholds.**This should only be done on day 1 of training** , before any stimuli or shocks etc. For subsequent days, you should use the threshold determined from day 1 (either included early on in the third chunk on here.) 

```{r}
BL<-subset(mydata,StimPeriod %in% 'BL') #this makes a subset of the data only using baseline values (if coded as BL)
M2SD<-mean(BL$Velocity,na.rm=T)+(2*sd(BL$Velocity,na.rm=T))
M3SD<-mean(BL$Velocity,na.rm=T)+(3*sd(BL$Velocity,na.rm=T))
ile95<-as.numeric(quantile(BL$Velocity,.95,na.rm=T))
ile975<-as.numeric(quantile(BL$Velocity,.975,na.rm=T))
ile99<-as.numeric(quantile(BL$Velocity,.99,na.rm=T))
ile995<-as.numeric(quantile(BL$Velocity,.995,na.rm=T))
ile999<-as.numeric(quantile(BL$Velocity,.999,na.rm=T))
M2SD;M3SD;ile95;ile975;ile99;ile995;ile999
```

Here is a way to make a quick data frame of these thresholds based on this experiment only and save it to your directory that you can pull up for later days of this experiment (in addition to/instead of analyzing using overall baseline data to get thresholds for darting). 
***Again, this should only be done on day 1***
```{r}
CURRENTthreshDF<-data.frame(M2SD,M3SD,ile95,ile975,ile99,ile995,ile999)
CURRENTthreshDF
save(CURRENTthreshDF, file="currentThresh.Rda")
```

You can choose your data threshold and only continue with that if wanted. If you want to objectively check all of them, that's a good idea if you are worried about a threshold being too low/too high in terms of missing actual darts, calling non-darts darts, and/or having too high of a baseline dart level. If you desire to check if these values are correlated with freezing for example, you can continue wiht all (with copy/pasting and adjusting a few bits) and run your correlations to see what level to choose.
*In this experiment, all thresholds except the 99th percentile (ile99; p = .08) lead to darting values which are highly correlated (during CS periods) with freezing. Therefore, this example will continue primarily using the ile99 threshold.

Below is calculating the total number of frames above threshold. This is not darting episodes and will be adjusted in the chunk.
The chunk two below is for using thresholds calcuated based on day 1 BL data for this experiment only. **If you are on another day of the experiment, just run the chunk immediately below to load in your current experiment's threshold data. ONLY run this chunk after Day 1. You will need to add a {r} after the ``` for it to run.**
```
load("currentThresh.Rda")
M2SD<-CURRENTthreshDF[1,1]
M3SD<-CURRENTthreshDF[1,2]
ile95<-CURRENTthreshDF[1,3]
ile975<-CURRENTthreshDF[1,4]
ile99<-CURRENTthreshDF[1,5]
ile995<-CURRENTthreshDF[1,6]
ile999<-CURRENTthreshDF[1,7]
M2SD;M3SD;ile95;ile975;ile99;ile995;ile999
```
```
M2SD_frames<-as.numeric(mydata$Velocity>M2SD)
M3SD_frames<-as.numeric(mydata$Velocity>M3SD)
ile95_frames<-as.numeric(mydata$Velocity>ile95)
ile975_frames<-as.numeric(mydata$Velocity>ile975)
ile99_frames<-as.numeric(mydata$Velocity>ile99)
ile995_frames<-as.numeric(mydata$Velocity>ile995)
ile999_frames<-as.numeric(mydata$Velocity>ile999)
```



This chunk is for using thresholds calculated based on all previously analyzed baseline data across experiments of the species you're using.
```{r}
#overM2SD_frames<-as.numeric(mydata$Velocity>threshDF$M2SD)
#overM3SD_frames<-as.numeric(mydata$Velocity>threshDF$M3SD)
#overile95_frames<-as.numeric(mydata$Velocity>threshDF$ile95)
#overile975_frames<-as.numeric(mydata$Velocity>threshDF$ile975)
#overile99_frames<-as.numeric(mydata$Velocity>threshDF$ile99)
overile995_frames<-as.numeric(mydata$Velocity>threshDF$ile995)
#overile999_frames<-as.numeric(mydata$Velocity>threshDF$ile999)
```
```{r}
#mydata<-cbind(mydata,M2SD_frames,M3SD_frames,ile95_frames,ile975_frames,ile99_frames,ile995_frames,ile999_frames,overM2SD_frames,overM3SD_frames,overile95_frames,overile975_frames,overile99_frames,overile995_frames,overile999_frames)
mydata<-cbind(mydata,overile995_frames)
head(mydata)
```


The code for determining darting episodes is below for all example thresholds used above, but only the 99th percentile (ile99) will run unless you add {r} after the upper ``` above the first line of each chunk. Each chunk will store a numeric string of darting episodes (at 3.75 Hz) and will produce a readout of the total number of darting episodes for the supplied threshold.
```
M2SD_frames[is.na(M2SD_frames)]<-0
M2SD_darts<-as.numeric(c(rep(0,length(M2SD_frames))))

if (M2SD_frames[1]==1){M2SD_darts[1]<-1}
if (M2SD_frames[2]==1 & M2SD_frames[1]==0){M2SD_darts[2]==1}
if (M2SD_frames[3]==1 & M2SD_frames[2]==0 & M2SD_frames[1]==0){M2SD_darts[3]==1}
irange<-4:length(M2SD_frames)
for (i in irange) {
  if (M2SD_frames[i] == 1 & M2SD_frames[i-1] == 0 & M2SD_frames[i-2]==0 & M2SD_frames[i-3]==0){ 
    M2SD_darts[i]<-1 
  }
}
sum(M2SD_darts)
```
```
M3SD_frames[is.na(M3SD_frames)]<-0
M3SD_darts<-as.numeric(c(rep(0,length(M3SD_frames))))

if (M3SD_frames[1]==1){M3SD_darts[1]<-1}
if (M3SD_frames[2]==1 & M3SD_frames[1]==0){M3SD_darts[2]==1}
if (M3SD_frames[3]==1 & M3SD_frames[2]==0 & M3SD_frames[1]==0){M3SD_darts[3]==1}
irange<-4:length(M3SD_frames)
for (i in irange) {
  if (M3SD_frames[i] == 1 & M3SD_frames[i-1] == 0 & M3SD_frames[i-2]==0 & M3SD_frames[i-3]==0){ 
    M3SD_darts[i]<-1 
  }
}
sum(M3SD_darts)
```
```
ile95_frames[is.na(ile95_frames)]<-0
ile95_darts<-as.numeric(c(rep(0,length(ile95_frames))))

if (ile95_frames[1]==1){ile95_darts[1]<-1}
if (ile95_frames[2]==1 & ile95_frames[1]==0){ile95_darts[2]==1}
if (ile95_frames[3]==1 & ile95_frames[2]==0 & ile95_frames[1]==0){ile95_darts[3]==1}
irange<-4:length(ile95_frames)
for (i in irange) {
  if (ile95_frames[i] == 1 & ile95_frames[i-1] == 0 & ile95_frames[i-2]==0 & ile95_frames[i-3]==0){ 
    ile95_darts[i]<-1 
  }
}
sum(ile95_darts)
```
```
ile975_frames[is.na(ile975_frames)]<-0
ile975_darts<-as.numeric(c(rep(0,length(ile975_frames))))

if (ile975_frames[1]==1){ile975_darts[1]<-1}
if (ile975_frames[2]==1 & ile975_frames[1]==0){ile975_darts[2]==1}
if (ile975_frames[3]==1 & ile975_frames[2]==0 & ile975_frames[1]==0){ile975_darts[3]==1}
irange<-4:length(ile975_frames)
for (i in irange) {
  if (ile975_frames[i] == 1 & ile975_frames[i-1] == 0 & ile975_frames[i-2]==0 & ile975_frames[i-3]==0){ 
    ile975_darts[i]<-1 
  }
}
sum(ile975_darts)
```
```
ile99_frames[is.na(ile99_frames)]<-0
ile99_darts<-as.numeric(c(rep(0,length(ile99_frames))))

if (ile99_frames[1]==1){ile99_darts[1]<-1}
if (ile99_frames[2]==1 & ile99_frames[1]==0){ile99_darts[2]==1}
if (ile99_frames[3]==1 & ile99_frames[2]==0 & ile99_frames[1]==0){ile99_darts[3]==1}
irange<-4:length(ile99_frames)
for (i in irange) {
  if (ile99_frames[i] == 1 & ile99_frames[i-1] == 0 & ile99_frames[i-2]==0 & ile99_frames[i-3]==0){ 
    ile99_darts[i]<-1 
  }
}
sum(ile99_darts)
```
```
ile995_frames[is.na(ile995_frames)]<-0
ile995_darts<-as.numeric(c(rep(0,length(ile995_frames))))

if (ile995_frames[1]==1){ile995_darts[1]<-1}
if (ile995_frames[2]==1 & ile995_frames[1]==0){ile995_darts[2]==1}
if (ile995_frames[3]==1 & ile995_frames[2]==0 & ile995_frames[1]==0){ile995_darts[3]==1}
irange<-4:length(ile995_frames)
for (i in irange) {
  if (ile995_frames[i] == 1 & ile995_frames[i-1] == 0 & ile995_frames[i-2]==0 & ile995_frames[i-3]==0){ 
    ile995_darts[i]<-1 
  }
}
sum(ile995_darts)
```
```
ile999_frames[is.na(ile999_frames)]<-0
ile999_darts<-as.numeric(c(rep(0,length(ile999_frames))))

if (ile999_frames[1]==1){ile999_darts[1]<-1}
if (ile999_frames[2]==1 & ile999_frames[1]==0){ile999_darts[2]==1}
if (ile999_frames[3]==1 & ile999_frames[2]==0 & ile999_frames[1]==0){ile999_darts[3]==1}
irange<-4:length(ile999_frames)
for (i in irange) {
  if (ile999_frames[i] == 1 & ile999_frames[i-1] == 0 & ile999_frames[i-2]==0 & ile999_frames[i-3]==0){ 
    ile999_darts[i]<-1 
  }
}
sum(ile999_darts)
```
```
overM2SD_frames[is.na(overM2SD_frames)]<-0
overM2SD_darts<-as.numeric(c(rep(0,length(overM2SD_frames))))

if (overM2SD_frames[1]==1){overM2SD_darts[1]<-1}
if (overM2SD_frames[2]==1 & overM2SD_frames[1]==0){overM2SD_darts[2]==1}
if (overM2SD_frames[3]==1 & overM2SD_frames[2]==0 & overM2SD_frames[1]==0){overM2SD_darts[3]==1}
irange<-4:length(overM2SD_frames)
for (i in irange) {
  if (overM2SD_frames[i] == 1 & overM2SD_frames[i-1] == 0 & overM2SD_frames[i-2]==0 & overM2SD_frames[i-3]==0){ 
    overM2SD_darts[i]<-1 
  }
}
sum(overM2SD_darts)
```
```
overM3SD_frames[is.na(overM3SD_frames)]<-0
overM3SD_darts<-as.numeric(c(rep(0,length(overM3SD_frames))))

if (overM3SD_frames[1]==1){overM3SD_darts[1]<-1}
if (overM3SD_frames[2]==1 & overM3SD_frames[1]==0){overM3SD_darts[2]==1}
if (overM3SD_frames[3]==1 & overM3SD_frames[2]==0 & overM3SD_frames[1]==0){overM3SD_darts[3]==1}
irange<-4:length(overM3SD_frames)
for (i in irange) {
  if (overM3SD_frames[i] == 1 & overM3SD_frames[i-1] == 0 & overM3SD_frames[i-2]==0 & overM3SD_frames[i-3]==0){ 
    overM3SD_darts[i]<-1 
  }
}
sum(overM3SD_darts)
```
```
overile95_frames[is.na(overile95_frames)]<-0
overile95_darts<-as.numeric(c(rep(0,length(overile95_frames))))

if (overile95_frames[1]==1){overile95_darts[1]<-1}
if (overile95_frames[2]==1 & overile95_frames[1]==0){overile95_darts[2]==1}
if (overile95_frames[3]==1 & overile95_frames[2]==0 & overile95_frames[1]==0){overile95_darts[3]==1}
irange<-4:length(overile95_frames)
for (i in irange) {
  if (overile95_frames[i] == 1 & overile95_frames[i-1] == 0 & overile95_frames[i-2]==0 & overile95_frames[i-3]==0){ 
    overile95_darts[i]<-1 
  }
}
sum(overile95_darts)
```
```
overile975_frames[is.na(overile975_frames)]<-0
overile975_darts<-as.numeric(c(rep(0,length(overile975_frames))))

if (overile975_frames[1]==1){overile975_darts[1]<-1}
if (overile975_frames[2]==1 & overile975_frames[1]==0){overile975_darts[2]==1}
if (overile975_frames[3]==1 & overile975_frames[2]==0 & overile975_frames[1]==0){overile975_darts[3]==1}
irange<-4:length(overile975_frames)
for (i in irange) {
  if (overile975_frames[i] == 1 & overile975_frames[i-1] == 0 & overile975_frames[i-2]==0 & overile975_frames[i-3]==0){ 
    overile975_darts[i]<-1 
  }
}
sum(overile975_darts)
```
```
overile99_frames[is.na(overile99_frames)]<-0
overile99_darts<-as.numeric(c(rep(0,length(overile99_frames))))

if (overile99_frames[1]==1){overile99_darts[1]<-1}
if (overile99_frames[2]==1 & overile99_frames[1]==0){overile99_darts[2]==1}
if (overile99_frames[3]==1 & overile99_frames[2]==0 & overile99_frames[1]==0){overile99_darts[3]==1}
irange<-4:length(overile99_frames)
for (i in irange) {
  if (overile99_frames[i] == 1 & overile99_frames[i-1] == 0 & overile99_frames[i-2]==0 & overile99_frames[i-3]==0){ 
    overile99_darts[i]<-1 
  }
}
sum(overile99_darts)
```
```{r}
overile995_frames[is.na(overile995_frames)]<-0
overile995_darts<-as.numeric(c(rep(0,length(overile995_frames))))

if (overile995_frames[1]==1){overile995_darts[1]<-1}
if (overile995_frames[2]==1 & overile995_frames[1]==0){overile995_darts[2]==1}
if (overile995_frames[3]==1 & overile995_frames[2]==0 & overile995_frames[1]==0){overile995_darts[3]==1}
irange<-4:length(overile995_frames)
for (i in irange) {
  if (overile995_frames[i] == 1 & overile995_frames[i-1] == 0 & overile995_frames[i-2]==0 & overile995_frames[i-3]==0){ 
    overile995_darts[i]<-1 
  }
}
sum(overile995_darts)
```
```
overile999_frames[is.na(overile999_frames)]<-0
overile999_darts<-as.numeric(c(rep(0,length(overile999_frames))))

if (overile999_frames[1]==1){overile999_darts[1]<-1}
if (overile999_frames[2]==1 & overile999_frames[1]==0){overile999_darts[2]==1}
if (overile999_frames[3]==1 & overile999_frames[2]==0 & overile999_frames[1]==0){overile999_darts[3]==1}
irange<-4:length(overile999_frames)
for (i in irange) {
  if (overile999_frames[i] == 1 & overile999_frames[i-1] == 0 & overile999_frames[i-2]==0 & overile999_frames[i-3]==0){ 
    overile999_darts[i]<-1 
  }
}
sum(overile999_darts)
```


Let's add the darting episodes to our dataframe.
```{r}
#mydata<-cbind(mydata,M2SD_darts,M3SD_darts,ile95_darts,ile975_darts,ile99_darts,ile995_darts,ile999_darts,overM2SD_darts,overM3SD_darts,overile95_darts,overile975_darts,overile99_darts,overile995_darts,overile999_darts)  # run this line if you have calculated darting episodes based on all potential thresholds. #,M2SD_darts,M3SD_darts,ile95_darts,ile975_darts,ile99_darts,ile995_darts,ile999_darts,overM2SD_darts,overM3SD_darts,overile95_darts,overile975_darts,overile99_darts,overile995_darts,overile999_darts)
mydata<-cbind(mydata,overile995_darts)      #If you only have one, or a few thresholds you are using, add the name of the number-of-darts column after the comma in this line. For example, here would just be adding the darts using the 99.5% threshold based on all available previous baseline data.
head(mydata)
```


Finally, let's count up the darting episodes in each Stimulus period for each animal. You can run just the lines relevant to the data you want.
```{r}
#Darts_M2SD<-aggregate(mydata$M2SD_darts,by=list(mydata$Animal,mydata$Group,mydata$Sex,mydata$StimPeriod),FUN="sum")
#Darts_M3SD<-aggregate(mydata$M3SD_darts,by=list(mydata$Animal,mydata$Group,mydata$Sex,mydata$StimPeriod),FUN="sum")
#Darts_95ile<-aggregate(mydata$ile95_darts,by=list(mydata$Animal,mydata$Group,mydata$Sex,mydata$StimPeriod),FUN="sum")
#Darts_975ile<-aggregate(mydata$ile975_darts,by=list(mydata$Animal,mydata$Group,mydata$Sex,mydata$StimPeriod),FUN="sum")
#Darts_99ile<-aggregate(mydata$ile99_darts,by=list(mydata$Animal,mydata$Group,mydata$Sex,mydata$StimPeriod),FUN="sum")
#Darts_995ile<-aggregate(mydata$ile995_darts,by=list(mydata$Animal,mydata$Group,mydata$Sex,mydata$StimPeriod),FUN="sum")
#Darts_999ile<-aggregate(mydata$ile999_darts,by=list(mydata$Animal,mydata$Group,mydata$Sex,mydata$StimPeriod),FUN="sum")
#Darts_overM2SD<-aggregate(mydata$overM2SD_darts,by=list(mydata$Animal,mydata$Group,mydata$Sex,mydata$StimPeriod),FUN="sum")
#Darts_overM3SD<-aggregate(mydata$overM3SD_darts,by=list(mydata$Animal,mydata$Group,mydata$Sex,mydata$StimPeriod),FUN="sum")
#Darts_over95ile<-aggregate(mydata$overile95_darts,by=list(mydata$Animal,mydata$Group,mydata$Sex,mydata$StimPeriod),FUN="sum")
#Darts_over975ile<-aggregate(mydata$overile975_darts,by=list(mydata$Animal,mydata$Group,mydata$Sex,mydata$StimPeriod),FUN="sum")
#Darts_over99ile<-aggregate(mydata$overile99_darts,by=list(mydata$Animal,mydata$Group,mydata$Sex,mydata$StimPeriod),FUN="sum")
Darts_over995ile<-aggregate(mydata$overile995_darts,by=list(mydata$Animal,mydata$Group,mydata$Sex,mydata$StimPeriod),FUN="sum")
#Darts_over999ile<-aggregate(mydata$overile999_darts,by=list(mydata$Animal,mydata$Group,mydata$Sex,mydata$StimPeriod),FUN="sum")
```

```{r}
#dartcount<-cbind(Darts_M2SD,Darts_M3SD$x,Darts_95ile$x,Darts_975ile$x,Darts_99ile$x,Darts_995ile$x,Darts_999ile$x,Darts_overM2SD$x,Darts_overM3SD$x,Darts_over95ile$x,Darts_over975ile$x,Darts_over99ile$x,Darts_over995ile$x,Darts_over999ile$x)   #Use this line if you want to use all threshold values.
#names(dartcount)<-c("Animal","Group","Sex","StimPeriod","Darts_M2SD","Darts_M3SD","Darts_95ile","Darts_975ile","Darts_99ile","Darts_995ile","Darts_999ile","Darts_overM2SD","Darts_overM3SD","Darts_over95ile","Darts_over975ile","Darts_over99ile","Darts_over995ile","Darts_over999ile")
dartcount<-Darts_over995ile  # Use this line if you just want to use one threshold value.
names(dartcount)<-c("Animal","Group","Sex","StimPeriod","Darts_over995ile")
dartcount
write.csv(dartcount,"DartCount_Exp1Day1.csv") #Rename your file here
```
And Voila, we Have dart counts that match with what is scored manually. The only difference/caveat is that there is one dart for each activity burst to the shock. This would be trivial to fix by just subtracting 1 from each ITI and POST stimPeriod, but it would be quite difficult and overly tedious to build into the code in any general way. If this is something we ever care about, you should check visually with the graph if there was a dart/activity burst and manually subtract one.(It would be sufficient to just scan the graphs to look for any which don't have activity bursts). This will also put a csv file (openable in Excel) in your working directory which can be easily sorted and rearranged for any analysis you might want/need to do with Excel and/or SPSS.

More stuff below to help visualize representative animals in pub-friendly images.
```{r}
minrate<-.266/60
end<-nObs*minrate
timeseq<-seq(from = minrate, to = end, by = minrate)
```     
```{r}
tpre1<-timeseq[PRE1];tpost1<-timeseq[POST1]
tpre2<-timeseq[PRE2];tpost2<-timeseq[POST2]
tpre3<-timeseq[PRE3];tpost3<-timeseq[POST3]
tpre4<-timeseq[PRE4];tpost4<-timeseq[POST4]
tpre5<-timeseq[PRE5];tpost5<-timeseq[POST5]
```

Below will give you a snapshot of the whole session for one animal.
```{r}
ggallfunc<-function(n){
ggn<-ggplot(mydata[which(mydata$Animal==n),], aes(timeseq,Velocity,color=StimPeriod)) + geom_line()+facet_wrap(Animal~Group,scales='free',labeller = labeller(.cols = label_parsed, .multi_line = FALSE))+geom_hline(yintercept=threshDF$ile995, linetype="dashed", color = "red")+theme_classic()+xlab("Time (min)")
print(ggn)}
```
```{r}
#ggallfunc() #This function will print the whole session for whatever animal number you use.
```

Below will give you a trial-by-trial plot of the whole session for one animal.
```{r}
ggarrangefunc<-function(n){
mydatan<-subset(mydata[which(mydata$Animal==n),])
ggn_T1<-ggplot(subset(mydatan,nObsString %in% c(PRE1:POST1)), aes(nObsString,Velocity,color=StimPeriod)) +
  geom_line()+
  geom_hline(yintercept=threshDF$ile995, linetype="dashed", color = "black")+
  geom_vline(xintercept = nObsString[mydata$StimPeriod=="CS1a"][1][],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nObsString[mydata$StimPeriod=="CS1b"][1][],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nObsString[mydata$StimPeriod=="IT1"][1][],linetype="dashed", color = "grey")+
  theme_classic()+
  ylab("Velocity (cm/s)")+
  scale_color_manual(name="Stimulus Period",labels=c("PreCS","Tone","White Noise","Shock"),values = c("darkgreen","blue","black","red"))+
  scale_x_continuous(name="Time (s)",limits = c(PRE1,POST1),labels =c("-10","0","10","20","30"),breaks=c(PRE1,nObsString[mydata$StimPeriod=="CS1a"][1][],nObsString[mydata$StimPeriod=="CS1b"][1][],nObsString[mydata$StimPeriod=="IT1"][1][],POST1))+
  geom_path(aes(group=1)) +
  scale_y_continuous(limits = c(0,125))+
  ggtitle(paste("Animal: ", n, mydatan$Group[which(mydatan$Animal==n)]))
#ggn_T1

ggn_T2<-ggplot(subset(mydatan,nObsString %in% c(PRE2:POST2)), aes(nObsString,Velocity,color=StimPeriod)) +
  geom_line()+
  geom_hline(yintercept=threshDF$ile995, linetype="dashed", color = "black")+
  geom_vline(xintercept = nObsString[mydata$StimPeriod=="CS2a"][1][],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nObsString[mydata$StimPeriod=="CS2b"][1][],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nObsString[mydata$StimPeriod=="IT2"][1][],linetype="dashed", color = "grey")+
  theme_classic()+
  ylab("Velocity (cm/s)")+
  scale_color_manual(name="Stimulus Period",labels=c("PreCS","Tone","White Noise","Shock"),values = c("blue","black","darkgreen","red"))+
  theme(legend.position = "none")+
  scale_x_continuous(name="Time (s)",limits = c(PRE2,POST2),labels =c("-10","0","10","20","30"),breaks=c(PRE2,nObsString[mydata$StimPeriod=="CS2a"][1][],nObsString[mydata$StimPeriod=="CS2b"][1][],nObsString[mydata$StimPeriod=="IT2"][1][],POST2))+
  geom_path(aes(group=1))+
  scale_y_continuous(limits = c(0,125))
#ggn_T2

ggn_T3<-ggplot(subset(mydatan,nObsString %in% c(PRE3:POST3)), aes(nObsString,Velocity,color=StimPeriod)) +
  geom_line()+
  geom_hline(yintercept=threshDF$ile995, linetype="dashed", color = "black")+
  geom_vline(xintercept = nObsString[mydata$StimPeriod=="CS3a"][1][],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nObsString[mydata$StimPeriod=="CS3b"][1][],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nObsString[mydata$StimPeriod=="IT3"][1][],linetype="dashed", color = "grey")+
  theme_classic()+
  ylab("Velocity (cm/s)")+
  scale_color_manual(name="Stimulus Period",labels=c("PreCS","Tone","White Noise","Shock"),values = c("blue","black","darkgreen","red"))+
  theme(legend.position = "none")+
  scale_x_continuous(name="Time (s)",limits = c(PRE3,POST3),labels =c("-10","0","10","20","30"),breaks=c(PRE3,nObsString[mydata$StimPeriod=="CS3a"][1][],nObsString[mydata$StimPeriod=="CS3b"][1][],nObsString[mydata$StimPeriod=="IT3"][1][],POST3))+
  geom_path(aes(group=1))+
  scale_y_continuous(limits = c(0,125))
#ggn_T3

ggn_T4<-ggplot(subset(mydatan,nObsString %in% c(PRE4:POST4)), aes(nObsString,Velocity,color=StimPeriod)) +
  geom_line()+
  geom_hline(yintercept=threshDF$ile995, linetype="dashed", color = "black")+
  geom_vline(xintercept = nObsString[mydata$StimPeriod=="CS4a"][1][],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nObsString[mydata$StimPeriod=="CS4b"][1][],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nObsString[mydata$StimPeriod=="IT4"][1][],linetype="dashed", color = "grey")+
  theme_classic()+
  ylab("Velocity (cm/s)")+
  scale_color_manual(name="Stimulus Period",labels=c("PreCS","Tone","White Noise","Shock"),values = c("blue","black","darkgreen","red"))+
  theme(legend.position = "none")+
  scale_x_continuous(name="Time (s)",limits = c(PRE4,POST4),labels =c("-10","0","10","20","30"),breaks=c(PRE4,nObsString[mydata$StimPeriod=="CS4a"][1][],nObsString[mydata$StimPeriod=="CS4b"][1][],nObsString[mydata$StimPeriod=="IT4"][1][],POST4))+
  geom_path(aes(group=1))+
  scale_y_continuous(limits = c(0,125))
#ggn_T4

ggn_T5<-ggplot(subset(mydatan,nObsString %in% c(PRE5:POST5)), aes(nObsString,Velocity,color=StimPeriod)) +
  geom_line()+
  geom_hline(yintercept=threshDF$ile995, linetype="dashed", color = "black")+
  geom_vline(xintercept = nObsString[mydata$StimPeriod=="CS5a"][1][],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nObsString[mydata$StimPeriod=="CS5b"][1][],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nObsString[mydata$StimPeriod=="POST"][1][],linetype="dashed", color = "grey")+
  theme_classic()+
  ylab("Velocity (cm/s)")+
  scale_color_manual(name="Stimulus Period",labels=c("Baseline","Tone","White Noise","Shock"),values = c("blue","black","darkgreen","red"))+
  theme(legend.position = "none")+
  scale_x_continuous(name="Time (s)",limits = c(PRE5,POST5),labels =c("-10","0","10","20","30"),breaks=c(PRE5,nObsString[mydata$StimPeriod=="CS5a"][1][],nObsString[mydata$StimPeriod=="CS5b"][1][],nObsString[mydata$StimPeriod=="POST"][1][],POST5))+
  geom_path(aes(group=1))+
  scale_y_continuous(limits = c(0,125))
#ggn_T5

figure<-ggarrange(ggn_T1, ggn_T2, ggn_T3, ggn_T4, ggn_T5,
                    labels = c("Trial 1", "Trial 2", "Trial 3","Trial 4","Trial 5"),label.y=0.2,label.x=.01,
                    ncol = 1, nrow = 5,
          common.legend = TRUE, legend = "top")
annotate_figure(figure,top = text_grob("Visualizing Tooth Growth", color = "red", face = "bold", size = 14))
print(figure)
}
```
```{r}
#ggarrangefunc()   #This function will print all 5 trials for whatever animal number you use.
```



##To get nice traces, run the chunks below
```{r}
mydata<-read.csv("Statistics-10-18_Experiment3_Day3Ext.csv")
nObs<-as.numeric(data.table(mydata)[,.N,by=Animal][1,2])
samplen<-length(mydata$Trial)/nObs
nObsString<-rep(c(1:nObs),samplen)
mydata<-cbind(mydata,nObsString)
mydata$Distance_cm<-as.numeric(as.character(mydata$Distance_cm))
mydata$Velocity<-as.numeric(as.character(mydata$Velocity))
mydata$Animal<-as.factor(mydata$Animal)
mydata[is.na(mydata)]<-0
mydata$Group = factor(mydata$Group, levels=c("TWN_T","TWN_WN","WNT_T","WNT_WN"),labels = c('`Tone-Noise Tone Test`','`Tone-Noise Noise Test`','`Noise-Tone Tone Test`','`Noise-Tone Noise Test`'))
mydata$StimPeriod<-factor(mydata$StimPeriod, levels = c("BL","CS1a","CS1b","IT1","CS2a","CS2b","ITI2","CS3a","CS3b","IT3","CS4a","CS4b","IT4","CS5a","CS5b","IT5","CS6a","CS6b","IT6","CS7a","CS7b","IT7","CS8a","CS8b","IT8","CS9a","CS9b","IT9","CS10a","CS10b","IT10","CS11a","CS11b","IT11","CS12a","CS12b","IT12","CS13a","CS13b","IT13","CS14a","CS14b","IT14","CS15a","CS15b","IT15","CS16a","CS16b","POST"),labels = c("BL","CS1a","CS1b","IT1","CS2a","CS2b","IT2","CS3a","CS3b","IT3","CS4a","CS4b","IT4","CS5a","CS5b","IT5","CS6a","CS6b","IT6","CS7a","CS7b","IT7","CS8a","CS8b","IT8","CS9a","CS9b","IT9","CS10a","CS10b","IT10","CS11a","CS11b","IT11","CS12a","CS12b","IT12","CS13a","CS13b","IT13","CS14a","CS14b","IT14","CS15a","CS15b","IT15","CS16a","CS16b","POST"))
PRE1<-mydata$nObsString[min(which(mydata$StimPeriod == "CS1a"))-38]
POST1<-mydata$nObsString[min(which(mydata$StimPeriod == "CS1b"))+76]
PRE2<-mydata$nObsString[min(which(mydata$StimPeriod == "CS2a"))-38]
POST2<-mydata$nObsString[min(which(mydata$StimPeriod == "CS2b"))+76]
PRE3<-mydata$nObsString[min(which(mydata$StimPeriod == "CS3a"))-38]
POST3<-mydata$nObsString[min(which(mydata$StimPeriod == "CS3b"))+76]
PRE4<-mydata$nObsString[min(which(mydata$StimPeriod == "CS4a"))-38]
POST4<-mydata$nObsString[min(which(mydata$StimPeriod == "CS4b"))+76]
PRE5<-mydata$nObsString[min(which(mydata$StimPeriod == "CS5a"))-38]
POST5<-mydata$nObsString[min(which(mydata$StimPeriod == "CS5b"))+76]
PRE6<-mydata$nObsString[min(which(mydata$StimPeriod == "CS6a"))-38]
POST6<-mydata$nObsString[min(which(mydata$StimPeriod == "CS6b"))+76]
PRE7<-mydata$nObsString[min(which(mydata$StimPeriod == "CS7a"))-38]
POST7<-mydata$nObsString[min(which(mydata$StimPeriod == "CS7b"))+76]
PRE8<-mydata$nObsString[min(which(mydata$StimPeriod == "CS8a"))-38]
POST8<-mydata$nObsString[min(which(mydata$StimPeriod == "CS8b"))+76]
PRE9<-mydata$nObsString[min(which(mydata$StimPeriod == "CS9a"))-38]
POST9<-mydata$nObsString[min(which(mydata$StimPeriod == "CS9b"))+76]
PRE10<-mydata$nObsString[min(which(mydata$StimPeriod == "CS10a"))-38]
POST10<-mydata$nObsString[min(which(mydata$StimPeriod == "CS10b"))+76]
PRE11<-mydata$nObsString[min(which(mydata$StimPeriod == "CS11a"))-38]
POST11<-mydata$nObsString[min(which(mydata$StimPeriod == "CS11b"))+76]
PRE12<-mydata$nObsString[min(which(mydata$StimPeriod == "CS12a"))-38]
POST12<-mydata$nObsString[min(which(mydata$StimPeriod == "CS12b"))+76]
PRE13<-mydata$nObsString[min(which(mydata$StimPeriod == "CS13a"))-38]
POST13<-mydata$nObsString[min(which(mydata$StimPeriod == "CS13b"))+76]
PRE14<-mydata$nObsString[min(which(mydata$StimPeriod == "CS14a"))-38]
POST14<-mydata$nObsString[min(which(mydata$StimPeriod == "CS14b"))+76]
PRE15<-mydata$nObsString[min(which(mydata$StimPeriod == "CS15a"))-38]
POST15<-mydata$nObsString[min(which(mydata$StimPeriod == "CS15b"))+76]
PRE16<-mydata$nObsString[min(which(mydata$StimPeriod == "CS16a"))-38]
POST16<-mydata$nObsString[min(which(mydata$StimPeriod == "CS16b"))+76]
```
```{r}
mydatan<-mydata
mydatan$StimPeriod<-factor(mydatan$StimPeriod,labels = c("preCS","Tone","Noise","postCS","CS2a","CS2b","IT2","CS3a","CS3b","IT3","CS4a","CS4b","IT4","CS5a","CS5b","IT5","CS6a","CS6b","IT6","CS7a","CS7b","IT7","CS8a","CS8b","IT8","CS9a","CS9b","IT9","CS10a","CS10b","IT10","CS11a","CS11b","IT11","CS12a","CS12b","IT12","CS13a","CS13b","IT13","CS14a","CS14b","IT14","CS15a","CS15b","IT15","CS16a","CS16b","POST"))
mydata1<-mydatan[which(nObsString %in% c(PRE1:POST1)),]
nObsTrial<-length(mydata1$Trial)/samplen
nTrialM<-rep(c(1:nObsTrial),samplen)
mydata1<-cbind(mydata1,nTrialM)

mydatan<-mydata
mydatan$StimPeriod<-factor(mydatan$StimPeriod,labels = c("BL","CS1a","CS1b","preCS","Tone","Noise","postCS","CS3a","CS3b","IT3","CS4a","CS4b","IT4","CS5a","CS5b","IT5","CS6a","CS6b","IT6","CS7a","CS7b","IT7","CS8a","CS8b","IT8","CS9a","CS9b","IT9","CS10a","CS10b","IT10","CS11a","CS11b","IT11","CS12a","CS12b","IT12","CS13a","CS13b","IT13","CS14a","CS14b","IT14","CS15a","CS15b","IT15","CS16a","CS16b","POST"))
mydata2<-mydatan[which(nObsString %in% c(PRE2:POST2)),]
nObsTrial<-length(mydata2$Trial)/samplen
nTrialM<-rep(c(1:nObsTrial),samplen)
mydata2<-cbind(mydata2,nTrialM)

mydatan<-mydata
mydatan$StimPeriod<-factor(mydatan$StimPeriod,labels = c("BL","CS1a","CS1b","IT1","CS2a","CS2b","preCS","Tone","Noise","postCS","CS4a","CS4b","IT4","CS5a","CS5b","IT5","CS6a","CS6b","IT6","CS7a","CS7b","IT7","CS8a","CS8b","IT8","CS9a","CS9b","IT9","CS10a","CS10b","IT10","CS11a","CS11b","IT11","CS12a","CS12b","IT12","CS13a","CS13b","IT13","CS14a","CS14b","IT14","CS15a","CS15b","IT15","CS16a","CS16b","POST"))
mydata3<-mydatan[which(nObsString %in% c(PRE3:POST3)),]
nObsTrial<-length(mydata3$Trial)/samplen
nTrialM<-rep(c(1:nObsTrial),samplen)
mydata3<-cbind(mydata3,nTrialM)

mydatan<-mydata
mydatan$StimPeriod<-factor(mydatan$StimPeriod,labels = c("BL","CS1a","CS1b","IT1","CS2a","CS2b","IT2","CS3a","CS3b","preCS","Tone","Noise","postCS","CS5a","CS5b","IT5","CS6a","CS6b","IT6","CS7a","CS7b","IT7","CS8a","CS8b","IT8","CS9a","CS9b","IT9","CS10a","CS10b","IT10","CS11a","CS11b","IT11","CS12a","CS12b","IT12","CS13a","CS13b","IT13","CS14a","CS14b","IT14","CS15a","CS15b","IT15","CS16a","CS16b","POST"))
mydata4<-mydatan[which(nObsString %in% c(PRE4:POST4)),]
nObsTrial<-length(mydata4$Trial)/samplen
nTrialM<-rep(c(1:nObsTrial),samplen)
mydata4<-cbind(mydata4,nTrialM)

mydatan<-mydata
mydatan$StimPeriod<-factor(mydatan$StimPeriod,labels = c("BL","CS1a","CS1b","IT1","CS2a","CS2b","IT2","CS3a","CS3b","IT3","CS4a","CS4b","preCS","Tone","Noise","postCS","CS6a","CS6b","IT6","CS7a","CS7b","IT7","CS8a","CS8b","IT8","CS9a","CS9b","IT9","CS10a","CS10b","IT10","CS11a","CS11b","IT11","CS12a","CS12b","IT12","CS13a","CS13b","IT13","CS14a","CS14b","IT14","CS15a","CS15b","IT15","CS16a","CS16b","POST"))
mydata5<-mydatan[which(nObsString %in% c(PRE5:POST5)),]
nObsTrial<-length(mydata5$Trial)/samplen
nTrialM<-rep(c(1:nObsTrial),samplen)
mydata5<-cbind(mydata5,nTrialM)

mydatan<-mydata
mydatan$StimPeriod<-factor(mydatan$StimPeriod,labels = c("BL","CS1a","CS1b","IT1","CS2a","CS2b","IT2","CS3a","CS3b","IT3","CS4a","CS4b","IT4","CS5a","CS5b","preCS","Tone","Noise","postCS","CS7a","CS7b","IT7","CS8a","CS8b","IT8","CS9a","CS9b","IT9","CS10a","CS10b","IT10","CS11a","CS11b","IT11","CS12a","CS12b","IT12","CS13a","CS13b","IT13","CS14a","CS14b","IT14","CS15a","CS15b","IT15","CS16a","CS16b","POST"))
mydata6<-mydatan[which(nObsString %in% c(PRE6:POST6)),]
nObsTrial<-length(mydata6$Trial)/samplen
nTrialM<-rep(c(1:nObsTrial),samplen)
mydata6<-cbind(mydata6,nTrialM)

mydatan<-mydata
mydatan$StimPeriod<-factor(mydatan$StimPeriod,labels = c("BL","CS1a","CS1b","IT1","CS2a","CS2b","IT2","CS3a","CS3b","IT3","CS4a","CS4b","IT4","CS5a","CS5b","IT5","CS6a","CS6b","preCS","Tone","Noise","postCS","CS8a","CS8b","IT8","CS9a","CS9b","IT9","CS10a","CS10b","IT10","CS11a","CS11b","IT11","CS12a","CS12b","IT12","CS13a","CS13b","IT13","CS14a","CS14b","IT14","CS15a","CS15b","IT15","CS16a","CS16b","POST"))
mydata7<-mydatan[which(nObsString %in% c(PRE7:POST7)),]
nObsTrial<-length(mydata7$Trial)/samplen
nTrialM<-rep(c(1:nObsTrial),samplen)
mydata7<-cbind(mydata7,nTrialM)

mydatan<-mydata
mydatan$StimPeriod<-factor(mydatan$StimPeriod,labels = c("BL","CS1a","CS1b","IT1","CS2a","CS2b","IT2","CS3a","CS3b","IT3","CS4a","CS4b","IT4","CS5a","CS5b","IT5","CS6a","CS6b","IT6","CS7a","CS7b","preCS","Tone","Noise","postCS","CS9a","CS9b","IT9","CS10a","CS10b","IT10","CS11a","CS11b","IT11","CS12a","CS12b","IT12","CS13a","CS13b","IT13","CS14a","CS14b","IT14","CS15a","CS15b","IT15","CS16a","CS16b","POST"))
mydata8<-mydatan[which(nObsString %in% c(PRE8:POST8)),]
nObsTrial<-length(mydata8$Trial)/samplen
nTrialM<-rep(c(1:nObsTrial),samplen)
mydata8<-cbind(mydata8,nTrialM)

mydatan<-mydata
mydatan$StimPeriod<-factor(mydatan$StimPeriod,labels = c("BL","CS1a","CS1b","IT1","CS2a","CS2b","IT2","CS3a","CS3b","IT3","CS4a","CS4b","IT4","CS5a","CS5b","IT5","CS6a","CS6b","IT6","CS7a","CS7b","IT7","CS8a","CS8b","preCS","Tone","Noise","postCS","CS10a","CS10b","IT10","CS11a","CS11b","IT11","CS12a","CS12b","IT12","CS13a","CS13b","IT13","CS14a","CS14b","IT14","CS15a","CS15b","IT15","CS16a","CS16b","POST"))
mydata9<-mydatan[which(nObsString %in% c(PRE9:POST9)),]
nObsTrial<-length(mydata9$Trial)/samplen
nTrialM<-rep(c(1:nObsTrial),samplen)
mydata9<-cbind(mydata9,nTrialM)

mydatan<-mydata
mydatan$StimPeriod<-factor(mydatan$StimPeriod,labels = c("BL","CS1a","CS1b","IT1","CS2a","CS2b","IT2","CS3a","CS3b","IT3","CS4a","CS4b","IT4","CS5a","CS5b","IT5","CS6a","CS6b","IT6","CS7a","CS7b","IT7","CS8a","CS8b","IT8","CS9a","CS9b","preCS","Tone","Noise","postCS","CS11a","CS11b","IT11","CS12a","CS12b","IT12","CS13a","CS13b","IT13","CS14a","CS14b","IT14","CS15a","CS15b","IT15","CS16a","CS16b","POST"))
mydata10<-mydatan[which(nObsString %in% c(PRE10:POST10)),]
nObsTrial<-length(mydata10$Trial)/samplen
nTrialM<-rep(c(1:nObsTrial),samplen)
mydata10<-cbind(mydata10,nTrialM)

mydatan<-mydata
mydatan$StimPeriod<-factor(mydatan$StimPeriod,labels = c("BL","CS1a","CS1b","IT1","CS2a","CS2b","IT2","CS3a","CS3b","IT3","CS4a","CS4b","IT4","CS5a","CS5b","IT5","CS6a","CS6b","IT6","CS7a","CS7b","IT7","CS8a","CS8b","IT8","CS9a","CS9b","IT9","CS10a","CS10b","preCS","Tone","Noise","postCS","CS12a","CS12b","IT12","CS13a","CS13b","IT13","CS14a","CS14b","IT14","CS15a","CS15b","IT15","CS16a","CS16b","POST"))
mydata11<-mydatan[which(nObsString %in% c(PRE11:POST11)),]
nObsTrial<-length(mydata11$Trial)/samplen
nTrialM<-rep(c(1:nObsTrial),samplen)
mydata11<-cbind(mydata11,nTrialM)

mydatan<-mydata
mydatan$StimPeriod<-factor(mydatan$StimPeriod,labels = c("BL","CS1a","CS1b","IT1","CS2a","CS2b","IT2","CS3a","CS3b","IT3","CS4a","CS4b","IT4","CS5a","CS5b","IT5","CS6a","CS6b","IT6","CS7a","CS7b","IT7","CS8a","CS8b","IT8","CS9a","CS9b","IT9","CS10a","CS10b","IT10","CS11a","CS11b","preCS","Tone","Noise","postCS","CS13a","CS13b","IT13","CS14a","CS14b","IT14","CS15a","CS15b","IT15","CS16a","CS16b","POST"))
mydata12<-mydatan[which(nObsString %in% c(PRE12:POST12)),]
nObsTrial<-length(mydata12$Trial)/samplen
nTrialM<-rep(c(1:nObsTrial),samplen)
mydata12<-cbind(mydata12,nTrialM)

mydatan<-mydata
mydatan$StimPeriod<-factor(mydatan$StimPeriod,labels = c("BL","CS1a","CS1b","IT1","CS2a","CS2b","IT2","CS3a","CS3b","IT3","CS4a","CS4b","IT4","CS5a","CS5b","IT5","CS6a","CS6b","IT6","CS7a","CS7b","IT7","CS8a","CS8b","IT8","CS9a","CS9b","IT9","CS10a","CS10b","IT10","CS11a","CS11b","IT11","CS12a","CS12b","preCS","Tone","Noise","postCS","CS14a","CS14b","IT14","CS15a","CS15b","IT15","CS16a","CS16b","POST"))
mydata13<-mydatan[which(nObsString %in% c(PRE13:POST13)),]
nObsTrial<-length(mydata13$Trial)/samplen
nTrialM<-rep(c(1:nObsTrial),samplen)
mydata13<-cbind(mydata13,nTrialM)

mydatan<-mydata
mydatan$StimPeriod<-factor(mydatan$StimPeriod,labels = c("BL","CS1a","CS1b","IT1","CS2a","CS2b","IT2","CS3a","CS3b","IT3","CS4a","CS4b","IT4","CS5a","CS5b","IT5","CS6a","CS6b","IT6","CS7a","CS7b","IT7","CS8a","CS8b","IT8","CS9a","CS9b","IT9","CS10a","CS10b","IT10","CS11a","CS11b","IT11","CS12a","CS12b","IT12","CS13a","CS13b","preCS","Tone","Noise","postCS","CS15a","CS15b","IT15","CS16a","CS16b","POST"))
mydata14<-mydatan[which(nObsString %in% c(PRE14:POST14)),]
nObsTrial<-length(mydata14$Trial)/samplen
nTrialM<-rep(c(1:nObsTrial),samplen)
mydata14<-cbind(mydata14,nTrialM)

mydatan<-mydata
mydatan$StimPeriod<-factor(mydatan$StimPeriod,labels = c("BL","CS1a","CS1b","IT1","CS2a","CS2b","IT2","CS3a","CS3b","IT3","CS4a","CS4b","IT4","CS5a","CS5b","IT5","CS6a","CS6b","IT6","CS7a","CS7b","IT7","CS8a","CS8b","IT8","CS9a","CS9b","IT9","CS10a","CS10b","IT10","CS11a","CS11b","IT11","CS12a","CS12b","IT12","CS13a","CS13b","IT13","CS14a","CS14b","preCS","Tone","Noise","postCS","CS16a","CS16b","POST"))
mydata15<-mydatan[which(nObsString %in% c(PRE15:POST15)),]
nObsTrial<-length(mydata15$Trial)/samplen
nTrialM<-rep(c(1:nObsTrial),samplen)
mydata15<-cbind(mydata15,nTrialM)

mydatan<-mydata
mydatan$StimPeriod<-factor(mydatan$StimPeriod,labels = c("BL","CS1a","CS1b","IT1","CS2a","CS2b","IT2","CS3a","CS3b","IT3","CS4a","CS4b","IT4","CS5a","CS5b","IT5","CS6a","CS6b","IT6","CS7a","CS7b","IT7","CS8a","CS8b","IT8","CS9a","CS9b","IT9","CS10a","CS10b","IT10","CS11a","CS11b","IT11","CS12a","CS12b","IT12","CS13a","CS13b","IT13","CS14a","CS14b","IT14","CS15a","CS15b","preCS","Tone","Noise","postCS"))
mydata16<-mydatan[which(nObsString %in% c(PRE16:POST16)),]
nObsTrial<-length(mydata16$Trial)/samplen
nTrialM<-rep(c(1:nObsTrial),samplen)
mydata16<-cbind(mydata16,nTrialM)
```


```{r}
##taken from Exp2 Day 3Ext, but did not split into 16 different data sets, going to try to find another way without doing that
#mydataall16<-rbind(mydata1,mydata2,mydata3,mydata4,mydata5,mydata6,mydata7,mydata8,mydata9,mydata10,mydata11,mydata12,mydata13,mydata14,mydata15,mydata16)
# Never mind, I like the flexibility of beign able to group different trials together. I'm going to do that. (above this chunk)
mydataall16<-rbind(mydata1,mydata2,mydata3,mydata4,mydata5,mydata6,mydata7,mydata8,mydata9,mydata10,mydata11,mydata12,mydata13,mydata14,mydata15,mydata16)
mydataall8_1<-rbind(mydata1,mydata2,mydata3,mydata4,mydata5,mydata6,mydata7,mydata8)
mydataall8_2<-rbind(mydata9,mydata10,mydata11,mydata12,mydata13,mydata14,mydata15,mydata16)
mydataall4_1<-rbind(mydata1,mydata2,mydata3,mydata4)
mydataall4_2<-rbind(mydata5,mydata6,mydata7,mydata8)
mydataall4_3<-rbind(mydata9,mydata10,mydata11,mydata12)
mydataall4_4<-rbind(mydata13,mydata14,mydata15,mydata16)
```
```{r}
mydataall16$StimPeriod<-as.character(mydataall16$StimPeriod)
mydataall16$StimPeriod[mydataall16$StimPeriod=='Tone'] <- 'preCS'
mydataall16$StimPeriod[mydataall16$StimPeriod=='Noise'] <- 'CS'
mydataall16$StimPeriod<-as.factor(mydataall16$StimPeriod)
mydataall16$StimPeriod = factor(mydataall16$StimPeriod, levels=c('preCS','CS','postCS'))
```
```{r}
ggplot(mydataall16, aes(nTrialM,Velocity,color=StimPeriod)) +
  geom_smooth(method="loess",span=.5,method.args = list(degree=2))+
  facet_grid(.~Group,labeller = labeller(.cols = label_parsed, .multi_line = FALSE))+
  theme_classic()+
  coord_cartesian(ylim=c(0,13))+
  theme(legend.position = "bottom")+
   ylab("Velocity (cm/s)")+
  labs(color="Stimulus Period")+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())
```
```{r}
#tiff("test.tiff", units="cm", width=36, height=14, res=300)

ggplot(mydataall16, aes(nTrialM,Velocity,color=StimPeriod))+# +linetype=StimPeriod, 
  geom_smooth(method="loess",span=.5,method.args = list(degree=2))+
  facet_grid(.~Group,labeller = labeller(.cols = label_parsed, .multi_line = FALSE))+
  theme_classic()+
  coord_cartesian(ylim=c(0,13))+
  theme(legend.position = "none")+
  theme(panel.spacing.x = unit(1, "lines"))+
  theme(text = element_text(size = 16))+    
  ylab("Velocity (cm/s)")+
  scale_x_continuous(name=element_blank(),limits = c(1,153),labels =c("-10s","0s","10s","20s","30s"),breaks=c(1,39,77,114,153))+
  scale_color_manual(values=c("black","black","black","black"))+
  geom_vline(xintercept = nTrialM[39],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nTrialM[77],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nTrialM[114],linetype="dashed", color = "grey")


#dev.off()
```
```{r}
#tiff("test.tiff", units="cm", width=36, height=14, res=300)

ggplot(mydataall16, aes(nTrialM,Velocity,color=StimPeriod))+# +linetype=StimPeriod, 
  geom_smooth(method="loess",span=.5,method.args = list(degree=2))+
  facet_grid(.~Group,labeller = labeller(.cols = label_parsed, .multi_line = FALSE))+
  theme_classic()+
  coord_cartesian(ylim=c(0,13))+
  theme(legend.position = "none")+
  theme(panel.spacing.x = unit(1, "lines"))+
  theme(text = element_text(size = 16))+   
  ylab("Velocity (cm/s)")+
  scale_x_continuous(name=element_blank(),limits = c(1,153),labels =c("0s","10s","20s"),breaks=c(39,77,114))+
  scale_color_manual(values=c("black","black","black","black"))+
  geom_vline(xintercept = nTrialM[39],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nTrialM[77],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nTrialM[114],linetype="dashed", color = "grey")


#dev.off()
```


```{r}
mydataall4_1$StimPeriod<-as.character(mydataall4_1$StimPeriod)
mydataall4_1$StimPeriod[mydataall4_1$StimPeriod=='Tone'] <- 'preCS'
mydataall4_1$StimPeriod[mydataall4_1$StimPeriod=='Noise'] <- 'CS'
mydataall4_1$StimPeriod<-as.factor(mydataall4_1$StimPeriod)
mydataall4_1$StimPeriod = factor(mydataall4_1$StimPeriod, levels=c('preCS','CS','postCS'))
```
```{r}
ggplot(mydataall4_1, aes(nTrialM,Velocity,color=StimPeriod)) +
  geom_smooth(method="loess",span=.5,method.args = list(degree=2))+
  facet_grid(.~Group,labeller = labeller(.cols = label_parsed, .multi_line = FALSE))+
  theme_classic()+
  coord_cartesian(ylim=c(0,17.5))+
  theme(legend.position = "bottom")+
  ylab("Velocity (cm/s)")+
  labs(color="Stimulus Period")+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())
```
```{r}
#tiff("test.tiff", units="cm", width=36, height=14, res=300)

ggplot(mydataall4_1, aes(nTrialM,Velocity,color=StimPeriod)) +
  geom_smooth(method="loess",span=.5,method.args = list(degree=2))+
  facet_grid(.~Group,labeller = labeller(.cols = label_parsed, .multi_line = FALSE))+
  theme_classic()+
  coord_cartesian(ylim=c(0,18))+
  theme(legend.position = "none")+
  theme(panel.spacing.x = unit(1, "lines"))+
  theme(text = element_text(size=16))+
  ylab("Velocity (cm/s)")+
  scale_x_continuous(name=element_blank(),limits = c(1,153),labels =c("-10s","0s","10s","20s","30s"),breaks=c(1,39,77,114,153))+
  scale_color_manual(values=c("black","black","black","black"))+
  geom_vline(xintercept = nTrialM[39],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nTrialM[77],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nTrialM[114],linetype="dashed", color = "grey")

#dev.off()
```

```{r}
#tiff("test.tiff", units="cm", width=36, height=14, res=300)

ggplot(mydataall4_1, aes(nTrialM,Velocity,color=StimPeriod)) +
  geom_smooth(method="loess",span=.5,method.args = list(degree=2))+
  facet_grid(.~Group,labeller = labeller(.cols = label_parsed, .multi_line = FALSE))+
  theme_classic()+
  coord_cartesian(ylim=c(0,18))+
  theme(legend.position = "none")+
  theme(panel.spacing.x = unit(1, "lines"))+
  theme(text = element_text(size=16))+
  ylab("Velocity (cm/s)")+
  scale_x_continuous(name=element_blank(),limits = c(1,153),labels =c("0s","10s","20s"),breaks=c(39,77,114))+
  scale_color_manual(values=c("black","black","black","black"))+
  geom_vline(xintercept = nTrialM[39],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nTrialM[77],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nTrialM[114],linetype="dashed", color = "grey")

#dev.off()
```


#OK so what I want is some way to compare amplitudes of velocity at different time bins for each curve.
#



```{r}
sum16<-summarySEwithin(data=mydataall16,measurevar="Velocity",betweenvars = "Group",withinvars = "nTrialM",idvar = "Animal",na.rm = F,conf.interval = .95)
head(sum16)
sum4<-summarySEwithin(data=mydataall4_1,measurevar="Velocity",betweenvars = "Group",withinvars = "nTrialM",idvar = "Animal",na.rm = F,conf.interval = .95)
```
```{r}
sum16$nTrialM<-as.numeric(sum16$nTrialM)
ggplot(sum16, aes(x=nTrialM, y=Velocity, group=Group)) +
    geom_line() +
    geom_errorbar(width=.1, aes(ymin=Velocity-se, ymax=Velocity+se)) +
   facet_grid(.~Group,labeller = labeller(.cols = label_parsed, .multi_line = FALSE))+
  theme_classic()+
  coord_cartesian(ylim=c(0,20))+
  theme(legend.position = "bottom")+
  theme(panel.spacing.x = unit(1, "lines"))+
   ylab("Velocity (cm/s)")+
  scale_x_continuous(name="Time (s)",limits = c(1,153),labels =c("-10","0","10","20","30"),breaks=c(1,39,77,114,153))+
  #scale_linetype_manual(values=c("dotted","longdash","solid","twodash"))+
  #scale_color_manual(values=c("black","black","black","black"))+
  geom_vline(xintercept = nTrialM[39],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nTrialM[77],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nTrialM[114],linetype="dashed", color = "grey")


```
```{r}
sum4$nTrialM<-as.numeric(sum4$nTrialM)
ggplot(sum4, aes(x=nTrialM, y=Velocity, group=Group)) +
    geom_line() +
    geom_errorbar(width=.1, aes(ymin=Velocity-se, ymax=Velocity+se)) +
   facet_grid(.~Group,labeller = labeller(.cols = label_parsed, .multi_line = FALSE))+
  theme_classic()+
  coord_cartesian(ylim=c(0,26))+
  theme(legend.position = "bottom")+
  theme(panel.spacing.x = unit(1, "lines"))+
   ylab("Velocity (cm/s)")+
  scale_x_continuous(name="Time (s)",limits = c(1,153),labels =c("-10","0","10","20","30"),breaks=c(1,39,77,114,153))+
  #scale_linetype_manual(values=c("dotted","longdash","solid","twodash"))+
  #scale_color_manual(values=c("black","black","black","black"))+
  geom_vline(xintercept = nTrialM[39],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nTrialM[77],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nTrialM[114],linetype="dashed", color = "grey")


```



#now to try to get ANOVA'S/mixed models for this data
#https://www.datanovia.com/en/lessons/mixed-anova-in-r/
```{r}
library(tidyverse)
library(ggpubr)
library(rstatix)
library(WRS2)
```
```{r}
mydatall16agg<-aggregate(mydataall16$Velocity,by=list(mydataall16$Animal,mydataall16$Group,mydataall16$nTrialM),FUN="mean")
names(mydatall16agg)<-c("Animal","Group","nTrialM","Velocity")
library(xlsx)
#write.xlsx(mydatall16agg, file = "mytest.xlsx", 
         #  sheetName="all16agg", append=F)
```

##This below is a way to collapse every two rows, binning into .533s blocks
```{r}
sortmydatall16agg<-mydatall16agg[order(mydatall16agg$Animal),]
sortmydatall16agg<-sortmydatall16agg %>% filter(row_number() %% 153 != 0)

df_summary <- sortmydatall16agg %>% as_tibble(rownames = "names") %>% 
      group_by(ceiling(1:n() / 2)) %>% 
      summarise(names = paste(names, collapse = "_"),
                Animal=Animal,
                Group=Group,
                Velocity = mean(Velocity)) %>% 
      select(-1)

toDelete <- seq(1, nrow(df_summary), 2)
sortmydatall16agg<-df_summary[ toDelete ,]
sortmydatall16agg$names<-NULL
names(sortmydatall16agg)<-c("nTrialM","Animal","Group","Velocity")

sortmydatall16agg$nTrialM<-rep(seq(1,76,1),times=32)
```
```{r}
sum16bin<-summarySEwithin(data=sortmydatall16agg,measurevar="Velocity",betweenvars = "Group",withinvars = "nTrialM",idvar = "Animal",na.rm = F,conf.interval = .95)
```
```{r}
#tiff("test.tiff", units="cm", width=36, height=14, res=300)

sum16bin$nTrialM<-as.numeric(sum16bin$nTrialM)
ggplot(sum16bin, aes(x=nTrialM, y=Velocity, group=Group)) +
  geom_line() +
  geom_errorbar(width=.1, aes(ymin=Velocity-se, ymax=Velocity+se)) +
  facet_grid(.~Group,labeller = labeller(.cols = label_parsed, .multi_line = FALSE))+
  theme_classic()+
  theme(text = element_text(size=16))+
  coord_cartesian(ylim=c(0,7.5))+
  theme(legend.position = "bottom")+
  theme(panel.spacing.x = unit(1, "lines"))+
  ylab("Velocity (cm/s)")+
  scale_x_continuous(name=element_blank(),limits = c(1,76),labels =c("-10s","0s","10s","20s","30s"),breaks=c(1,20,39,58,76))+
  geom_vline(xintercept = nTrialM[20],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nTrialM[39],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nTrialM[58],linetype="dashed", color = "grey")

#dev.off()
```
```{r}
#tiff("test.tiff", units="cm", width=36, height=14, res=300)

sum16bin$nTrialM<-as.numeric(sum16bin$nTrialM)
ggplot(sum16bin, aes(x=nTrialM, y=Velocity, group=Group)) +
  geom_line() +
  geom_errorbar(width=.1, aes(ymin=Velocity-se, ymax=Velocity+se)) +
  facet_grid(.~Group,labeller = labeller(.cols = label_parsed, .multi_line = FALSE))+
  theme_classic()+
  theme(text = element_text(size=16))+
  coord_cartesian(ylim=c(0,7.5))+
  theme(legend.position = "bottom")+
  theme(panel.spacing.x = unit(1, "lines"))+
  ylab("Velocity (cm/s)")+
  scale_x_continuous(name=element_blank(),limits = c(1,76),labels =c("0s","10s","20s"),breaks=c(20,39,58))+
  geom_vline(xintercept = nTrialM[20],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nTrialM[39],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nTrialM[58],linetype="dashed", color = "grey")

#dev.off()
```


```{r}
mydatall4_1agg<-aggregate(mydataall4_1$Velocity,by=list(mydataall4_1$Animal,mydataall4_1$Group,mydataall4_1$nTrialM),FUN="mean")
names(mydatall4_1agg)<-c("Animal","Group","nTrialM","Velocity")
#library(xlsx)
#write.xlsx(mydatall4_1agg, file = "mytest.xlsx", 
           #sheetName="all4_1agg", append=F)
```

##This below is a way to collapse every two rows, binning into .533s blocks
```{r}
sortmydatall4_1agg<-mydatall4_1agg[order(mydatall4_1agg$Animal),]
sortmydatall4_1agg<-sortmydatall4_1agg %>% filter(row_number() %% 153 != 0)

df_summary4 <- sortmydatall4_1agg %>% as_tibble(rownames = "names") %>% 
      group_by(ceiling(1:n() / 2)) %>% 
      summarise(names = paste(names, collapse = "_"),
                Animal=Animal,
                Group=Group,
                Velocity = mean(Velocity)) %>% 
      select(-1)

toDelete <- seq(1, nrow(df_summary4), 2)
sortmydatall4_1agg<-df_summary4[ toDelete ,]
sortmydatall4_1agg$names<-NULL
names(sortmydatall4_1agg)<-c("nTrialM","Animal","Group","Velocity")

sortmydatall4_1agg$nTrialM<-rep(seq(1,76,1),times=32)
```
```{r}
sum4_1bin<-summarySEwithin(data=sortmydatall4_1agg,measurevar="Velocity",betweenvars = "Group",withinvars = "nTrialM",idvar = "Animal",na.rm = F,conf.interval = .95)
```
```{r}
#tiff("test.tiff", units="cm", width=36, height=14, res=300)

sum4_1bin$nTrialM<-as.numeric(sum4_1bin$nTrialM)
ggplot(sum4_1bin, aes(x=nTrialM, y=Velocity, group=Group)) +
  geom_line() +
  geom_errorbar(width=.1, aes(ymin=Velocity-se, ymax=Velocity+se)) +
  facet_grid(.~Group,labeller = labeller(.cols = label_parsed, .multi_line = FALSE))+
  theme_classic()+
  theme(text = element_text(size=16))+
  coord_cartesian(ylim=c(0,18))+
  theme(legend.position = "bottom")+
  theme(panel.spacing.x = unit(1, "lines"))+
  ylab("Velocity (cm/s)")+
  scale_x_continuous(name=element_blank(),limits = c(1,76),labels =c("-10s","0s","10s","20s","30s"),breaks=c(1,20,39,58,76))+
  geom_vline(xintercept = nTrialM[20],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nTrialM[39],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nTrialM[58],linetype="dashed", color = "grey")

#dev.off()
```
```{r}
#tiff("test.tiff", units="cm", width=36, height=14, res=300)

sum4_1bin$nTrialM<-as.numeric(sum4_1bin$nTrialM)
ggplot(sum4_1bin, aes(x=nTrialM, y=Velocity, group=Group)) +
  geom_line() +
  geom_errorbar(width=.1, aes(ymin=Velocity-se, ymax=Velocity+se)) +
  facet_grid(.~Group,labeller = labeller(.cols = label_parsed, .multi_line = FALSE))+
  theme_classic()+
  theme(text = element_text(size=16))+
  coord_cartesian(ylim=c(0,18))+
  theme(legend.position = "bottom")+
  theme(panel.spacing.x = unit(1, "lines"))+
  ylab("Velocity (cm/s)")+
  scale_x_continuous(name=element_blank(),limits = c(1,76),labels =c("0s","10s","20s"),breaks=c(20,39,58))+
  geom_vline(xintercept = nTrialM[20],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nTrialM[39],linetype="dashed", color = "grey")+
  geom_vline(xintercept = nTrialM[58],linetype="dashed", color = "grey")

#dev.off()
```

```{r}
sortmydatall16agg %>%
  group_by(nTrialM, Group) %>%
  get_summary_stats(Velocity, type = "mean_sd")
```

```{r}
bxp <- ggboxplot(
  sortmydatall16agg, x = "nTrialM", y = "Velocity",
  color = "Group", palette = "jco"
  )
bxp

```
```{r}
sortmydatall16agg %>%
  group_by(nTrialM, Group) %>%
  identify_outliers(Velocity)

```
##so there are tons of outliers. I probably need to use WRS2 package to analyze which is robust to outliers.
```{r}
sortmydatall16agg %>%
  group_by(nTrialM, Group) %>%
  shapiro_test(Velocity)

```
#just about nothing is normal/everything is normal when collapsed across all trials within an animal
```{r}
ggqqplot(mydatall4_1agg, "Velocity", ggtheme = theme_bw()) +
  facet_grid(nTrialM ~ Group)
```
##honestly too many cells here
#I separately want to think about how to bin this into larger bins, like 1 second bins. Maybe just combine two bins to get .5s bins
```{r}
sortmydatall16agg %>%
  group_by(nTrialM) %>%
  levene_test(Velocity~Group)

```
#does not pass homogeneity of variance by any standard/at any point, but I may want to bin these together. Does pass when all 16 trials collapsed within animal
```{r}
box_m(sortmydatall16agg[, "Velocity", drop = FALSE], sortmydatall16agg$Group)
```
#does not pass homogeneity of covariances by Box's test of equality of covariance matrices collapsed or not on all 16 trials

#need to ungroup these binned data to get R to run ANOVA
```{r}
sortmydatall16agg2<-ungroup(sortmydatall16agg)
sortmydatall4_1agg2<-ungroup(sortmydatall4_1agg)
sortmydatall16agg3<-sortmydatall16agg2[which(sortmydatall16agg$nTrialM > 38 & sortmydatall16agg$nTrialM<59), ]
sortmydatall4_1agg3<-sortmydatall4_1agg2[which(sortmydatall4_1agg$nTrialM > 38 & sortmydatall4_1agg$nTrialM<59), ]

```

# Two-way mixed ANOVA test
```{r}
res.aov <- anova_test(
  data = sortmydatall4_1agg3, dv = Velocity, wid=Animal,
  between = Group, within = nTrialM
  )
get_anova_table(res.aov)
#res.aov
```
#post-hoc analysis
```{r}
# Effect of group at each time point
one.way <- sortmydatall4_1agg2 %>%
  group_by(nTrialM) %>%
  anova_test(dv = Velocity, wid = Animal, between = Group) %>%
  get_anova_table() %>%
  adjust_pvalue(method = "bonferroni")
one.way
```
```{r}
# Pairwise comparisons between group levels
pwc <- sortmydatall4_1agg3 %>%
  group_by(nTrialM) %>%
  pairwise_t_test(Velocity ~ Group, p.adjust.method = "bonferroni")
pwc
```
```{r}
# Effect of time at each level of group
one.way2 <- mydatall16agg %>%
  group_by(Group) %>%
  anova_test(dv = Velocity, wid = Animal, within = nTrialM) %>%
  get_anova_table() %>%
  adjust_pvalue(method = "bonferroni")
one.way2
```
```{r}
# Pairwise comparisons between time points at each group levels
# Paired t-test is used because we have repeated measures by time
##too long/too many compariosn with this data-set
pwc2 <- mydatall16agg %>%
  group_by(Group) %>%
  pairwise_t_test(
    Velocity ~ nTrialM, paired = TRUE, 
    p.adjust.method = "bonferroni"
    ) %>%
  select(-df, -statistic, -p) # Remove details
pwc2
```

#If non-signif interaction, analyze any main effects:
```{r}
timemain<-mydatall16agg %>%
  pairwise_t_test(
    Velocity ~ nTrialM, paired = TRUE, 
    p.adjust.method = "bonferroni"
  )
```

```{r}
groupmain<-mydatall16agg %>%
  pairwise_t_test(
    Velocity ~ Group, 
    p.adjust.method = "bonferroni"
  )

```