# Κεφάλαιο 15
# Υπολογισμός της διπαραγοντικής ανάλυσης διακύμανσης ανεξάρτητων δειγμάτων στο R

# Φορτώνουμε τα ακόλουθα πακέτα
library(foreign)
library(car)
library(ez)
library(lsr)
library(heplots)
library(multcomp)
library(effectsize)

# Διαβάζουμε τα δεδομένα του ερευνητικού παραδείγματος από το αρχείο και τα αποδίδουμε στο πλαίσιο δεδομένων ch15.1
ch15.1<-read.spss('chapter15_1.sav', to.data.frame=T)

# Βασικός έλεγχος του πλαισίου δεδομένων
str(ch15.1)
summary(ch15.1)

# Αλλαγή ονομασίας των επιπέδων για τη γραφική απεικόνιση
levels(ch15.1$instructions) <- c("Σ Ελέγχου","Πειραματική Σ 1","Πειραματική Σ 2")

# Διάγραμμα αλληλεπίδρασης των δύο μεταβλητών
with(ch15.1,interaction.plot(instructions,age,mistakes))
# Πιο "φροντισμένο"…
with(ch15.1,interaction.plot(instructions,age,mistakes, type="b",pch=c(21,21),lty=c(1,1),lwd=c(2,2), col=c("red","blue"),cex.axis=0.8,las=T,ylim=c(0,45), ylab="Αριθμός λαθών",xlab="",trace.label="",leg.bty="o"))

# Έλεγχος ομοιογένειας διακύμανσης κατά Levene προσαρμοσμένος στο SPSS
leveneTest(mistakes~instructions*age,ch15.1,center=mean)

# Ανάλυση διακύμανσης
a1 <- aov(mistakes~instructions*age,ch15.1)
summary(a1)

# Εναλλακτικός τύπος που δίνει ίδια αποτελέσματα με το SPSS και σε άνισες ομάδες
Anova(lm(mistakes~instructions*age,ch15.1),type=c("III"))

# Έλεγχος προσαρμογής
par(mfrow=c(2,2)) # θα φανούν και τα τέσσερα διαγράμματα μαζί
plot(a1)
par(mfrow=c(1,1)) 

# Έλεγχος υπολοίπων
shapiro.test(a1$residuals)

# Μέγεθος επίδρασης (η2, η2p, ω2, ω2p) από διάφορες βιβλιοθήκες
etasq(a1, partial=F) 		# υπολογίζεται η τιμή η2
etasq(a1) 			# υπολογίζεται η τιμή η2p
etaSquared(a1)
omega_squared(a1, partial=F) 	# υπολογίζεται η τιμή ω2
omega_squared(a1) 		# υπολογίζεται η τιμή ω2p

# Εναλλακτικά, ανάλυση διακύμανσης με έλεγχο ομοιογένειας και διάγραμμα αλληλεπίδρασης με διαστήματα εμπιστοσύνης από τη βιβλιοθήκη ez
ch15.1$subject <- paste("s",rownames(ch15.1),sep="") 	# χρειάζεται μεταβλητή ομαδοποίησης
ezANOVA(data=ch15.1,dv=.(mistakes),between=.(instructions,age),wid=.(subject))
ezPlot(data=ch15.1,dv=.(mistakes),between=.(instructions,age),wid=.(subject),x=.(instructions),split=.(age))

# Ανάλυση απλών επιδράσεων
s1 <- summary(a1)[[1]]
#  :κρατάμε τον παρονομαστή από την παραγοντική ανάλυση
MS2 <- s1["Residuals","Mean Sq"]
df2 <- s1["Residuals","Df"]
# Αναλύουμε την επίδραση της μιας μεταβλητής σε κάθε επίπεδο της άλλης
for (i in levels(ch15.1$instructions)) {
	a1i <- aov(mistakes~age,ch15.1,subset=(instructions==i)) # oneway ANOVA
	s1i <- summary(a1i)[[1]]
# :χρησιμοποιούμε τον αριθμητή από την απλή ανάλυση
MS1 <- s1i["age","Mean Sq"]
df1 <- s1i["age","Df"]
F <- MS1/MS2
p <- 1.0-pf(F,df1,df2)
cat(sprintf("%15s : F = %6.3f  df1=%d  df2=%d  p=%.3f",i,F,df1,df2,p),fill=T)}

# Επαναλαμβάνουμε για την άλλη μεταβλητή
for (a in levels(ch15.1$age)) {
a1a <- aov(mistakes~instructions,ch15.1,subset=(age==a)) # oneway ANOVA
s1a <- summary(a1a)[[1]]
# :χρησιμοποιούμε τον αριθμητή από την απλή ανάλυση
MS1 <- s1a["instructions","Mean Sq"]
df1 <- s1a["instructions","Df"]
F <- MS1/MS2
p <- 1.0-pf(F,df1,df2)
cat(sprintf("%15s : F = %6.3f df1=%d df2=%d p=%.3f",a,F,df1,df2,p),fill=T)}

# Εκ των υστέρων πολλαπλές συγκρίσεις 
TukeyHSD(a1)
TukeyHSD(a1,"instructions")

# Έλεγχος αντιθετικών γραμμικών συνδυασμών με τη βιβλιοθήκη multcomp
# Πίνακας αντιθετικών συνδυασμών Helmert
K <- rbind("Level 1 vs later  " =c(1,-0.5,-0.5),"Level 2 vs Level 3" =c(0,1,-1))

