library(tidyverse)
library(ggplot2)
library(ordinal)
library(RVAideMemoire)
library(lme4)
library(lmerTest)
library(gridExtra)
library(emmeans)
library(glmmTMB)
library(DHARMa)
library(car)

setwd("/Users/louisetosetto/My Drive/Neuro_Lab/Damsels/Behaviour/DFTrial_Third/DFBehaviourStats")

dfbehave <- read.csv(file = "dfbehaviourdata.csv", header = TRUE)

#Then can subtract file name into new columns
dfbehave <- dfbehave%>%
  mutate(df = substr(file_name, 7, 9)) %>%
  mutate(grating= substr(file_name, 1, 4)) %>%
  mutate(orient= substr(file_name, 5, 5)) 

#create variable of time
dfbehave$time <-rep(c(1:20),times=60)

#clean data
dfbehave$df <- as.factor(dfbehave$df)
dfbehave$grating <- as.factor(dfbehave$grating)
dfbehave$orient <- as.factor(dfbehave$orient)
dfbehave$distance <- as.numeric(dfbehave$distance)
dfbehave$trial <- as.factor(dfbehave$trial)
dfbehave$time <- as.numeric(dfbehave$time)

#DISTANCE FROM GRATING

#Because there is an upper limit of distance from grating due to the tank size. These data can be transformmed 
#from 0 - 1  - becomes proportional data and we have a beta distribution
dfbehave$distance_norm <-dfbehave$distance/max(dfbehave$distance)
hist(dfbehave$distance_norm)

#to run model the y values must be 0 < y < 1
dfbehave$distance_norm[dfbehave$distance_norm==1] <- 0.9999
max(dfbehave$distance_norm)

#DISTANCE FROM GRATINGS 

#Run glmm for beta distributed for distance from grating over time
beta1 <- glmmTMB(distance_norm ~ grating * orient * time + (1|trial/df), data = dfbehave, family = beta_family())
Anova(beta1, type="III")

#Distance from Grating  - (drop interaction with time)
beta2 <- glmmTMB(distance_norm ~ grating * orient + time + (1|trial/df), data = dfbehave, family = beta_family())
Anova(beta2, type = "III")

#Run diagnostic plots (library DHARMa)
r = simulateResiduals(beta2, plot = TRUE)

#pairwise comparisons
emmeans(beta2, specs = pairwise ~ grating | orient) #within group comparisons type
emmeans(beta2, specs = pairwise ~ orient | grating)

#DISTANCE MOVED IN RESPONSE TO DIFFERENT GRATINGS

hist(dfbehave$travelsum) #left skewed
hist(log(dfbehave$travelsum)) #normal distribution

#Run lmm for distance moved every 15 seconds over time
m0 <- lmer(log(travelsum) ~ grating * orient * time + (1|trial/df), data = dfbehave)
Anova(m0, type = "III")

#Distance Moved - (drop interaction with time)
m1 <- lmer(log(travelsum) ~ grating * orient + time + (1|trial/df), data = dfbehave) 
Anova(m1, type="III") #assess main effects of no interaction

#check assumptions of model
hist(residuals(m1)) #check normality of residuals
plot(m1) #check linearity of residuals
qqnorm(residuals(m1)) 
qqline(residuals(m1))

#pairwise comparisons
emmeans(m1, specs = pairwise ~ grating | orient, type = "response")
emmeans(m1, specs = pairwise ~ orient | grating, type = "response")