#######################################################
# Load appropriate libraries 
#######################################################
library(dplyr)
library(lme4)
library(lmerTest)
library(multcomp)
library(rms)
library(nlme)
library(ggplot2)
library(viridis)
theme_set(theme_bw())

setwd('C:\\Users\\morrile2\\Documents\\MULTIS\\LayerEffectonStiffness')

#######################################################
# Functions
#######################################################

run_model_noRand <- function(equation, ex_data){
  ML.model <- lm(equation, data=ex_data)
  return(ML.model)
}

run_model <- function(equation, ex_data){
  ML.model <- lmer(equation, data=ex_data)
  return(ML.model)
}


#######################################################
# Load data
#######################################################
# Read in data 
ind.data <- read.csv("dat/001_MasterList_indentation_orig.csv", header=TRUE, sep=",", na.strings="NA", dec=".")

ind.data['norm_muscle_disp'] <- ind.data['muscle_disp']/ind.data['MaxForce']
ind.data['norm_fat_disp'] <- ind.data['fat_disp']/ind.data['MaxForce']

#######################################################
# Data summary plots
#######################################################
ggplot(ind.data, aes(x=Muscle, y=norm_muscle_disp, color=Location)) + geom_point(alpha=.5, size=1)+
  theme(text = element_text(size=6), legend.key.size = unit(0.3, "cm"))+
  labs(x=bquote("Muscle Thickness"), y=bquote("Normalized muscle displacement (mm/N)"))+
  guides(colour = guide_legend(override.aes = list(size=2)))
ggsave('doc/Figures/mdisp_mthick.png', height = 5, width = 8, units='cm')

ggplot(ind.data, aes(x=Muscle, y=norm_fat_disp, color=Location)) + geom_point(alpha=.5, size=1)+
  theme(text = element_text(size=6), legend.key.size = unit(0.3, "cm"))+
  labs(x=bquote("Muscle Thickness"), y=bquote("Normalized fat displacement (mm/N)"))+
  guides(colour = guide_legend(override.aes = list(size=2)))
ggsave('doc/Figures/fdisp_mthick.png', height = 5, width = 8, units='cm')

ggplot(ind.data, aes(x=Fat, y=norm_muscle_disp, color=Location)) + geom_point(alpha=.5, size=1)+
  theme(text = element_text(size=6), legend.key.size = unit(0.3, "cm"))+
  labs(x=bquote("Fat Thickness"), y=bquote("Normalized muscle displacement (mm/N)"))+
  guides(colour = guide_legend(override.aes = list(size=2)))
ggsave('doc/Figures/mdisp_fthick.png', height = 5, width = 8, units='cm')

ggplot(ind.data, aes(x=Fat, y=norm_fat_disp, color=Location)) + geom_point(alpha=.5, size=1)+
  theme(text = element_text(size=6), legend.key.size = unit(0.3, "cm"))+
  labs(x=bquote("Fat Thickness"), y=bquote("Normalized fat displacement (mm/N)"))+
  guides(colour = guide_legend(override.aes = list(size=2)))
ggsave('doc/Figures/fdisp_fthick.png', height = 5, width = 8, units='cm')


#######################################################
# Data summary plots - split by location
#######################################################

ggplot(ind.data, aes(x=Muscle, y=norm_muscle_disp, color=Location)) + geom_point(alpha=.5, size=1)+
  theme(text = element_text(size=6), legend.key.size = unit(0.3, "cm"))+
  facet_wrap(vars(Location), nrow=2, scales='free')+
  labs(x=bquote("Muscle Thickness"), y=bquote("Normalized muscle displacement (mm/N)"))+
  guides(colour = guide_legend(override.aes = list(size=2)))
ggsave('doc/Figures/mdisp_mthick_splitloc.png', height = 9, width = 16.5, units='cm')

ggplot(ind.data, aes(x=Fat, y=norm_fat_disp, color=Location)) + geom_point(alpha=.5, size=1)+
  theme(text = element_text(size=6), legend.key.size = unit(0.3, "cm"))+
  facet_wrap(vars(Location), nrow=2, scales='free')+
  labs(x=bquote("Fat Thickness"), y=bquote("Normalized fat displacement (mm/N)"))+
  guides(colour = guide_legend(override.aes = list(size=2)))
ggsave('doc/Figures/fdisp_fthick_splitloc.png', height = 9, width = 16.5, units='cm')

#######################################################
# Statistical Analysis
#######################################################

# Stats model fit to all locations together 
M1 = run_model(muscle_disp ~ Muscle + Fat + (1|SubID) + (1|Location), ind.data)
summary(M1)

M2 = run_model(fat_disp ~ Muscle + Fat + (1|SubID) + (1|Location), ind.data)
summary(M2)