Έκδοση: 07 / 06 / 2024

Προετοιμασία περιβάλλοντος R
Οι παρακάτω εντολές, ελέγχουν αν έχουν εγκατασταθεί τα απαιτούμενα πακέτα για την εκτέλεση του συνόλου του κώδικα της ενότητας.

list.of.packages <- c("onewaytests.", "HH", "Rmisc", "gplots", "agricolae", "lsr", "car", "MASS", "moments", "afex", "plyr", "heplots", "magrittr", "dunn.test", "conover.test")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

Όλες οι συναρτήσεις που χρησιμοποιούνται στην ενότητα αυτή βρίσκονται στο αρχείο της γλώσσας R: MyRFunctions.R.

Για να γίνουν διαθέσιμες για χρήση αρκεί να φορτωθούν στο περιβάλλον της R, εκτελώντας την εντολή

source("https://utopia.duth.gr/epdiaman/files/kedivim/MyRFunctions.R")

(Στην περίπτωση όπου γίνει λήψη τοπικά το αρχείο MyRFunctions.R, πρέπει να αλλαξει ανάλογα και η διεύθυνση του αρχείου)

1 Εισαγωγή

Με τον όρο Ανάλυση Διακύμανσης (Analysis of Variation ή ANOVA) περιγράφεται η μεγάλη οικογένεια στατιστικών ελέγχων που έχει στόχο τον εντοπισμό διαφοροποιήσεων μεταξύ των μέσων τιμών μίας ή περισσοτέρων εξαρτημένων συνεχών μεταβλητών ανάμεσα σε ανεξάρτητες ή εξαρτημένες ομάδες που ορίζονται από κάποιον ή κάποιους ανεξάρτητους ποιοτικούς παράγοντες, ελέγχοντας παράλληλα για την επίδραση κάποιων συνεχών παραμέτρων.

2 Ανάλυση Διακύμανσης με έναν παράγοντα (1-way ANOVA).

Χρησιμοποιούμε ανάλυση διακύμανσης με έναν παράγοντα (1 - way ANOVA) όταν θέλουμε να συγκρίνουμε μέσες τιμές από δύο ή περισσότερες ομάδες. Αποτελεί τη γενίκευση του t - test για δύο ανεξάρτητα δείγματα.

Παράδειγμα
Μια επιχείρηση επιχειρώντας να βρει την απαραίτητη διάρκεια εκπαίδευσης των πωλητών της για την προώθηση ενός νέου προϊόντος επιλέγει με τυχαίο τρόπο 31 πωλητές τους οποίους χωρίζει σε τρεις ομάδες (Α, n = 10), (Β, n = 11), (Γ, n = 10), τους εκπαιδεύει για 1, 2 και 3 ημέρες αντίστοιχα και παρακολουθεί τις πωλήσεις τους τον επόμενο μήνα. Τα δεδομένα κατσχωρούνται στο dataframe anova.df, με τη μεταβλητή group («Ομάδα εκπαίδευσης») και τιμές 1 ( “Α”), 2 ( “Β”), 3 (“Γ”) και τη μεταβλητή perform ( «Πωλήσεις») όπου καταχωρούνται οι πωλήσεις του κάθε πωλητή.

perform.A = c(63.3, 68.3, 86.7, 52.8, 75, 58, 69.5, 32.7, 60.9, 58.2)
perform.B = c(72.8, 88.2, 80.8, 71.3, 81.5, 47.6, 81, 81.4, 83, 76, 74.5)
perform.C = c(82.3, 89.7, 81, 85.1, 74.1, 75.9, 74.7, 81.1, 76.4, 81.8)
sales = c(perform.A, perform.B, perform.C)
group = c(rep(1, 10), rep(2, 11), rep(3, 10))
group = factor(group, levels = c(1, 2, 3), labels = c("Α", "Β", "Γ"))

anova.df = data.frame(group, sales)

Hmisc::label(anova.df$group) = 'Ομάδα'
Hmisc::label(anova.df$sales) = 'Πωλήσεις'

print(anova.df)
##    group sales
## 1      Α  63.3
## 2      Α  68.3
## 3      Α  86.7
## 4      Α  52.8
## 5      Α  75.0
## 6      Α  58.0
## 7      Α  69.5
## 8      Α  32.7
## 9      Α  60.9
## 10     Α  58.2
## 11     Β  72.8
## 12     Β  88.2
## 13     Β  80.8
## 14     Β  71.3
## 15     Β  81.5
## 16     Β  47.6
## 17     Β  81.0
## 18     Β  81.4
## 19     Β  83.0
## 20     Β  76.0
## 21     Β  74.5
## 22     Γ  82.3
## 23     Γ  89.7
## 24     Γ  81.0
## 25     Γ  85.1
## 26     Γ  74.1
## 27     Γ  75.9
## 28     Γ  74.7
## 29     Γ  81.1
## 30     Γ  76.4
## 31     Γ  81.8

Βρίσκουμε τις μέσες τιμές των πωλήσεων ανά ομάδα, αξιοποιώντας τη συνάρτηση my_explore_vars_by_group:

my_explore_vars_by_group(anova.df, 'sales', "group")
Descriptive Statistics of Πωλήσεις among levels of Ομάδα
Group N M (SD) 95% C.I.
1 10 62.5 (14.3) 52.3 - 72.8
2 11 76.2 (10.7) 69 - 83.4
3 10 80.2 (4.98) 76.6 - 83.8
Total 31 73.1 (12.8) 68.4 - 77.8
Tukeys HSD Post Hoc Test for variable: Πωλήσεις
The Tukey HSD is not a Post Hoc Test per se. It may provide valuable results as an independent test as well.
Difference Lower Upper p
2-1 13.651 2.061 25.24 0.018
3-1 17.67 5.808 29.532 0.003
3-2 4.019 -7.57 15.609 0.671

Αν μΑ, μΒ και μΓ είναι οι μέσες πωλήσεις των πωλητών της Α, Β και Γ ομάδας αντίστοιχα, αναζητούμε την απάντηση στην ερώτηση σχετικά με το αν απορρίπτεται ή όχι η στατιστική υπόθεση:

Η0: μΑ = μΒ = μΓ,

έναντι της ερευνητικής υπόθεσης:

Η1: όχι η Η0.

2.1 Υλοποίηση ανάλυσης διακύμανσης

Έχοντας επιβεβαιώσει τις απαραίτητες προϋποθέσεις της κανονικότητας και της ομοιογένειας, εφαρμόζουμε τη συνάρτηση my_ANOVA.

all.ANOVA.output = my_ANOVA(data = anova.df, model = sales ~ group)
all.ANOVA.output$htmlTable
ANOVA Results for Πωλήσεις
SS df F p
(Intercept) 39112.516 1 340.354 <0.001
group 1725.74 2 7.509 0.002
Residuals 3217.682 28
LSD Test for variable: Ομάδα
Fisher’s LSD is a series of pairwise t-tests, with each test using the mean squared error from the significant ANOVA as its pooled variance estimate (and naturally taking the associated degrees of freedom). It is a valid test, in case where the ANOVA is significant.
N Πωλήσεις Groups
Α 10 62.54 b
Β 11 76.191 a
Γ 10 80.21 a
Tukeys HSD Post Hoc Test for variable: Ομάδα
The Tukey HSD is not a Post Hoc Test per se. It may provide valuable results even if ANOVA is not significant.
Difference Lower Upper p
Β-Α 13.651 2.061 25.24 0.018
Γ-Α 17.67 5.808 29.532 0.003
Γ-Β 4.019 -7.57 15.609 0.671
Standard Deviation of Πωλήσεις among levels of all levels of group.
Level Α Β Γ
Group Size 10 11 10
SD 14.30 10.70 4.98
Levene test of variance equality (homogeneity) of Πωλήσεις over all levels of group.
F(2, 28) = 1.489, p = 0.243. Homogeneity hypothesis is confirmed.
Bartlett’s test of variance equality (homogeneity) of Πωλήσεις over all levels of group.
If you have strong evidence that your data do in fact come from a normal, or nearly normal, distribution, then Bartlett’s test has better performance.
c2(2) = 8.214, p = 0.016. Homogeneity hypothesis is not supported!
Normality test for model residuals.
If the main goal of an ANOVA is to see whether or not certain effects are significant, then the assumption of normality of the residuals is only required for small samples, thanks to the central limit theorem. With sample sizes of a few hundred participants even extreme violations of the normality assumptions are unproblematic. So mild violations of this assumptions are usually no problem with sample sizes exceeding 30.
Statistic Description
Skewness & Kurtosis
  Skewness -0.965 Normal distribution has skew = 0. For a unimodal distribution, negative skew commonly indicates that the tail is on the left side of the distribution, and positive skew indicates that the tail is on the right.
  Kurtosis 5.62 Normal distribution has kurtosis = 3. Values over 3 indicates a platykurtic distribution and values less than 3 indicates a leptokurtic distribution.
Normality Tests. Η0: This sample is from a normal distribution vs Η1: Not the Η0.
  Shapiro–Wilk W = 0.884, p = 0.003 This test is more appropriate method for small sample sizes (<50 samples) although it can also be handling on larger sample size.
  Lilliefors D = 0.181, p = 0.011 The Lilliefors test uses the same calculations as the Kolmogorov-Smirnov test, but it is more conservative in the sense that the Lilliefors Test is less likely to show that data is normally distributed.

Παρατηρούμε πως η στατιστική σημαντικότητα του ελέγχου (στήλη p) είναι ίση με p = 0,002 μικρότερο από το όριο του 0,05 άρα η στατιστική υπόθεση:

Η0: μΑ = μΒ = μΓ,

απορρίπτεται, δηλαδή οι τρεις ομάδες έχουν μέση τιμή με στατιστικά σημαντική διαφορά (F(2, 28) = 7,509, p = 0,002).

Προϋποθέσεις
Απαραίτητες θεωρητικές προϋποθέσεις για την ανάλυση διακύμανσης είναι:

• η κανονικότητα της κατανομής των τιμών της συνεχούς μεταβλητής σε κάθε μία ομάδα, και

• η ομοιογένεια των ομάδων (η τυπική απόκλιση δεν διαφέρει σημαντικά μεταξύ των ομάδων).

Στην περίπτωση όπου η κανονικότητα των υπολοίπων δεν γίνεται αποδεκτή, τότε μειώνεται η εγκυρότητα του μοντέλου, δηλαδή του συνόλου των επεξηγηματικών μεταβλητών που επιλέχθηκαν για να εξηγήσουν τη μεταβλητότητα της εξαρτημένης μεταβλητής.

Στην περίπτωση όπου η ομοιογένεια των ομάδων δεν γίνεται αποδεκτή, τότε αυξάνεται το σφάλμα τύπου ΙΙ (απόρριψη της Η0 και ανίχνευση στατιστικών διαφορών ενώ στην πραγματικότητα δεν υπάρχουν).

Οι δύο αυτές προϋποθέσεις ελέγχονται στο output της συνάρτησης my_ANOVA που χρησιμοποιήσαμε. Πιο συγκεκριμένα, παρατηρούμε ότι:

Η προϋπόθεση της κανονικότητας των υπολοίπων δεν τηρείται (Δοκιμασία Lilliefors, D = 0,181, p = 0,011), γεγονος που υποδεικνύει πως μόνο η ομάδα δεν εξηγεί αποτελεσματικά τη μεταβλητότητα των Πωλήσεων.

Η προϋπόθεση της ομοιογένειας δεν είναι σαφώς αποδεκτή καθώς η δοκιμασία Levene δεν την απορρίπτει (F(2, 28) = 1.489, p = 0.243) αλλά απορρίπτεται από τη δοκιμασία Bartlett (χ2(2) = 8.214, p = 0.016). Ωστόσο, παρατηρούμε πως η σχέση της μεγαλύτερης με τη μικρότερη διακύμανση δεν ξεπερνάει το 4. Καθώς: “A rule of thumb for heteroscedasticity is that the maximum group variance can be as much as 4 times the minimum group variance without too much damage to your analysis.” (https://stats.stackexchange.com/questions/135232/bartletts-test-vs-levenes-test), αποδεχόμαστε την υπόθεση της ομοιογένειας.

Συμπερασματικά, μπορούμε να δεχθούμε την ύπαρξη διαφορών μεταξύ των τριών ομάδων, με την επισήμανση πως υπάρχουν και άλλοι παράγοντες που επεξηγούν (και ίσως εξαφανίζουν) αυτήν την διαφορά.

Κάποια επιπλέον διαγράμματα με τα οποία αξιολογείται η ποιότητα του γραμμικού μοντέλου δίνονται με την εντολή plot(fit) η οποία παράγει 4 διαγράμματα:

layout(matrix(c(1,2,3,4),2,2))
plot(all.ANOVA.output$res.aov)

Εκ των υστέρων (post - hoc) συγκρίσεις.
Το είδος της διαφοροποίησης διευκρινίζεται από τους εκ υστέρων ελέγχους LSD (Least Significance Difference) και Tukey’s, σύμφωνα με τους οποίους:

• Οι ομάδες Β, Γ είναι διαφορετικές από την ομάδα Α.

• Η ομάδα Β και Γ είναι συγγενείς.

Η διαφοποίηση μεταξύ των τριών ομάδων αναπαριστάται γραφικά με ένα διάγραμμα error bar:

library(gplots)
plotmeans(sales~group,xlab="Ομάδα",ylab="Πωλήσεις", main="Πωλήσεις ανά ομάδα (95% C.I.)")

Σχόλιο
Η επιλογή του post hoc ελέγχου είναι σημαντική απόφαση καθώς κάθε μία από τις διαθέσιμες μεθόδους έχει τις δικές της ιδιαιτερότητες. Ο ερευνητής μπορεί να γνωρίζει πως στην περίπτωση ομοιογενών ομάδων χρησιμοποιείται:
- Η Fisher’s LSD μέθοδος όταν ελέγχονται 3 ομάδες.
- Η Tukey’s b όταν υπάρχουν 3 έως 6 ομάδες.
- Η Tukey’s HSD (Honestly Significant Difference) όταν υπάρχουν περισσότερες από 6 ανεξάρτητες ομάδες.
…ενώ αν η ομοιογένεια των ομάδων δεν είναι δεδομένη τότε μπορεί να χρησιμοποιήσει (α) τη δοκιμασία Tamhane T2 αν η υπόθεση της κανονιότητας δεν απορρίπτεται.
(β) τη δοκιμασία Games-Howell.

Μία γραφική σύνοψη των διαθέσιμων εκ των υστέρων ελέγχων δίνεται στο επόμενο σχήμα:

Διάγραμμα
Διάγραμμα

Σημειώσεις
1. Μία στατιστική δοκιμασία χαρακτηρίζεται ως συντηρητική (conservative) όταν τείνει να έχει μεγαλύτερη τιμή στατιστικής σημαντικότητας p από την πραγματική, δυσχεραίνονται με τον τρόπο αυτό την απόρριψη της Η0.
2. Πηγή διαγράμματος: https://www.ijntse.com/upload/1447070311130.pdf

2.2 Ανάλυση διακύμανσης όταν δεν τηρούνται οι προϋποθέσεις της ομοιογένειας και της κανονικότητας των υπολοίπων

Όταν μία από τις δύο απαραίτητες προϋποθέσεις της παραμετρικής ανάλυσης διακύμανσης (κανονικότητα και ομοιογένεια) δεν είναι αποδεκτή τότε ο ερευνητής πρέπει να επιλέξει μία διαφορετική μέθοδο για να αναλύσει τα δεδομένα του. Παρακάτω παρουσιάζονται δύο διαδικασίες που μπορεί να χρησιμοποιηθούν αντί της παραμετρικής ANOVA.

2.2.1 Ανάλυση διακύμανσης σε ανομοιογενείς ομάδες

Στην περίπτωση όπου η κανονικότητα των υπολοίπων είναι δεδομένη αλλά η ομοιογένεια δεν γίνεται αποδεκτή τότε προτείνεται η εφαρμογή της δοκιμασίας Welch, με την εντολή:

library(onewaytests)
welch.test(sales ~ group, anova.df)
## 
##   Welch's Heteroscedastic F Test (alpha = 0.05) 
## ------------------------------------------------------------- 
##   data : sales and group 
## 
##   statistic  : 6.678925 
##   num df     : 2 
##   denom df   : 15.86982 
##   p.value    : 0.007861148 
## 
##   Result     : Difference is statistically significant. 
## -------------------------------------------------------------

2.2.2 Μη κανονική κατανομή και ομοιογένεια. Δοκιμασία Kruskal Wallis

Στην περίπτωση όπου η κανονικότητα των ανεξάρτητων ομάδων δεν είναι δεδομένη αλλά είναι γνωστό πως τα δύο δείγματα προέρχονται από κατανομή με όμοια μορφή τότε η λύση είναι η εφαρμογή της μη παραμετρικής μεθόδου Kruskal Wallis.

Αποφασίζουμε τη χρήση της δοκιμασίας Kruskal Wallis στις εξής περιπτώσεις:
• Η μεταβλητή δεν είναι συνεχής αλλά διατακτική (π.χ. παίρνει τις τιμές 1, 2, 3, 4, 5 όπως σε ένα ερωτηματολόγιο) και έχουμε μικρό δείγμα.
• Υπάρχει μεγάλη απόκλιση από την κανονική κατανομή.

Η δοκιμασία Kruskal Wallis πρακτικά είναι εφαρμόσιμη σε κάθε περίπτωση καθώς δεν προϋποθέτει όμοιο μέγεθος ομάδων ή οποιαδήποτε υπόθεση για την κατανομή της εξαρτημένης μεταβλητής.

Ωστόσο, αν είναι όμοιες οι κατανομές τότε ελέγχεται η υπόθεση πως συνολικά οι κατανομές διαφέρουν:

Η0: Οι 3 ανεξάρτητες ομάδες έχουν κατανομή που δεν διαφέρει σημαντικά,

έναντι της:

Η1: όχι η Η0.

Αν είναι ανόμοιες οι κατανομές τότε ελέγχεται η υπόθεση ότι οι διάμεσοι διαφέρουν.

Η0: δΑ = δΒ = δΓ.

έναντι της:

Η1: όχι η Η0.

Παραδείγματα κατανομών με ίδιο ή διαφορετικό σχήμα παρουσιάζονται στο συγκριτικό διάγραμμα που ακολουθεί:

Κατανομές με ίδιο και διαφορετικό σχήμα Πηγή: https://statistics.laerd.com/spss-tutorials/kruskal-wallis-h-test-using-spss-statistics.php

Στην περίπτωσή μας, ένα συγκριτικό διάγραμμα των κατανομών μπορεί να προκύψει με τη συνάρτηση my_hist_multiple_density_sidebyside. Η εφαρμογή της προσφέρει το διάγραμμα:

my_hist_multiple_density_sidebyside(anova.df$sales, anova.df$group, xlab = 'Πωλήσεις')

Για να εφαρμόσουμε τη δοκιμασία Kruskal Wallis αρκεί να εκτελέσουμε την εντολή:

kruskal.test(sales ~ group, anova.df)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  sales by group
## Kruskal-Wallis chi-squared = 10.245, df = 2, p-value = 0.005962

Παρατηρούμε ότι στην περίπτωσή μας οι κατανομές είναι ανόμοιες άρα, ελέγχεται η υπόθεση

Η0: δΑ = δΒ = δΓ.

έναντι της:

Η1: όχι η Η0.

Η υπόθεση Η0 απορρίπτεται όταν το p είναι μικρότερο από το 0,05 (ή άλλο όριο που έχει τεθεί από τον ερευνητή πριν τη συλλογή των δεδομένων).

Παρατηρούμε πως για τα δεδομένα του παραδείγματος η μη παραμετρική δοκιμασία υποστηρίζει πως οι 3 ομάδες έχουν στατιστικά διαφορετική διάμεσο (χ2(2) = 10,245, p = 0,006).

2.3 Εντοπισμός ομοιογενών ομάδων με μη παραμετρική μεθοδολογία

Ο τρόπος που προτείνεται είναι η εφαρμογή της δοκιμασίας Dunn:

library(dunn.test)
dunn.test(sales, group)
##   Kruskal-Wallis rank sum test
## 
## data: sales and group
## Kruskal-Wallis chi-squared = 10.2448, df = 2, p-value = 0.01
## 
## 
##                          Comparison of sales by group                          
##                                 (No adjustment)                                
## Col Mean-|
## Row Mean |          Α          Β
## ---------+----------------------
##        Β |  -2.306923
##          |    0.0105*
##          |
##        Γ |  -3.086791  -0.852508
##          |    0.0010*     0.1970
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2

Εναλλακτικά, μπορεί να χρησιμοποιηθεί η δοκιμασία Conover, η οποία θεωρείται περισσότερο αξιόπιστη από τη δοκιμασία Dunn για μικρά δείγματα (n < 30):

library(conover.test)
conover.test(sales, group)
##   Kruskal-Wallis rank sum test
## 
## data: sales and group
## Kruskal-Wallis chi-squared = 10.2448, df = 2, p-value = 0.01
## 
## 
##                          Comparison of sales by group                          
##                                 (No adjustment)                                
## Col Mean-|
## Row Mean |          Α          Β
## ---------+----------------------
##        Β |  -2.746445
##          |    0.0052*
##          |
##        Γ |  -3.674897  -1.014931
##          |    0.0005*     0.1594
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2

Μία ακόμα επιλογή, είναι ο υπολογισμός της τάξης των παρατηρήσεων και η καταχώρησή τους στην μεταβλητή rank.sales, με την εντολή

anova.df$rank.sales = rank(anova.df$sales)

και η εφαρμογή της δοκιμασίας Tukey HSD στη νέα μεταβλητή rank.sales:

rank.fit = aov(rank.sales ~ group, data = anova.df)
TukeyHSD(rank.fit, conf.level=.95)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = rank.sales ~ group, data = anova.df)
## 
## $group
##          diff        lwr      upr     p adj
## Β-Α  9.163636  0.9078579 17.41941 0.0272569
## Γ-Α 12.550000  4.0999413 21.00006 0.0027805
## Γ-Β  3.386364 -4.8694149 11.64214 0.5737891

