## Replication material for:
## Permutation tests for regression, ANOVA and comparison of signals :
## the permuco package
## December 2019

suppressWarnings(RNGversion("3.5.0"))

##=========================================
## 6. Tutorial for the permuco package

install.packages("permuco")
library("permuco")
##=========================================


##=========================================
## 6.1. Fixed effects model

emergencycost$LOSc <- scale(emergencycost$LOS, scale = FALSE)

set.seed(42)
mod_cost_0 <- aovperm(cost ~ LOSc*sex*insurance, 
                      data = emergencycost, np = 100000)
mod_cost_0
##=========================================
emergencycost$LOS14 <- emergencycost$LOS - 14
mod_cost_14 <- aovperm(cost ~ LOS14*sex*insurance, 
                       data = emergencycost, 
                       np = 100000)
mod_cost_14
##=========================================

contrasts(emergencycost$insurance) <- contr.sum
contrasts(emergencycost$sex) <- contr.sum

contrasts(emergencycost$sex)

modlm_cost_14 <- lmperm(cost ~ LOS14*sex*insurance,
                        data = emergencycost, 
                        np = 100000)
modlm_cost_14
##=========================================

contrasts(emergencycost$insurance) <- contr.treatment
emergencycost$insurance <- relevel(emergencycost$insurance,
                                   ref = "public")

contrasts(emergencycost$insurance)

contrasts(emergencycost$sex) <- contr.sum
contrasts(emergencycost$sex)

mod_cost_se <- aovperm(cost ~ LOSc*sex*insurance, 
                       data = emergencycost,
                       np = 100000, coding_sum = FALSE)
mod_cost_se
##=========================================


##=========================================
## 6.2. Repeated measures ANCOVA
jpah2016$bmic <- scale(jpah2016$bmi, scale = FALSE)

mod_jpah2016 <- aovperm(iapa ~ bmic*condition*time + Error(id/(time)), 
                        data = jpah2016, method = "Rd_kheradPajouh_renaud")
mod_jpah2016

par(pin = c(1.5,1.5))
plot(mod_jpah2016, effect = c("bmic", "condition", "bmic:condition"))
##=========================================

##=========================================
## 6.3. EEG experiment in attention shifting
electrod_O1 <-
  clusterlm(attentionshifting_signal ~ visibility*emotion*direction
            + Error(id/(visibility*emotion*direction)),
            data = attentionshifting_design)

plot(electrod_O1)

summary(electrod_O1)$visibility


full_electrod_O1 <-
  clusterlm(attentionshifting_signal ~ visibility*emotion*direction
            + Error(id/(visibility*emotion*direction)),
            data = attentionshifting_design, P = electrod_O1$P,
            method = "Rde_kheradPajouh_renaud",
            multcomp = c("troendle", "tfce", "clustermass",
                         "bonferroni", "holm", "benjaminin_hochberg"))
plot(full_electrod_O1, multcomp = "tfce", enhanced_stat = TRUE)

plot(full_electrod_O1, multcomp = "troendle")

summary(full_electrod_O1, multcomp = "troendle")$visibility
##=========================================

##=========================================
## A.1. ANOVA and ANCOVA

library("permuco")
library("lmPerm")
library("flip")
library("GFD")
set.seed(42)

emergencycost$LOSc <- scale(emergencycost$LOS, scale = FALSE)
contrasts(emergencycost$sex) <- contr.sum 
contrasts(emergencycost$insurance) <- contr.sum 

X <- model.matrix( ~ sex+insurance, data = emergencycost)[, -1]
colnames(X) <- c("sex_num", "insurance_num")
emergencycost <- data.frame(emergencycost,X)

anova_permuco <- aovperm(cost ~ sex*insurance, data = emergencycost)
anova_GFD <- GFD(cost ~ sex*insurance, data = emergencycost,
                 CI.method = "perm", nperm  = 5000)

ancova_permuco <- aovperm(cost ~ LOSc*sex*insurance, data = emergencycost,
                          method = "huh_jhun")
ancova_flip <- flip(cost ~1, X = ~ sex_num, Z = ~ LOSc*insurance_num*sex_num
                    - sex_num, data = emergencycost, statTest = "ANOVA",
                    perms = 5000)
ancova_lmPerm <- aovp(cost ~ LOS*sex*insurance, data = emergencycost, 
                      seqs = FALSE, nCycle = 1)

anova_permuco
anova_GFD
ancova_permuco
ancova_flip
summary(ancova_lmPerm)
##=========================================


##=========================================
## A.2. Repeated Measures ANOVA

jpah2016$id <- as.factor(jpah2016$id)
jpah2016$bmic <- scale(jpah2016$bmi,scale = FALSE)
contrasts(jpah2016$time) <- contr.sum
contrasts(jpah2016$condition) <- contr.sum

rancova_permuco <- aovperm(iapa ~ bmic*condition*time + Error(id/(time)),
                           data = jpah2016)
rancova_lmPerm <- aovp(iapa ~ bmic*condition*time + Error(id/(time)),
                       data = jpah2016, nCycle = 1, seqs = FALSE)

rancova_permuco
summary(rancova_lmPerm)