#!/usr/bin/env Rscript3.5.1
rm(list=ls());
#----library
# core
library(tidyverse);
library(broom); library(broom.mixed);
# analysis
library(lme4); library(lmerTest);
# plot
library(ggplot2);
# estimated marginal means
library(emmeans); library(multcomp); #----linear mixed model
lmer_ <- function(){ #----load data path <- "/Users/mdl-admin/Desktop/mdl-R33-analysis/output/analysis/html/model/lmer/csv/" df <- read.csv(file=file.path(path, 'onset_data.csv'), header=TRUE) #----drop data # os #samples too small df <- df[(!df$os=="cos"),] #----normalize trial to [0,1] (recommended by Jason) df$TrialNum <- lapply(df$TrialNum, function(x){((x - 0)/(197 - 0))}) #----set type # set as factor df$os <- factor(df$os) df$trialType <- factor(df$trialType) df$participant <- factor(df$participant) # set trial as numeric (recommended by Jason) df$TrialNum <- as.numeric(df$TrialNum) #----run model model <- lmerTest::lmer(sqrt(diff_stim) ~ os + trialType + TrialNum + (1+TrialNum|participant), data=df) #----getting x,y coordinates for qqplot # create plot gg <- ggplot(model) + ggplot2::stat_qq(aes(sample = .resid, colour = factor(trialType))) + ggplot2::geom_abline(linetype = "dotted") + theme_bw() # convert to tibble gg <- ggplot_build(gg)[["data"]][[1]] %>% dplyr::select(sample, theoretical) #----creating residuals tibbles from model # including raw data, residuals vs fitted # .resid=residuals, .fitted=predicted values .estimate=estimate of fixed effect residuals <- broom::augment(model) %>% dplyr::select(participant, sqrt.diff_stim., os, trialType, TrialNum, .resid, .fitted) #----merge residuals and qq tibbles residuals <- merge(x=residuals, y=gg, by.x=".resid", by.y="sample", all.x = TRUE) #----prepare data for export # convert output to tibble output <- broom.mixed::tidy(model) # convert summary to tibble summary <- broom.mixed::glance(model) #----return model output, model, residuals, model summary return(list(output, model, residuals, summary))
}