Από το μέρος $groups του output είναι φανερό πως τα αποτελέσματα της ομαδοποίησης για τα δεδομένα του παραδείγματος που μελετούμε είναι παρόμοια με αυτά που προκύπτουν από την εφαρμογή της διαδικασίας στις ίδιες τις τιμές της μεταβλητής.

2.4 Καταγραφή ανάλυσης διακύμανσης με έναν παράγοντα.

Στοιχεία που πρέπει να αναφερθούν
• τα στοιχεία των δύο μετρήσεων (πλήθος παρατηρήσεων, η μέση τιμή και η τυπική απόκλιση για κάθε μία από αυτές),
• η τιμή του στατιστικού F (στήλη t στο output) ως απόλυτη τιμή,
• οι 2 βαθμοί ελευθερίας (στήλη df, γραμμή Between Groups και γραμμή Within Groups με αυτήν τη σειρά),
• η στατιστική σημαντικότητα p της διαφοροποίησης από τη στατιστική υπόθεση Η0 (από τη στήλη Sig.).

Μη σημαντική διαφοροποίηση
Η ανάλυση διακύμανσης με έναν παράγοντα χρησιμοποιήθηκε για να ανιχνεύσει τη στατιστική σημαντικότητα των διαφοροποιήσεων μεταξύ των μέσων τιμών των τριών ανεξάρτητων ομάδων. Βρέθηκε πως δεν υπάρχει σημαντική διαφοροποίηση των μέσων τιμών των πωλήσεων μεταξύ των διαφορετικών ομάδων, F(2, 28) = 2,772, p = 0,08.

Σημαντική διαφοροποίηση.
Για την ανίχνευση της διαφοροποίησης των μέσων τιμών των πωλήσεων μεταξύ των τριών ανεξάρτητων ομάδων χρησιμοποιήθηκε η ανάλυση διακύμανσης με έναν παράγοντα. Βρέθηκε πως υπάρχει σημαντική διαφοροποίηση των μέσων τιμών των πωλήσεων μεταξύ των τριών διαφορετικών ομάδων, F(2, 28) = 7,512, p = 0,002.
Η σύγκριση των ομάδων ανά ζεύγη με τη μέθοδο LSD του Fisher κατέδειξε πως η ομάδα Α (M = 62,5, 95% CI [52,3, 72,8]) έχει σημαντικά μικρότερη μέση τιμή (σε επίπεδο σημαντικότητας p < 0,05) από τις ομάδες Β (M = 76,2, 95% CI [69,0, 83,4]) και Γ (M = 80,2, 95% CI [76,7, 83,8]) οι οποίες μεταξύ τους δεν διαφέρουν σημαντικά.

3 Ανάλυση Διακύμανσης με δύο παράγοντες (2-way ANOVA).

Χρησιμοποιούμε ανάλυση διακύμανσης με δύο παράγοντες (2 - way ANOVA) όταν θέλουμε να ελέγξουμε τη σημαντικότητα της επίδρασης δύο παραγόντων στις τιμές μίας συνεχούς μεταβλητής. Αποτελεί τη γενίκευση της ανάλυσης διακύμανσης με ένα παράγοντα. Απαραίτητες προϋποθέσεις είναι (α) η κανονικότητα της κατανομής των τιμών της συνεχούς μεταβλητής σε κάθε ένα συνδυασμό των δύο παραγόντων, (β) η ομοιογένεια των ομάδων (η τυπική απόκλιση δεν διαφέρει σημαντικά μεταξύ όλων των συνδυασμών των δύο παραγόντων), προϋπόθεση που είναι ιδιαίτερα σημαντική όταν δεν υπάρχει ίδιο πλήθος παρατηρήσεων σε κάθε έναν συνδυασμό των παραγόντων (unbalanced design).

Παράδειγμα
Μια επιχείρηση επιχειρώντας να βρει την απαραίτητη διάρκεια εκπαίδευσης των πωλητών της για την προώθηση ενός νέου προϊόντος επιλέγει με τυχαίο τρόπο 31 πωλητές τους οποίους χωρίζει σε τρεις ομάδες (Α, n = 10), (Β, n = 11), (Γ, n = 10), τους εκπαιδεύει για 1, 2 και 3 ημέρες αντίστοιχα και παρακολουθεί τις πωλήσεις τους τον επόμενο μήνα. Στην R, τα δεδομένα καταχωρούνται σε τρεις μεταβλητές, τη μεταβλητή group («Ομάδα εκπαίδευσης») και τιμές 1 (ομάδα Α), 2 (ομάδα Β), 3 (ομάδα Γ), τη μεταβλητή gender όπου καταχωρείται το φύλο του πωλητή με τιμές 0 (Γυναίκα) και 1 (Άνδρας) και τη μεταβλητή perform (Πωλήσεις) όπου καταχωρούνται οι πωλήσεις του κάθε πωλητή. Οι παραπάνω εντολές ορίζουν το πλαίσιο δεδομένων που απαιτείται για την ανάλυση.

perform.A = c(63.3, 68.3, 86.7, 52.8, 75, 58, 69.5, 32.7, 60.9, 58.2)
A.gender = c(1, 1, 1, 0, 0, 0, 1, 0, 0, 0)
perform.B = c(72.8, 88.2, 80.8, 71.3, 81.5, 47.6, 81, 81.4, 83, 76, 74.5)
B.gender = c(1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1)
perform.C = c(82.3, 89.7, 81, 85.1, 74.1, 75.9, 74.7, 81.1, 76.4, 81.8)
C.gender = c(0, 0, 1, 0, 1, 1, 1, 1, 1, 0)

sales = c(perform.A, perform.B, perform.C)
gender = c(A.gender, B.gender, C.gender)
group = c(rep(1, 10), rep(2, 11), rep(3, 10))
gender = factor(gender, levels = c(0, 1), labels = c("Γυναίκα", "Άνδρας"))
group = factor(group, levels = c(1, 2, 3), labels = c("Α", "Β", "Γ"))

anova.df = data.frame(group, gender, sales)

Hmisc::label(anova.df$group) = 'Ομάδα'
Hmisc::label(anova.df$gender) = 'Φύλο'
Hmisc::label(anova.df$sales) = 'Πωλήσεις'

print(anova.df)
##    group  gender sales
## 1      Α  Άνδρας  63.3
## 2      Α  Άνδρας  68.3
## 3      Α  Άνδρας  86.7
## 4      Α Γυναίκα  52.8
## 5      Α Γυναίκα  75.0
## 6      Α Γυναίκα  58.0
## 7      Α  Άνδρας  69.5
## 8      Α Γυναίκα  32.7
## 9      Α Γυναίκα  60.9
## 10     Α Γυναίκα  58.2
## 11     Β  Άνδρας  72.8
## 12     Β Γυναίκα  88.2
## 13     Β Γυναίκα  80.8
## 14     Β  Άνδρας  71.3
## 15     Β Γυναίκα  81.5
## 16     Β  Άνδρας  47.6
## 17     Β Γυναίκα  81.0
## 18     Β Γυναίκα  81.4
## 19     Β Γυναίκα  83.0
## 20     Β Γυναίκα  76.0
## 21     Β  Άνδρας  74.5
## 22     Γ Γυναίκα  82.3
## 23     Γ Γυναίκα  89.7
## 24     Γ  Άνδρας  81.0
## 25     Γ Γυναίκα  85.1
## 26     Γ  Άνδρας  74.1
## 27     Γ  Άνδρας  75.9
## 28     Γ  Άνδρας  74.7
## 29     Γ  Άνδρας  81.1
## 30     Γ  Άνδρας  76.4
## 31     Γ Γυναίκα  81.8

Στο σημείο αυτό απαιτούνται οι εντολές:

contrasts(anova.df$gender) = contr.helmert
contrasts(anova.df$group) = contr.helmert

που ορίζουν τις κατάλληλες συγκρίσεις για τους παράγοντες «Φύλο» και «Ομάδα». Ο λόγος για τον οποίον απαιτούνται οι παραπάνω εντολές είναι πως η ανάλυση διακύμανσης είναι μία διαδικασία που βασίζεται στη γραμμική παλινδρόμηση.

Στη γραμμική παλινδρόμηση το γραμμικό μοντέλο που χρησιμοποιείται μπορεί να έχει διάφορες μορφές και είναι σημαντικό για τον ερευνητή να διαμορφωθεί αυτό κατάλληλα, ώστε οι υπολογισμοί που θα ακολουθήσουν στα πλαίσια της ανάλυσης διακύμανσης να είναι αξιόπιστοι.

Ελέγχουμε τις στατιστικές υποθέσεις:

Η0,1: Ο παράγοντας ομάδα δεν επιδρά σημαντικά στις τιμές των πωλήσεων.

Η0,2: Ο παράγοντας φύλο δεν επιδρά σημαντικά στις τιμές των πωλήσεων.

Η0,3: Ο συνδυασμός των παραγόντων φύλου και ομάδας δεν επιδρά σημαντικά στις τιμές των πωλήσεων.

3.1 Υλοποίηση δοκιμασίας

Για να την υλοποιήσουμε τον έλεγχο διακύμανσης αρκεί να εκτελέσουμε τις εντολές:

my.ANOVA.output = my_ANOVA(data = anova.df, model = sales ~ group + gender + group:gender, show.eta.square = TRUE)
my.ANOVA.output$htmlTable
ANOVA Results for Πωλήσεις
SS df F p η2 η2p
(Intercept) 156735.236 1 2054.498 <0.001
group 1380.314 2 9.047 0.001 0.279 0.42
gender 39.866 1 0.523 0.476 0.008 0.02
group:gender 1261.4 2 8.267 0.002 0.255 0.398
Residuals 1907.221 25
LSD Test for variable: Ομάδα
Fisher’s LSD is a series of pairwise t-tests, with each test using the mean squared error from the significant ANOVA as its pooled variance estimate (and naturally taking the associated degrees of freedom). It is a valid test, in case where the ANOVA is significant.
N Πωλήσεις Groups
Α 10 62.54 b
Β 11 76.191 a
Γ 10 80.21 a
Tukeys HSD Post Hoc Test for variable: Ομάδα
The Tukey HSD is not a Post Hoc Test per se. It may provide valuable results even if ANOVA is not significant.
Difference Lower Upper p
Β-Α 13.651 2.061 25.24 0.018
Γ-Α 17.67 5.808 29.532 0.003
Γ-Β 4.019 -7.57 15.609 0.671
LSD Test for variable: Φύλο
Fisher’s LSD is a series of pairwise t-tests, with each test using the mean squared error from the significant ANOVA as its pooled variance estimate (and naturally taking the associated degrees of freedom). It is a valid test, in case where the ANOVA is significant.
N Πωλήσεις Groups
Άνδρας 14 72.657 a
Γυναίκα 17 73.435 a
Tukeys HSD Post Hoc Test for variable: Φύλο
The Tukey HSD is not a Post Hoc Test per se. It may provide valuable results even if ANOVA is not significant.
Difference Lower Upper p
Άνδρας-Γυναίκα -0.778 -10.411 8.854 0.87
Standard Deviation of Πωλήσεις among levels of all level combinations of group, gender.
Level Α.Γυναίκα Β.Γυναίκα Γ.Γυναίκα Α.Άνδρας Β.Άνδρας Γ.Άνδρας
Group Size 6 7 4 4 4 6
SD 13.80 3.60 3.62 10.20 12.70 3.09
Levene test of variance equality (homogeneity) of Πωλήσεις over all level combinations of group, gender.
F(5, 25) = 0.865, p = 0.518. Homogeneity hypothesis is confirmed.
Bartlett’s test of variance equality (homogeneity) of Πωλήσεις over all level combinations of group, gender.
If you have strong evidence that your data do in fact come from a normal, or nearly normal, distribution, then Bartlett’s test has better performance.
c2(5) = 16.859, p = 0.005. Homogeneity hypothesis is not supported!
Normality test for model residuals.
If the main goal of an ANOVA is to see whether or not certain effects are significant, then the assumption of normality of the residuals is only required for small samples, thanks to the central limit theorem. With sample sizes of a few hundred participants even extreme violations of the normality assumptions are unproblematic. So mild violations of this assumptions are usually no problem with sample sizes exceeding 30.
Statistic Description
Skewness & Kurtosis
  Skewness -0.645 Normal distribution has skew = 0. For a unimodal distribution, negative skew commonly indicates that the tail is on the left side of the distribution, and positive skew indicates that the tail is on the right.
  Kurtosis 5.322 Normal distribution has kurtosis = 3. Values over 3 indicates a platykurtic distribution and values less than 3 indicates a leptokurtic distribution.
Normality Tests. Η0: This sample is from a normal distribution vs Η1: Not the Η0.
  Shapiro–Wilk W = 0.904, p = 0.009 This test is more appropriate method for small sample sizes (<50 samples) although it can also be handling on larger sample size.
  Lilliefors D = 0.195, p = 0.004 The Lilliefors test uses the same calculations as the Kolmogorov-Smirnov test, but it is more conservative in the sense that the Lilliefors Test is less likely to show that data is normally distributed.

Παρατηρούμε ότι στο output υπάρχουν οι δείκτες η2 που δείχνουν την επιρροή των παραγόντων «Φύλο» και «Ομάδα», όπως και του συνδυασμού τους στην εξαρτημένη μεταβλητή των πωλήσεων.

Έλεγχος προϋποθέσεων
Πριν την ανάγνωση του output πρέπει να τεκμηριωθεί η αξιοπιστία των στατιστικών που υπολογίστηκαν. Η κανονικότητα της κατανομής των υπολοίπων απορρίπτεται (W = 0,904, p = 0,009), ωστόσο η υπόθεση πως το δείγμα μας προέρχεται από ομοιογενή πληθυσμό δεν απορρίπτεται F(5, 25) = 0,866, p = 0,518.

Αποτελέσματα
Στον πίνακα “ANOVA Results for Πωλήσεις” κρύβονται τα βασικά αποτελέσματα της δοκιμασίας. Πιο συγκεκριμένα, αναφορικά με το παράδειγμά μας, παρατηρούμε πως:
• Ως προς τον παράγοντα ομάδα (group), είναι F(2, 25) = 9,047, p = 0,001.
• Ως προς τον παράγοντα φύλο (gender), είναι F(1, 25) = 0,523, p = 0,476.
• Ως προς το συνδυασμό ομάδας και φύλου (group*gender) είναι F(2, 25) = 8,267, p = 0,002.

Από το output της 2 - way ANOVA με δύο παράγοντες Α και Β ερμηνεύουμε πρώτα τη γραμμή για τη επιρροή του συνδυασμού των παραγόντων \(Α*Β\) στην εξαρτημένη μεταβλητή. Αν υπάρχει στατιστικά σημαντική διαφοροποίηση ως προς το συνδυασμό \(A*B\) τότε η επιρροή μίας εκ των δύο μεμονωμένων διαφοροποιήσεων ως προς τον παράγοντα Α ή τον παράγοντα Β αρμόζει να καταγραφούν ως σημαντική επιρροή του παράγοντα στο μέρος του δείγματος που ορίζεται από την τιμή 0 του άλλου παράγοντα! Δηλαδή, αν \(Α*Β\) είναι σημαντικός παράγοντας διαφοροποίησης και ο παράγοντας Α καταγράφεται ως σημαντικός από τα αποτελέσματα της ANOVA τότε αυτή η στατιστικώς σημαντική επιρροή είναι ορθό να καταγραφεί ως μόνο για τις παρατηρήσεις με Β = 0.

Οι παραπάνω στατιστικοί υπολογισμοί ερμηνεύονται ως εξής: Η στατιστική υπόθεση Η0,3 απορρίπτεται, F(2, 25) = 8,267, p = 0,002, η2μερικό = 0,398, δηλαδή οι πωλήσεις των τριών ομάδων διαφοροποιούνται στις γυναίκες με διαφορετικό τρόπο από ότι στους άνδρες,

Ως προς την υπόθεση Η0,1 μπορούμε να καταγράψουμε πως ο παράγοντας ομάδα (group) είναι σημαντικός παράγοντας επιρροής της μεταβλητής πωλήσεις στο μέρος του δείγματος που ορίζεται από τη συνθήκη gender = 0, δηλαδή για τις γυναίκες!

Ως προς την υπόθεση Η0,2 δεν καταγράφεται στατιστικά σημαντική επίδραση του φύλου στις τιμές των πωλήσεων.

Απομένει να διευκρινιστεί το είδος της διαφορετικότητας των τριών ομάδων μεταξύ ανδρών και γυναικών. Για το σκοπό αυτό, αρκεί η εκτέλεση της συνάρτησης my_explore_vars_by_group:

my_explore_vars_by_group(anova.df, 'sales', "group")
Descriptive Statistics of Πωλήσεις among levels of Ομάδα
Group N M (SD) 95% C.I.
1 10 62.5 (14.3) 52.3 - 72.8
2 11 76.2 (10.7) 69 - 83.4
3 10 80.2 (4.98) 76.6 - 83.8
Total 31 73.1 (12.8) 68.4 - 77.8
Tukeys HSD Post Hoc Test for variable: Πωλήσεις
The Tukey HSD is not a Post Hoc Test per se. It may provide valuable results as an independent test as well.
Difference Lower Upper p
2-1 13.651 2.061 25.24 0.018
3-1 17.67 5.808 29.532 0.003
3-2 4.019 -7.57 15.609 0.671
htmlTable::htmlTable(by(anova.df$sales, anova.df[,c("group", "gender")], my_mean_sd), caption = '<b>Descriptive Statistics of Πωλήσεις among all level combinations of Ομάδα and Φύλο.</b>')
Descriptive Statistics of Πωλήσεις among all level combinations of Ομάδα and Φύλο.
gender
Γυναίκα Άνδρας
group
  Α M = 56.3 (13.8), 95% C.I. 45.3 - 67.3 M = 72.0 (10.2), 95% C.I. 62.0 - 81.9
  Β M = 81.7 (3.60), 95% C.I. 79.0 - 84.4 M = 66.6 (12.7), 95% C.I. 54.1 - 79.0
  Γ M = 84.7 (3.62), 95% C.I. 81.2 - 88.3 M = 77.2 (3.09), 95% C.I. 74.7 - 79.7

Ο παραπάνω κώδικας εφαρμόζει τη συνάρτηση my_mean_sd σε κάθε συνδυασμό των παραγόντων «Ομάδα» και «Φύλο»

Παρατηρούμε πως οι πωλήσεις των γυναικών αυξάνονται με το πλήθος των ημερών ενώ οι πωλήσεις των ανδρών μειώνονται τη 2η ημέρα για να αυξηθούν ξανά την 3η.

Ως προς τη διαφοροποίηση μεταξύ ανδρών και γυναικών απομένει ο εντοπισμός των κρίσιμων διαφορών στις ημέρες εκπαίδευσης για κάθε ένα φύλο ξεχωριστά.

Ο κανόνας με τον οποίο γίνεται η καταγραφή των επιμέρους αποτελεσμάτων είναι απλός: αν δύο συνδυασμοί έχουν 95% διαστήματα εμπιστοσύνης που δεν επικαλύπτονται τότε υπάρχει στατιστικά σημαντική διαφοροποίηση μεταξύ των δύο μέσων τιμών.

Καθώς από τον πίνακα προκύπτει πως για κάθε ένα φύλο ξεχωριστά, το 95% διάστημα εμπιστοσύνης για τις ομάδες Α, Β και Γ είναι αντίστοιχα:
• Γυναίκες: Α(45,3 έως 67), Β(79 έως 84) και Γ(81,2 έως 88)
• Άνδρες: Α(62 έως 82), Β(54,1 έως 79) και Γ(74,7 έως 80)
Συνάγουμε πως:
• Στις γυναίκες, η ομάδα Α έχει σημαντικά μικρότερες πωλήσεις από τις ομάδες Β και Γ οι οποίες δεν διαφέρουν σημαντικά μεταξύ τους.
• Στους άνδρες οι 3 ομάδες δεν διαφέρουν μεταξύ τους.

Μπορούμε να πάρουμε και ένα συγκριτικό διάγραμμα των παραπάνω τιμών αξιοποιώντας τη συνάρτηση my_interaction_plot. Η εφαρμογή της προσφέρει το παρακάτω διάγραμμα:

my_interaction_plot(anova.df$sales, anova.df$group, anova.df$gender, legend.label = 'Φύλο', xlab="Ομάδα",  ylab="Πωλήσεις",  main="Αλληλεπίδραση φύλου και ομάδας στις Πωλήσεις")

