#!/usr/bin/env Rscript3.5.1
rm(list=ls());
#----library
# core
library(tidyverse);
library(broom);
# analysis
library(lme4); library(lmerTest);
# plot
library(ggplot2);
# estimated marginal means
library(emmeans); library(multcomp);
    
#----repeated measures anova model
anova_ <- function(){
    #----load data
    path <- "/Users/mdl-admin/Desktop/mdl-R33-analysis/output/analysis/html/model/anova/csv/"
    df <- read.csv(file=file.path(path, 'gaze_bias.csv'), header=TRUE)
    
    #----set type
    df$trialType <- factor(df$trialType)
    df$cesd_group <- factor(df$cesd_group)
    df$participant <- factor(df$participant)
    
    #----run model
    lmer_ <- lmerTest::lmer(gaze_bias ~ cesd_group + trialType + (1|participant), data=df)
    aov_ <- stats::anova(lmer_)
    
    #----Estimated Marginal Means (Least-Squares Means) of all factor levels
    pair_ <- emmeans::emmeans(lmer_, c('cesd_group','trialType'),
                              type = "response", adjust = "tukey")
    # Comparison of Differences Between Levels of Factor
    lsmeans_ <- summary(emmeans::as.glht(pairs(pair_)), test=adjusted("free")) %>% 
       broom::tidy()
    
    #----getting x,y coordinates for qqplot
    # create plot
    gg <- ggplot(lmer_) + 
      ggplot2::stat_qq(aes(sample = .resid, colour = factor(trialType))) + 
      ggplot2::geom_abline(linetype = "dotted")
    
    # convert to tibble
    gg <- ggplot_build(gg)[["data"]][[1]] %>% 
      dplyr::select(sample, theoretical)
    
    #----creating residuals tibble from model
    # including raw data, residuals vs fitted)
    # .resid=residuals, .fitted=predicted values .estimate=estimate of fixed effect
    residuals <- broom::augment(lmer_) %>% 
       dplyr::select(participant, gaze_bias, cesd_group, trialType, .resid, .fitted)
    
    #----merge residuals and qq tibble
    residuals <- as.data.frame(merge(x=residuals, y=gg, by.x=".resid", by.y="sample", all.x = TRUE))
    
    #----prepare data for export
    # convert model output to tibble
    output <- broom::tidy(aov_)
    # convert model summary to tibble
    summary <- broom::glance(lmer_)
    
    #----return model output, model, residuals, lsmeans, model summary
    return(list(output, lmer_, residuals, lsmeans_, summary))
}