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

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

# Απόδοση των δεδομένων του ερευνητικού παραδείγματος στο πλαίσιο δεδομένων ch14a
ch14a<-read.spss('chapter14_1.sav', to.data.frame=T)

# Υπολογισμός περιγραφικών στατιστικών μέτρων για τις 4 συνθήκες του πειράματος
describeBy(ch14a$performance,ch14a$reinforcer,mat=T,digits=3)

# Υπολογισμός του κριτηρίου Levene για τον έλεγχο της ισότητας των διακυμάνσεων
leveneTest(performance~reinforcer,ch14a,center=mean)

# Υπολογισμός της Ανάλυσης Διακύμανσης
anova <- aov(performance~reinforcer,ch14a)
anova
summary(anova)

# Υπολογισμός των μεγεθών επίδρασης (ω τετράγωνο και η τετράγωνο)
omega_sq(anova)
eta_sq(anova)

# Υπολογισμός όλων των συγκρίσεων κατά ζεύγη με τη χρήση t-test (διόρθωση Bonferroni)
pairwise.t.test(ch14a$performance,ch14a$reinforcer,pool.sd=T,p.adjust.method="bonf")

# Υπολογισμός των αντιθετικών γραμμικών συνδυασμών
cont <- rbind("Όλα έναντι σιωπής"=c(1,1,1,-3),
                   "Θετική - Μομφή"   =c(1,1,-2,0),
                   "Έπαινος-Γλύκισμα" =c(1,-1,0,0))
contrasts <- glht(anova,linfct=mcp(reinforcer=cont))
summary(contrasts,test=adjusted(type="none"))

# Υπολογισμός των εκ των υστέρων συγκρίσεων με τη χρήση του Tukey HSD
TukeyHSD(anova)

# Κατασκευή διαγράμματος διαστημάτων εμπιστοσύνης των μέσων διαφορών των παραπάνω συγκρίσεων
# Το όρισμα «mar» του «par» ορίζει το πλάτος των περιθωρίων με τη σειρά: κάτω, αριστερά, πάνω, δεξιά. Η προεπιλογή είναι να ορίσουμε «αριστερά» σε 4, εδώ το έχουμε αλλάξει σε 12 ώστε να χωρέσουν τα ονόματα των συγκρίσεων. 
par(mar=c(4,12,4,1)+.1)     # Δοκιμάστε να σχεδιάσετε το διάγραμμα χωρίς να τρέξετε την εντολή αυτή…
plot(TukeyHSD(anova))

# Κατασκευή διαγράμματος διαστημάτων εμπιστοσύνης των μέσων όρων των τεσσάρων συνθηκών
describeBy(ch14a$performance,ch14a$reinforcer,mat=T, digits=3)->stats
rownames(stats) <- levels(ch14a$reinforcer)
par(las=T)
par(cex.axis=0.7)
error.bars(stats=stats,eyes=F,sd=T,xlab="Ενισχυτής", ylab="95% ΔΕ Επίδοση",main="",ylim=c(30,90))

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

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

# Απόδοση των δεδομένων του ερευνητικού παραδείγματος στο πλαίσιο δεδομένων ch14b
ch14b<-read.spss('chapter14_2.sav', to.data.frame=T)

# Για τη συγκεκριμένη ανάλυση τα δεδομένα πρέπει να μετατραπούν σε μορφή long (μία μέτρηση σε κάθε σειρά)
longch14b<-melt(ch14b, id="participant", measured = c("anx1","anx2","anx3"))

# Οργάνωση των δεδομένων ώστε να εμφανίζονται στη σειρά οι μετρήσεις κάθε συμμετέχοντα
longch14b<-longch14b[order(longch14b$participant),]

# Ο απλούστερος υπολογισμός με τη χρήση του πακέτου ez
anova1<-ezANOVA(data=longch14b, dv=.(value), wid=.(participant), within=.(variable), detailed=FALSE, type=2)
anova1

# Κατασκευή του διαγράμματος διαστημάτων εμπιστοσύνης
ezPlot(data=longch14b, dv=.(value), wid=.(participant),within=.(variable), x=.(variable),do_lines=FALSE)

# Επιλογή των 3 μεταβλητών με τις τιμές των 3 συνθηκών της ανεξάρτητης μεταβλητής
myvars <- c("anx1", "anx2", "anx3")
ch14b <- ch14b[myvars]

# Οργάνωση των δεδομένων ώστε να εμφανίζονται στη σειρά οι μετρήσεις κάθε συμμετέχοντα – Μέθοδος 2
attr(ch14b,"codepage")<-NULL
write.table(ch14b,"chapter14_2a.Rdata")
longch14b <- reshape(ch14b,direction="long",varying=c(1:3),
			v.names="anxiety",times=names(ch14b),
			timevar="time",idvar="subject")
longch14b$time <- as.factor(longch14b$time)
longch14b$subject <- as.factor(paste("s",longch14b$subject,sep=""))
attr(longch14b,"reshapeLong")<-NULL
write.table(longch14b,"chapter14_2b.Rdata")