Συνοπτικά, από το σύνολο των αποτελεσμάτων του παραδείγματος φαίνεται πως η επιχειρηματική απόφαση σχετικά με την εκπαίδευση των στελεχών της επιχείρησης που θα συνδυάζει την αποτελεσματικότητα με την οικονομία στους πόρους, είναι η εκπαίδευση 2 ημερών για τις γυναίκες και η εκπαίδευση μίας ημέρας για τους άνδρες.

3.2 Συγκρίσεις (contrasts) ενός γραμμικού μοντέλου

Η παρατήρηση αυτή έχει σκοπό να διευκρινίσει την αναγκαιότητα των εντολών

contrasts(anova.df$gender) = contr.helmert
contrasts(anova.df$group) = contr.helmert
που απαιτήθηκαν για την ορθή εφαρμογή της ανάλυσης διακύμανσης με δύο παράγοντες. Οι παραπάνω εντολές ορίζουν τις συγκρίσεις (contrasts ) μεταξύ των επιπέδων των παραγόντων «Φύλο (gender)» και «Ομάδα (group)» στην ανάλυση διακύμανσης, δηλαδή τον τρόπο με τον οποίο θα αναπαρασταθούν στην εξίσωση που προσαρμόζεται στα δεδομένα, τα επίπεδα των δύο παραγόντων. Η προκαθορισμένη επιλογή των συγκρίσεων για την R είναι η contr.treatment η οποία δημιουργεί απλές ψευδομεταβλητές (dummy variables) για κάθε μία μεταβλητή, δηλαδή τοποθετεί στο γραμμικό μοντέλο τους παράγοντες ως τους γραμμικούς συνδυασμούς:
Φύλο = Χ1
Ομάδα = Υ1 + Υ2

όπου οι Χ1, Υ1, Υ2, λαμβάνουν τις τιμές 0 και 1 και ορίζουν τις τιμές των παραγόντων σύμφωνα με τον παρακάτω πίνακα συγκρίσεων:

Κωδικοποίηση παραγόντων σύμφωνα με contr.treatment

Παράγοντας Φύλο Μεταβλητή X1
Άνδρας 1
Γυναίκα 0
Παράγοντας Ομάδα Μεταβλητή Y1 Μεταβλητή Y2
Α 1 0
Β 0 1
Γ 0 0

Με την παραπάνω κωδικοποίηση, γίνεται προσπάθεια προσαρμογής στα δεδομένα, του μοντέλου

sales = Χ1 + Υ1 + Υ2

και επιτυγχάνεται η σύγκριση μεταξύ Ανδρών (Χ1 = 1) και Γυναικών (Χ1 = 0) όπως και η σύγκριση των ομάδων Α (Υ1 = 1, Υ2 = 0) και Β (Υ1 = 0, Υ2 = 1) με την ομάδα Γ (Υ1 = 0, Υ2 = 0).

Ενδεικτικά, στα δεδομένα του παραδείγματος, μπορούμε να δούμε τις συγκρίσεις και το αποτέλεσμα της γραμμικής παλινδρόμησης, με τις εντολές:

contrasts(anova.df$gender) = contr.treatment
print(contrasts(anova.df$gender))
##         2
## Γυναίκα 0
## Άνδρας  1
contrasts(anova.df$group) = contr.treatment
print(contrasts(anova.df$group))
##   2 3
## Α 0 0
## Β 1 0
## Γ 0 1
lm.model = lm(sales ~ group + gender, anova.df)
print(lm.model)
## 
## Call:
## lm(formula = sales ~ group + gender, data = anova.df)
## 
## Coefficients:
## (Intercept)       group2       group3      gender2  
##      63.574       13.557       18.187       -2.584

Η παραπάνω προκαθορισμένη κωδικοποίηση είναι η πιο απλή και είναι κατάλληλη για χρήση σε γραμμικά μοντέλα σε εφαρμογές γραμμικής παλινδρόμησης. Ωστόσο, στην ανάλυση διακύμανσης απαιτείται οι παραπάνω στήλες - διανύσματα να είναι ένα ορθογώνιο σύνολο διανυσμάτων, δηλαδή να έχουν τις ιδιότητες:

(α) το άθροισμα των στοιχείων τους να είναι ίσο με 0.
(β) το εσωτερικό γινόμενο οποιωνδήποτε δύο από αυτές να είναι ίσο με 0.

Οι παραπάνω δύο ιδιότητες είναι απαραίτητες γιατί, στην περίπτωση αυτή, το συνολικό άθροισμα τετραγώνων που θα αποδοθεί σε κάθε παράγοντα είναι δυνατό να διαχωριστεί στα μέρη που ορίζονται από κάθε μία σύγκριση και έτσι να ελεγχθούν οι επιμέρους υποθέσεις που περιγράφονται από τις συγκρίσεις.

Πιο συγκεκριμένα, αν και υπάρχει δυνατότητα εισαγωγής μη-ορθογωνίων αντιθέσεων, αυτές έχουν μικρότερη ισχύ από τις ορθογώνιες υπό την έννοια πως έχουν μεγαλύτερη ικανότητα ανίχνευσης σημαντικών διαφορών μεταξύ των μέσων των ομάδων, με μικρότερη πιθανότητα ψευδώς θετικών (λάθη τύπου Ι) και ψευδώς αρνητικών (λάθη τύπου II) αποτελεσμάτων.

Στην γλώσσα R, ο τρόπος για να επιβάλουμε τις παραπάνω συνθήκες είναι απλά να ορίσουμε μία κατάλληλη διαδικασία υπολογισμού συγκρίσεων (contrasts) και ο απλούστερος τρόπος για να το πετύχουμε αυτό είναι με τις συγκρίσεις κατά Helmert, οι οποίες επιτυγχάνουν τη σύγκριση ενός επιπέδου με τη μέση τιμή όλων των προηγούμενων επιπέδων (2ο επίπεδο με το 1ο και το 3ο επίπεδο με τη μέση τιμή του 1ου και του 2ου).

Ενδεικτικά, η προσαρμογή γραμμικού μοντέλου με αυτές τις συγκρίσεις θα οδηγήσει στα εξής αποτελέσματα:

contrasts(anova.df$gender) = contr.helmert
contrasts(anova.df$group) = contr.helmert
lm.model = lm(sales ~ group + gender, anova.df)
print(contrasts(anova.df$gender))
##         [,1]
## Γυναίκα   -1
## Άνδρας     1
print(contrasts(anova.df$group))
##   [,1] [,2]
## Α   -1   -1
## Β    1   -1
## Γ    0    2
print(lm.model)
## 
## Call:
## lm(formula = sales ~ group + gender, data = anova.df)
## 
## Coefficients:
## (Intercept)       group1       group2      gender1  
##      72.863        6.778        3.803       -1.292

Στην R, οι συγκρίσεις (contrasts) καταχωρούνται είτε σε κάθε μία μεταβλητή ξεχωριστά, όπως για παράδειγμα

contrasts(gender) = matrix(c(0.7071, -0.7071, -0.7071, 0.7071), 2, 2) # ή contrasts(gender) = contr.poly

είτε ως μεταβλητή περιβάλλοντος για χρήση σε κάθε μία εφαρμογή όπου γίνεται προσαρμογή μοντέλων όπως στις συναρτήσεις aov ή lm,

options(contrasts = c("contr. poly ", "contr.poly"))

Στη τελευταία περίπτωση, το πρώτο όρισμα αφορά τις συγκρίσεις που θα ορίζονται στις ποιοτικές μεταβλητές και το δεύτερο τις συγκρίσεις που θα ορίζονται για τις κατηγορικές μεταβλητές. Στον παρακάτω πίνακα παρουσιάζονται οι διαθέσιμες προκαθορισμένες συγκρίσεις τόσο στην R όσο και στο SPSS σε αντιπαράθεση μεταξύ τους.

R SPSS
contr.SAS SIMPLE (Κάθε επίπεδο του παράγοντα συγκρίνεται με το τελευταίο)
contr.treatment Ίδιο με την επιλογή SIMPLE μόνο που κάθε επίπεδο του παράγοντα συγκρίνεται με το πρώτο και όχι με το τελευταίο.
contr.sum DEVIATION (Απόκλιση από το γενικό μέσο όρο)
contr.poly POLYNOMIAL (Πολυωνυμικές συγκρίσεις) ή DIFFERENCE (Διαφορές των αντίστροφων συγκρίσεων Helmert)
contr.helmert HELMERT (Συγκρίσεις Helmert) ή REPEATED (Σύγκριση κάθε επιπέδου με τη μέση τιμή των προηγουμένων επιπέδων)
contrasts() SPECIAL (Συγκρίσεις ορισμένες από το χρήστη)

Χρήσιμες αναφορές:
1.https://www.ibm.com/support/knowledgecenter/fr/SSLVMB_20.0.0/com.ibm.spss.statistics.help/syn_unianova_contrast.htm
2. http://rstudio-pubs-static.s3.amazonaws.com/65059_586f394d8eb84f84b1baaf56ffb6b47f.html
3. http://www.unc.edu/courses/2006spring/ecol/145/001/docs/lectures/lecture26.htm

3.3 Παρατηρήσεις

3.3.1 Οι τρεις τρόποι υπολογισμού των τετραγωνικών αθροισμάτων

Στην ανάλυση διακύμανσης υπάρχουν τρεις μεθόδοι για την κατανομή του συνολικού αθροίσματος τετραγώνων σε συνιστώσες που σχετίζονται με τις διαφορετικές πηγές διακύμανσης των δεδομένων. Αυτές οι τρεις μέθοδοι οδηγούν σε όμοια συμπεράσματα όταν οι υποομάδες έχουν ανάλογο μέγεθος (balanced design) αλλά δύναται να διαφέρουν πολύ όταν υπάρχουν σημαντικά διαφορετικές ομάδες ως προς το μέγεθος (unbalanced design)

Άθροισμα τύπου I
Τα αθροίσματα τύπου Ι ακολυθούν ιεραρχική προσέγγιση, όπου κάθε μεταβλητή πρόβλεψης προστίθεται στο μοντέλο με προκαθορισμένη σειρά και το άθροισμα των τετραγώνων που σχετίζονται με κάθε μεταβλητή υπολογίζεται με βάση τη μοναδική συνεισφορά αυτής της μεταβλητής στο μοντέλο, αφού ληφθούν υπόψη όλες οι μεταβλητές που έχουν είδη εισαχθεί στο μοντέλο. Τα αθροίσματα τύπου Ι μπορούν να χρησιμοποιηθούν όταν η σειρά με την οποία εισάγονται οι μεταβλητές πρόβλεψης στο μοντέλο έχει νόημα και μπορεί να ερμηνευθεί βάσει κάποιας θεωρητικής προσέγγισης.

Άθροισμα τύπου II
Για τον υπολογισμό αθροισμάτων τύπου ΙΙ, κάθε επεξηγηματική μεταβλητή προστίθεται στο μοντέλο με συγκεκριμένη σειρά βάσει ενός κριτηρίου όπως η στατιστική σημασία. Το άθροισμα των τετραγώνων που σχετίζονται με κάθε μεταβλητή υπολογίζεται με βάση τη μοναδική συνεισφορά αυτής της μεταβλητής στο μοντέλο, αφού ληφθούν υπόψη όλες οι άλλες μεταβλητές στο μοντέλο.

Άθροισμα τύπου III
Με τα αθροίσματα τύπου III υπολογίζονται τα αθροίσματα τετραγώνων με βάση τη μοναδική συνεισφορά της κάθε μίας μεταβλητής στο μοντέλο, αφού ληφθούν υπόψη όλες οι άλλες μεταβλητές στο μοντέλο. Τα αθροίσματα τύπου III είναι κατάλληλα για μοντέλα με μη ισορροπημένα δεδομένα και θεωρούνται πιο αξιόπιστα όταν υπάρχουν αποκλίσεις απο τις προϋποθέσεις εφαρμογής της ανάλυσης διακύμανσης. Ειδικότερα, θεωρούνται πιο κατάλληλα αν εξετάζεται επιπλέον και αλληλεπίδραση μεταξύ των παραγόντων.

3.3.2 Ερμηνεία της εξάρτησης δύο παραγόντων

Στην δοκιμασία 2-way ANOVA, αν Α και Β οι δύο παράγοντες και ο συνδυασμός \(Α*Β\) δεν προκύπτει στατιστικά σημαντικός παράγοντας επιρροής των τιμών της εξαρτημένης μεταβλητής τότε η ερμηνεία του output ως προς τους παράγοντες Α και Β γίνεται με τον ίδιο ακριβώς τρόπο με αυτόν που χρησιμοποιείται στη 1-way ANOVA: Αν η σημαντικότητα p είναι μικρότερη από το 0,05 (ή όποιο άλλο όριο απόρριψης έχει τεθεί) τότε ο παράγοντας θεωρείται πως επηρεάζει τις τιμές της εξαρτημένης μεταβλητής και στη συνέχεια ο ερευνητής μπορεί να χρησιμοποιήσει τους διαθέσιμους post - hoc ελέγχους για τη διευκρίνιση της ομοιογένειας των ανεξάρτητων ομάδων. Αν αντίθετα, όπως στο παράδειγμά μας, ο συνδυασμός \(Α*Β\) είναι στατιστικά σημαντικός παράγοντας επιρροής των τιμών της εξαρτημένης μεταβλητής τότε η μόνη ορθή ερμηνεία των υπολοίπων αποτελεσμάτων είναι αυτή που περιγράφεται στο κείμενο του παραδείγματος. Πράγματι, αν:

Yijk = μ + ai + bj + abij + eijk

είναι το μοντέλο πρόβλεψης της k -οστής τιμής της εξαρτημένης μεταβλητής Y στο συνδυασμό του i επιπέδου του παράγοντα Α και του j επιπέδου του παράγοντα Β (μ η γενική μέση τιμή της Υ) τότε η ύπαρξη στατιστικά σημαντικής επίδρασης τους συνδυασμού Α*Β στην τιμή του Υ αντιστοιχεί στο γεγονός πως οι παράμετροι abij είναι (εν γένει) διάφοροι του μηδενός και άρα πως για τις τιμές της Y που βρίσκονται στο συνδυασμό των επιπέδων i και j των παραγόντων Α και Β δεν είναι σαφές αν υπάρχει και τι μέγεθος έχει η μεμονωμένη επίδραση του κάθε ενός παράγοντα ξεχωριστά. Εξαίρεση σε αυτή την ασάφεια είναι η περίπτωση όπου ο ένας παράγοντας είναι 0 περίπτωση στην οποία το μοντέλο απλοποιείται σε

Yijk = μ + ai + eijk ή Yijk = μ + bj + eijk,

περιέχει μόνο τον μη μηδενικό παράγοντα και αυτός είναι ο λόγος για τον ιδιαίτερο τρόπο ερμηνείας των επιρροών των μεμονωμένων παραγόντων στην περίπτωση στατιστικά σημαντικής επίδρασης του συνδυασμού τους.

3.3.3 Eta Square (η2) και Partial Eta Square (μερικό η2)

Στα πλαίσια οποιαδήποτε διαδικασίας ANOVA, υπάρχει δυνατότητα να υπολογι­στεί και να αναφερθεί για κάθε έναν από τους παράγοντες ή συνδυασμούς πα­ραγόντων ο δείκτης η2 (Eta Square) και ο μερικός δείκτης η2 (Partial Eta Square) ως μέτρα της επιρροής του παράγοντα ή του συνδυασμού παραγόντων στην εξαρτημένη μεταβλητή. Φανερά, έχει νόημα η αναφορά σε ένα τέτοιο στατιστικό μόνο όταν ο παράγοντας ή ο συνδυασμός παραγόντων προκύπτει στατιστικά σημαντικός. Ο δείκτης η2 (Eta Square) ορίζεται ως ο λόγος της διακύμανσης που εξηγείται από τον παράγοντα ή το συνδυασμό παραγόντων προς τη συνολική διακύμανσης του δείγματος. (Tabachnick & Fidell 2001, σελ. 54 - 55, Thompson 2006, σελ. 317 - 319).

Ο μερικός δείκτης η2 (Partial Eta Square) ορίζεται ως ο λόγος της διακύμανσης που εξηγείται από τον παράγοντα ή το συνδυασμό παραγόντων προς τη διακύ­μανση αυτή αυξημένη προς τη διακύμανση του σφάλματος του μοντέλου.

Καθώς στη δοκιμασία 1 - way ANOVA υπάρχει μόνο ένας παράγοντας και ισχύει SSπαράγοντα + SSσφάλμα = SSσύνολο, συμπεραίνουμε πως σε αυτήν την περίπτωση ελέγ­χου το Eta Square (η2) είναι ίσο με το Partial Eta Square (μερικό η2). Το ίδιο δεν συμβαίνει σε περισσότερο πολύπλοκες αναλύσεις όπου το Partial Eta Square (μερικό η2) αναμένεται να είναι μεγαλύτερο από το Eta Square (η2), ενδεχομένως και δραματικά μεγαλύτερο. Επιπλέον, το άθροισμα των Eta Square (μερικό η2) για όλους τους παράγοντες πρέπει να έχει άθροισμα 1 ενώ δεν υπάρχει ανάλογος περιορισμός για το Partial Eta Square (μερικό η2). Τέλος, ως μία αναφορά για το πιθανό μέγεθος που μπορεί να έχει ο δείκτης Eta Square (η2) αναφέρουμε πως στις κοινωνικές επιστήμες περιμένουμε να είναι από 0,01 έως 0,09 ενώ η γενικό­τερη ερμηνεία του είναι (Cohen, 1988):
• 0,01 ~ μικρό,
• 0,06 ~ μεσαίο,
• 0,14 ~ μεγάλο.

Ο δείκτης η2 (Eta Square) αρμόζει να υπολογιστεί και να αναφερθεί σε δοκιμασίες ANOVA που είναι
(α) ισορροπημένες (balanced) δηλαδή έχουν το ίδιο πλήθος παρατηρήσεων σε κάθε έναν από τους συνδυασμούς των παραγόντων και
(β) οι παρατηρήσεις είναι ανεξάρτητες.
Ο μερικός δείκτης η2 αρμόζει να αναφερθεί σε δοκιμασίες ANOVA που δεν έχουν ανεξάρτητες παρατηρήσεις, όπως στις επαναλαμβανόμενες μετρήσεις - repeated measures ή ακόμα και στη διαδικασία MANOVA όπου οι εξαρτημένες μεταβλητές πρέπει να έχουν κάποια συσχέτιση.

4 Πολλαπλές μετρήσεις. Αξιολογώντας την αποτελεσματικότητα μίας δίαιτας με το πέρασμα του χρόνου.

Η ανάλυση των πολλαπλών μετρήσεων (repeated measures) χρησιμοποιείται όταν θέλουμε να ελέγξουμε τη στατιστική σημαντικότητα της μεταβολής της μέσης τιμής μεταξύ τριών ή περισσοτέρων μετρήσεων σε μία συνεχή παράμετρο ενός πληθυσμού ενώ υπάρχει επιπλέον η δυνατότητα ελέγχου της διαφοροποίησης των τιμών μεταξύ δύο ή περισσοτέρων ανεξάρτητων ομάδων. Απαραίτητες προϋποθέσεις είναι
(α) η κανονικότητα της κατανομής των τιμών της μεταβλητής σε κάθε μία από τις μετρήσεις (και μεταξύ των ανεξάρτητων ομάδων αν υπάρχουν),
(β) ομοιογένεια μεταξύ των ανεξάρτητων ομάδων (αν υπάρχουν)
(γ) απουσία ιδιαζόντων ή/και ακραίων τιμών (outliers - extreme values)
(δ) η τυπική απόκλιση μεταξύ όλων των συνδυασμών των διαφορών των μετρήσεων και των ομάδων (αν υπάρχουν) να μη διαφέρει στατιστικά (σφαιρικότητα των παρατηρήσεων).

4.1 Παράδειγμα 1ο (γενίκευση του παραδείγματος του paired samples t - test)

Με σκοπό την αξιολόγηση μίας δίαιτας, επιλέχθηκαν 16 ενήλικες εθελοντές οι οποίοι την ακολούθησαν πιστά για 12 μήνες. Το βάρος μετρήθηκε στην αρχή, στους 6 και στους 12 μήνες και τα αποτελέσματα καταχωρήθηκαν στις μεταβλητές weight1, weight2 και weight3. Επιπλέον, για κάθε έναν εθελοντή καταγράφηκε το φύλο του στη μεταβλητή gender με την κωδικοποίηση 0 = Γυναίκα, 1 = Άνδρας.

gender = c(0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1)
gender = factor(gender, levels = c(0, 1), labels = c("Γυναίκα", "Άνδρας"))
weight1 = c(81.2, 76.7, 75.7, 81.2, 71.7, 71.2, 68.5, 89.8, 107.5, 105.7, 99.3, 100.7, 90.3, 105.7, 98.0, 116.6)
weight2 = c(78.0, 73.0, 73.0, 78.5, 69.9, 64.9, 63.5, 87.1, 102.1, 102.5, 97.1, 95.3, 87.5, 102.5, 93.4, 112.9)
weight3 = c(78.4, 72.1, 73.7, 78.5, 69.7, 65.4, 63.3, 85.5, 101.3, 102.2, 95.7, 93.9, 86.5, 102.5, 93.9, 113.9)
dataF = data.frame(gender, weight1, weight2, weight3)
Hmisc::label(dataF$gender) = 'Φύλο'
Hmisc::label(dataF$weight1) = '1η μέτρηση'
Hmisc::label(dataF$weight2) = '2η μέτρηση'
Hmisc::label(dataF$weight3) = '3η μέτρηση'

