
USES R X 64 4.0; DOES NOT WORK IN R 3.6

install.packages(c('reshape','car'))
library(car)
library(reshape)

#there is a handy library called schoRsch that gives you partial eta squared from ezANOVA
#http://finzi.psych.upenn.edu/library/schoRsch/html/anova_out.html


zz <- file.path("U:","My Documents")

File.Location <- file.path(zz,"rep_meas_formattednogroup.csv")
y <- read.csv(File.Location)
y <- data.frame(y)

# obtain mean and sd for y1 in each group

library(data.table)
#dt <- data.table(y)
#y$group = rep(1:2,each=9)
#dt[,list(mean=mean(y$y1),sd=sd(y$y1)),by=y$group]

subject <- c(1:18)

# create a data frame which contains only the response and subject indicator
isr <- data.frame(y,subject)

# group means for y1; note: could gave created group within R
y$group = rep(1:2,each=9)
tapply(y$y1,y$group,mean)
tapply(y$y1,y$group,sd)

# can look at distributions of repeated measures (normaliity assumptions about the distribution of differences)
# see later for  a method for assessing residuals
boxplot(y$y2-y$y1)
boxplot(y$y2-y$y1,subset=y$group==1)

# we will just put the response into long format so removing group from the data frame; will generate this again after we have put the response
# into long format

# we need to get rid of the between subjects factor which we will recreate later for the formatting to work

isr$group <- NULL

library(reshape)
mdata = melt(isr, id=c("subject"))
colnames(mdata) = c("subject", "var", "score")
mdata <- data.frame(mdata)
#mdata <- arrange(mdata,subject)
head(mdata)


# more generally could construct the long formatted group  variable from the group in the repmeas_formatted.csv input file
group2 <- rep(y$group,6)

mood <- rep(c(1:3),each=18)
mood <- rep(mood, 2)
time <- rep(c(1:2),each=54)


#simple effect comparing mood (Within subjects factor) in group1 (Between subjects factor)

# Paired t-tests for three mood groups using all data
mdata <- data.frame(mdata,mood,time)
attach(mdata)

#install.packages("ez")
library(ez)

# there are several ways to do mixed ANOVA in R; Field recommends ezANOVA which is in the ez library
# this produces the Greenhouse-Geisser corrected p-values

# We have to make sure that ezANOVA will recognize the group inputs as factors which we can do using recode
# otherwise we get incorrect dfs and also would not obtain the Greenhouse-Geisser output

# factor usually works at constructing factors but doesn't seem to work with ezANOVA

#mdata$mood2 <- factor(mdata$mood, labels = c('A','B','C'))
#mdata$time <- factor(mdata$time, labels=c('base','follow-up'))
#mdata$group2 <- factor(mdata$group, labels=c('controls','patients'))

#  so we use recode which does

mdata$mood2 <- recode(mdata$mood, "1='A'; 2='B'; 3='C'")
mdata$time <- recode(mdata$time, "1='base';2='follow-up'")
#mdata$group2 <- recode(mdata$group, "1='controls';2='patients'")
mdata$group2 <- recode(mdata$group2, "1='controls';2='patients'")


# can fit WW model
out <- ezANOVA(data=mdata, dv=.(score), wid=.(subject), within=.(mood2,time),detailed=T,type=3)

# REPEAT BUT WITH GROUP

install.packages(c('reshape','car'))
library(car)
library(reshape)

zz <- file.path("U:","My Documents")

File.Location <- file.path(zz,"rep_meas_formattednogroup.csv")
y <- read.csv(File.Location)
y <- data.frame(y)

# obtain mean and sd for y1 in each group

library(data.table)
dt <- data.table(y)
dt[,list(mean=mean(y1),sd=sd(y1)),by=group]

subject <- c(1:18)

# create a data frame which contains only the response and subject indicator
isr <- data.frame(y,subject)

library(reshape)
mdata = melt(isr, id=c("subject"))
colnames(mdata) = c("subject", "var", "score")
mdata <- data.frame(mdata)
#mdata <- arrange(mdata,subject)
head(mdata)


# more generally could construct the long formatted group  variable from the group in the repmeas_formatted.csv input file
group2 <- rep(c(1,2),each=9)

mood <- rep(c(1:3),each=18)
mood <- rep(mood, 2)
time <- rep(c(1:2),each=54)


#simple effect comparing mood (Within subjects factor) in group1 (Between subjects factor)

# Paired t-tests for three mood groups using all data
mdata <- data.frame(mdata,mood,time,group2)
attach(mdata)

#install.packages("ez")
library(ez)

# there are several ways to do mixed ANOVA in R; Field recommends ezANOVA which is in the ez library
# this produces the Greenhouse-Geisser corrected p-values

# We have to make sure that ezANOVA will recognize the group inputs as factors which we can do using recode
# otherwise we get incorrect dfs and also would not obtain the Greenhouse-Geisser output

# factor usually works at constructing factors but doesn't seem to work with ezANOVA

#mdata$mood2 <- factor(mdata$mood, labels = c('A','B','C'))
#mdata$time <- factor(mdata$time, labels=c('base','follow-up'))
#mdata$group2 <- factor(mdata$group, labels=c('controls','patients'))

#  so we use recode which does

mdata$mood2 <- recode(mdata$mood, "1='A'; 2='B'; 3='C'")
mdata$time <- recode(mdata$time, "1='base';2='follow-up'")
#mdata$group2 <- recode(mdata$group, "1='controls';2='patients'")
mdata$group2 <- recode(mdata$group2, "1='controls';2='patients'")

out <- ezANOVA(data=mdata, dv=.(score), wid=.(subject), within=.(mood2,time),between=.(group2),detailed=T,type=3)

# extracting residuals see https://stackoverflow.com/questions/30238213/ggplot2-residuals-with-ezanova

#install.packages("afex")
library(afex)
out <- ezANOVA(data=mdata, dv=.(score), wid=.(subject), within=.(mood2,time),between=.(group2),detailed=T,type=3)
#install.packages("afex")
library(afex)
out <- ezANOVA(data=mdata, dv=.(score), wid=.(subject), within=.(mood2,time),between=.(group2),detailed=T,type=3)
anova_2 <- aov_ez("subject","score",mdata,within=c("mood2","time"),between=c("group2"))
nice(anova_2)
residuals(anova_2)
#rstandard(anova_2)
#anova_1 <- ez.glm("subject","score",mdata,within=c("mood2","time"),between=c("group2"),return="lm")
#anova_2 <- ez.glm("subject","score",mdata,within=c("mood2","time"),between=c("group2"),return="full")
#nice.anova(anova_2$Anova)
#residuals(anova_1)
#rstandard(anova_1)


#install.packages("magrittr")
#install.packages("tidyr")
#install.packages("mlr3misc")
library(magrittr)
library(tidyr)
library(mlr3misc)
anova_1 <- ezANOVA(data=mdata, dv=.(score), wid=.(subject), within=.(mood2,time),between=.(group2),detailed=T,return_aov=T)
ezanova_residuals <- purrr::map(anova_1$aov, residuals)
ezanova_residuals[[2]]
hist(ezanova_residuals[[2]])
shapiro.test(ezanova_residuals[[2]]) 