# Απλός υπολογισμός – Μέθοδος 2 (δεν περιλαμβάνει έλεγχο σφαιρικότητας ούτε δυνατότητα συγκρίσεων)
anova2 <- aov(anxiety~time+Error(subject),longch14b)
summary(anova2)

# Υπολογισμός μεγεθών επίδρασης (η2, η2p, generalized η2) μέσω του πακέτου DescTools
EtaSq(anova2, type=1, anova=T)

# Υπολογισμός με τη βοήθεια του πακέτου car – Μέθοδος 3
# (Περιλαμβάνεται έλεγχος σφαιρικότητας και διόρθωση Greenhouse-Geisser, Huynh-Feld)
ma1 <- as.matrix(ch14b)
lm1 <- lm(ma1~1)
kb <- factor(colnames(ma1))
anova3 <- Anova(lm1,idata=data.frame(kb), idesign=~kb,type="III")
summary(anova3,multi=F)

# Πραγματοποίηση του ελέγχου της σφαιρικότητας
mauchly.test(lm1,X=~1)

# Πραγματοποίηση πολλαπλών συγκρίσεων με τα πακέτα nlme και multcomp
anova4 <- lme(anxiety~time,random=~1|subject, data=longch14b)
anova(anova4)

# Πραγματοποίηση όλων των συγκρίσεων 
summary(glht(anova4,linfct=mcp(time="Tukey")))

# Πραγματοποίηση διαδοχικών συγκρίσεων (ανάλογες με την επιλογή του SPSS "Repeated")
summary(glht(anova4,linfct=mcp(time="Sequen"))) 

# Υπολογισμός του κριτηρίου Kruskal-Wallis στο R

# Απόδοση των δεδομένων του ερευνητικού παραδείγματος στο πλαίσιο δεδομένων ch14c
ch14c<-read.spss('chapter14_3.sav', to.data.frame=T)

# Υπολογισμός του κριτηρίου
kruskal.test(bullying~grade,ch14c)

# Υπολογισμός του κριτηρίου Friedman στο R

# Απόδοση των δεδομένων του ερευνητικού παραδείγματος στο πλαίσιο δεδομένων ch14d
ch14d<-read.spss('chapter14_4.sav', to.data.frame=T)

# Υπολογισμός του κριτηρίου
friedman.test(as.matrix(ch14d))

# Οργάνωση των δεδομένων ώστε να εμφανίζονται στη σειρά οι μετρήσεις κάθε συμμετέχοντα – Μέθοδος 2
longch14d <- reshape(ch14d,direction="long", varying=c(1:3),v.names="satisfaction", timevar="Time",times=names(ch14d))
longch14d$Time <- as.factor(longch14d$Time)
friedman.test(satisfaction~Time|id,longch14d)

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

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

# Απόδοση των δεδομένων του αρχείου chapter14_5.sav στο πλαίσιο δεδομένων ch14e
ch14e<-read.spss('chapter14_5.sav', to.data.frame=T)

# Έλεγχος προϋποθέσεων 1: Γραμμικότητα
# Κατασκευή διαγράμματος διασποράς μεταξύ της συμμεταβλητής και της εξαρτημένης μεταβλητής
# Προσθήκη γραμμών παλινδρόμησης
ggscatter(ch14e, x = "commitment", y = "satisfaction", color = "org_type", add = "reg.line" )+stat_regline_equation( aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = org_type))

# Έλεγχος προϋποθέσεων 2: Ομοιογένεια των κλίσεων παλινδρόμησης
ch14e %>% anova_test(satisfaction ~ org_type*commitment)

# Έλεγχος προϋποθέσεων 3: Κανονικότητα σφαλμάτων
model <- lm(satisfaction ~ commitment + org_type, data = ch14e)
model.metrics <- augment(model) %>% select(-.hat, -.sigma, -.fitted, -.se.fit)
# Remove details
# Έλεγχος κανονικότητας σφαλμάτων με το κριτήριο shapiro wilk
shapiro_test(model.metrics$.resid)

# Έλεγχος προϋποθέσεων 4: Ομοιογένεια των διακυμάνσεων
model.metrics %>% levene_test(.resid ~ org_type)

# Έλεγχος προϋποθέσεων 5: Έκτοπες/ακραίες τιμές
model.metrics %>% filter(abs(.std.resid) > 3) %>% as.data.frame()

# Ανάλυση συνδιακύμανσης
res.aov <- ch14e %>% anova_test(satisfaction ~ commitment + org_type) 
get_anova_table(res.aov)

# Εκ των υστέρων πολλαπλές συγκρίσεις με τον έλεγχο Bonferroni
pwc <- ch14e %>% emmeans_test( satisfaction ~ org_type, covariate = commitment, p.adjust.method = "bonferroni")
pwc

# Εκτιμώμενοι οριακοί μέσοι όροι
get_emmeans(pwc)

# Κατασκευή διαγράμματος 
pwc <- pwc %>% add_xy_position(x = "org_type", fun = "mean_se")
ggline(get_emmeans(pwc), x = "org_type", y = "emmean") + geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
	width = 0.2) + stat_pvalue_manual(pwc, hide.ns = TRUE, tip.length = FALSE) + labs( subtitle = get_test_label
	(res.aov, detailed = TRUE), caption = get_pwc_label(pwc))