Επιθυμούμε να βρούμε:
• Αν η δίαιτα είχε αποτέλεσμα στη μείωση του βάρους των εθελοντών.
• Τη διάρκεια (6 ή 12 μήνες) που είναι περισσότερο αποτελεσματική.
• Αν υπάρχει αλληλεπίδραση του φύλου και της αποτελεσματικότητας της δίαιτας.

4.2 Περιγραφή δεδομένων

Για την περιγραφή των δεδομένων μπορούμε να χρησιμοποιήσουμε τις συναρτήσεις explore_vars και my_explore_vars_by_group explore_vars και my_explore_vars_by_group.

Εκτελούμε τις εντολές:

my_explore_vars(dataF, c('weight1', 'weight2', 'weight3'))
Descriptive statistics of 1η μέτρηση, 2η μέτρηση, 3η μέτρηση
Variable N Missing Valid M (SD) CI95
1 1η μέτρηση 16 0 16 90.0 (15.2) 81.9 - 98.1
2 2η μέτρηση 16 0 16 86.3 (15.2) 78.2 - 94.4
3 3η μέτρηση 16 0 16 86.0 (15.1) 78 - 94.1
my_explore_vars_by_group(dataF, c('weight1', 'weight2', 'weight3'), 'gender', Tukeys.test = FALSE)
Descriptive Statistics of 1η μέτρηση among levels of Φύλο
Group N M (SD) 95% C.I.
1 7 75.2 (4.96) 70.6 - 79.8
2 9 102 (8.51) 95 - 108
Total 16 90 (15.2) 81.9 - 98.1
Descriptive Statistics of 2η μέτρηση among levels of Φύλο
Group N M (SD) 95% C.I.
1 7 71.5 (5.86) 66.1 - 77
2 9 97.8 (8.2) 91.5 - 104
Total 16 86.3 (15.2) 78.2 - 94.4
Descriptive Statistics of 3η μέτρηση among levels of Φύλο
Group N M (SD) 95% C.I.
1 7 71.6 (5.91) 66.1 - 77
2 9 97.3 (8.83) 90.5 - 104
Total 16 86 (15.1) 78 - 94.1

Μία απλή επισκόπηση των μέσων τιμών καταδεικνύει πως:
• οι άνδρες έχουν μεγαλύτερο βάρος από τις γυναίκες.
• η δίαιτα είναι ιδιαίτερα αποτελεσματική τους πρώτους 6 μήνες.
• η δίαιτα είναι λιγότερο αποτελεσματική τους μήνες 7 - 12 και αυτό μόνο για τους άνδρες.

Το ερώτημα που πρέπει να απαντήσουμε είναι αν οι παραπάνω διαφοροποιήσεις στο δείγμα είναι τόσο μεγάλες ώστε να τεκμηριώνουν ανάλογες διαφοροποιήσεις στον πληθυσμό, δηλαδή σε όλους όσους θα ακολουθήσουν αυτήν την δίαιτα.

4.3 Υλοποίηση δοκιμασίας

Εκτελούμε τις εντολές

repeated.ANOVA = my_Repeated_ANOVA(dataF, betweencols = 'gender', withincols = c('weight1', 'weight2', 'weight3'))
repeated.ANOVA$anova.table
Repeated ANOVA table: Greenhouse-Geisser correction applied.
Effect df MSE F p.value
1 gender 1, 14 163.04 49.36 *** <.001
2 time 1.57, 22.05 0.97 99.35 *** <.001
3 gender:time 1.57, 22.05 0.97 0.69 .480
1η μέτρηση
Normality statistics
Statistic Description
Skewness & Kurtosis
  Skewness 0.157 Normal distribution has skew = 0. For a unimodal distribution, negative skew commonly indicates that the tail is on the left side of the distribution, and positive skew indicates that the tail is on the right.
  Kurtosis 2.832 Normal distribution has kurtosis = 3. Values over 3 indicates a platykurtic distribution and values less than 3 indicates a leptokurtic distribution.
Normality Tests. Η0: This sample is from a normal distribution vs Η1: Not the Η0.
  Shapiro–Wilk W = 0.965, p = 0.756 This test is more appropriate method for small sample sizes (<50 samples) although it can also be handling on larger sample size.
  Lilliefors D = 0.131, p = 0.657 The Lilliefors test uses the same calculations as the Kolmogorov-Smirnov test, but it is more conservative in the sense that the Lilliefors Test is less likely to show that data is normally distributed.
2η μέτρηση
Normality statistics
Statistic Description
Skewness & Kurtosis
  Skewness 0.211 Normal distribution has skew = 0. For a unimodal distribution, negative skew commonly indicates that the tail is on the left side of the distribution, and positive skew indicates that the tail is on the right.
  Kurtosis 2.555 Normal distribution has kurtosis = 3. Values over 3 indicates a platykurtic distribution and values less than 3 indicates a leptokurtic distribution.
Normality Tests. Η0: This sample is from a normal distribution vs Η1: Not the Η0.
  Shapiro–Wilk W = 0.967, p = 0.787 This test is more appropriate method for small sample sizes (<50 samples) although it can also be handling on larger sample size.
  Lilliefors D = 0.103, p = 0.919 The Lilliefors test uses the same calculations as the Kolmogorov-Smirnov test, but it is more conservative in the sense that the Lilliefors Test is less likely to show that data is normally distributed.
3η μέτρηση
Normality statistics
Statistic Description
Skewness & Kurtosis
  Skewness 0.302 Normal distribution has skew = 0. For a unimodal distribution, negative skew commonly indicates that the tail is on the left side of the distribution, and positive skew indicates that the tail is on the right.
  Kurtosis 2.802 Normal distribution has kurtosis = 3. Values over 3 indicates a platykurtic distribution and values less than 3 indicates a leptokurtic distribution.
Normality Tests. Η0: This sample is from a normal distribution vs Η1: Not the Η0.
  Shapiro–Wilk W = 0.969, p = 0.821 This test is more appropriate method for small sample sizes (<50 samples) although it can also be handling on larger sample size.
  Lilliefors D = 0.114, p = 0.833 The Lilliefors test uses the same calculations as the Kolmogorov-Smirnov test, but it is more conservative in the sense that the Lilliefors Test is less likely to show that data is normally distributed.
Mauchly’s Test for Sphericity. Tests the hypothesis that all pairs of measurements have equal variance in their difference.
Mauchly’s W p
time 0.621 0.045
gender:time 0.621 0.045

Στο παραπάνω output κρύβονται οι απαντήσεις σε όλα τα ερωτήματα που μπορούν να τεθούν για τα δεδομένα που μελετούμε.

Αρχικά παρατηρούμε πως η απαραίτητη προϋπόθεση της κανονικότητας της κατανομής των υπολοίπων για κάθε μία από τις μετρήσεις είναι αποδεκτή.

Στη συνέχεια, παρατηρούμε πως η προϋπόθεση της σφαιρικότητας απορρίπτεται. Πιο συγκεκριμένα, είναι p = 0.045 < 0,05 άρα η υπόθεση Η0: σ13,Α = σ23,Α = σ12,Α = σ13,Γ = σ23,Γ = σ12,Γ, της σφαιρικότητας για τα δεδομένα μας απορρίπτεται (σij,A και σij,Γ η τυπική απόκλιση των διαφορών μεταξύ της i και j μέτρησης για τους άνδρες και τις γυναίκες αντίστοιχα).

Αυτό έχει ως αποτέλεσμα να εφαρμόζεται η διόρθωση Greenhouse-Geisser για τον υπολογισμό των πιθανοτήτων, βάσει των οποίων απορρίπτονται οι επιμέρους υποθέσεις.

Ειδικότερα, στη συνάρτηση my.Repeated.ANOVA αυτό γίνεται σε κάθε περίπτωση καθώς πρόκειται για περισσοτερό συντηρητική επιλογή που οδηγεί σε πιο ασφαλή συμπεράσματα.

Περισσότερες πληροφορίες για την επιλογή αυτή:
https://afex.singmann.science/forums/topic/default-sphericity-correction-method-gg
https://stats.stackexchange.com/q/2824/442
https://stats.stackexchange.com/q/2492/442

4.3.1 Αποτελέσματα

Παρατηρούμε πως υπάρχει σημαντική επίδραση του χρονικού παράγοντα στις τιμές του βάρους (F(1,57; 22,05) = 99,35, p < 0.001), ενώ δεν υπάρχει σημαντική αλληλεπίδραση του χρονικού παράγοντα και του φύλου ως προς τις τιμές του βάρους (F(1,57; 22.05) = 0,69, p = 0,48). Επιπλέον, συμπεραίνουμε πως η επίδραση του φύλου στις τιμές του βάρους είναι στατιστικά σημαντική, μία παρατήρηση αναμενόμενη λόγω της μεγάλης διαφοράς στο μέσο βάρος των ανδρών και των γυναικών F(1, 14) = 49,36, p < 0,001.

Τα παραπάνω αποτελέσματα αναπαριστώνται στο κατάλληλο ευθειόγραμμα, το οποίο παράγεται με τον εξής κώδικα:

data.long = repeated.ANOVA$data.long
# Αλλαγή των labels από weight1, weight2, weight3 σε 1η, 2η, 3η
data.long$time=plyr::revalue(data.long$time, c("weight1"="1η", "weight2"="2η", "weight3"="3η"))
# Δημιουργία του διαγράμματος
par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
interaction.plot(data.long$time, data.long$gender, 
                                    data.long$DependentScore, type="b", col=c("red","blue"), 
                                     legend=F, lty=c(1,2), lwd=2, pch=c(18,24), 
                                     xlab="Μέτρηση",  ylab="Βάρος",
                                     main="Εξέλιξη βάρους ανά φύλο", 
                                     ylim = c(65, 105))

# Προσθήκη υπομνήματος
legend("topright",c("Γυναίκα","Άνδρας"), bty="n", lty=c(1,2), lwd=2, pch=c(18,24), col=c("red","blue"), title="Φύλο",inset=c(-0.3,0))

Συμπερασματικά για το 1ο παράδειγμα
Οι άνδρες έχουν στατιστικά μεγαλύτερο βάρος από τις γυναίκες, υπάρχει μεταβολή του βάρους μεταξύ των τριών μετρήσεων και δεν παρουσιάζεται διαφορετική συμπεριφορά στην αλλαγή αυτή μεταξύ των δύο φύλων. Οι παραπάνω παρατηρήσεις επαληθεύονται και από την παρατήρηση του διαγράμματος των μέσων τιμών όπου οι δύο γραμμές (μπλε: άνδρες, κόκκινη: γυναίκες) εμφανίζουν αρνητική κλίση και είναι περίπου παράλληλες.

4.4 Παράδειγμα 2ο

Υποθέτουμε πως για τον έλεγχο της δίαιτας του 2ου παραδείγματος συλλέξαμε τα δεδομένα των παρακάτω μεταβλητών. Καταχωρούμε τα δεδομένα και εφαρμόζουμε με τον ίδιο ακριβώς τρόπο τη δοκιμασία repeated measures όπως και στο 1ο παράδειγμα.

gender = c(1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0)
gender = factor(gender, levels = c(0, 1), labels = c("Γυναίκα", "Άνδρας"))
weight1 = c(89.8, 107.5, 105.7, 81.2, 99.3, 76.7, 100.7, 75.7, 90.3, 105.7, 81.2, 71.7, 71.2, 98, 116.6, 68.5)
weight2 = c(90.8, 108.5, 106.7, 80.2, 100.3, 75.7, 101.7, 74.8, 91.3, 106.7, 80.2, 70.7, 70.2, 99, 117.6, 67.5)
weight3 = c(91.8, 109.5, 107.7, 79.2, 101.3, 74.7, 102.7, 73.7, 92.3, 107.7, 79.2, 69.7, 69.2, 100, 118.6, 66.5)
dataF = data.frame(gender, weight1, weight2, weight3)
Hmisc::label(dataF$gender) = 'Φύλο'
Hmisc::label(dataF$weight1) = '1η μέτρηση'
Hmisc::label(dataF$weight2) = '2η μέτρηση'
Hmisc::label(dataF$weight3) = '3η μέτρηση'

4.4.1 Περιγραφή δεδομένων

Εκτελούμε τις εντολές:

my_explore_vars(dataF, c('weight1', 'weight2', 'weight3'))
Descriptive statistics of 1η μέτρηση, 2η μέτρηση, 3η μέτρηση
Variable N Missing Valid M (SD) CI95
1 1η μέτρηση 16 0 16 90.0 (15.2) 81.9 - 98.1
2 2η μέτρηση 16 0 16 90.1 (16.1) 81.5 - 98.7
3 3η μέτρηση 16 0 16 90.2 (17.0) 81.2 - 99.3
my_explore_vars_by_group(dataF, c('weight1', 'weight2', 'weight3'), 'gender', Tukeys.test = FALSE)
Descriptive Statistics of 1η μέτρηση among levels of Φύλο
Group N M (SD) 95% C.I.
1 7 75.2 (4.96) 70.6 - 79.8
2 9 102 (8.51) 95 - 108
Total 16 90 (15.2) 81.9 - 98.1
Descriptive Statistics of 2η μέτρηση among levels of Φύλο
Group N M (SD) 95% C.I.
1 7 74.2 (4.96) 69.6 - 78.8
2 9 103 (8.51) 96 - 109
Total 16 90.1 (16.1) 81.5 - 98.7
Descriptive Statistics of 3η μέτρηση among levels of Φύλο
Group N M (SD) 95% C.I.
1 7 73.2 (4.96) 68.6 - 77.8
2 9 104 (8.51) 97 - 110
Total 16 90.2 (17) 81.2 - 99.3

Μία απλή επισκόπηση των μέσων τιμών καταδεικνύει πως:
• οι άνδρες έχουν μεγαλύτερο βάρος από τις γυναίκες.
• υπάρχει μικρή μείωση βάρους για τις γυναίκες.
• υπάρχει αύξηση βάρους για τους άνδρες.

4.4.2 Υλοποίηση δοκιμασίας

Εκτελούμε, τη συνάρτηση

repeated.ANOVA = my_Repeated_ANOVA(dataF, betweencols = c('gender'), withincols =  c('weight1', 'weight2', 'weight3'))
data.long = repeated.ANOVA$data.long
repeated.ANOVA$anova.table
Repeated ANOVA table: Greenhouse-Geisser correction applied.
Effect df MSE F p.value
1 gender 1, 14 155.79 60.88 *** <.001
2 time 2, 28 0.00 1.31 .285
3 gender:time 2, 28 0.00 77176.31 *** <.001
1η μέτρηση
Normality statistics
Statistic Description
Skewness & Kurtosis
  Skewness 0.157 Normal distribution has skew = 0. For a unimodal distribution, negative skew commonly indicates that the tail is on the left side of the distribution, and positive skew indicates that the tail is on the right.
  Kurtosis 2.832 Normal distribution has kurtosis = 3. Values over 3 indicates a platykurtic distribution and values less than 3 indicates a leptokurtic distribution.
Normality Tests. Η0: This sample is from a normal distribution vs Η1: Not the Η0.
  Shapiro–Wilk W = 0.965, p = 0.756 This test is more appropriate method for small sample sizes (<50 samples) although it can also be handling on larger sample size.
  Lilliefors D = 0.131, p = 0.657 The Lilliefors test uses the same calculations as the Kolmogorov-Smirnov test, but it is more conservative in the sense that the Lilliefors Test is less likely to show that data is normally distributed.
2η μέτρηση
Normality statistics
Statistic Description
Skewness & Kurtosis
  Skewness 0.155 Normal distribution has skew = 0. For a unimodal distribution, negative skew commonly indicates that the tail is on the left side of the distribution, and positive skew indicates that the tail is on the right.
  Kurtosis 2.831 Normal distribution has kurtosis = 3. Values over 3 indicates a platykurtic distribution and values less than 3 indicates a leptokurtic distribution.
Normality Tests. Η0: This sample is from a normal distribution vs Η1: Not the Η0.
  Shapiro–Wilk W = 0.965, p = 0.753 This test is more appropriate method for small sample sizes (<50 samples) although it can also be handling on larger sample size.
  Lilliefors D = 0.131, p = 0.651 The Lilliefors test uses the same calculations as the Kolmogorov-Smirnov test, but it is more conservative in the sense that the Lilliefors Test is less likely to show that data is normally distributed.
3η μέτρηση
Normality statistics
Statistic Description
Skewness & Kurtosis
  Skewness 0.157 Normal distribution has skew = 0. For a unimodal distribution, negative skew commonly indicates that the tail is on the left side of the distribution, and positive skew indicates that the tail is on the right.
  Kurtosis 2.832 Normal distribution has kurtosis = 3. Values over 3 indicates a platykurtic distribution and values less than 3 indicates a leptokurtic distribution.
Normality Tests. Η0: This sample is from a normal distribution vs Η1: Not the Η0.
  Shapiro–Wilk W = 0.965, p = 0.756 This test is more appropriate method for small sample sizes (<50 samples) although it can also be handling on larger sample size.
  Lilliefors D = 0.131, p = 0.657 The Lilliefors test uses the same calculations as the Kolmogorov-Smirnov test, but it is more conservative in the sense that the Lilliefors Test is less likely to show that data is normally distributed.
Mauchly’s Test for Sphericity. Tests the hypothesis that all pairs of measurements have equal variance in their difference.
Mauchly’s W p

Η υπόθεση της σφαιρικότητας δεν γίνεται αποδεκτή, συνεπώς έχει εφαρμοστεί η διόρθωση κατά Hyunh-Feldt. Βρίσκουμε πως δεν υπάρχει σημαντική επίδραση του χρονικού παράγοντα στις τιμές του βάρους (F(2, 28) = 1,313, p = 0,285) ωστόσο υπάρχει στατιστική αλληλεπίδραση μεταξύ του χρονικού παράγοντα και του φύλου (F(2, 28) = 77176,3, p < 0,001) και σημαντική επίδραση του φύλου στις τιμές του βάρους (F(1, 14) = 60,871, p < 0,001).

4.4.3 Αποτελέσματα

Συμπερασματικά για το 2ο παράδειγμα
Οι άνδρες έχουν στατιστικά μεγαλύτερο βάρος από τις γυναίκες, δεν υπάρχει μεταβολή του βάρους μεταξύ των τριών μετρήσεων, ωστόσο υπάρχει διαφοροποίηση στη χρονική μεταβολή του βάρους μεταξύ ανδρών και γυναικών. Τα παραπάνω αποτελέσματα αναδεικνύονται περαιτέρω στο διάγραμμα των μέσων τιμών: υπάρχει αύξηση του βάρους για τους άνδρες με το πέρασμα του χρόνου ενώ αντίθετα στις γυναίκες το βάρος μειώνεται.

data.long = repeated.ANOVA$data.long
# Αλλαγή των labels από weight1, weight2, weight3 σε 1η, 2η, 3η
data.long$time=plyr::revalue(data.long$time, c("weight1"="1η", "weight2"="2η", "weight3"="3η"))
# Δημιουργία του διαγράμματος
par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
interaction.plot(data.long$time, data.long$gender, 
                                    data.long$DependentScore, type="b", col=c("red","blue"), 
                                     legend=F, lty=c(1,2), lwd=2, pch=c(18,24), 
                                     xlab="Μέτρηση",  ylab="Βάρος",
                                     main="Εξέλιξη βάρους ανά φύλο", 
                                     ylim = c(65, 105))

# Προσθήκη υπομνήματος
legend("topright",c("Γυναίκα","Άνδρας"), bty="n", lty=c(1,2), lwd=2, pch=c(18,24), col=c("red","blue"), title="Φύλο",inset=c(-0.3,0))

4.5 Παράδειγμα 3ο

Υποθέτουμε πως για τον έλεγχο της δίαιτας του 3ου παραδείγματος συλλέξαμε τα δεδομένα των παρακάτω μεταβλητών. Καταχωρούμε τα δεδομένα και εφαρμόζουμε με τον ίδιο ακριβώς τρόπο τη δοκιμασία repeated measures όπως και στο 1ο παράδειγμα.

gender = c(1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0)
gender = factor(gender, levels = c(0, 1), labels = c("Γυναίκα", "Άνδρας"))
weight1 = c(89.8, 107.5, 105.7, 81.2, 99.3, 76.7, 100.7, 75.7, 90.3, 105.7, 81.2, 71.7, 71.2, 98, 116.6, 68.5)
weight2 = c(88.5, 106.3, 104.1, 76.7, 97.3, 72.3, 99.3, 71.2, 88.7, 104.1, 77.1, 66.9, 67.1, 96.7, 115.2, 64.3)
weight3 = c(83, 101.9, 100, 75.4, 92.9, 70.5, 94.7, 69.4, 84, 99.3, 75.2, 65.2, 66.1, 92.4, 111, 62.9)
dataF = data.frame(gender, weight1, weight2, weight3)
Hmisc::label(dataF$gender) = 'Φύλο'
Hmisc::label(dataF$weight1) = '1η μέτρηση'
Hmisc::label(dataF$weight2) = '2η μέτρηση'
Hmisc::label(dataF$weight3) = '3η μέτρηση'