# Στήσιμο του πίνακα αντιθέσεων βάσει των οδηγιών και παραδειγμάτων του multcomp
tmp <- expand.grid(instructions=unique(ch15.1$instructions),age=unique(ch15.1$age))
X <- model.matrix(~instructions*age,data=tmp)
KK <- cbind(K,K)/2
summary(glht(a1, linfct = KK %*% X),test=adjusted(type="none"))

# Υπολογισμός της διπαραγοντικής ανάλυσης διακύμανσης εξαρτημένων δειγμάτων στο R

# Φορτώνουμε τα ακόλουθα πακέτα
library(foreign)

# Διαβάζουμε τα δεδομένα του ερευνητικού παραδείγματος από το αρχείο και τα αποδίδουμε στο πλαίσιο δεδομένων ch15.2w
ch15.2w<-read.spss('chapter15_2.sav', to.data.frame=T)

# Βασικός έλεγχος του πλαισίου δεδομένων
str(ch15.2w)
summary(ch15.2w)

# Μετατροπή των δεδομένων σε μορφή long (μια σειρά ανά μέτρηση, όλες οι επαναλήψεις σε μία μεταβλητή)
ch15.2l <- reshape(ch15.2w,direction="long",varying=c(1:6),timevar="Exercise:Day",times=names(ch15.2w),v.names="Fatigue")
ch15.2l$Exercise <- as.factor(substr(ch15.2l$"Exercise:Day",1,5))
ch15.2l$Day <- as.factor(substr(ch15.2l$"Exercise:Day",7,10))
ch15.2l$id <- as.factor(paste("S",ch15.2l$id,sep=""))
rownames(ch15.2l) <- NULL
ch15.2l$"Exercise:Day" <- NULL

# Ανάλυση και διάγραμμα αλληλεπίδρασης με διαστήματα εμπιστοσύνης, με τη βιβλιοθήκη ez 
library(ez)
ezANOVA(data=ch15.2l,dv=.(Fatigue),wid=.(id),within=.(Exercise,Day))
ezPlot(data=ch15.2l,dv=.(Fatigue),wid=.(id),within=.(Exercise,Day),x=.(Day),split=.(Exercise),y_lab="Κόπωση",x_lab="Ημέρα",split_lab="Άσκηση",
	levels=list(Day=list(new_names=c("Ημέρα 1","Ημέρα 2","Ημέρα 3")),Exercise=list(new_names=c("Βαριά","Ελαφριά"))))

# Για πολλαπλές συγκρίσεις, συμπεριλαμβάνουμε κάθε φορά μόνο τις υπό σύγκριση συνθήκες
bc <- length(levels(ch15.2l$Day)) # πλήθος συγκρίσεων, για διόρθωση Bonferroni
for (d in levels(ch15.2l$Day)) {
	data <- droplevels(subset(ch15.2l,Day!=d))
	dlev <- levels(data$Day)
	aac1 <- aov(Fatigue~Day*Exercise+Error(id/Day*Exercise),data)
	p<-summary(aac1)[["Error: id:Day"]][[1]]["Day","Pr(>F)"]*bc
	if (p > 1.0) { p <- 1.0 }
	cat(paste(paste(dlev,collapse="-"),": p ",round(p,3)),fill=T)
 }

# Υπολογισμός της διπαραγοντικής ανάλυσης διακύμανσης εξαρτημένων δειγμάτων στο R

# Φορτώνουμε τα ακόλουθα πακέτα
library(foreign)
library(ez)
library(emmeans)

# Διαβάζουμε τα δεδομένα του ερευνητικού παραδείγματος από το αρχείο και τα αποδίδουμε στο πλαίσιο δεδομένων ch15.3w
ch15.3w<-read.spss('chapter15_3.sav', to.data.frame=T)

# Βασικός έλεγχος του πλαισίου δεδομένων
str(ch15.3w)
summary(ch15.3w)

# Μετατροπή των δεδομένων σε μορφή long (μια σειρά ανά μέτρηση, όλες οι επαναλήψεις σε μία μεταβλητή)
Participants <- as.factor(c(ch15.3w$Participant_ID,ch15.3w$Participant_ID))
mem_perf <- c(ch15.3w$auditory_task, ch15.3w$visual_task)
memory_task<-as.factor(rep(c("auditory","visual"), each=297))
condition <- as.factor(c(ch15.3w$condition,ch15.3w$condition))
ch15.3l <- data.frame(Participants,mem_perf,memory_task,condition)

# Ανάλυση (μέθοδος 1)
MixedANOVA <- aov(mem_perf~memory_task*condition+Error(Participants/memory_task),ch15.3l)
summary(MixedANOVA)

# Ανάλυση και διάγραμμα αλληλεπίδρασης με διαστήματα εμπιστοσύνης με το πακέτο ez (ίδια αποτελέσματα με το SPSS)
ezANOVA(data=ch15.3l,dv=.(mem_perf),wid=.(Participants),within=.(memory_task),between=.(condition),type=3, detailed=T)
ezPlot(data=ch15.3l,dv=.(mem_perf),wid=.(Participants),within=.(memory_task),between=.(condition),x=.(memory_task),split=.(condition), 
	y_lab="Επίδοση στις μνημονικές δοκιμασίες", x_lab="", split_lab="Ομάδα", levels=list(memory_task=list(new_names=c("Auditory","Visual"))))

# Υπολογισμός του ελέγχου των απλών επιδράσεων του παράγοντα condition στον παράγοντα task με το πακέτο emmeans
Simple1<-emmeans(MixedANOVA,~memory_task|condition)
test(pairs(Simple1), joint=TRUE)