Στο αρχείο αυτό θα προσπαθήσουμε να παρουσιάσουμε κάποια βασικά σημεία της χρήσης της R που συνήθως δημιουργούν προβλήματα, με έμφαση στις μετατροπές μεταξύ τύπων δεδομένων. Τα θέματα συζητούνται αποσπασματικά γι’ αυτό και το αρχείο δεν έχει ιδιαίτερη δομή. Μπορείτε να πλοηγηθείτε στην αντίστοιχη ενότητα χωρίς υποχρεωτικά να έχετε διαβάσει το σύνολό του.
Οι παράγοντες είναι στοιχεία των dataframes που αντιστοιχούν σε κατηγορικές μεταβλητές. Οι παράγοντες είναι στην ουσία διανύσματα χαρακτήρων που όμως μεταφράζονται πάντα σε κατηγορικά δεδομένα όταν χειριζόμαστε πλαίσια δεδομένων. Στο παρακάτω παράδειγμα βλέπουμε τη λειτουργία των παραγόντων στο ενσωματωμένο στην R πλαίσιο δεδομένων iris:
str(iris)## 'data.frame': 150 obs. of 5 variables:
## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...is.factor(iris$Species)## [1] TRUEΣτην τελευταία γραμμή του output της \(str()\) βλέπουμε ότι το διάνυσμα Species που περιέχει τα ονόματα των ειδών ίριδος αναγνωρίζεται ως factor με 3 επίπεδα (levels) τα οποία αντιστοιχούν στις τρείς διαφορετικές τιμές που παίρνει το διάνυσμα Species. Επί της ουσίας όλα τα διανύσματα που είναι factors είναι διανύσματα χαρακτήρων ωστόσο εφόσον χαρακτηριστούν ως factors η συνάρτηση \(is.character()\) δεν επαληθεύεται:
is.character(iris$Species)## [1] FALSEΓια ποιο λόγο η R μετατρέπει τα διανύσματα χαρακτήρων σε factors; Αυτό συμβαίνει επειδή οι παράγοντες είναι ιδιαίτερα χρήσιμοι για μια σειρά από χειρισμούς σε πλαίσια δεδομένων.
Παρά το γεγονός ότι η R θα αποδώσει την ιδιότητα παράγοντα σε κάθε διάνυσμα χαρακτήρων μέσα σε ένα dataframe, αυτό έχει νόημα μόνο στην περίπτωση που ο αριθμός των μοναδικών τιμών μέσα στο διάνυσμα είναι περιορισμένος ώστε να επιτρέπει την ομαδοποίηση. Το πλήθος των μοναδικών τιμών μπορεί κανείς να δει με τη χρήση της συνάρτησης \(levels()\) η οποία τις επιστρέφει σε ένα διάνυσμα χαρακτήρων:
levels(iris$Species)## [1] "setosa" "versicolor" "virginica"Στην πράξη οποιοδήποτε διάνυσμα χαρακτήρων μπορεί να μετατραπεί σε παράγοντα με τη χρήση της συνάρτησης \(factor()\):
a_vector<-c(rep("X",5),rep("Y",3), rep("Z",1))
a_factor<-factor(a_vector)
levels(a_factor)## [1] "X" "Y" "Z"Όταν παίρνουμε τα επίπεδα ενός παράγοντα, αυτά οργανώνονται πάντοτε με αλφαβητική σειρά. Στο παράδειγμα που ακολουθεί δημιουργούμε ένα πλαίσιο δεδομένων με τρία διανύσματα (εισόδημα, ηλικία και επίπεδο εκπαίδευσης), από τα οποία τα δύο πρώτα είναι αριθμητικές τιμές και το τελευταίο κατηγορική μεταβλητή που μετατρέπεται αυτόματα σε παράγοντα κατά τη δημιουργία του dataframe.
# Data in vectors
income <- c(32,51,12,23,26,27,24,26,25,18,30,22,24,21,25,27,18)
age <- c(48,59,26,23,37,32,20,25,45,55,44,42,51,30,35,29,19)
education <- c("middle","high","basic","high","high","middle","high","high","middle","basic", "high", "basic", "basic", "middle", "middle", "middle", "basic")
ageclass<-c("higher","higher","lower","lower","higher","lower","lower", "lower","higher","higher","higher","higher","higher","lower","higher","lower", "lower")
# Create the dataframe.
income_data <- data.frame(income,age,education,ageclass)
levels(income_data$education)## [1] "basic" "high" "middle"Σε αρκετές περιπτώσεις θέλουμε μια διαφορετική σειρά στα επίπεδα του παράγοντα. Αυτό γίνεται με την αντικατάσταση του παράγοντα από έναν παράγοντα όπου τα επίπεδα ορίζονται αναλυτικά όπως στο παράδειγμα:
income_data$education<-factor(income_data$education, levels=c("basic", "middle", "high"))
levels(income_data$education)## [1] "basic" "middle" "high"Μια πρώτη συνάρτηση που δρα σε παράγοντες είναι η \(table()\) που επενεργεί αποκλειστικά σε διανύσματα παραγόντων για να μας δώσει το πλήθος του κάθε level στο διάνυσμα:
table(income_data$education)##
## basic middle high
## 5 6 6Η R επιτρέπει την κατηγοριοποίηση μέσω πολλαπλών παραγόντων με τη χρήση της συνάρτησης \(table()\). Η συγκεκριμένη συνάρτηση μπορεί να δεχτεί ως είσοδο δύο ή περισσότερα διανύσματα παραγόντων και να δημιουργήσει έναν πίνακα “σύμπτωσης” (contingency table) που περιέχει τον αριθμό των στοιχείων των συνδυασμών τους.
table(income_data$education, income_data$ageclass)##
## higher lower
## basic 3 2
## middle 3 3
## high 3 3Η βασική χρησιμότητα των παραγόντων είναι η κατηγοριοποίηση των δεδομένων σε ένα dataframe. Η R αναγνωρίζει τις μεταβλητές που μπορούν να χρησιμοποιηθούν ως παράγοντες και μπορεί να κατατάξει τις υπόλοιπες μεταβλητές με βάση τα επίπεδά τους πριν κάνει υπολογισμούς σε διανυσματικό επίπεδο. Χαρακτηριστικές συναρτήσεις που επιτρέπουν υπολογισμούς αυτού του είδους είναι η \(by()\) και η \(aggregate()\). Ας δούμε πώς μπορούμε να εκμεταλλευτούμε τους παράγοντες με αυτές στο παράδειγμα του income_data που δημιουργήσαμε προηγουμένως:
by(income_data$income, income_data$education, mean)## income_data$education: basic
## [1] 18.8
## ------------------------------------------------------------
## income_data$education: middle
## [1] 26.16667
## ------------------------------------------------------------
## income_data$education: high
## [1] 30Στο πρώτο παράδειγμα η \(by(a,b,c)\) ορίζει αρχικά το διάνυσμα στο οποίο θέλει να πραγματοποιήσει τον υπολογισμό, στη δεύτερη θέση το διάνυσμα παραγόντων με το οποίο θα κάνει την κατηγοριοποίηση και στην τρίτη θέση τη συνάρτηση που θα εφαρμοστεί για τους υπολογισμούς. Το αποτέλεσμα εδώ είναι η μέση τιμή του εισοδήματος ανά μορφωτικό επίπεδο. Αντίστοιχη λειτουργία πραγματοποιεί η \(aggregate()\) με τη βασική διαφορά ότι τα διανύσματα παραγόντων που θα χρησιμοποιηθούν δίνονται με τη μορφή λίστας. Στο παρακάτω παράδειγμα υπολογίζουμε τη μέση τιμή των ηλικιών ανά μορφωτικό επίπεδο:
aggregate(income_data$age, by=list(Education=education), mean)## Education x
## 1 basic 38.60000
## 2 high 34.66667
## 3 middle 36.50000Ένα σημαντικό πλεονέκτημα της \(aggregate()\) είναι ότι μπορεί να εφαρμοστεί σε συνδυασμούς παραγόντων που οδηγούν σε πιο πολύπλοκες κατηγοριοποιήσεις. Στη συνέχεια βλέπουμε πώς μπορούμε να χρησιμοποιήσουμε την \(aggregate()\) για να υπολογίσουμε χαρακτηριστικές τιμές για αυτούς τους συνδυασμούς, κάτι που δεν μπορούμε να κάνουμε με την \(by()\).
aggregate(income_data$income, by=list(Education=education, Age=ageclass), mean)## Education Age x
## 1 basic higher 21.33333
## 2 high higher 35.66667
## 3 middle higher 27.33333
## 4 basic lower 15.00000
## 5 high lower 24.33333
## 6 middle lower 25.00000Από το οποίο προκύπτει ότι υψηλότερο μέσο εισόδημα έχουν τα μεγαλύτερης ηλικίας και μορφωτικού επιπέδου άτομα.
Οι λίστες (lists) είναι η πιο πολύπλοκη κατηγορία δεδομένων στην R. Αποτελούνται από σύνολα διανυσμάτων που δεν είναι υποχρεωτικό να είναι του ίδιου τύπου αλλά ούτε και του ίδιου μήκους. Στην πιο διευρυμένη μορφή τους μπορούν να περιέχουν ακόμα και πίνακες ή και συναρτήσεις ως στοιχεία. Μια λίστα δημιουργείται με τη χρήση της ομώνυμης συνάρτησης, η οποία δέχεται μια παράθεση στοιχείων, όπως στο παράδειγμα:
l<-list(c(1,3,5), c("Steven", "Robbie", "Kenny", "Ian"), TRUE, c(0.9, 0.8))
l## [[1]]
## [1] 1 3 5
##
## [[2]]
## [1] "Steven" "Robbie" "Kenny" "Ian"
##
## [[3]]
## [1] TRUE
##
## [[4]]
## [1] 0.9 0.8Βλέπουμε ότι η λίστα απαρτίζεται από τις παραμέτρους που περάσαμε στην \(list()\) με τη σειρά που το κάναμε. Για να δημιουργήσουμε μια πιο καλά οργανωμένη εκδοχή της θα δώσουμε ονόματα στα στοιχεία της \(l\) με τη χρήση της \(names\):
names(l)<-c("odds", "names", "logical", "decimals")
l## $odds
## [1] 1 3 5
##
## $names
## [1] "Steven" "Robbie" "Kenny" "Ian"
##
## $logical
## [1] TRUE
##
## $decimals
## [1] 0.9 0.8Προσέξτε ότι τώρα τα ονόματα των στοιχείων της λίστας βρίσκονται στη θέση που προηγουμένων υπήρχε αρίθμηση σε διπλές αγκύλες. Ίσως το πιο σημαντικό σημείο σε ό,τι αφορά τις λίστες είναι η χρήση αυτών των διπλών αγκυλών. Είτε φέρουν ονόματα, είτε όχι, η θέση των στοιχείων μιας λίστας δίνεται από έναν αριθμό μέσα σε διπλές αγκύλες αντί για μονές όπως ισχύει για τα διανύσματα και τους πίνακες. Η διαφορά στη χρήση μονών και διπλών αγκυλών είναι αρκετά λεπτή και μπορεί να διαφύγει σε πρώτο επίπεδο. Προσέξτε τη διαφορά μεταξύ:
l[1]## $odds
## [1] 1 3 5και:
l[[1]]## [1] 1 3 5Πρακτικά δε φαίνεται να υπάρχει καμία. Δείτε όμως τι είδους αντικείμενο κρύβεται πίσω από τις δύο αυτές δομές:
class(l[1])## [1] "list"class(l[[1]])## [1] "numeric"Χρήση δηλαδή των μονών αγκυλών επιστρέφει ένα αντικείμενο τύπου λίστας ενώ οι διπλές αγκύλες επιστρέφουν το διάνυσμα του πρώτου στοιχείου. Αυτό είναι πολύ σημαντικό καθώς όπως έχουμε συζητήσει και παραπάνω συγκεκριμένες συναρτήσεις μπορούν να εφαρμοστούν μόνο σε συγκεκριμένους τύπους αντικειμένων. Πώς μπορούμε τώρα να εντοπίσουμε ένα στοιχείο εντός των διανυσμάτων που απαρτίζουν μια λίστα; Ακολουθώντας τη λογική που έχουμε δει ως τώρα, αρκεί να προσδιορίσουμε αρχικά το διάνυσμα που θέλουμε με διπλές αγκύλες και στη συνέχεια να ορίσουμε σε αυτό την επιθυμητή θέση με μονές. Έτσι π.χ. το τρίτο όνομα στο στοιχείο \(names\) της \(l\) δίνεται ως εξής:
l[[2]][3]## [1] "Kenny"Εναλλακτικά, η πρόσβαση σε αυτό μπορεί να γίνει και με τη χρήση του ονόματος του διανύσματος όπως είδαμε στα dataframes:
l$names[3]## [1] "Kenny"Στην περίπτωση που θέλουμε να πάρουμε πίσω όλη τη λίστα σε ένα διάνυσμα μπορούμε να την “καταρρεύσουμε” σε διάνυσμα με την συνάρτηση \(unlist()\):
unlist(l)## odds1 odds2 odds3 names1 names2 names3 names4 logical
## "1" "3" "5" "Steven" "Robbie" "Kenny" "Ian" "TRUE"
## decimals1 decimals2
## "0.9" "0.8"Όπως φαίνεται, στην περίπτωση ονοματισμένης λίστας, τα ονόματα διατηρούνται και στο διάνυσμα, κάτι που είναι αρκετά βολικό.
Οι χειρισμοί των λιστών μοιάζουν πολύ με τους αντίστοιχους των διανυσμάτων. Έτσι μπορούμε να ενώσουμε δύο ή περισσότερες λίστες με την συνάρτηση \(c()\):
l1<-list(c("A","B"), c(1,2,3))
l2<-list(c("C","D","E"), c(T,F,F,F))
l3<-c(l1,l2)
l3## [[1]]
## [1] "A" "B"
##
## [[2]]
## [1] 1 2 3
##
## [[3]]
## [1] "C" "D" "E"
##
## [[4]]
## [1] TRUE FALSE FALSE FALSEΗ βασικότερη χρησιμότητα των λιστών είναι η δυνατότητα εφαρμογής συναρτήσεων στα στοιχεία τους με τη χρήση των συναρτήσεων \(apply()\) που λειτουργούν καλώντας άλλες συναρτήσεις και εφαρμόζοντάς τις απευθείας, με αποτέλεσμα την ταχύτατη εκτέλεσή τους. Απαραίτητη προϋπόθεση είναι τα στοιχεία της λίστας να μπορούν να αποτελέσουν όρισμα της συνάρτησης. Στην περίπτωση των λιστών, το πλεονέκτημα που έχουμε σε σχέση με τους πίνακες είναι ότι τα διανύσματα που μπορούμε να αναλύσουμε δεν είναι υποχρεωτικά ίδιου μήκους. Έστω για παράδειγμα ότι έχουμε μια σειρά από διανύσματα των οποίων θέλουμε να υπολογίσουμε τις μέγιστες τιμές. Μπορούμε να αποφύγουμε μια επαναληπτική κλήση της αντίστοιχης συνάρτησης \(max()\) ενσωματώνοντας τα διανύσματα σε μια λίστα και καλώντας την \(max()\) σε μία κλήση της \(sapply()\):
x<-c(1, 8, 9, 12, 5)
y<-c(6, 7, 0.8, 2, 5.1)
z<-c(-10, 5, 2, 3.8)
my_list<-list(x,y,z)
sapply(my_list, FUN=max)## [1] 12 7 5Το αποτέλεσμα είναι ένα διάνυσμα τριών τιμών που αντιστοιχεί στις μέγιστες τιμές των \(x\), \(y\) και \(z\). Το ίδιο αποτέλεσμα καταχωρημένο σε λίστα προκύπτει αν αντί για την \(sapply()\) καλέσουμε την \(lapply()\):
lapply(my_list, FUN=max)## [[1]]
## [1] 12
##
## [[2]]
## [1] 7
##
## [[3]]
## [1] 5Με τον όρο Ανάλυση Διακύμανσης ή Ανάλυση Διασποράς (Analysis of Variance) αναφερόμαστε σε μια δέσμη από στατιστικές αναλύσεις μοντελοποίησης που σκοπό έχουν την εκτίμηση της διαφοράς των μέσων τιμών μεταξύ ομάδων σε ευρύτερα δείγματα. Όπως και πολλά άλλα μεθοδολογικά εργαλεία, η ANOVA αναπτύχθηκε από τον Ronald Fisher και στηρίζεται στην διαισθητικά απλή αλλά τεχνικά μάλλον πολύπλοκη βασική αρχή της ανάλυσης της παρατηρούμενης διακύμανσης μιας μεταβλητής σε διαφορετικές συνιστώσες που σχετίζονται με εγγενείς τάσεις και τυχαία σφάλματα. Η ANOVA είναι η μεθοδολογική προσέγγιση επιλογής για την σύγκριση περισσότερων από δύο δείγματα. Εννοιολογικά δε διαφέρει πολύ από μια προσέγγιση πολλαπλών ζευγαρωτών t-test, ωστόσο είναι γενικά πιο αυστηρή και έτσι οδηγεί σε πιο συντηρητικά αποτελέσματα με περιορισμένα σφάλματα τύπου Ι, οδηγεί δηλαδή σπανιότερα στην απόρριψη μηδενικών υποθέσεων που είναι αληθείς.
Η βάση της ANOVA είναι η εκτίμηση της διαφοράς στη μέση τιμή μιας παραμετρικής μεταβλητής σε σχέση με μια ή περισσότερες κατηγορικές μεταβλητές. Ανάλογα με το σχεδιασμό του πειράματος, των αριθμών των μεταβλητών και των ομάδων δειγμάτων υπάρχουν διάφορες παραλλαγές της ANOVA. Έτσι υπάρχει μονο- (one-way) και πολυ-παραγοντική ανάλυση (n-way ANOVA) ανάλογα με το αν η διακύμανση μελετάται σε σχέση με μία ή περισσότερες κατηγορικές μεταβλητές, η πολυμεταβλητή (mulitvariate ANOVA ή MANOVA) όταν η μετρούμενες μεταβλητές απόκρισεις είναι παραπάνω από μία, ενώ ανάλογα με το αν υπάρχουν επιπλέον παραμετρικές μεταβλητές που μπορούν να επηρεάσουν την μετρούμενη (target variable) χρησιμοποιείται ένας συνδυασμός ΑNOVA και παλινδρόμησης (regression) που ονομάζεται Analysis of Co-Variance.
Σε πρακτικό επίπεδο, παρότι χρησιμοποιείται στον έλεγχο υποθέσεων, η ANOVA αποτελεί στην ουσία μια κατηγορία μοντέλων που ονομάζονται γενικευμένα γραμμικά μοντέλα (generalized linear models). Σημείο εκκίνησης είναι η μηδενική υπόθεση ότι οι τιμές όλων των ομάδων/κατηγοριών που εξετάζονται είναι τυχαία δείγματα που προέρχονται από την ίδια κατανομή. Η διαφορά με άλλους (ζευγαρωτούς) ελέγχους υποθέσεων είναι ότι η ANOVA λειτουργεί προσπαθώντας να μοντελοποιήσει την σχέση μεταξύ της μεταβλητής που μετρούμε, που ονομάζεται και μεταβλητή απόκρισης (response variable) και των κατηγορικών μεταβλητών που ονομάζονται και επεξηγηματικές (explanatory variables). Ο στόχος του μοντέλου είναι να εκτιμηθεί η σχέση αυτή σε επίπεδο διασποράς.
H R διαθέτει διάφορες συναρτήσεις για την διενέργεια της ANOVA. Σε αυτή την φάση θα εξετάσουμε την \(aov()\), η οποία ακολουθεί την γενική σύνταξη των γραμμικών μοντέλων που θα δούμε πιο αναλυτικά σε επόμενα κεφάλαια. Συνοπτικά η σύνταξή της είναι:
aov(ResponseVariable ~ A + B + C …, data=data.frame)
όπου ResponseVariable είναι η παραμετρική μεταβλητή που μας ενδιαφέρει να ελέγξουμε και Α, Β, C, κλπ είναι οι κατηγορικές μεταβλητές στη βάση των οποίων θέλουμε να κάνουμε τον έλεγχο. Σε μια προσέγγιση, στην οποία οι κατηγορικές μεταβλητές είναι ανεξάρτητες, αυτές δίνονται με τον τελεστή “+” ανάμεσά τους, ενώ σε περίπτωση που έχουμε συσχετιζόμενες κατηγορικές μεταβλητές (factorial ANOVA) μπορούμε να το δηλώσουμε με "*".
aov(ResponseVariable ~ A * B * C …, data=data.frame)
Ο πιο σημαντικός τελεστής και στις δύο περιπτώσεις είναι ο τελεστής σύνδεσης (~) ο οποίος στην R υποκαθιστά τη σχέση f() των μαθηματικών μοντέλων, υποδηλώνει δηλαδή την σχέση εξαρτημένων (αριστερά) και ανεξάρτητων (δεξιά) μεταβλητών.
H εφαρμογή των συναρτήσεων που θα παρουσιάσουμε στη συνέχεια δεν διαφέρει μεταξύ των παραπάνω περιπτώσεων, είναι ωστόσο σημαντικό να έχουμε μια καλή εποπτεία των δεδομένων μας και των παραμέτρων που θέλουμε να εξετάσουμε.
Πριν προχωρήσουμε στην πρακτική εφαρμογή να επισημάνουμε κάτι σημαντικό σχετικά με την δομή των δεδομένων πριν τη διενέργεια ANOVA στην R. Τα δεδομένα θα πρέπει να έχουν τη μορφή dataframe με τις δύο βασικές μεταβλητές α) τη συγκρινόμενη αριθμητική μεταβλητή και β) την κατηγορική μη-αριθμητική μεταβλητή που ορίζει τα groups. Παρακάτω βλέπετε ένα έτοιμο παράδειγμα που διαθέτει η R:
chickwts## weight feed
## 1 179 horsebean
## 2 160 horsebean
## 3 136 horsebean
## 4 227 horsebean
## 5 217 horsebean
## 6 168 horsebean
## 7 108 horsebean
## 8 124 horsebean
## 9 143 horsebean
## 10 140 horsebean
## 11 309 linseed
## 12 229 linseed
## 13 181 linseed
## 14 141 linseed
## 15 260 linseed
## 16 203 linseed
## 17 148 linseed
## 18 169 linseed
## 19 213 linseed
## 20 257 linseed
## 21 244 linseed
## 22 271 linseed
## 23 243 soybean
## 24 230 soybean
## 25 248 soybean
## 26 327 soybean
## 27 329 soybean
## 28 250 soybean
## 29 193 soybean
## 30 271 soybean
## 31 316 soybean
## 32 267 soybean
## 33 199 soybean
## 34 171 soybean
## 35 158 soybean
## 36 248 soybean
## 37 423 sunflower
## 38 340 sunflower
## 39 392 sunflower
## 40 339 sunflower
## 41 341 sunflower
## 42 226 sunflower
## 43 320 sunflower
## 44 295 sunflower
## 45 334 sunflower
## 46 322 sunflower
## 47 297 sunflower
## 48 318 sunflower
## 49 325 meatmeal
## 50 257 meatmeal
## 51 303 meatmeal
## 52 315 meatmeal
## 53 380 meatmeal
## 54 153 meatmeal
## 55 263 meatmeal
## 56 242 meatmeal
## 57 206 meatmeal
## 58 344 meatmeal
## 59 258 meatmeal
## 60 368 casein
## 61 390 casein
## 62 379 casein
## 63 260 casein
## 64 404 casein
## 65 318 casein
## 66 352 casein
## 67 359 casein
## 68 216 casein
## 69 222 casein
## 70 283 casein
## 71 332 caseinΕδώ βλέπουμε την αριθμητική μεταβλητή (weight) πουλερικών που έχουν λάβει διαφορετικό συμπλήρωμα διατροφής (feed). Το feed εδώ θα λειτουργήσει σαν κατηγορία για να δούμε με ποιο διατροφικό συμπλήρωμα αυξάνεται περισσότερο το βάρος των πουλιών.
Στην περίπτωση που θα δούμε παρακάτω μελετάμε τη μάζα ζώων (πουλερικών) στη βάση διαφορετικών διατροφικών συμπληρωμάτων. Εδώ έχουμε την απλούστερη περίπτωση μιας μεταβλητής απόκρισης (βάρος) και μιας κατηγορικής μεταβλητής (διατροφικό συμπλήρωμα). Με βάση τα παραπάνω η κλήση της συνάρτησης \(aov()\) θα είναι:
aov(weight~feed, data=chickwts)## Call:
## aov(formula = weight ~ feed, data = chickwts)
##
## Terms:
## feed Residuals
## Sum of Squares 231129.2 195556.0
## Deg. of Freedom 5 65
##
## Residual standard error: 54.85029
## Estimated effects may be unbalancedΤο output από μόνο του δε φαίνεται να μας λέει πολλά καθώς αναγράφει μόνο μέρος των στατιστικών που έχει υπολογίσει η ANOVA χωρς να γίνεται λόγος για στατιστική σημαντικότητα. Τις πληροφορίες αυτές παίρνουμε αν αποθηκεύσουμε το αποτέλεσμα της ανάλυσης σε μια μεταβλητή και εκτελέσουμε πάνω της την \(summary()\):
aov(weight~feed, data=chickwts)->fit
summary(fit)## Df Sum Sq Mean Sq F value Pr(>F)
## feed 5 231129 46226 15.37 5.94e-10 ***
## Residuals 65 195556 3009
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1To αποτέλεσμα σε αυτήν την περίπτωση είναι αρκετά πιο πλούσιο σε πληροφορία. Οι βαθμοί ελευθερίας υπήρχαν και πριν και είναι 5 για την κατηγορική μεταβλητή \(feed\) (N-1, με N=6 διατροφικά συμπληρώματα) ενώ η αντίστοιχη τιμή df για τα σφάλματα είναι 65, όπως προκύπτει αν αφαιρέσουμε από το συνολικό πλήθος των τιμών (Ν-1, για Ν=71) τους βαθμούς ελευθερίας που έχουν ενσωματωθεί στην κατηγορική μεταβλητή (Ν=5).
Η τιμή που κυρίως μας ενδιαφέρει εδώ είναι η F και η πιθανότητα μια τέτοια η μεγαλύτερη τιμή να έχει προκύψει δεδομένου ότι οι κατανομές έχουν την ίδια διασπορά.
Ακόμα κι έτσι όμως δεν έχουμε τα δεδομένα που θέλουμε γιατί μας λείπουν οι ποσοτικές μεταβολές των μέσων τιμών βάρους ανά κατηγορία. Οι τιμές αυτές εμπεριέχονται στο αποτέλεσμα της ANOVA (fit) και μπορούμε να τις λάβουμε εφαρμόζοντας μια επιπλέον συνάρτηση:
TukeyHSD(fit)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = weight ~ feed, data = chickwts)
##
## $feed
## diff lwr upr p adj
## horsebean-casein -163.383333 -232.346876 -94.41979 0.0000000
## linseed-casein -104.833333 -170.587491 -39.07918 0.0002100
## meatmeal-casein -46.674242 -113.906207 20.55772 0.3324584
## soybean-casein -77.154762 -140.517054 -13.79247 0.0083653
## sunflower-casein 5.333333 -60.420825 71.08749 0.9998902
## linseed-horsebean 58.550000 -10.413543 127.51354 0.1413329
## meatmeal-horsebean 116.709091 46.335105 187.08308 0.0001062
## soybean-horsebean 86.228571 19.541684 152.91546 0.0042167
## sunflower-horsebean 168.716667 99.753124 237.68021 0.0000000
## meatmeal-linseed 58.159091 -9.072873 125.39106 0.1276965
## soybean-linseed 27.678571 -35.683721 91.04086 0.7932853
## sunflower-linseed 110.166667 44.412509 175.92082 0.0000884
## soybean-meatmeal -30.480519 -95.375109 34.41407 0.7391356
## sunflower-meatmeal 52.007576 -15.224388 119.23954 0.2206962
## sunflower-soybean 82.488095 19.125803 145.85039 0.0038845Η εφαρμογή της \(TukeyHSD()\) μας επιστρέφει έναν πίνακα όπου οι συγκρίσεις είναι στην πρώτη στήλη ακολουθούμενες από τέσσερις αριθμητικές τιμές που είναι α) η διαφορά των μέσων τιμών για τη δεδομένη σύγκριση β-γ) τα κάτω-πάνω όρια του διαστήματος εμπιστοσύνης και δ) η διορθωμένη τιμή p-value λαμβάνοντας υπ-οψιν τις πολλαπλές συγκρίσεις.
Αν π.χ. μας ενδιαφέρει το πόσο καλύτεροι είναι οι ηλίοσποροι από το λιναρόσπορο η διαφορά μέσων τιμών είναι αυτή που αντιστοιχεί στο sunflower-linseed=110.166667.
Tελευταία επισήμανση θα είναι το πως απομωνώνουμε τις τιμές από αυτόν τον πίνακα. Αυτό γίνεται με τη χρήση δεικτών που θα ακολουθούν τη δομή του πίνακα αποτελεσμάτων ανά στηλη. Ετσι οι παρακάτω εντολές:
results<-TukeyHSD(fit)
results[[1]][1:2]## [1] -163.3833 -104.8333Δίνουν τις πρώτες δύο τιμές diff του πίνακα. Αντίστοιχα αν κάποιος θέλει την τιμή diff και την τιμή p-value για τη σύγκριση sunflower-linseed θα πρέπει να υπολογίσει τους δείκτες στους οποίους βρίσκονται αυτές. Π.χ. η diff είναι η:
results[[1]][12]## [1] 110.1667ενώ για την αντίστοιχη p-value θα πρέπει να προσθέσει τρία μήκη του πίνακα ώστε να περάσει από τις στήλες lwr, upr στην τέταρτη 12+ 15*3 = 57
results[[1]][57]## [1] 8.843233e-05