4.5.1 Περιγραφή δεδομένων

Εκτελούμε τις εντολές:

my_explore_vars(dataF, c('weight1', 'weight2', 'weight3'))
Descriptive statistics of 1η μέτρηση, 2η μέτρηση, 3η μέτρηση
Variable N Missing Valid M (SD) CI95
1 1η μέτρηση 16 0 16 90.0 (15.2) 81.9 - 98.1
2 2η μέτρηση 16 0 16 87.2 (16.5) 78.4 - 96
3 3η μέτρηση 16 0 16 84.0 (15.2) 75.9 - 92.1
my_explore_vars_by_group(dataF, c('weight1', 'weight2', 'weight3'), 'gender', Tukeys.test = FALSE)
Descriptive Statistics of 1η μέτρηση among levels of Φύλο
Group N M (SD) 95% C.I.
1 7 75.2 (4.96) 70.6 - 79.8
2 9 102 (8.51) 95 - 108
Total 16 90 (15.2) 81.9 - 98.1
Descriptive Statistics of 2η μέτρηση among levels of Φύλο
Group N M (SD) 95% C.I.
1 7 70.8 (4.97) 66.2 - 75.4
2 9 100 (8.54) 93.5 - 107
Total 16 87.2 (16.5) 78.4 - 96
Descriptive Statistics of 3η μέτρηση among levels of Φύλο
Group N M (SD) 95% C.I.
1 7 69.2 (4.86) 64.8 - 73.7
2 9 95.5 (8.81) 88.7 - 102
Total 16 84 (15.2) 75.9 - 92.1

Μία απλή επισκόπηση των μέσων τιμών καταδεικνύει πως:
• οι άνδρες έχουν μεγαλύτερο βάρος από τις γυναίκες.
• υπάρχει μείωση βάρους τόσο για τους άνδρες όσο και για τις γυναίκες.
• Οι γυναίκες επωφελούνται περισσότερο το πρώτο εξάμηνο της δίαιτας ενώ οι άνδρες το δεύτερο.

4.5.2 Υλοποίηση δοκιμασίας

Εκτελούμε, τη συνάρτηση

repeated.ANOVA = my_Repeated_ANOVA(dataF, betweencols = c('gender'), withincols =  c('weight1', 'weight2', 'weight3'))
data.long = repeated.ANOVA$data.long
repeated.ANOVA$anova.table
Repeated ANOVA table: Greenhouse-Geisser correction applied.
Effect df MSE F p.value
1 gender 1, 14 158.49 55.39 *** <.001
2 time 1.57, 22.00 0.09 2009.69 *** <.001
3 gender:time 1.57, 22.00 0.09 161.78 *** <.001
1η μέτρηση
Normality statistics
Statistic Description
Skewness & Kurtosis
  Skewness 0.157 Normal distribution has skew = 0. For a unimodal distribution, negative skew commonly indicates that the tail is on the left side of the distribution, and positive skew indicates that the tail is on the right.
  Kurtosis 2.832 Normal distribution has kurtosis = 3. Values over 3 indicates a platykurtic distribution and values less than 3 indicates a leptokurtic distribution.
Normality Tests. Η0: This sample is from a normal distribution vs Η1: Not the Η0.
  Shapiro–Wilk W = 0.965, p = 0.756 This test is more appropriate method for small sample sizes (<50 samples) although it can also be handling on larger sample size.
  Lilliefors D = 0.131, p = 0.657 The Lilliefors test uses the same calculations as the Kolmogorov-Smirnov test, but it is more conservative in the sense that the Lilliefors Test is less likely to show that data is normally distributed.
2η μέτρηση
Normality statistics
Statistic Description
Skewness & Kurtosis
  Skewness 0.185 Normal distribution has skew = 0. For a unimodal distribution, negative skew commonly indicates that the tail is on the left side of the distribution, and positive skew indicates that the tail is on the right.
  Kurtosis 2.827 Normal distribution has kurtosis = 3. Values over 3 indicates a platykurtic distribution and values less than 3 indicates a leptokurtic distribution.
Normality Tests. Η0: This sample is from a normal distribution vs Η1: Not the Η0.
  Shapiro–Wilk W = 0.965, p = 0.756 This test is more appropriate method for small sample sizes (<50 samples) although it can also be handling on larger sample size.
  Lilliefors D = 0.121, p = 0.769 The Lilliefors test uses the same calculations as the Kolmogorov-Smirnov test, but it is more conservative in the sense that the Lilliefors Test is less likely to show that data is normally distributed.
3η μέτρηση
Normality statistics
Statistic Description
Skewness & Kurtosis
  Skewness 0.146 Normal distribution has skew = 0. For a unimodal distribution, negative skew commonly indicates that the tail is on the left side of the distribution, and positive skew indicates that the tail is on the right.
  Kurtosis 2.937 Normal distribution has kurtosis = 3. Values over 3 indicates a platykurtic distribution and values less than 3 indicates a leptokurtic distribution.
Normality Tests. Η0: This sample is from a normal distribution vs Η1: Not the Η0.
  Shapiro–Wilk W = 0.969, p = 0.821 This test is more appropriate method for small sample sizes (<50 samples) although it can also be handling on larger sample size.
  Lilliefors D = 0.121, p = 0.768 The Lilliefors test uses the same calculations as the Kolmogorov-Smirnov test, but it is more conservative in the sense that the Lilliefors Test is less likely to show that data is normally distributed.
Mauchly’s Test for Sphericity. Tests the hypothesis that all pairs of measurements have equal variance in their difference.
Mauchly’s W p
time 0.619 0.044
gender:time 0.619 0.044

Η υπόθεση της σφαιρικότητας δεν γίνεται αποδεκτή, συνεπώς έχει εφαρμοστεί η διόρθωση κατά Hyunh-Feldt. Βρίσκουμε πως υπάρχει σημαντική επίδραση του χρονικού παράγοντα στις τιμές του βάρους (F(1.57, 22) = 2009,7, p < 0,001), στατιστική αλληλεπίδραση μεταξύ του χρονικού παράγοντα και του φύλου (F(1.57, 22) = 161,8, p < 0,001) αλλά και σημαντική επίδραση του φύλου στις τιμές του βάρους (F(1, 14) = 55,39, p < 0,001).

4.5.3 Αποτελέσματα

Συμπερασματικά για το 3ο παράδειγμα
Οι άνδρες έχουν στατιστικά μεγαλύτερο βάρος από τις γυναίκες, υπάρχει μεταβολή του βάρους μεταξύ των τριών μετρήσεων, ενώ υπάρχει περαιτέρω διαφοροποίηση στη χρονική μεταβολή του βάρους μεταξύ ανδρών και γυναικών. Τα παραπάνω αποτελέσματα παρουσιάζονται στο διάγραμμα των μέσων τιμών: υπάρχει μείωση του βάρους τόσο για τους άνδρες όσο και για τις γυναίκες αλλά αυτή η μείωση έχει διαφορετικό χαρακτήρα στα δύο φύλα, γεγονός που αντικατοπτρίζεται στο διάγραμμα ως μη παραλληλία των δύο γραμμών.

# Αλλαγή των labels από weight1, weight2, weight3 σε 1η, 2η, 3η
data.long$time=plyr::revalue(data.long$time, c("weight1"="1η", "weight2"="2η", "weight3"="3η"))
# Δημιουργία του διαγράμματος
par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
interaction.plot(data.long$time, data.long$gender, 
                                    data.long$DependentScore, type="b", col=c("red","blue"), 
                                     legend=F, lty=c(1,2), lwd=2, pch=c(18,24), 
                                     xlab="Μέτρηση",  ylab="Βάρος",
                                     main="Εξέλιξη βάρους ανά φύλο", 
                                     ylim = c(65, 105))

# Προσθήκη υπομνήματος
legend("topright",c("Γυναίκα","Άνδρας"), bty="n", lty=c(1,2), lwd=2, pch=c(18,24), col=c("red","blue"), title="Φύλο",inset=c(-0.3,0))

4.6 Παράδειγμα 4ο

Υποθέτουμε πως για τον έλεγχο της δίαιτας του 4ου παραδείγματος συλλέξαμε τα δεδομένα των παρακάτω μεταβλητών. Καταχωρούμε τα δεδομένα και εφαρμόζουμε με τον ίδιο ακριβώς τρόπο τη δοκιμασία repeated measures όπως και στο 1ο παράδειγμα.

gender = c(1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0)
gender = factor(gender, levels = c(0, 1), labels = c("Γυναίκα", "Άνδρας"))
weight1 = c(89.8, 107.5, 105.7, 81.2, 99.3, 76.7, 100.7, 75.7, 90.3, 105.7, 81.2, 71.7, 71.2, 98, 116.6, 68.5)
weight2 = c(88.4, 107.2, 106.1, 80.4, 98, 77.5, 100.1, 75.5, 88.5, 104.1, 79.8, 69.8, 71.7, 96.6, 117.4, 70.2)
weight3 = c(89.6, 108.8, 105.4, 79.6, 96.5, 78.1, 100.1, 74.7, 89.7, 103, 81.5, 69.7, 71.1, 94.8, 118.2, 70.1)
dataF = data.frame(gender, weight1, weight2, weight3)
Hmisc::label(dataF$gender) = 'Φύλο'
Hmisc::label(dataF$weight1) = '1η μέτρηση'
Hmisc::label(dataF$weight2) = '2η μέτρηση'
Hmisc::label(dataF$weight3) = '3η μέτρηση'

4.6.1 Περιγραφή δεδομένων

Εκτελούμε τις εντολές:

my_explore_vars(dataF, c('weight1', 'weight2', 'weight3'))
Descriptive statistics of 1η μέτρηση, 2η μέτρηση, 3η μέτρηση
Variable N Missing Valid M (SD) CI95
1 1η μέτρηση 16 0 16 90.0 (15.2) 81.9 - 98.1
2 2η μέτρηση 16 0 16 89.5 (15.1) 81.4 - 97.5
3 3η μέτρηση 16 0 16 89.4 (15.1) 81.4 - 97.5
my_explore_vars_by_group(dataF, c('weight1', 'weight2', 'weight3'), 'gender', Tukeys.test = FALSE)
Descriptive Statistics of 1η μέτρηση among levels of Φύλο
Group N M (SD) 95% C.I.
1 7 75.2 (4.96) 70.6 - 79.8
2 9 102 (8.51) 95 - 108
Total 16 90 (15.2) 81.9 - 98.1
Descriptive Statistics of 2η μέτρηση among levels of Φύλο
Group N M (SD) 95% C.I.
1 7 75 (4.47) 70.9 - 79.1
2 9 101 (9.26) 93.6 - 108
Total 16 89.5 (15.1) 81.4 - 97.5
Descriptive Statistics of 3η μέτρηση among levels of Φύλο
Group N M (SD) 95% C.I.
1 7 75 (4.84) 70.5 - 79.4
2 9 101 (9.33) 93.5 - 108
Total 16 89.4 (15.1) 81.4 - 97.5

Μία απλή επισκόπηση των μέσων τιμών καταδεικνύει πως:
• οι άνδρες έχουν μεγαλύτερο βάρος από τις γυναίκες.
• υπάρχει πολύ μικρή μείωση βάρους τόσο για τους άνδρες όσο και για τις γυναίκες.

4.6.2 Υλοποίηση δοκιμασίας

Εκτελούμε, τη συνάρτηση

repeated.ANOVA = my_Repeated_ANOVA(dataF, betweencols = c('gender'), withincols =  c('weight1', 'weight2', 'weight3'))
data.long = repeated.ANOVA$data.long
repeated.ANOVA$anova.table
Repeated ANOVA table: Greenhouse-Geisser correction applied.
Effect df MSE F p.value
1 gender 1, 14 167.58 47.37 *** <.001
2 time 1.73, 24.16 0.98 1.59 .226
3 gender:time 1.73, 24.16 0.98 0.60 .531
1η μέτρηση
Normality statistics
Statistic Description
Skewness & Kurtosis
  Skewness 0.157 Normal distribution has skew = 0. For a unimodal distribution, negative skew commonly indicates that the tail is on the left side of the distribution, and positive skew indicates that the tail is on the right.
  Kurtosis 2.832 Normal distribution has kurtosis = 3. Values over 3 indicates a platykurtic distribution and values less than 3 indicates a leptokurtic distribution.
Normality Tests. Η0: This sample is from a normal distribution vs Η1: Not the Η0.
  Shapiro–Wilk W = 0.965, p = 0.756 This test is more appropriate method for small sample sizes (<50 samples) although it can also be handling on larger sample size.
  Lilliefors D = 0.131, p = 0.657 The Lilliefors test uses the same calculations as the Kolmogorov-Smirnov test, but it is more conservative in the sense that the Lilliefors Test is less likely to show that data is normally distributed.
2η μέτρηση
Normality statistics
Statistic Description
Skewness & Kurtosis
  Skewness 0.235 Normal distribution has skew = 0. For a unimodal distribution, negative skew commonly indicates that the tail is on the left side of the distribution, and positive skew indicates that the tail is on the right.
  Kurtosis 3.168 Normal distribution has kurtosis = 3. Values over 3 indicates a platykurtic distribution and values less than 3 indicates a leptokurtic distribution.
Normality Tests. Η0: This sample is from a normal distribution vs Η1: Not the Η0.
  Shapiro–Wilk W = 0.956, p = 0.587 This test is more appropriate method for small sample sizes (<50 samples) although it can also be handling on larger sample size.
  Lilliefors D = 0.126, p = 0.717 The Lilliefors test uses the same calculations as the Kolmogorov-Smirnov test, but it is more conservative in the sense that the Lilliefors Test is less likely to show that data is normally distributed.
3η μέτρηση
Normality statistics
Statistic Description
Skewness & Kurtosis
  Skewness 0.523 Normal distribution has skew = 0. For a unimodal distribution, negative skew commonly indicates that the tail is on the left side of the distribution, and positive skew indicates that the tail is on the right.
  Kurtosis 3.066 Normal distribution has kurtosis = 3. Values over 3 indicates a platykurtic distribution and values less than 3 indicates a leptokurtic distribution.
Normality Tests. Η0: This sample is from a normal distribution vs Η1: Not the Η0.
  Shapiro–Wilk W = 0.957, p = 0.606 This test is more appropriate method for small sample sizes (<50 samples) although it can also be handling on larger sample size.
  Lilliefors D = 0.135, p = 0.603 The Lilliefors test uses the same calculations as the Kolmogorov-Smirnov test, but it is more conservative in the sense that the Lilliefors Test is less likely to show that data is normally distributed.
Mauchly’s Test for Sphericity. Tests the hypothesis that all pairs of measurements have equal variance in their difference.
Mauchly’s W p
time 0.721 0.12
gender:time 0.721 0.12

Στην περίπτωση αυτή, υποστηρίζεται η σφαιρικότητα των παρατηρήσεων (p = 0,102).Βρίσκουμε πως δεν υπάρχει σημαντική επίδραση του χρονικού παράγοντα στις τιμές του βάρους (F(2, 28) = 1,586, p = 0.223), αλλά ούτε σημαντική αλληλεπίδραση μεταξύ του χρονικού παράγοντα και του φύλου (F(2, 28) = 0,605, p = 0,553). Υπάρχει ωστόσο σημαντική επίδραση του φύλου στις τιμές του βάρους (F(1, 14) = 47,372, p < 0,001).

4.6.3 Αποτελέσματα

Συμπερασματικά για το 4ο παράδειγμα
Οι άνδρες έχουν στατιστικά μεγαλύτερο βάρος από τις γυναίκες, ωστόσο δεν υπάρχει μεταβολή του βάρους μεταξύ των τριών μετρήσεων, ενώ επιπλέον δεν υπάρχει διαφοροποίηση στη χρονική μεταβολή του βάρους μεταξύ ανδρών και γυναικών.

Τα παραπάνω αποτελέσματα παρουσιάζονται στο διάγραμμα των μέσων τιμών: υπάρχει σχετική σταθερότητα του βάρους μεταξύ των τριών μετρήσεων ενώ δεν αναδεικνύεται κάποια διαφοροποίηση μεταξύ ανδρών και γυναικών.

# Αλλαγή των labels από weight1, weight2, weight3 σε 1η, 2η, 3η
data.long$time=plyr::revalue(data.long$time, c("weight1"="1η", "weight2"="2η", "weight3"="3η"))
# Δημιουργία του διαγράμματος
par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
interaction.plot(data.long$time, data.long$gender, 
                                    data.long$DependentScore, type="b", col=c("red","blue"), 
                                     legend=F, lty=c(1,2), lwd=2, pch=c(18,24), 
                                     xlab="Μέτρηση",  ylab="Βάρος",
                                     main="Εξέλιξη βάρους ανά φύλο", 
                                     ylim = c(65, 105))

# Προσθήκη υπομνήματος
legend("topright",c("Γυναίκα","Άνδρας"), bty="n", lty=c(1,2), lwd=2, pch=c(18,24), col=c("red","blue"), title="Φύλο",inset=c(-0.3,0))

4.7 Ταυτότητα της ανάλυσης πολλαπλών μετρήσεων (repeated measures ANOVA)

Στατιστικές υποθέσεις

Μηδενικές υποθέσεις Η0
Η0,1: Δεν υπάρχει επίδραση του χρονικού παράγοντα στη μέση τιμή του βάρους.
Η0,2: Δεν υπάρχει αλληλεπίδραση μεταξύ του χρονικού παράγοντα και του φύλο ως προς την τιμή του βάρους.
Η0,3: Δεν υπάρχει στατιστική διαφορά στο βάρος μεταξύ των διαφορετικών επιπέδων του φύλου.

Εναλλακτικές υποθέσεις Η1
Η1,1: όχι η Η0,1.
Η1,2: όχι η Η0,2.
Η1,3: όχι η Η0,3.

Ερευνητικός στόχος
Ανίχνευση της επίδρασης του χρονικού παράγοντα και του φύλου στη μέση τιμή του βάρους.

Προϋποθέσεις
(α) σφαιρικότητα των παρατηρήσεων (ελέγχουμε όταν υπάρχουν περισσότερες από 2 επαναλαμβανόμενες μετρήσεις).
(β) ομοιογένεια ως προς τις ανεξάρτητες ομάδες του φύλου.
(γ) κανονικότητα των υπολοίπων (προκύπτει από την με την κανονικότητα των επαναλαμβανόμενων μετρήσεων σε όλες τις κατηγορίες της ανεξάρτητης μεταβλητής)
(δ) απουσία ακραίων τιμών (outliers - extreme values)

Σημείωση: ο παραπάνω πίνακας αναφέρεται στο σενάριο του παραδείγματος για α = 0,05.

4.8 Ενδεικτική καταγραφή δοκιμασίας ανάλυσης πολλαπλών μετρήσεων (repeated measures ANOVA)

Στοιχεία που πρέπει να αναφερθούν
• τα στοιχεία των δύο μετρήσεων (πλήθος παρατηρήσεων, η μέση τιμή και η τυπική απόκλιση για κάθε μία από αυτές),
• η τιμή του στατιστικού F (στήλη t στο output) ως απόλυτη τιμή,
• οι 2 βαθμοί ελευθερίας (στήλη df, γραμμή Between Groups και γραμμή Within Groups με αυτήν τη σειρά),
• η στατιστική σημαντικότητα p της διαφοροποίησης από τη στατιστική υπόθεση Η0 (από τη στήλη Sig.).

Υπόδειγμα 1
Βρέθηκε πως υπάρχει σημαντική επίδραση του χρονικού παράγοντα στις τιμές του βάρους (F(1,690; 23,665) = 99,352, p < 0.001), ενώ δεν υπάρχει σημαντική αλληλεπίδραση του χρονικού παράγοντα και του φύλου ως προς τις τιμές του βάρους (F(1,690; 23,665) = 0,687, p = 0,489).

Υπόδειγμα 2
Βρέθηκε πως δεν υπάρχει σημαντική επίδραση του χρονικού παράγοντα στις τιμές του βάρους (F(2, 28) = 1,109, p = 0,344) ωστόσο υπάρχει στατιστική αλληλεπίδραση μεταξύ του χρονικού παράγοντα και του φύλου (F(2, 28) = 288312,7, p < 0,001) και σημαντική επίδραση του φύλου στις τιμές του βάρους (F(1, 14) = 60,949, p < 0,001).

