#######################################################
# Load appropriate libraries 
#######################################################
library(lme4)
library(lmerTest)
library(multcomp)

#######################################################
# Functions
#######################################################
run_model <- function(equation, ex_data){
  ML.model <- lmer(equation, data=ex_data)
  ML.model
  
  # Check normality of residuals using qq plot and Shapiro-Wilk Test
  qqnorm(residuals(ML.model))
  qqline(residuals(ML.model))
  return(ML.model)
}

#######################################################
# Mixed-Effect model
#######################################################
# Read in data 
ind.data <- read.csv("/home/morrile2/Documents/Multis/studies/Indentation/dat/002_MasterList_indentation.csv", header=TRUE, sep=",", na.strings="NA", dec=".")

# Log transform of stiffness data
ind.data$Total_Stiff_log <- log(ind.data$Total_Stiff)

# Create mixed model with all variables
M1 = run_model(Total_Stiff_log ~ Location + Age + Gender + BMI + ActivityLevel + (1|SubID), ind.data)
summary(M1)
anova(M1)
summary(glht(M1, linfct=mcp(Location="Tukey")), test=adjusted("bonferroni"))
summary(glht(M1, linfct=mcp(ActivityLevel="Tukey")), test=adjusted("bonferroni"))

# Log transform of modulus data
ind.data$Total_EMod_log <- log(ind.data$Total_EMod)

# Create mixed model with all variables
M1 = run_model(Total_EMod_log ~ Location + Age + Gender + BMI + ActivityLevel + (1|SubID), ind.data)
summary(M1)
anova(M1)
summary(glht(M1, linfct=mcp(Location="Tukey")), test=adjusted("bonferroni"))
summary(glht(M1, linfct=mcp(ActivityLevel="Tukey")), test=adjusted("bonferroni"))