#-----------------------------------------------------------------------------#
# Title: The French press: A repeatable and high-throughput approach to 
# exercising zebrafish
# Authors: Usui, T. et al.
# Description: Script for data analysis
# Date created: 25/04/2017
# Date updated: 11/11/2017
#-----------------------------------------------------------------------------#


# 1. Packages and data loading
#-----------------------------------------------------------------------------#

# Remove any old objects from the R environment
rm(list=ls())

# Load packages
library(car)
library(lme4)
library(lmerTest)
library(rptR)

# Read in data object.
MaxSpeedObj <- readRDS("./DataObjects.rds")


# 2. Data prep and checks
#-----------------------------------------------------------------------------#

# Prep data type
VarNum <- c("speed", "length", "weight", "temp")
VarFac <- c("sex", "treatment", "measure", "fish_ID")
MaxSpeedObj[VarFac] <- lapply(MaxSpeedObj[VarFac], factor)
MaxSpeedObj[VarNum] <- lapply(MaxSpeedObj[VarNum], as.character)
MaxSpeedObj[VarNum] <- lapply(MaxSpeedObj[VarNum], as.numeric)

# Mean-centering
MaxSpeedObj$temp.centred <- scale(MaxSpeedObj$temp, center=TRUE, scale=FALSE)
MaxSpeedObj$length.centred <- scale(MaxSpeedObj$length, center=TRUE, scale=FALSE)


# 3. LMER for investigating maximum speed and weight
#-----------------------------------------------------------------------------#

# 3.1 lmer model for speed
SpLmer <- lmer(speed ~ measure*treatment + sex + temp.centred + length.centred + (1|fish_ID), data=MaxSpeedObj, REML=T)

# Contrast model with sham-control group as reference
MaxSpeedObj$treatment <- relevel(MaxSpeedObj$treatment, ref = "sham")
SpLmerCont <- lmer(speed ~ measure*treatment + sex + temp.centred + length.centred + (1|fish_ID), data=MaxSpeedObj, REML=T)
  
# 3.2 lmer model for weight
WeLmer <- lmer(weight  ~ measure*treatment + sex + (1 | fish_ID), data=MaxSpeedObj, REML=T)


# 4. Repeatability analysis
#-----------------------------------------------------------------------------#

# 4.1 Unadjusted repeatability
RepUnAj <- rpt(speed ~ (1|fish_ID), grname="fish_ID", data=MaxSpeedObj, datatype="Gaussian", nboot=10000, npermut=10000, parallel=TRUE)

# 4.2 Adjusted repeatability
RepAj <- rpt(speed ~ sex + length + measure*treatment + (1|fish_ID), grname=c("fish_ID"), data=MaxSpeedObj, datatype="Gaussian", nboot=10000, npermut=10000, parallel=TRUE)
  
# 4.3 Enhanced agreement repeatability  
RepEnAj <- rpt(speed ~ sex + length + measure*treatment + (1|fish_ID), grname=c("fish_ID"), data=MaxSpeedObj, datatype="Gaussian", nboot=10000, npermut=10000, parallel=TRUE, adjusted=FALSE)