Υπόδειγμα 3
Βρέθηκε πως υπάρχει σημαντική επίδραση του χρονικού παράγοντα στις τιμές του βάρους ((F(1.522, 21.304) = 2378,864, p < 0,001), στατιστική αλληλεπίδραση μεταξύ του χρονικού παράγοντα και του φύλου (F(1.522, 21.304) = 184,630, p < 0,001) αλλά και σημαντική επίδραση του φύλου στις τιμές του βάρους (F(1, 14) = 55,046, p < 0,001).

Υπόδειγμα 4
Βρέθηκε πως δεν υπάρχει σημαντική επίδραση του χρονικού παράγοντα στις τιμές του βάρους ((F(2, 28) = 1,569, p = 0.226), αλλά ούτε σημαντική αλληλεπίδραση μεταξύ του χρονικού παράγοντα και του φύλου (F(2, 28) = 0,584, p = 0,564). Υπάρχει ωστόσο σημαντική επίδραση του φύλου στις τιμές του βάρους (F(1, 14) = 47,372, p < 0,001).

5 Παρατηρήσεις

5.1 Eta Square (η2) και Partial Eta Square (μερικό η2)

Οι δείκτες Eta Square και Partial Eta Square είναι δυνατό να υπολογιστούν και να αναφερθούν, διαδικασία που παραλείπεται για λόγους απλότητας της παρουσίασης. Ο ενδιαφερόμενος αναγνώστης παραπέμπεται στην δεύτερη παρατήρηση της παραγράφου της δοκιμασίας 2-way ANOVA.

5.2 Σφαιρικότητα

Αν υπάρχει μεγάλη απόκλιση από την προϋπόθεση της σφαιρικότητας τότε ίσως να πρέπει ο ερευνητής εφαρμόσει MANOVA που δεν απαιτεί αυτή την προϋπόθεση χωρίς όμως να θεωρεί τις παρατηρήσεις ως δείγματα από τα ίδια υποκείμενα. Στην περίπτωση αυτή, λόγω της διαφορετικής διαχείρισης των δεδομένων από τη MANOVA απαιτείται μεγαλύτερο δείγμα από τη REPEATED MEASURES για το ίδιο επίπεδο αξιοπιστίας του τελικού αποτελέσματος. (Field 1998, Keselman κ.α. 2001).

5.3 Post-hoc συγκρίσεις

Δεν αναφερθήκαμε στην εκ των υστέρων σύγκριση των μετρήσεων για την ανίχνευση στατιστικώς σημαντικών διαφοροποιήσεων μεταξύ τους. Ο ερευνητής για να βρει στατιστικά όμοιες μετρήσεις μπορεί να εφαρμόσει πολλαπλά t-test για ζευγαρωτές παρατηρήσεις μεταξύ όλων των συνδυασμών των μετρήσεων, διορθώνοντας ωστόσο το όριο στατιστικής σημαντικότητας της διαφοροποίησης από p = 0.05 σε p = 0.05/3 = 0,016. Η εντολή που απαιτείται είναι η εξής:

my_t_test_paired_samples(dataF, c('weight1', 'weight2', 'weight3'))
1η μέτρηση (90 ± 15.2) vs 2η μέτρηση (89.5 ± 15.1), N = 16 pairs. H0: μ1η μέτρηση = μ2η μέτρηση vs H1: μ1η μέτρηση ≠ μ2η μέτρηση. H0 is not rejected (t(15) = 1.93, p = 0.073).

1η μέτρηση (90 ± 15.2) vs 3η μέτρηση (89.4 ± 15.1), N = 16 pairs. H0: μ1η μέτρηση = μ3η μέτρηση vs H1: μ1η μέτρηση ≠ μ3η μέτρηση. H0 is not rejected (t(15) = 1.4, p = 0.181).

2η μέτρηση (89.5 ± 15.1) vs 3η μέτρηση (89.4 ± 15.1), N = 16 pairs. H0: μ2η μέτρηση = μ3η μέτρηση vs H1: μ2η μέτρηση ≠ μ3η μέτρηση. H0 is not rejected (t(15) = 0.0909, p = 0.929).

Από το output των παραπάνω εντολών συνάγεται πως δεν υπάρχει σημαντική διαφοροποίηση μεταξύ των τριών μετρήσεων, επιβεβαιώνοντας το αποτέλεσμα της Repeated ANOVA.

6 Ανάλυση Διακύμανσης με παράλληλο έλεγχο ως προς συσχετισμένες μεταβλητές (ANCOVA).

Χρησιμοποιούμε ανάλυση διακύμανσης με παράλληλο έλεγχο ως προς συσχετισμένες μεταβλητές (ANCOVA) όταν θέλουμε να ελέγξουμε τη σημαντικότητα της επίδρασης ενός ή περισσοτέρων παραγόντων στις τιμές μίας συνεχής μεταβλητής ελέγχοντας παράλληλα ως προς τις τιμές μίας ή περισσοτέρων συνεχών μεταβλητών οι οποίες συσχετίζονται με αυτήν και ενδέχεται να αποτελούν μέρη της επιρροής των παραγόντων στις τιμές της.
Προϋποθέσεις της πολλαπλής ανάλυσης διακύμανσης είναι οι εξής:
• Τα υπόλοιπα (ή σφάλματα - residuals) είναι ανεξάρτητα από τις μεταβλητές που συμμετέχουν στο μοντέλο.
• Ομοιογένεια των υπολοίπων σε όλες τις κατηγορίες που ορίζονται από τους παράγοντες.
• Τα υπόλοιπα ακολουθούν κανονική κατανομή με μέση τιμή 0.

6.1 Παράδειγμα

Ένας ιδιοκτήτης φροντιστηρίου ξένων γλωσσών ενδιαφέρεται να εξετάσει αν υπάρχει διαφοροποίηση μεταξύ ανδρών και γυναικών ως προς την επίδοση στις εξετάσεις Proficiency. Ωστόσο, καθώς έχει τη πεποίθηση πως η ηλικία του υποψηφίου επηρεάζει την επίδοση θέλει να ανιχνεύσει τη διαφοροποίηση της επίδοσης ως προς το φύλο παίρνοντας υπόψη του και την επιρροή της ηλικίας σε αυτήν. Για το σκοπό αυτό συνέλλεξε στοιχεία από 12 σπουδαστές του φροντιστηρίου και τα καταχώρησε στις παρακάτω μεταβλητές.

exam = c(65, 65, 60, 70, 55, 80, 40, 90, 50, 100, 30, 95)
gender = c(0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0)
gender = factor(gender, levels = c(0, 1), labels = c("Γυναίκα", "Άνδρας")) 
age = c(29, 28, 22, 19, 18, 28, 16, 25, 17, 35, 15, 32)
ancova.df = data.frame(exam, gender, age)

Hmisc::label(ancova.df$exam) = 'Εξετάσεις'
Hmisc::label(ancova.df$age) = 'Ηλικία'
Hmisc::label(ancova.df$gender) = 'Φύλο'

Απαιτήθηκαν τρεις μεταβλητές, η μεταβλητή exam (Επίδοση) στην οποία καταχωρήθηκε η τελική επίδοση στις εξετάσεις, η μεταβλητή gender (Φύλο) και τιμές 0: Γυναίκα και 1: Άνδρας και η μεταβλητή age όπου καταχωρείται η ηλικία του υποψηφίου. Αρχικά, με την εντολή:

my_explore_vars_by_group(ancova.df, 'exam', 'gender')
Descriptive Statistics of Εξετάσεις among levels of Φύλο
Group N M (SD) 95% C.I.
1 7 80.7 (14.6) 67.3 - 94.2
2 5 47 (12) 32 - 62
Total 12 66.7 (21.7) 52.9 - 80.4
Tukeys HSD Post Hoc Test for variable: Εξετάσεις
The Tukey HSD is not a Post Hoc Test per se. It may provide valuable results as an independent test as well.
Difference Lower Upper p
2-1 -33.714 -51.466 -15.962 0.002

βρίσκουμε πως οι γυναίκες έχουν μέση επίδοση 80,7 μονάδες (SD = 14,6) έναντι 47 μονάδων των ανδρών συνυποψήφιών τους (SD = 12). Επιπλέον, εκτελώντας την εντολή

my_t_test_independent_samples(ancova.df, 'exam', 'gender')
Εξετάσεις over Φύλο H0: μΓυναίκα = μΆνδρας vs H1: μΓυναίκα ≠ μΆνδρας. H0 is rejected.
Group Γυναίκα (N = 7): 80.7 ± 14.6 vs Group Άνδρας (N = 5): 47 ± 12, t(10) = 4.232, p = 0.002.
Equality of variances: F(6, 4) = 1.461, p = 0.744.

εφαρμόζουμε τη δοκιμασία t - test και βρίσκουμε πως η ισότητα των δύο μέσων τιμών απορρίπτεται (t(10) = 4,232, p = 0,002, δίπλευρος έλεγχος) άρα φαίνεται πως πράγματι υπάρχει διαφορά στην επίδοση μεταξύ ανδρών και γυναικών. Ωστόσο, όπως σωστά υποψιάζεται ο ιδιοκτήτης του φροντιστηρίου υπάρχει και συσχέτιση της επίδοσης με την ηλικία του υποψηφίου. Εκτελώντας την εντολή

my_cor_table(ancova.df, c('exam', 'age'))
Pearson Correlation Coefficients (*:p<.05, **:p<.01, ***:p<.001)
Εξετάσεις Ηλικία
Εξετάσεις - .859***
Ηλικία .859*** -

βλέπουμε πως υπάρχει σαφής θετική γραμμική σχέση μεταξύ της ηλικίας και της επίδοσης (r(12) = 0,859, p < 0,001).

Άρα, ανακύπτει με φυσικό τρόπο το ερώτημα:

Είναι η διαφοροποίηση στην επίδοση μεταξύ ανδρών και γυναικών ανεξάρτητη από την ηλικία του υποψηφίου ή μήπως προκύπτει συμπτωματικά ενώ η ουσιαστική επιρροή στην επίδοση είναι αυτή της ηλικίας;

Η απάντηση σε αυτό το ερώτημα μπορεί να προκύψει με εφαρμογή της διαδικασίας ANCOVA. Εκτελούμε τις εντολές

options(contrasts = c("contr.treatment", "contr.poly"))

my.ancova = my_ANOVA(ancova.df, exam ~ age + gender, show.eta.square = TRUE)
my.ancova$htmlTable
ANOVA Results for Εξετάσεις
SS df F p η2 η2p
(Intercept) 159.721 1 1.244 0.294
age 695.902 1 5.42 0.045 0.135 0.376
gender 198.218 1 1.544 0.245 0.038 0.146
Residuals 1155.527 9
LSD Test for variable: Φύλο
Fisher’s LSD is a series of pairwise t-tests, with each test using the mean squared error from the significant ANOVA as its pooled variance estimate (and naturally taking the associated degrees of freedom). It is a valid test, in case where the ANOVA is significant.
N Εξετάσεις Groups
Άνδρας 5 47 b
Γυναίκα 7 80.714 a
Tukeys HSD Post Hoc Test for variable: Φύλο
The Tukey HSD is not a Post Hoc Test per se. It may provide valuable results even if ANOVA is not significant.
Difference Lower Upper p
Άνδρας-Γυναίκα -33.714 -51.466 -15.962 0.002
Standard Deviation of Εξετάσεις among levels of all levels of gender.
Level Γυναίκα Άνδρας
Group Size 7 5
SD 14.6 12.0
Levene test of variance equality (homogeneity) of Εξετάσεις over all levels of gender.
F(1, 10) = 0.622, p = 0.448. Homogeneity hypothesis is confirmed.
Bartlett’s test of variance equality (homogeneity) of Εξετάσεις over all levels of gender.
If you have strong evidence that your data do in fact come from a normal, or nearly normal, distribution, then Bartlett’s test has better performance.
c2(1) = 0.152, p = 0.697. Homogeneity hypothesis is confirmed.
Pearson Correlation Coefficients (*:p<.05, **:p<.01, ***:p<.001)
Εξετάσεις Ηλικία
Εξετάσεις - .859***
Ηλικία .859*** -
Normality test for model residuals.
If the main goal of an ANOVA is to see whether or not certain effects are significant, then the assumption of normality of the residuals is only required for small samples, thanks to the central limit theorem. With sample sizes of a few hundred participants even extreme violations of the normality assumptions are unproblematic. So mild violations of this assumptions are usually no problem with sample sizes exceeding 30.
Statistic Description
Skewness & Kurtosis
  Skewness -0.564 Normal distribution has skew = 0. For a unimodal distribution, negative skew commonly indicates that the tail is on the left side of the distribution, and positive skew indicates that the tail is on the right.
  Kurtosis 2.149 Normal distribution has kurtosis = 3. Values over 3 indicates a platykurtic distribution and values less than 3 indicates a leptokurtic distribution.
Normality Tests. Η0: This sample is from a normal distribution vs Η1: Not the Η0.
  Shapiro–Wilk W = 0.898, p = 0.147 This test is more appropriate method for small sample sizes (<50 samples) although it can also be handling on larger sample size.
  Lilliefors D = 0.241, p = 0.053 The Lilliefors test uses the same calculations as the Kolmogorov-Smirnov test, but it is more conservative in the sense that the Lilliefors Test is less likely to show that data is normally distributed.

Έλεγχος προϋποθέσεων Πριν την ανάγνωση του output της διαδικασίας Anova, επαληθεύουμε από το output πως τηρείται η κανονικότητα της κατανομής των υπολοίπων του γραμμικού μοντέλου. Ειδικότερα:
• ο συντελεστής κυρτότητας είναι 2,149
• ο συντελεστής ασυμμετρίας είναι -0,564
και πως η υπόθεση πως αυτά προέρχονται από κανονική κατανομή δεν απορρίπτεται (Shapiro - Wilk W = 0,896, p = 0,147).

Περαιτέρω, ως προς την προϋπόθεση της ομοιογένειας, από τον αντίστοιχο output παρατηρούμε πως η υπόθεση πως το δείγμα μας προέρχεται από ομοιογενή πληθυσμό δεν απορρίπτεται F(1, 10) = 0,622, p = 0,448.

6.2 Αποτελέσματα

Συνεχίζουμε προχωρώντας στην ερμηνεία του output. Με τη διαδικασία που ορίσαμε ελέγχουμε τις στατιστικές υποθέσεις:

Η0,1: Η ηλικία του υποψηφίου δεν επιδρά σημαντικά στην επίδοση.
Η0,2: Ο παράγοντας φύλο δεν επιδρά σημαντικά στην επίδοση.

Παρατηρούμε πως:
• Ως προς τον παράγοντα φύλο (gender): F(1, 9) = 1,544, p = 0,245.
• Ως προς τη επιρροή της ηλικίας (age): F(1, 9) = 5,420, p = 0,045.

Οι παραπάνω στατιστικοί υπολογισμοί ερμηνεύονται ως εξής:
• Η στατιστική υπόθεση Η0,1 απορρίπτεται, F(1, 9) = 5,420, p = 0,045, η2μερικό = 0,376, δηλαδή η ηλικία ενός υποψηφίου επηρεάζει με στατιστικά σημαντικό τρόπο την επίδοση του στις εξετάσεις Proficiency.
• Ως προς την υπόθεση Η0,2 δεν καταγράφεται στατιστικά σημαντική επίδραση του φύλου στις τιμές των πωλήσεων.

Είναι φανερό πως το τελικό συμπέρασμα είναι διαφορετικό από αυτό που θα καταλήγαμε αν βασιζόμασταν αποκλειστικά στο t - test μεταξύ ανδρών και γυναικών για την εξαγωγή των συμπερασμάτων.

Στην παράγραφο αυτή, φανερώνεται ο βασικός προβληματισμός κάθε ερευνητή αλλά και ο περιορισμός της στατιστικής επιστήμης ως μέθοδος δημιουργίας νέας γνώσης. Παρατηρήσαμε πως το τελικό αποτέλεσμα είναι διαφορετικό από αυτό που θα καταλήγαμε αν δεν είχαμε συμπεριλάβει και την ηλικία ως παράγοντα μεταβολής της επίδοσης των μαθητών. Το ηθικό δίδαγμα είναι για την ορθή υλοποίηση μίας έρευνας που βασίζεται σε συλλογή δεδομένων είναι απαραίτητος ο καλός σχεδιασμός πριν τη συλλογή των δεδομένων, η επιλογή των κατάλληλων μεταβλητών και της κατάλληλης δοκιμασίας. Δυστυχώς, ο μόνος τρόπος για να ανιχνευθεί ο λάθος σχεδιασμός είναι η επανάληψη της διαδικασίας από ανεξάρτητο ερευνητή, με την συμπερίληψη επιπλέον μεταβλητών.

7 Παρατηρήσεις

7.1 Εφαρμογή γραμμικής παλινδρόμησης στη θέση της ANCOVA

Για τα δεδομένα του δείγματος είναι δυνατό να χρησιμοποιηθεί η απλή γραμμική παλινδρόμηση για την εκμαίευση του ίδιου αποτελέσματος. Πράγματι, αυτό είναι δυνατό, γιατί η επίδοση και η ηλικία είναι συνεχείς μεταβλητές ενώ το φύλο παίρνει μόνο δύο τιμές (0: Γυναίκα, 1: Άνδρας) άρα η αύξηση από 0 στο 1 ερμηνεύεται άμεσα ως η διαφοροποίηση μεταξύ ανδρών και γυναικών.

Εκτελούμε τις εντολές

my.model = my_lm(ancova.df, exam ~ age + gender)
my.model$htmlreport
Linear regression results for Εξετάσεις
Model: exam ~ const + age + gender
  95% C.I.
B SE t p   Lower Upper
Constant 26.438 23.704 1.115 0.294   -27.183 80.059
Ηλικία 1.938 0.833 2.328 0.045   0.055 3.822
Φύλο -13.554 10.909 -1.243 0.245   -38.232 11.123
Coefficient of determination
R2 = 0.776
Homoscedasticity Breusch–Pagan test
x2(1) = 1.991, p = 0.37. Homoscedasticity assumption is not rejected. Model seems to be valid.
Normality tests of models’ residuals
Statistic Description
Skewness & Kurtosis
  Skewness -0.564 Normal distribution has skew = 0. For a unimodal distribution, negative skew commonly indicates that the tail is on the left side of the distribution, and positive skew indicates that the tail is on the right.
  Kurtosis 2.149 Normal distribution has kurtosis = 3. Values over 3 indicates a platykurtic distribution and values less than 3 indicates a leptokurtic distribution.
Normality Tests. Η0: This sample is from a normal distribution vs Η1: Not the Η0.
  Shapiro–Wilk W = 0.898, p = 0.147 This test is more appropriate method for small sample sizes (<50 samples) although it can also be handling on larger sample size.
  Lilliefors D = 0.241, p = 0.053 The Lilliefors test uses the same calculations as the Kolmogorov-Smirnov test, but it is more conservative in the sense that the Lilliefors Test is less likely to show that data is normally distributed.

Από το παραπάνω output, προτείνεται το γραμμικό μοντέλο:

Επίδοση = 26,438 - 13,553* Φύλο + 1,938 * Ηλικία.

Από την ανάγνωση του output, παρατηρούμε πως η παραπάνω εξίσωση μπορεί να χρησιμοποιηθεί για την πρόβλεψη της επίδοσης με αξιοπρεπή επιτυχία (συντελεστής προσδιορισμού R2 = 0,776 = 77,6%). Από τις δύο ανεξάρτητες μεταβλητές το Φύλο προκύπτει μη σημαντικό μέρος για το γραμμικό μοντέλο (t(10) = 1,243, p = 0,245) ενώ η Ηλικία είναι σημαντικός παράγοντας μεταβολής της επίδοσης (t(10) = 2,328, p = 0,045).
Σημειώνεται πως αν στη θέση της μεταβλητής Φύλο είχαμε κάποια άλλη ποιοτική μεταβλητή με περισσότερες από 2 τιμές τότε η εφαρμογή της γραμμικής παλινδρόμησης αντί της ANCOVA δεν θα είχε νόημα, τουλάχιστον όχι δίχως κάποια αντικατάσταση με ψευδομεταβλητές (dummy variables).

8 Πολλαπλή Ανάλυση Διακύμανσης (MANOVA).

Χρησιμοποιούμε πολλαπλή ανάλυση διακύμανσης (MANOVA) όταν θέλουμε να ελέγξουμε τη σημαντικότητα της επίδρασης μίας ομάδας ανεξάρτητων παραγόντων στις τιμές δύο ή περισσοτέρων συνεχών μεταβλητών οι οποίες έχουν κοινή νοηματική βάση.

Προϋποθέσεις της πολλαπλής ανάλυσης διακύμανσης:

(α) Ιδανικό είναι να υπάρχει το ίδιο πλήθος παρατηρήσεων σε κάθε έναν συνδυασμό των ανεξάρτητων παραγόντων (balanced design), το οποίο πλήθος πρέπει να είναι μεγαλύτερο από το πλήθος των ανεξάρτητων παραγόντων. Επιπλέον, ο λόγος των συχνοτήτων δύο συνδυασμών δεν πρέπει να είναι μεγαλύτερος από 1:1,5. Αν σε κάθε έναν συνδυασμό των ανεξάρτητων παραγόντων υπάρχει πλήθος μεγαλύτερο από 30 τότε η δοκιμασία είναι έγκυρη ακόμα και όταν υπάρχει απόκλιση από την κανονικότητα ή την ομοιογένεια.

(β) Απαιτείται ανεξαρτησία μεταξύ των παρατηρήσεων (άρα η MANOVA δεν μπορεί να εφαρμοστεί για επαναλαμβανόμενες μετρήσεις - repeated measures).

(γ) Οι εξαρτημένες μεταβλητές πρέπει να έχουν γραμμική σχέση μεταξύ τους. Από την άλλη μεριά, αν υπάρχει άριστη γραμμική σχέση μεταξύ δύο εξαρτημένων μεταβλητών τότε δεν έχει νόημα η συμμετοχή και των δύο στην ανάλυση.

(δ) Κανονικότητα
Απαιτείται η πολυμεταβλητή κατανομή όλων των εξαρτημένων μεταβλητών να είναι κανονική. Αυτό ελέγχεται με τον υπολογισμό της απόστασης Mahalanobis (Mahalanobis distance), στατιστικό με το οποίο εντοπίζονται παρατηρήσεις στις εξαρτημένες μεταβλητές που επηρεάζουν αρνητικά την πολυμεταβλητή κανονικότητα. Για τον έλεγχο αυτό προτείνεται επίπεδο εμπιστοσύνης 0,001.

(ε) Ομοιογένεια
Αν το πλήθος των παρατηρήσεων σε κάθε έναν συνδυασμό των ανεξαρτήτων παραγόντων δεν είναι το ίδιο, τότε ο ερευνητής πρέπει να ελέγξει την ισότητα των πινάκων που περιέχουν τις συνδιακυμάνσεις του συνόλου των εξαρτημένων μεταβλητών στους συνδυασμούς των ανεξάρτητων. Αυτό μπορεί να γίνει με το στατιστικό M του Box (Box’s M test).

8.1 Παράδειγμα

Μια επιχείρηση επιθυμεί να βρει την απαραίτητη διάρκεια εκπαίδευσης των πωλητών της για τη μέγιστη αποτελεσματικότητά τους στην προώθηση ενός νέου προϊόντος στις πωλήσεις παράλληλα με την προώθηση ενός νέου είδους ασφάλειας ζωής για τους πελάτες της. Για να το πετύχει αυτό επιλέγει με τυχαίο τρόπο 31 πωλητές τους οποίους χωρίζει σε τρεις ομάδες (Α, n = 10), (Β, n = 11), (Γ, n = 10), τους εκπαιδεύει για 1, 2 και 3 ημέρες αντίστοιχα, παρακολουθεί τις πωλήσεις τους και καταγράφει τα συμβόλαια τους για τον επόμενο μήνα. Τα δεδομένα στην R, καταχωρούνται σε τέσσερις μεταβλητές, τη μεταβλητή group (Ομάδα εκπαίδευσης) και τιμές 1 (Α), 2 (Β), 3 (Γ), τη μεταβλητή gender όπου καταχωρείται το φύλο του πωλητή με τιμές 0 (Γυναίκα) και 1 (Άνδρας), τη μεταβλητή perform (Πωλήσεις) όπου καταχωρούνται οι πωλήσεις του κάθε πωλητή και την μεταβλητή contracts όπου καταχωρείται το πλήθος των συμβολαίων που σύναψε ο πωλητής.

Αρχικά, χρησιμοποιούνται οι παρακάτω εντολές για την καταχώρηση των δεδομένων

perform.A = c(63.3, 68.3, 86.7, 52.8, 75, 58, 69.5, 32.7, 60.9, 58.2)
contracts.A = c(24, 30, 33, 19, 30, 22, 28, 13, 17, 18)
A.gender = c(1, 1, 1, 0, 0, 0, 1, 0, 0, 0)

perform.B = c(72.8, 88.2, 80.8, 71.3, 81.5, 47.6, 81, 81.4, 83, 76, 74.5)
contracts.B = c(25, 26, 29, 22, 41, 15, 40, 37, 39, 29, 26)
B.gender = c(1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1)

perform.C = c(82.3, 89.7, 81, 85.1, 74.1, 75.9, 74.7, 81.1, 76.4, 81.8)
contracts.C = c(38, 42, 48, 42, 40, 45, 33, 31, 47, 38)
C.gender = c(0, 0, 1, 0, 1, 1, 1, 1, 1, 0)

sales = c(perform.A, perform.B, perform.C)
contracts = c(contracts.A , contracts.B , contracts.C)
gender = c(A.gender, B.gender, C.gender)
group = c(rep(1, 10), rep(2, 11), rep(3, 10))

gender =factor(gender, levels = c(0,1), labels=c("Γυναίκα", "Άνδρας")) 
group = factor(group, levels = c(1,2,3), labels = c("Α", "Β", "Γ"))

manova.df = data.frame(group, gender, sales, contracts)

Hmisc::label(manova.df$group) = 'Ομάδα εκπαίδευσης'
Hmisc::label(manova.df$gender) = 'Φύλο'
Hmisc::label(manova.df$sales) = 'Πωλήσεις'
Hmisc::label(manova.df$contracts) = 'Συμβόλαια'

Επιπλέον, ως προς τις μεταβλητές «Πωλήσεις» και «Συμβόλαια» υπολογίζουμε το συντελεστή συσχέτισης Pearson και δημιουργούμε το αντίστοιχο διάγραμμα διασποράς με τις εντολές

plot(contracts, sales, xlab = "Πωλήσεις", ylab = "Συμβόλαια", 
     main = "Συσχέτιση μεταξύ των πωλήσεων και των συμβολαίων", 
     cex = 1.3, cex.main = 1.3, cex.axis = 1.3, cex.lab = 1.3)

my_cor_table(manova.df, c('contracts', 'sales'))
Pearson Correlation Coefficients (*:p<.05, **:p<.01, ***:p<.001)
Συμβόλαια Πωλήσεις
Συμβόλαια - .763***
Πωλήσεις .763*** -

Παρατηρούμε πως οι μεταβλητές «Πωλήσεις» και «Συμβόλαια»:
• Είναι εννοιολογικά συμβατές καθώς και οι δύο περιγράφουν την αποτελεσματικότητα του πωλητή.
• Δεν έχουν μη γραμμική σχέση όπως φανερώνεται από το αντίστοιχο διάγραμμα διασποράς, άρα μπορούμε να τις χρησιμοποιήσουμε μαζί με πολυμεταβλητό γραμμικό μοντέλο για να τις προβλέψουμε.
• Δεν έχουν πλήρη γραμμική σχέση (Pearson’s r(31) = 0,763, p < .001) άρα έχει νόημα να τις θεωρήσουμε ως δύο ξεχωριστές αναπαραστάσεις της ικανότητας ενός πωλητή.

Στα πλαίσια της πολυμεταβλητής ανάλυσης διακύμανσης (MANOVA), ελέγχουμε τις παρακάτω στατιστικές υποθέσεις:

Η0,1: Ο παράγοντας ομάδα δεν επιδρά σημαντικά στην αποτελεσματικότητα στις πωλήσεις (όπως αυτή εκφράζεται από τις πωλήσεις και τα συμβόλαια).


Η0,2: Ο παράγοντας φύλο δεν επιδρά σημαντικά στην αποτελεσματικότητα στις πωλήσεις (όπως αυτή εκφράζεται από τις πωλήσεις και τα συμβόλαια).


Η0,3: Ο συνδυασμός των παραγόντων φύλου και ομάδας δεν επιδρά σημαντικά στην αποτελεσματικότητα στις πωλήσεις.


Δίνοντας την εντολή

my_contigency_table(manova.df, c('group', 'gender'))
Cross tabulation of Ομάδα εκπαίδευσης, Φύλο
gender  
Γυναίκα Άνδρας   Sum
group
  Α 6 4   10
  Β 7 4   11
  Γ 4 6   10
  Sum 17 14   31

παρατηρούμε πως δεν έχουμε το ίδιο πλήθος παρατηρήσεων σε κάθε ένα συνδυασμό των ανεξαρτήτων παραγόντων (unbalanced design), άρα θα πρέπει να είμαστε ιδιαίτερα προσεκτικοί στην επαλήθευση των προϋποθέσεων της πολυμεταβλητής κανονικότητας του ζεύγους των μεταβλητών (πωλήσεις, πλήθος συμβολαίων) όπως αυτό εκφράζεται από τη απόσταση Mahalanobis όπως και από την ομοιογένεια όπως αυτή εκφράζεται από το στατιστικό M του Box.

8.2 Έλεγχος προϋποθέσεων

Πριν την υλοποίηση της δοκιμασίας πρέπει να επαληθευτεί η προϋπόθεση της ομοιογένειας των τιμών της κάθε μίας εξαρτημένης μεταβλητής ως προς τις ανεξάρτητες μεταβλητές.

options(contrasts = c("contr.treatment", "contr.poly"))
my_check_homogeneity_of_column(manova.df, 'sales', c('group', 'gender'))
Standard Deviation of Πωλήσεις among levels of Interaction of Ομάδα εκπαίδευσης, Φύλο
Level Α.Γυναίκα Β.Γυναίκα Γ.Γυναίκα Α.Άνδρας Β.Άνδρας Γ.Άνδρας
Group Size 6 7 4 4 4 6
SD 13.80 3.60 3.62 10.20 12.70 3.09
Levene test of variance equality (homogeneity) of Πωλήσεις over Interaction of Ομάδα εκπαίδευσης, Φύλο
F(5, 25) = 0.865, p = 0.518. Homogeneity hypothesis is confirmed.
Bartlett’s test of variance equality (homogeneity) of Πωλήσεις over Interaction of Ομάδα εκπαίδευσης, Φύλο
If you have strong evidence that your data do in fact come from a normal, or nearly normal, distribution, then Bartlett’s test has better performance.
c2(5) = 16.859, p = 0.005. Homogeneity hypothesis is not supported!
my_check_homogeneity_of_column(manova.df, 'contracts', c('group', 'gender'))
Standard Deviation of Συμβόλαια among levels of Interaction of Ομάδα εκπαίδευσης, Φύλο
Level Α.Γυναίκα Β.Γυναίκα Γ.Γυναίκα Α.Άνδρας Β.Άνδρας Γ.Άνδρας
Group Size 6 7 4 4 4 6
SD 5.78 6.21 2.31 3.77 4.97 7.28
Levene test of variance equality (homogeneity) of Συμβόλαια over Interaction of Ομάδα εκπαίδευσης, Φύλο
F(5, 25) = 0.939, p = 0.473. Homogeneity hypothesis is confirmed.
Bartlett’s test of variance equality (homogeneity) of Συμβόλαια over Interaction of Ομάδα εκπαίδευσης, Φύλο
If you have strong evidence that your data do in fact come from a normal, or nearly normal, distribution, then Bartlett’s test has better performance.
c2(5) = 4.185, p = 0.523. Homogeneity hypothesis is confirmed.
leveneTest(sales ~ group * gender, manova.df, center = mean)
## Levene's Test for Homogeneity of Variance (center = mean)
##       Df F value Pr(>F)
## group  5  2.0796 0.1017
##       25
leveneTest(contracts ~ group * gender, manova.df, center = mean)
## Levene's Test for Homogeneity of Variance (center = mean)
##       Df F value Pr(>F)
## group  5   1.782 0.1532
##       25

Από το output των τελευταίων εντολών, συνάγουμε πως η ομοιογένεια των τιμών της κάθε μίας εξαρτημένης μεταβλητής ως προς τις ανεξάρτητες μεταβλητές δεν απορρίπτεται (Πωλήσεις: F(5, 25) = 2,08, p = 0,102 και Συμβόλαια: F(5, 25) = 1,78, p = 0,153).

Συνεχίζουμε, υπολογίζοντας το στατιστικό M του Box, με τη συνάρτηση boxM. Πριν την υλοποίηση της δοκιμασίας πρέπει να δημιουργηθεί μία νέα μεταβλητή όπου θα φιλοξενήσει την αλληλεπίδραση των δύο παραγόντων.

library(magrittr)
group1 = unclass(manova.df$group) %>% as.numeric 
gender1 = unclass(manova.df$gender) %>% as.numeric 
manova.df$interaction = factor(group1 * 10 + gender1)

Τώρα, αρκεί η εντολή

library(heplots)
boxM(manova.df[, 3:4], manova.df[, "interaction"])
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  manova.df[, 3:4]
## Chi-Sq (approx.) = 30.434, df = 15, p-value = 0.01045

Παρατηρούμε, πως, από το στατιστικό M του Box απορρίπτεται η υπόθεση πως οι πινάκες των συνδιακυμάνσεων των εξαρτημένων μεταβλητών είναι ίσοι μεταξύ των συνδυασμών της μεταβλητής ομάδα (χ2(15) = 30,434, p = 0,010), γεγονός που σύμφωνα με τη βιβλιογραφία, επιβάλει την ανάγνωση του στατιστικού Pillai Trace ως περισσότερο έγκυρης επιλογής έναντι των υπολοίπων, για τον έλεγχο των στατιστικών υποθέσεων Η0,1, Η0,2, Η0,3. (Tabachnick & Fidell, 1983)

8.3 Υλοποίηση

Εκτελούμε τις εντολές

library(car)
lm.manova = lm(cbind(sales,contracts) ~ group + gender + group*gender, manova.df,   contrasts=list(group=contr.sum, gender=contr.sum))
man.res = Manova(lm.manova, type="III", test ="Pillai")
summary(man.res, multivariate=TRUE)
## 
## Type III MANOVA Tests:
## 
## Sum of squares and products for error:
##               sales contracts
## sales     1907.2208  628.8167
## contracts  628.8167  796.6310
## 
## ------------------------------------------
##  
## Term: (Intercept) 
## 
## Sum of squares and products for the hypothesis:
##               sales contracts
## sales     156735.24  66384.42
## contracts  66384.42  28116.78
## 
## Multivariate Tests: (Intercept)
##                  Df test stat approx F num Df den Df     Pr(>F)    
## Pillai            1   0.98831 1014.265      2     24 < 2.22e-16 ***
## Wilks             1   0.01169 1014.265      2     24 < 2.22e-16 ***
## Hotelling-Lawley  1  84.52207 1014.265      2     24 < 2.22e-16 ***
## Roy               1  84.52207 1014.265      2     24 < 2.22e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------
##  
## Term: group 
## 
## Sum of squares and products for the hypothesis:
##              sales contracts
## sales     1380.314  1254.427
## contracts 1254.427  1346.956
## 
## Multivariate Tests: group
##                  Df test stat  approx F num Df den Df     Pr(>F)    
## Pillai            2 0.7605186  7.669726      4     50 6.7098e-05 ***
## Wilks             2 0.3210793  9.177519      4     48 1.3686e-05 ***
## Hotelling-Lawley  2 1.8603590 10.697064      4     46 3.2618e-06 ***
## Roy               2 1.7119068 21.398835      2     25 3.8377e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------
##  
## Term: gender 
## 
## Sum of squares and products for the hypothesis:
##              sales contracts
## sales     39.86608 16.223382
## contracts 16.22338  6.602057
## 
## Multivariate Tests: gender
##                  Df test stat  approx F num Df den Df  Pr(>F)
## Pillai            1 0.0208619 0.2556761      2     24 0.77648
## Wilks             1 0.9791381 0.2556761      2     24 0.77648
## Hotelling-Lawley  1 0.0213063 0.2556761      2     24 0.77648
## Roy               1 0.0213063 0.2556761      2     24 0.77648
## 
## ------------------------------------------
##  
## Term: group:gender 
## 
## Sum of squares and products for the hypothesis:
##               sales contracts
## sales     1261.3998  780.5544
## contracts  780.5544  574.9240
## 
## Multivariate Tests: group:gender
##                  Df test stat  approx F num Df den Df     Pr(>F)    
## Pillai            2 0.5728122  5.016966      4     50 0.00177001 ** 
## Wilks             2 0.4763247  5.387202      4     48 0.00115074 ** 
## Hotelling-Lawley  2 0.9962500  5.728438      4     46 0.00079302 ***
## Roy               2 0.8788744 10.985931      2     25 0.00037695 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

8.4 Αποτελέσματα

Από την ανάγνωση του μέρους του output με τίτλο Multivariate Tests στη γραμμή Pillai Trace, συμπεραίνουμε πως ως προς την αποτελεσματικότητα ενός πωλητή όπως αυτή εκφράζεται από τις προώθηση του νέου προϊόντος αλλά και το πλήθος των ασφαλιστικών συμβολαίων που αυτός συνάπτει με τους πελάτες προκύπτουν τα εξής συμπεράσματα:

  1. Η στατιστική υπόθεση Η0,1 απορρίπτεται, F(4, 50) = 7,670, p < 0,001, δηλαδή υπάρχει σημαντική διαφοροποίηση στην αποτελεσματικότητα στις πωλήσεις μεταξύ των τριών ομάδων.

  2. Η στατιστική υπόθεση Η0,2 δεν απορρίπτεται, F(2, 24) = 0,523, p = 0,776, δηλαδή δεν υπάρχει σημαντική διαφοροποίηση στην αποτελεσματικότητα στις πωλήσεις μεταξύ των δύο φύλων.

  3. Η στατιστική υπόθεση Η0,3 απορρίπτεται, F(4, 50) = 5,017, p = 0,002, δηλαδή η αποτελεσματικότητα στις πωλήσεις των τριών ομάδων δεν διαφοροποιούνται με τον ίδιο τρόπο στους άνδρες και στις γυναίκες.

Φανερά, τα παραπάνω αποτελέσματα είναι μόνο η αρχή της ανάλυσης. Τώρα πρέπει να ακολουθήσει η διευκρίνηση των διαφοροποιήσεων. Ο ερευνητής όμως πρέπει να γνωρίζει πως η συνέχιση της ερμηνείας πρέπει να γίνει με πολύ προσοχή καθώς η γραμμική εξάρτηση των δύο εξαρτημένων μεταβλητών αφαιρεί την αξιοπιστία από τους επιμέρους υπολογισμούς των σημαντικοτήτων των ελέγχων που παρουσιάζονται σε αυτόν τον πίνακα, δηλαδή με απλά λόγια δεν είναι σωστό να βγάλουμε συμπεράσματα απλά από την ανάγνωση των p των επιμέρους 2-way ANOVA για τις πωλήσεις και τα συμβόλαια ξεχωριστά.

Για την αντιμετώπιση αυτού του προβλήματος, στη βιβλιογραφία προτείνεται η stepdown analysis μέθοδος κατά την οποία πρώτα επιλέγεται η πιο σημαντική από τις εξαρτημένες μεταβλητές και ακολουθούν οι υπόλοιπες μία προς μία διορθώνοντας κάθε φορά την επίδραση των προηγουμένων πιο σημαντικών μεταβλητών.

Ο εντοπισμός της πλέον σημαντικής από τις δύο μεταβλητής γίνεται με ποιοτικά κριτήρια βάσει των πεποιθήσεων του ερευνητή ή ποσοτικά με την σύγκριση των λόγων F των επιμέρους ελέγχων: Ο μεγαλύτερος αντιστοιχεί στην πιο στατιστικά σημαντική μεταβλητή. Στην περίπτωσή μας επιλέγουμε τη δεύτερη μέθοδο.

Στο σημείο αυτό πρέπει να εκτελεστούν οι εντολές

lm.anova.sales = lm(sales ~ group + gender + group*gender, manova.df, 
                                         contrasts=list(group=contr.sum, gender=contr.sum))
Anova(lm.anova.sales, type="III") 
## Anova Table (Type III tests)
## 
## Response: sales
##              Sum Sq Df   F value    Pr(>F)    
## (Intercept)  156735  1 2054.4977 < 2.2e-16 ***
## group          1380  2    9.0466  0.001107 ** 
## gender           40  1    0.5226  0.476458    
## group:gender   1261  2    8.2673  0.001754 ** 
## Residuals      1907 25                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

και

lm.anova.contracts = lm(contracts ~ group + gender + group*gender, 
               manova.df, contrasts=list(group=contr.sum, gender=contr.sum))

Anova(lm.anova.contracts, type="III")
## Anova Table (Type III tests)
## 
## Response: contracts
##               Sum Sq Df  F value    Pr(>F)    
## (Intercept)  28116.8  1 882.3654 < 2.2e-16 ***
## group         1347.0  2  21.1352 4.231e-06 ***
## gender           6.6  1   0.2072  0.652910    
## group:gender   574.9  2   9.0212  0.001123 ** 
## Residuals      796.6 25                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Παρατηρούμε πως τόσο στην εξάρτηση από την ομάδα (group) όσο και στο εξάρτηση από το συνδυασμό ομάδας και φύλου (\(group*gender\)), η μεγαλύτερη τιμή του στατιστικού F αντιστοιχεί στα Συμβόλαια και ως εκ τούτου θεωρούμε αυτή τη μεταβλητή ως βασικότερη. Διαβάζοντας τα αντίστοιχα στατιστικά από τον πίνακα καταλήγουμε πως η ομάδα αποτελεί σημαντικό παράγοντα διαφοροποίησης των τιμών των Συμβολαίων F(2, 25) = 21,135, p < 0,001, ωστόσο το ίδιο δεν συμβαίνει για το φύλο F(1, 25) = 0,207, p = 0,653, ενώ ο συνδυασμός ομάδας και φύλου αποτελεί σημαντικό παράγοντα διαφοροποίησης των τιμών της, F(2, 25) = 9,021, p = 0,001.

Απομένει ο εντοπισμός της επιρροής της ομάδας και του φύλου στην άλλη εξαρτημένη μεταβλητή (πωλήσεις) διορθώνοντας ως προς την επίδραση των συμβολαίων, κάτι που μπορεί να επιτευχθεί με τη διαδικασία ANCOVA.

Αρκεί να εκτελεστούν οι εντολές

lm.post.hoc.ancova = lm(contracts ~ group + gender + group * gender + sales, manova.df, contrasts=list(group=contr.sum, gender=contr.sum))
Anova(lm.post.hoc.ancova, type="III")
## Anova Table (Type III tests)
## 
## Response: contracts
##              Sum Sq Df F value   Pr(>F)    
## (Intercept)   16.59  1  0.6758 0.419134    
## group        475.47  2  9.6820 0.000826 ***
## gender         0.23  1  0.0095 0.923209    
## sales        207.32  1  8.4434 0.007755 ** 
## group:gender 155.37  2  3.1638 0.060321 .  
## Residuals    589.31 24                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Από το τελευταίο output, παρατηρούμε πως ακόμα και ύστερα από τη διόρθωση των πωλήσεων (sales) στο μέγεθος των συμβολαίων (contracts) η επίδραση της ομάδας στις τιμές των συμβολαίων είναι στατιστικά σημαντική F(2, 24) = 9,682, p = 0,001, ενώ ο συνδυασμός group*gender δεν προκύπτει πλέον στατιστικά σημαντικός παράγοντας διαφοροποίησης των τιμών της F(2, 24) = 3,164, p = 0,060.

Συμπερασματικά, η ομάδα είναι σημαντικός παράγοντας διαφοροποίησης της αποτελεσματικότητας ενός πωλητή όπως αυτή εκφράζεται από τις συνολικές πωλήσεις αλλά και το πλήθος των ασφαλιστικών συμβολαίων. Επιπλέον, στα Συμβόλαια υπάρχει διαφορετική εξέλιξη των τιμών μεταξύ ανδρών και γυναικών στις τρεις ανεξάρτητες ομάδες, ενώ δεν προκύπτει το ίδιο συμπέρασμα για τις πωλήσεις.

Η περαιτέρω ανάλυση θα μπορούσε να συνεχίσει με τον εντοπισμό ομοιογενών ομάδων, ωστόσο αυτό είναι ένα θέμα που καλύφθηκε στην 2-way ANOVA παράγραφο.

8.5 Παρατηρήσεις

8.5.1 Η απόσταση Mahalanobis

Δεν αναφερθήκαμε στον τρόπο υπολογισμού της απόστασης Mahalanobis για κάθε μία από τις παρατηρήσεις μας, διαδικασία που μας επιτρέπει να εντοπίζουμε ιδιαίτερα ακραίες παρατηρήσεις ως προς τις τιμές του συνόλου των εξαρτημένων μεταβλητών που πρόκειται να συμμετάσχουν στη διαδικασία MANOVA, και οι οποίες μειώνουν την εγκυρότητα της μεθόδου, καθώς επηρεάζουν την πολυμεταβλητή κανονικότητα των εξαρτημένων μεταβλητών. Για τον υπολογισμό της απόστασης ο ερευνητής αρκεί να εκτελέσει τον παρακάτω κώδικα:

mahal.manova.df = as.numeric(as.matrix(manova.df[, -c(1, 2)]))
mahal.manova.df = matrix(mahal.manova.df,nrow = 31,ncol = 3)
mahal_dist = mahalanobis(mahal.manova.df,colMeans(mahal.manova.df), cov(mahal.manova.df))
manova.df$MD = round(mahal_dist, 1)

Με τον παραπάνω κώδικα, υπολογίζεται η απόσταση Mahalanobis και καταχωρείται στη νέα μεταβλητή MD του data frame manova.df. Μία γραφική επισκόπηση των τιμών της μπορεί να προκύψει ως εξής:

hist(manova.df$MD, main = "Απόσταση Mahalanobis", 
                       col = c("grey45"),  xlab = "Mahalanobis distance", 
                       ylab = "Συχνότητα",  ylim = c(0,25))

Από το ιστόγραμμα των τιμών παρατηρούμε πως υπάρχει κάποια τιμή αρκετά μεγαλύτερη από τις άλλες.

Το ερώτημα που ανακύπτει είναι: τι σημαίνει ιδιαίτερα μεγάλη απόσταση Mahanalobis;

Η απάντηση είναι: Μία απόσταση θεωρείται ιδιαίτερα μεγάλη όταν είναι μεγαλύτερη από την κρίσιμη τιμή της κατανομής χ2 με β..ε. όσους και το πλήθος των εξαρτημένων μεταβλητών σε επίπεδο εμπιστοσύνης 0,001.

Εκτελούμε την εντολή

qchisq(0.999, 2, lower.tail=TRUE)
## [1] 13.81551

και βρίσκουμε πως η κρίσιμη τιμή για την περίπτωσή μας είναι το 13,82. Τέλος, δημιουργούμε μία νέα μεταβλητή που στο data frame στην οποία θα καταχωρηθεί η πληροφορία για το αν η παρατήρηση είναι ιδιάζουσα τιμή ή όχι.

manova.df$outlier = "No"
manova.df$outlier[manova.df$MD > 13.82]
## character(0)

Παρατηρούμε πως καμία από τις παρατηρήσεις μας δεν έχει απόσταση Mahanalonis μεγαλύτερη από τον αριθμό 13,82, δηλαδή η προϋπόθεση της πολυμεταβλητής κανονικότητας στις εξαρτημένες μεταβλητές είναι αποδεκτή για τα δεδομένα του παραδείγματός μας.

9 Δραστηριότητες

9.1 Εφαρμογή στο data frame parents

Στα πλαίσια διπλωματικής εργασίας με τίτλο “Ψυχολογικά Χαρακτηριστικά Γονέων με Παιδιά με Αναπηρία” καταγράφηκαν ορισμένα δημογραφικά και ψυχομετρικά χαρακτηριστικά γονεων από 100 οικογένειες δύο νομων της χώρας. Τα δεδομένα είναι σε ένα αρχείο Excel διαθέσιμα για λήψη εδώ και θα αξιοποιηθούν στη συνέχεια σε δραστηριότητες εξοικείωσης.
1) Να κάνετε λήψη τα δεδομένα.
2) Να εισάγετε τα δεδομένα στο R Studio ως data.frame με όνομα parents
3) Εκτελέστε την εντολή head(parents)
4) Μετατρέψτε όσες στήλες πρέπει από αλφαριθμητικές (character) σε παράγοντες (factors) με τις εντολές:

