# Κεφάλαιο 7
# Έλεγχος κανονικότητας ερευνητικών δεδομένων στο R

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

# Διαβάζουμε τα δεδομένα του κεφαλαίου από το αρχείο και τα αποδίδουμε στο πλαίσιο δεδομένων ch7
ch7<-read.spss('chapter7_1.sav', to.data.frame=T)

## ΕΛΕΓΧΟΣ ΕΛΛΕΙΠΟΥΣΩΝ ΤΙΜΩΝ
# Βασικός έλεγχος του πλαισίου δεδομένων ch7
str(ch7)

# Πίνακες συχνοτήτων των δύο μεταβλητών
table(ch7$performance)
table(ch7$motivation)

# Αφαιρούμε τις δύο «αδύνατες» τιμές (>9) από τη μεταβλητή performance αντικαθιστώντας τις με ελλείπουσες τιμές (NA)
ch7$performance[ch7$performance>9]<-NA

# Έλεγχος του πλαισίου δεδομένων (εμφανίζονται οι ελλείπουσες τιμές ανά μεταβλητή)
summary(ch7)

# Πλήθος ελλειπουσών τιμών ανά μεταβλητή
apply(apply(ch7,2,is.na),2,sum)

# Ποσοστό (%) ελλειπουσών τιμών ανά μεταβλητή
round(100*apply(apply(ch7,2,is.na),2,sum)/nrow(ch7),1)

# Λίγο πιο απλά, ξεχωριστά η αναλογία για κάθε μεταβλητή
sum(is.na(ch7$performance))/length(ch7$performance)
sum(is.na(ch7$motivation))/length(ch7$motivation)

# Εξετάζουμε τις εγγραφές με ελλείπουσες τιμές
ch7[is.na(ch7$performance) | is.na(ch7$motivation),]

# Πλήθος ελλειπουσών τιμών ανά εγγραφή
ch7mv<-cbind(ch7,apply(apply(ch7,1,is.na),2,sum))
names(ch7mv)[3]<-"N missing"
ch7mv[ch7mv[["N missing"]]>0,]

# Το πακέτο mice παρέχει χρήσιμες συναρτήσεις για τον έλεγχο ελλειπουσών τιμών
md.pattern(ch7)

## ΙΣΤΟΓΡΑΜΜΑΤΑ
# Ιστογράμματα
hist(ch7$performance)
hist(ch7$motivation)

# Εναλλακτικά, μπορούμε να απεικονίσουμε  τους πίνακες συχνοτήτων
barplot(table(ch7$performance),las=1)
barplot(table(ch7$motivation),las=1)

# Ιστογράμματα με κανονική κατανομή αναφοράς 
# Καταρχάς αφαιρούμε τις ελλείπουσες τιμές
performance<-na.omit(ch7$performance)

# Ιστόγραμμα τιμών με τα κατάλληλα όρια
hist(performance,col="grey80",breaks=(min(performance)-1): max(performance),las=1,ylim=c(0,45),xlab="Performance")

# Διάνυσμα 100 τιμών για τις οποίες θα υπολογίσουμε την κανονική καμπύλη ώστε να φαίνεται συνεχής και ομαλή
x<-seq(min(performance)-1,max(performance),length.out=100)

# Εμφάνιση της καμπύλης με κατάλληλη προσαρμογή αξόνων
lines(x,dnorm(x+0.5,mean(performance),sd(performance))*length(performance))

# Το 0.5 μετακινεί οριζόντια την καμπύλη να συμφωνεί με το κέντρο κάθε μπάρας (αντί το σύνορο)
# Ο πολλαπλασιασμός με το πλήθος των τιμών μετατρέπει τη σχετική συχνότητα του dnorm σε αντίστοιχο πλήθος

# Επανάληψη για το κίνητρο
motivation<-na.omit(ch7$motivation)
hist(motivation,col="grey80",breaks=(min(motivation)-1):
max(motivation),las=1,xlab="Motivation")
x<-seq(min(motivation)-1,max(motivation),length.out=100)
lines(x,dnorm(x+0.5,mean(motivation),sd(motivation))*length(motivation))

## ΣΥΜΜΕΤΡΙΑ/ΛΟΞΟΤΗΤΑ ΚΑΙ ΚΥΡΤΩΣΗ
# Οι βασικοί τύποι
sk.p<- (sum((performance-mean(performance))^3)
	/length(performance)) / ((sum((performance-
	mean(performance))^2)/length(performance))^(3/2))
sk.m<- (sum((motivation-mean(motivation))^3)
	/length(motivation)) / ((sum((motivation-
	mean(motivation))^2)/length(motivation))^(3/2))
ku.p<- (sum((performance-mean(performance))^4)
	/length(performance)) / ((sum((performance-
	mean(performance))^2)/length(performance))^2) -3
ku.m<- (sum((motivation-mean(motivation))^4)
	/length(motivation)) / ((sum((motivation-
	mean(motivation))^2)/length(motivation))^2) -3

# τυπικό σφάλμα λοξότητας και κύρτωσης
Np<-length(performance)
Nm<-length(motivation)