cols <- c("age_category", "gender", "marital", "place_live", "grad_level_new", "job", "schico", "schico_drugs", "child_gender", "child_disease")
parents[cols] <- lapply(parents[cols], factor)  
  1. Εκτελέστε την εντολή summary(parents) για να εξοικειωθείτε με τα περιεχόμενα του αρχείου
    στ) Μεταβάλετε τις ετικέτες του παράγοντα child_gender από “αγορι” και “κορη” σε “Αγόρι” και “Κορίτσι¨αντίστοιχα, με τις εντολές
parents$child_gender = plyr::revalue(parents$child_gender, c("αγορι"="Αγόρι", "κορη"="Κορίτσι"))
  1. Μεταβάλετε τη σειρά των επιπέδων του παράγοντα child_disease με την εντολή:
parents$child_disease = ordered(parents$child_disease, levels = c('Σε αρχικά στάδια η νόσος', 2, 3, 4, 'Πολύ βαριά'))
  1. Εισάγετε ετικέτες σε ορισμένες από τις μεταβλητές που θα χρησιμοποιήσουμε στη συνέχεια.
library(Hmisc)
label(parents$gender) = 'Φύλο'
label(parents$marital) = 'Οικογενειακή κατάσταση'
label(parents$job) = 'Εργασιακή κατάσταση'
label(parents$age_category) = 'Ηλικιακή κατηγορία'
label(parents$schico) = 'Καταφυγή σε ψυχολόγο'
label(parents$schico_drugs) = 'Λήψη φαρμάκων'
label(parents$child_gender) = 'Φύλο Παιδιού'
label(parents$child_disease) = 'Εκτίμηση ασθένειας παιδιού'
label(parents$hdhq.total) = 'Επιθετικότητα (HDHQ)'
label(parents$scl_90_gen_deikt_sympt) = 'Ψυχοπαθολογία (GSI)'
label(parents$oas_total_score) = 'Εξωτερική Ντροπή'
label(parents$ess_total) = 'Εσωτερική Ντροπή'
label(parents$stai_present) = 'Άγχος ως κατάσταση'
label(parents$stai_general) = 'Άγχος ως χαρακτηριστικό της προσωπικότητας'
label(parents$grad_level_new) = 'Εκπαιδευτικό επίπεδο'

Θέλοντας να ελέγξουμε την υπόθεση πως η ψυχοπαθολογία των γονέων διαφοροποιείται σημαντικά μεταξύ των επιπέδων του εκπαιδευτικού επιπέδου.

  1. Ελέγξτε την κανονικότητα της κατανομής της ψυχοπαθολογίας scl_90_gen_deikt_sympt ανά επίπεδο του παράγοντα grad_level_new
my_check_normality_of_column_by_group(parents$scl_90_gen_deikt_sympt, parents$grad_level_new)
  1. Ελέγξτε την ομοιογένεια της κατανομής της ψυχοπαθολογίας ανά επίπεδο του παράγοντα grad_level_new.
car::leveneTest(parents$scl_90_gen_deikt_sympt, parents$grad_level_new, center = mean)
HH::hov(scl_90_gen_deikt_sympt ~ grad_level_new, data = parents)

Οπτικοποιήστε τις τυχόν διαφορές με την εντολή:

library(HH)
hovPlot(scl_90_gen_deikt_sympt ~ grad_level_new, data = parents)
  1. Eφαρμόστε την παραμετρική δοκιμασία ANOVA.
model1 = my_ANOVA(scl_90_gen_deikt_sympt ~ grad_level_new, data = parents)

Εφαρμόστε διαγνωστικά εργαλεία για να ελέγξετε την κατανομή των υπολοίπων. Οπτικοποιήστε την κατανομή των υπολοίπων ώστε να επαληθεύσετε ή να επιβεβαιώσετε οπτικά την απόσταση από την κανονικότητα. Αξιοποιήστε τη συνάρτηση my_histFrequency_with_normal_curve.

model.residuals = model1$res.aov$residuals
my_histFrequency_with_normal_curve(model.residuals)
  1. Αν οι προϋποθέσεις της κανονικότητας και της ομοιογένειας τηρούνται τότε προχωρήστε στην ανάγνωση του output.

Αν η ομοιογένεια δεν γίνεται αποδεκτή αλλά γίνεται αποδεκτή η κανονικότητα των υπολοίπων τότε εφαρμόστε τη δοκιμασία Welch:

library(onewaytests)
welch.test(scl_90_gen_deikt_sympt ~ grad_level_new, parents)

Αν ούτε η ομοιογένεια ούτε η κανονικότητα των υπολοίπων δεν γίνονται αποδεκτές τότε ελέγξτε την ομοιότητα των κατανομών με τη συνάρτηση my_hist_multiple_density_sidebyside.

my_hist_multiple_density_sidebyside(parents$scl_90_gen_deikt_sympt, parents$grad_level_new)

Προχωρήστε με τη δοκιμασία Kruskal Wallis:

library(onewaytests)
welch.test(scl_90_gen_deikt_sympt ~ grad_level_new, parents)
  1. Εκ των υστέρων έλεγχοι.

Ερμηνεύστε τα αποτελέσματα των ελέγχων Fisher’s LSD και Tukey’s HSD.

Μη παραμετρικές δοκιμασίες Dunn και Conover:

Για την καλύτερη παρουσίαση των αποτελεσμάτων, δημιουργούμε μία νέα μεταβλητή με συντομότερες ετικέτες:

parents$grad_level_new2 = plyr::revalue(parents$grad_level_new, c("Απόφοιτος Δημοτικού ή Αναλφάβητος"="Δημοτικό", 
                                                            "Απόφοιτος Γυμνασίου"="Γυμνάσιο", 
                                                            "Απόφοιτος Λυκείου"="Λυκείο", 
                                                            "Απόφοιτος Α.Ε.Ι. / Τ.Ε.Ι."="Α.Ε.Ι. / Τ.Ε.Ι."))

Δοκιμασία Dunn:

dunn.test::dunn.test(parents$scl_90_gen_deikt_sympt, parents$grad_level_new2)

Δοκιμασία Conover:

conover.test::conover.test(parents$scl_90_gen_deikt_sympt, parents$grad_level_new2)

Θέλουμε να ελέγξουμε την υπόθεση πως η ψυχοπαθολογία των γονέων διαφοροποιείται σημαντικά μεταξύ των επιπέδων του φύλου και του εκπαιδευτικού επιπέδου, με παράλληλο έλεγχο τυχόν αλληλεπίδρασης των δύο μεταβλητών.

Τόσο η υπόθεση της ομοιογένειας όσο και της κανονικότητας θα ελεγχθεί για τα υπόλοιπα του μοντέλου, μετά τους υπολογισμούς.

  1. Ελέγξτε αν πρόκειται για ισορροπημένο (balanced) ή μη ισορροπημένο σχεδιασμό (unbalanced design)
table(parents$grad_level_new, parents$gender)

Αν δεν προκύπτει ισορροπημένος σχεδιασμός, τότε μπορείτε να δημιουργήσετε μία νέα μεταβλητή η οποία θα έχει πιο ισορροπημένες κατηγορίες.

parents$grad_level_new2 = plyr::revalue(parents$grad_level_new, c("Απόφοιτος Δημοτικού ή Αναλφάβητος"="Α/θμια", "Απόφοιτος Γυμνασίου"="Β/θμια", "Απόφοιτος Λυκείου"="Β/θμια", "Απόφοιτος Α.Ε.Ι. / Τ.Ε.Ι."="Γ/θμια"))
  1. Eφαρμόστε την παραμετρική δοκιμασία ANOVA.
model2 = my_ANOVA(scl_90_gen_deikt_sympt ~ grad_level_new + gender + grad_level_new:gender, data = parents)

14α) Οπτικοποιήστε την κατανομή των υπολοίπων ώστε να επαληθεύσετε ή να επιβεβαιώσετε οπτικά την απόσταση από την κανονικότητα. Αξιοποιήστε τη συνάρτηση my_histFrequency_with_normal_curve: my_histFrequency_with_normal_curve.

my_histFrequency_with_normal_curve(model2$res.aov$residuals)

14β) Ελέγξτε την κυρτότητα και τη συμμετρία της κατανομής.

14γ) Ερμηνεύστε τους στατιστικούς ελέγχους κανονικότητας

  1. Για την ερμηνεία τυχόν αλληλεπίδρασης προχωρήστε στο αντίστοιχο διάγραμμα.
parents$grad_level_new = ordered(parents$grad_level_new, levels = c('Απόφοιτος Δημοτικού ή Αναλφάβητος', 
                                                                    'Απόφοιτος Γυμνασίου', 'Απόφοιτος Λυκείου',
                                                                    'Απόφοιτος Α.Ε.Ι. / Τ.Ε.Ι.'))
interaction.plot(parents$grad_level_new, parents$gender, 
                 parents$scl_90_gen_deikt_sympt, type="b", col=c("red","blue"), 
                 legend=T, trace.label = 'Φύλο', lty=c(1,2), lwd=2, pch=c(18,24), 
                 xlab="Εκπαιδευτικό Επίπεδο",  ylab="SCL - 90",
                 main="Αλληλεπίδραση φύλου και εκπαιδευτικού επιπέδου στην ψυχοπαθολογία", 
                 ylim = c(0, 2))
  1. Για την ερμηνεία τυχόν κύριας επιρροής προχωρήστε σε ένα από τα παρακάτω διαγράμματα
plot(parents$grad_level_new, parents$scl_90_gen_deikt_sympt, type="b", col=c("red","blue"), 
                 legend=T, trace.label = 'Φύλο', 
                 xlab="Εκπαιδευτικό Επίπεδο",  ylab="SCL - 90",
                 main="Επίδραση εκπαιδευτικού επιπέδου στην ψυχοπαθολογία", 
                 ylim = c(0, 2.5))

ή

gplots::plotmeans(scl_90_gen_deikt_sympt~grad_level_new2, data = parents, xlab="Εκπαιδευτικό Επίπεδο",ylab="SCL-90", main="Επίδραση εκπαιδευτικού επιπέδου στην ψυχοπαθολογία (95% C.I.)")
  1. Επαναλάβατε την μελέτη για τις κλίμακες του άγχους και της επιθετικότητας.

9.2 Εφαρμογή στο data frame acceptance

Στα πλαίσια διπλωματικής με αντικείμενο τον αξιολόγηση μίας νέας παρέμβασης με στόχο την ευαισθητοποίηση των μαθητών σε διάφορες μορφές αναπηρίας, χρησιμοποιήθηκε το ερωτηματολόγιο Acceptance Scale for Kindergartners—Revised (ASK-R). Στην έρευνα συμμετείχαν 32 μαθητές από τους οποίους οι 12 ήταν στην πειραματική ομάδα και οι 20 στην ομάδα ελέγχου. Για κάθε έναν μαθητή καταγράφηκε το φύλο (1 = Αγόρι, 2 = Κορίτσι), η ηλικία του (σε μήνες) και οι αποκρίσεις του πριν και μετά την παρέμβαση. Τα αποτελέσματα αποθηκεύτηκαν σε ένα αρχείο SPSS που είναι διαθέσιμα για λήψη εδώ και θα αξιοποιηθούν στη συνέχεια σε δραστηριότητες εξοικείωσης.
1) Να κάνετε λήψη τα δεδομένα.
2) Να εισάγετε τα δεδομένα στο R Studio ως data.frame με όνομα acceptance
3) Εκτελέστε την εντολή print(acceptance)
4) Μετατρέψτε όσες στήλες πρέπει από αλφαριθμητικές (character) σε παράγοντες (factors) με τις εντολές:

cols <- c("Group", "gender")
acceptance[cols] <- lapply(acceptance[cols], factor)  
  1. Εκτελέστε την εντολή summary(acceptance) για να εξοικειωθείτε με τα περιεχόμενα του αρχείου
  2. Προσθέστε ετικέτες στις μεταβλητές taxi, diglosia, math.diskolies, με τις εντολές
acceptance$Group = plyr::revalue(acceptance$Group, c("1"="Πειραματική", "2"="Ελέγχου"))
acceptance$gender = plyr::revalue(acceptance$gender, c("1"="Αγόρι", "2"="Κορίτσι"))

Θέλουμε να ελέγξουμε την υπόθεση πως η παρέμβαση βελτίωσε σημαντικά την αποδοχή των παιδιών απέναντι στους ανάπηρους συμμαθητές τους, με παράλληλο έλεγχο ως προς το φύλο και ως προς την ομάδα ελέγχου.

  1. Εφαρμόστε Repeated ANOVA ελέγχοντας μόνο ως προς την παρέμβαση.
my_Repeated_ANOVA(acceptance, withincols =  c('ASKpre', 'ASKpost'))
  1. Εφαρμόστε Repeated ANOVA ελέγχοντας μόνο ως προς την ομάδα. Τι παρατηρείτε;
my_Repeated_ANOVA(acceptance, betweencols = c('Group'), withincols =  c('ASKpre', 'ASKpost'))
  1. Εφαρμόστε Repeated ANOVA ελέγχοντας μόνο ως προς το φύλο. Τι παρατηρείτε;
my_Repeated_ANOVA(acceptance, betweencols = c('gender'), withincols =  c('ASKpre', 'ASKpost'))
  1. Εκτελέστε την εντολή
my_Repeated_ANOVA(acceptance, betweencols = c("gender", 'Group'), withincols =  c('ASKpre', 'ASKpost'))
  1. Εφαρμόστε ένα t - test για να ελέγξετε μόνο την διαφοροποίηση μεταξύ των δύο μετρήσεων. Τι παρατηρείτε;
my_t_test_paired_samples(acceptance, c('ASKpre', 'ASKpost'))

9.3 Εφαρμογή στο data frame students

Στα πλαίσια διπλωματικής με αντικείμενο τον τρόπο επεξεργασίας μαθηματικών προβλημάτων από μαθητές Λυκείου, καταγράφηκαν για 200 μαθητές, ορισμένα δημογραφικά χαρακτηριστικά τους όπως και ο τρόπος με τον οποίον αντιμετώπισαν 12 μαθηματικά προβλήματα. Τα αποτελέσματα αποθηκεύτηκαν σε ένα αρχείο SPSS που είναι διαθέσιμα για λήψη εδώ και θα αξιοποιηθούν στη συνέχεια σε δραστηριότητες εξοικείωσης.
1) Να κάνετε λήψη τα δεδομένα.
2) Να εισάγετε τα δεδομένα στο R Studio ως data.frame με όνομα students
3) Εκτελέστε την εντολή head(students)
4) Μετατρέψτε όσες στήλες πρέπει από αλφαριθμητικές (character) σε παράγοντες (factors) με τις εντολές:

cols <- c("taxi", "tmima", "diglosia", "math.diskolies")
students[cols] <- lapply(students[cols], factor)  
  1. Εκτελέστε την εντολή summary(students) για να εξοικειωθείτε με τα περιεχόμενα του αρχείου
  2. Προσθέστε ετικέτες στις μεταβλητές taxi, diglosia, math.diskolies, με τις εντολές
students$taxi = plyr::revalue(students$taxi, c("1"="A", "2"="B", "3"="Γ"))
students$diglosia = plyr::revalue(students$diglosia, c("1"="Όχι", "2"="Ναι"))
students$math.diskolies = plyr::revalue(students$math.diskolies, c("0"="Όχι", "1"="Ναι"))
  1. Αξιοποιήστε τις μεθόδους της 8ης ενότητας για να μελετήσετε τα δεδομένα.

  1. Επικοινωνία: Επαμεινώνδας Διαμαντοπουλος, Τμήμα ΗΜ/ΜΥ, Δ.Π.Θ. Email: . Τηλ: 25410 79757, 6944683327↩︎