# from Tabachnick & Fidell (6th ed., p. 79)
se.sk.p<-sqrt(6/Np)
se.ku.p<-sqrt(24/Np)
se.sk.m<-sqrt(6/Nm)
se.ku.m<-sqrt(24/Nm)

# from Wright & Herrington (2011)
se.sk<-sqrt((6*(N-2))/((N+1)*(N+3)))
se.ku<-sqrt((24*N*(N-2)*(N-3))/((N+3)*(N+5)*(N+1)^2))

# from http://www.tc3.edu/instruct/sbrown/stat/shape.htm
se.sk<-sqrt((6*N*(N-1))/((N-2)*(N+1)*(N+3)))
se.ku<-2*se.sk*sqrt((N^2-1)/((N-3)*(N+5)))

# Οι σχετικές συναρτήσεις υπάρχουν σε διάφορα πακέτα, π.χ. στο e1071
# (Οι συνηθισμένοι ορισμοί χρησιμοποιούνται με type=1, για το αποτέλεσμα του SPSS δίνουμε type=2)
skewness(performance,type=1)
skewness(motivation,type=1)
kurtosis(performance,type=1)
kurtosis(motivation,type=1)

# Το πακέτο psych (χρησιμοποιεί λίγο διαφορετικούς τύπους)
skew(performance)
skew(motivation)
kurtosi(performance)
kurtosi(motivation)

# Εναλλακτικά, από το πακέτο psych, μαζί με άλλα περιγραφικά
describe(ch7)

# Τυπικές τιμές ασυμμετρίας και κύρτωσης (διαιρούμε με το αντίστοιχο τυπικό σφάλμα)
z.sk.p<-sk.p/se.sk.p
z.sk.m<-sk.m/se.sk.m
z.ku.p<-ku.p/se.ku.p
z.ku.m<-ku.m/se.ku.m

# Αντίστοιχες πιθανότητες 
2*(1-pnorm(z.sk.p))
2*(1-pnorm(z.sk.m))
2*(1-pnorm(z.ku.p))
2*(1-pnorm(z.ku.m))

# Τα όρια πιθανότητας που αναφέρονται στο βιβλίο
round(2*(1-pnorm(2.58)),2) # .01
round(2*(1-pnorm(3.29)),3) # .001
# ...μπορούν να υπολογιστούν ως εξής (διαιρούμε δια 2 για αμφίπλευρο έλεγχο):
round(qnorm(1-(0.01/2)),2)
round(qnorm(1-(0.001/2)),2)


## ΔΙΑΓΡΑΜΜΑΤΑ ΚΑΙ ΕΛΕΓΧΟΣ ΚΑΝΟΝΙΚΟΤΗΤΑΣ
# Q-Q plot
qqnorm(performance)
qqline(performance)
qqnorm(motivation)
qqline(motivation)

# Το SPSS αντιστρέφει τους άξονες, θέτοντας τις παρατηρούμενες τιμές στον οριζόντιο 
qqnorm(performance,datax=T)
qqline(performance,datax=T)
qqnorm(motivation,datax=T)
qqline(motivation,datax=T)

# «Διορθωμένο» κανονικό γράφημα
as.data.frame(qqnorm(performance,datax=T,plot.it=F))->dfqp
if(nrow(dfqp)%%4 %in% c(1,3)){dfqp$y[dfqp$y==max(dfqp$y)]<-NA} 
p2<-dfqp[dfqp$y %in% fivenum(dfqp$y)[c(2,4)],] 
slope<-(diff(p2$y)/diff(p2$x))
dfqp$detrend<-scale(dfqp$x*slope-dfqp$y,scale=F)
plot(dfqp$x,dfqp$detrend)
abline(h=0)

# Εναλλακτικός υπολογισμός διόρθωσης
intercept<-p2$y[1]-p2$x[1]*slope
dfqp$detrend<-dfqp$y-(intercept+dfqp$x*slope)
as.data.frame(qqnorm(motivation,datax=T,plot.it=F))->dfqm
if(nrow(dfqm)%%4 %in% c(1,3)){dfqm$y[dfqm$y==max(dfqm$y)]<-NA}
m2<-dfqm[dfqm$y %in% fivenum(dfqm$y)[c(2,4)],] 
slope<-(diff(m2$y)/diff(m2$x))
dfqm$detrend<-scale(dfqm$x*slope-dfqm$y,scale=F)
plot(dfqm$x,dfqm$detrend)
abline(h=0)

# Θηκογράμματα
boxplot(performance)
boxplot(motivation)

# Στατιστικοί έλεγχοι κανονικότητας
# Kolmogorov-Smirnov (όχι αυτό που υπολογίζεται από το SPSS!)
ks.test(performance,pnorm)
ks.test(motivation,pnorm)

# Lilliefors (Kolmogorov-Smirnov), από το πακέτο nortest - αυτό είναι το κριτήριο που υπολογίζει το SPSS
lillie.test(performance)
lillie.test(motivation)

# Shapiro-Wilk
shapiro.test(performance)
shapiro.test(motivation)