Έκδοση: 03 / 05 / 2024

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

list.of.packages <- c("pROC", "ggplot2", "OptimalCutpoints", "plyr", "htmlTable", "ppcor", "psych", "car", "rms", "lmtest", "broom", "see", "performance")
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 Ευθεία

Η ευθεία είναι το γεωμετρικό αντικείμενο που προκύπτει αν ενώσουμε δύο σημεία και προεκτείνουμε απεριόριστα και ως προς τις δύο κατευθύνσεις.

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

y = α x + β.

Το α ονομάζεται κλίση της ευθείας και δηλώνει τη μεταβολή στην τιμή του y όταν το x αυξηθεί κατά 1 μονάδα.

Το β ονομάζεται σταθερός όρος και αντιστοιχεί στην τιμή που παίρνει το y όταν το x είναι 0.

2 Διάγραμμα Διασποράς

Ένας καθηγητής θέλει να μελετήσει τη σχέση μεταξύ της ενδιάμεσης έκθεσης προόδου και της τελικής επίδοσης των φοιτητών του τμήματός του. Για το σκοπό αυτό θα αξιοποιήσει ιστορικά στοιχεία για 14 φοιτητές για τους οποίους κατέγραψε την επίδοση στην ενδιάμεση πρόοδο (μεταβλητή test) και την τελική επίδοση του φοιτητή (μεταβλητή exam).

test = c(5.5, 6, 7, 5, 8, 2, 9, 10, 2, 3, 4, 6.5, 8.5, 1)
exam = c(6.5, 6, 8, 7.5, 7, 4, 8, 10, 1, 5, 5, 6, 9, 5)
my.df = data.frame(test, exam)
Hmisc::label(my.df$test) = 'Πρόοδος'
Hmisc::label(my.df$exam) = 'Τελικές εξετάσεις'

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

sjPlot::plot_scatter(my.df, test, exam)

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

Είναι εμφανές πως τα σημεία βρίσκονται «κοντά» σε κάποια ευθεία, γεγονός που ερμηνεύεται ως υψηλή αναλογικότητα μεταξύ των Χ και Y.

3 Συντελεστής συσχέτισης Pearson

Ο συντελεστής συσχέτισης r δύο μεταβλητών Χ και Υ είναι ένας αριθμός μεταξύ -1 και +1 ο οποίος δείχνει το είδος και την ένταση της γραμμικότητας δύο συνε­χών μεταβλητών, δηλαδή μας ενημερώνει για το πόσο ρεαλιστική είναι η υπόθεση πως Υ = α*Χ για κάποιο πραγματικό αριθμό α. Ο συντελεστής συσχέτισης αξιολογείται σύμφωνα με τον παρακάτω πίνακα:

Μέτρο συντελεστή r Είδος σχέσης
\[0 ≤ |r| < 0,3\] Μη γραμμική σχέση
\[0,3 ≤ |r| < 0,7\] Ασθενής γραμμική σχέση \[0,3 ≤ r < 0,7\] Ασθενής θετική γραμμική σχέση
\[ \] \[-0,7 ≤ r < -0,3\] Ασθενής αρνητική γραμμική σχέση
\[0,7 ≤ |r| ≤ 1\] Ισχυρή γραμμική σχέση \[0,7 ≤ r ≤ 1\] Ισχυρή θετική γραμμική σχέση
\[-1 ≤ r ≤ -0,7\] Ισχυρή αρνητική γραμμική σχέση
Ενδεικτικά παραδείγματα, παρουσιάζονται στην παρακάτω εικόνα (πηγή: Wikipedia)
Παραδείγματα συντελεστών συσχέτισης
Παραδείγματα συντελεστών συσχέτισης
Στην παρακάτω εικόνα παρουσιάζεται η μεγάλη ποικιλία που μπορεί να λάβει ο συντελεστής συσχέτισης. Προσέξτε πως η μη γραμμική σχέση δύο μεταβλητών, δεν αντικατοπτρίζεται σε αυξημένο συντελεστσή.
Παραδείγματα συντελεστών συσχέτισης
Παραδείγματα συντελεστών συσχέτισης

Υπολογισμός
Για τον υπολογισμό του συντελεστή συσχέτισης αρκεί να εφαρμόσουμε τη συνάρτηση my.cor.table.

my_cor_table(my.df, c('test', 'exam'))
Pearson Correlation Coefficients (*:p<.05, **:p<.01, ***:p<.001)
Πρόοδος Τελικές εξετάσεις
Πρόοδος - .853***
Τελικές εξετάσεις .853*** -

Παρατηρούμε πως r(14) = 0,853, p < 0,001, άρα επαληθεύεται η ισχυρή θετική γραμμική σχέση μεταξύ του βαθμού προόδου και του βαθμού στην τελική εξέταση.

4 Απλή Γραμμική Παλινδρόμηση (Simple Linear Regression - SLR).

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

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

Για τον υπολογισμό του γραμμικού μοντέλου αρκεί να εφαρμόσουμε τη συνάρτηση my.lm.

my.model = my_lm(my.df, exam ~ test)
my.model$htmlreport
Linear regression results for Τελικές εξετάσεις
Model: exam ~ const + test
  95% C.I.
B SE t p   Lower Upper
Constant 2.501 0.747 3.348 0.006   0.873 4.129
Πρόοδος 0.684 0.121 5.651 <0.001   0.42 0.947
Coefficient of determination
R2 = 0.727
Homoscedasticity Breusch–Pagan test
x2(1) = 2.905, p = 0.088. Homoscedasticity assumption is not rejected. Model seems to be valid.
Normality tests of models’ residuals
Statistic Description
Skewness & Kurtosis
  Skewness -0.704 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.619 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.943, p = 0.455 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.136, p = 0.69 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.

Οι συντελεστές της εξίσωσης είναι στην στήλη B. Παρατηρούμε, ότι η ευθεία γραμμικής παλινδρόμησης έχει εξίσωση:

Εξετάσεις = 2,501 + 0,684 * Πρόοδος

Για την γραφική αναπαράσταση της ευθείας γραμμικής παλινδρόμησης μπορούμε να εφαρμόσουμε ξανά τη συνάρτηση plot_scatter:

sjPlot::plot_scatter(my.df, test, exam, fit.line = 'lm')

Αξιοποίηση της γραμμικής εξίσωσης πρόβλεψης
Η εξίσωση γραμμικής παλινδρόμησης αξιοποιείται με τρεις τρόπους.

Α) Ο πρώτος είναι για την πρόβλεψη των τιμών της εξαρτημένης από την/τις ανεξάρτητες. Για παράδειγμα, η πρόβλεψη για την τελική επίδοση στις εξετάσεις για έναν φοιτητή που είχε επίδοση 7 στην Πρόοδο είναι:

  Εξετάσεις = 2,501 + 0,684 * 7 = 7,289 ≈ 7,3

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

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

Γ) Ο τρίτος τρόπος είναι η ερμηνεία του σταθερού όρου ο οποίος αντιστοιχεί στην τιμή της εξαρτημένης μεταβλητής όταν η τιμή της ανεξάρτητης είναι ίση με μηδέν.

Στο παράδειγμά μας, προκύπτει πως ο φοιτητής που θα γράψει 0 στην Πρόοδο περιμένουμε να γράψει 2,501 στις εξετάσεις.

Συντελεστής Προσδιορισμού
Η συνολική αποτελεσματικότητα του γραμμικού μοντέλου όπως αυτή εκφράζεται με το συντελεστή προσδιορισμού R2 (Multiple R square) είναι ίση με 0,727 = 72,7% δηλαδή το 72,7% της συνολικής μεταβλητότητας της τελικής επίδοσης ενός φοιτητή ερμηνεύεται από την επίδοσή του στην Πρόοδο. Συμπεραίνουμε ανάλογα πως το υπόλοιπο 0,273 = 27,3% οφείλεται σε άλλους παράγοντες που δεν καταγράφονται στην έρευνά μας.

Στήλες t, p
Το στατιστικό t = B / Std.Error χρησιμοποιείται για να ελεγχθεί για κάθε μία από τις ανεξάρτητες μεταβλητές που συμμετέχουν στο γραμμικό μοντέλο, η στατιστική υπόθεση:

Η0: Β = 0,

έναντι της:

Η1: Β ≠ 0.

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

Η στατιστική σημαντικότητα που αντιστοιχεί σε αυτόν τον (δίπλευρο) έλεγχο βρίσκεται στη στήλη p. Οι βαθμοί ελευθερίας για τον έλεγχο αυτό είναι df = n - 2 = 14 - 2 = 12.

Για τα δεδομένα του παραδείγματός μας παρατηρούμε πως t(12) =5, 651, p < 0,001, άρα η Πρόοδος είναι σημαντική μεταβλητή στο γραμμικό μοντέλο. Επιπλέον, παρατηρούμε πως η σταθερά (Constant) του μοντέλου είναι στατιστικά μη μηδενική ποσότητα στο γραμμικό μοντέλο που μελετούμε (t(12) = 3,348, p = 0,005).

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

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

Για τα δεδομένα του παραδείγματός μας καταγράφεται ότι αύξηση μίας μονάδας στην Πρόοδο αντιστοχεί σε μέση αύξηση από 0,42 έως 0,947 μονάδες την επίδοση στις τελικές εξετάσεις, ενώ επίδοση στην Πρόοδο ίση με 0 αντιστοχεί σε μέση επίδοση στις τελικές εξετάσεις από 0,873 έως 4,129 μονάδες.

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

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

Δοκιμασία ελέγχου κανονικότητας των υπολοίπων
Η υπόθεση πως τα υπόλοιπα (ypredicted - yobserved) έχουν κανονική κατανομή ελέγχεται με τη δοκιμασία Lilliefors η οποία είναι μία περισσότερο έγκυρη εκδοχή της δοκιμασίας Kolmogorov-Smirnov. Από το Output παρατηρούμε πως η υπόθεση αυτή δεν απορρίπτεται.

Δοκιμασία ελέγχου ομοσκεδαστικότητας
Η υπόθεση της ομοσκεδαστικότητας ελέγχεται με τη δοκιμασία Breusch-Pagan. Από το Output παρατηρούμε πως η υπόθεση της ομοσκεδαστικότητας δεν απορρίπτεται.

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

fit = my.model$fit
round(fit$residuals, 3)
## Τελικές εξετάσεις 
##      1      2      3      4      5      6      7      8      9     10     11 
##  0.239 -0.603  0.713  1.581 -0.971  0.132 -0.654  0.662 -2.868  0.448 -0.236 
##     12     13     14 
## -0.945  0.688  1.815

Ένα επιπλέον χρήσιμο διάγραμμα αναπαράστασης είναι το εξής:

sjPlot::plot_residuals(fit)

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

my_histFrequency_with_normal_curve(fit$residuals)

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

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

4.1.1 Τυποποιημένοι Συντελεστές (Standardized Coefficients)

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

lm(scale(exam) ~ scale(test))
## 
## Call:
## lm(formula = scale(exam) ~ scale(test))
## 
## Coefficients:
## (Intercept)  scale(test)  
##  -5.844e-17    8.526e-01

4.1.2 Γραμμικό μοντέλο χωρίς σταθερό όρο

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

lm(exam ~ test - 1)
## 
## Call:
## lm(formula = exam ~ test - 1)
## 
## Coefficients:
##  test  
## 1.047

4.1.3 Η προϋπόθεση της κανονικότητας

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

Παράδειγμα
Εξαρτημένη μεταβλητή Υ: δίτιμη (0 - 1) ποιοτική μεταβλητή, ανεξάρτητη μεταβλητή Χ = Υ + έναν τυχαίο αριθμό από τη Ν(0,1) κατανομή.

Φανερά, τόσο η Υ όσο και η Χ δεν ακολουθούν την κανονική κατανομή αλλά τα υπόλοιπα του γραμμικού μοντέλου Υ = α*Χ + β την ακολουθούν εκ κατασκευής. Η αναπαράσταση του παραπάνω παραδείγματος μπορεί να γίνει με τον εξής κώδικα R:

library(psych)
x = rbinom(100, 1, 0.3)
y = rnorm(length(x), 5 + x * 5, 1)
scatter.hist(x, y, correl=TRUE, density=TRUE, ellipse=FALSE, xlab="x", ylab="y")

fit = lm(y~x)
summary(fit)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.34784 -0.64030 -0.04472  0.41942  2.89136 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.8523     0.1102   44.04   <2e-16 ***
## x             5.0454     0.1918   26.30   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9019 on 98 degrees of freedom
## Multiple R-squared:  0.8759, Adjusted R-squared:  0.8747 
## F-statistic:   692 on 1 and 98 DF,  p-value: < 2.2e-16

4.1.4 Standardized και Studentized Residual Value

Ως υπόλοιπο (residual) μίας παρατήρησης ονομάζεται η διαφορά της πραγματικής τιμής της εξαρτημένης μεταβλητής από την προβλεπόμενη βάσει του γραμμικού μοντέλου ενώ το τυποποιημένο υπόλοιπο (standardized residual) είναι το αποτέλεσμα της διαίρεσης του υπολοίπου με την τυπική απόκλιση όλων των υπολοίπων.
Το τυποποιημένο κατά student υπόλοιπο (studentized residual) είναι το αποτέλεσμα της διαίρεσης του υπολοίπου με το τυπικό σφάλμα του.

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

5 Πολλαπλή Γραμμική Παλινδρόμηση (Multiple Linear Regression - MLR).

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

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

Παράδειγμα
Ένας καθηγητής θέλει να δημιουργήσει γραμμικό μοντέλο πρόβλεψης της τελικής επίδοσης των φοιτητών του τμήματός του από την ενδιάμεση εξέταση προόδου, το φύλο και τις απουσίες από το μάθημα. Για το σκοπό αυτό θα αξιοποιήσει ιστορικά στοιχεία για 14 φοιτητές για τους οποίους κατέγραψε την επίδοση στην Πρόοδο (μεταβλητή test), την τελική επίδοση του φοιτητή (μεταβλητή exam), το φύλο του φοιτητή (μεταβλητή gender) και τις απουσίες του φοιτητή από τις διαλέξεις του μαθήματος (μεταβλητή absences).

test = c(5.5, 6, 7, 5, 8, 2, 9, 10, 2, 3, 4, 6.5, 8.5, 1)
exam = c(6.5, 6, 8, 7.5, 7, 4, 8, 10, 1, 5, 5, 6, 9, 5)
gender = c(0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1)
gender = factor(gender, levels = c(0, 1), labels = c("Γυναίκα", "Άνδρας"))
absences = c(2, 3, 1, 1, 2, 5, 0, 0, 7, 5, 4, 3, 0, 4)
exam.df = data.frame(test, exam, gender, absences)
Hmisc::label(exam.df$test) = 'Βαθμός προόδου'
Hmisc::label(exam.df$exam) = 'Βαθμός εξετάσεων'
Hmisc::label(exam.df$gender) = 'Φύλο'
Hmisc::label(exam.df$absences) = 'Απουσίες'

Με τις παρακάτω εντολές υπολογίζουμε τους συντελεστές συσχέτισης Pearson μεταξύ των τριών συνεχών μεταβλητών test, exam και absences:

my_cor_table(exam.df, c('exam', 'test', 'absences'))
Pearson Correlation Coefficients (*:p<.05, **:p<.01, ***:p<.001)
Βαθμός εξετάσεων Βαθμός προόδου Απουσίες
Βαθμός εξετάσεων - .853*** -.967***
Βαθμός προόδου .853*** - -.866***
Απουσίες -.967*** -.866*** -

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

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

• Συντελεστής συσχέτισης Εξετάσεις με Πρόοδο r = 0,853,
• Συντελεστής συσχέτισης Εξετάσεις με Απουσίες r = - 0,967,

Η αποθαρρυντική παρατήρηση είναι πως ο συντελεστής συσχέτισης μεταξύ των δύο ανεξάρτητων παρατηρήσεων Πρόοδο και Απουσίες που σκοπεύουμε να εισάγουμε στη γραμμική εξίσωση είναι ίσος με r = -0,866, δηλαδή οι δύο ανεξάρτητες παρατηρήσεις έχουν ισχυρή γραμμική σχέση μεταξύ τους, γεγονός που ονομάζεται συγγραμικότητα (collinearity) και δημιουργεί σοβαρό πρόβλημα στον υπολογισμό των συντελεστών του μοντέλου.

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

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

my.model = my_lm(exam.df, exam ~ test + gender + absences)
my.model$htmlreport
Linear regression results for Βαθμός εξετάσεων
Model: exam ~ const + test + gender + absences
  95% C.I.
B SE t p   Lower Upper
Constant 8.466 1.101 7.688 <0.001   6.012 10.919
Βαθμός προόδου 0.052 0.125 0.417 0.686   -0.226 0.33
Φύλο 0.389 0.47 0.827 0.428   -0.659 1.437
Απουσίες -1.018 0.179 -5.696 <0.001   -1.416 -0.62
Coefficient of determination
R2 = 0.94
Homoscedasticity Breusch–Pagan test
x2(1) = 2.197, p = 0.533. Homoscedasticity assumption is not rejected. Model seems to be valid.
Normality tests of models’ residuals
Statistic Description
Skewness & Kurtosis
  Skewness 0.398 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.099 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.91, p = 0.158 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.225, p = 0.052 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.

Η επιρροή των παραγόντων αναπαριστάται διαγραμματικά από το ακόλουθο διάγραμμα

library(sjPlot)    
plot_model(my.model$fit, axis.labels=c("Απουσίες", "Φύλο", 'Πρόοδος'), show.values=TRUE, show.p=TRUE, title="Παράγοντες μοντέλου πρόβλεψης τελικού βαθμού")

Ερμηνεία του output

H συνολική αποτελεσματικότητα του γραμμικού μοντέλου περιγράφεται από το συντελεστή προσδιορισμού R square. Παρατηρούμε πως R2 = 0,940 = 94% δηλαδή το 94% της συνολικής μεταβλητότητας της τελικής επίδοσης ενός φοιτητή ερμηνεύεται από την επίδοσή του στην Πρόοδο, το φύλο του και τις απουσίες του από το μάθημα στη διάρκεια του εξαμήνου. Ανάλογα, το υπόλοιπο 0,06 = 6% οφείλεται σε άλλους παράγοντες που δεν καταγράφονται στην έρευνά μας.

Η δοκιμασία της ανάλυσης διακύμανσης (ANOVA) μας διαβεβαιώνει πως το σύνολο της μεταβλητότητας του τελικού βαθμού εξέτασης που ερμηνεύεται από το γραμμικό μοντέλο που ορίσαμε είναι στατιστικά σημαντικός F(3, 10) = 51,8, p < 0,001.

Το γραμμικό μοντέλο πρόβλεψης παρουσιάζεται στον τμήμα του output «Coefficients».

Η εξίσωσης πρόβλεψης και κάποια σχόλια…
Στη στήλη Estimate παρουσιάζονται οι συντελεστές του μοντέλου. Η εξίσωση πρόβλεψης που προκύπτει από τα δεδομένα μας είναι:

Εξετάσεις = 8,466 + 0,052 * Πρόοδος + 0,389 * Φύλο - 1,018 *Απουσίες

Καθώς κάποια/ες από τις ανεξάρτητες μεταβλητές που συμμετέχουν στο γραμμικό μοντέλο δεν έχουν σημαντικά διαφορετικούς από το μηδέν συντελεστές (στον πίνακα Coefficients, το αντίστοιχο p είναι μεγαλύτερο από 0,05), το γραμμικό μοντέλο μπορεί να χρησιμοποιηθεί αποκλειστικά για την πρόβλεψη των τιμών της εξαρτημένης μεταβλητής. Με απλά λόγια, η διαδικασία προβλέπει με επιτυχία τις τιμές της εξαρτημένης αλλά δεν μπορεί να εξηγήσει τη σχέση της εξαρτημένης με κάθε μία από τις ανεξάρτητες.

Για παράδειγμα, η πρόβλεψη για την τελική επίδοση στις εξετάσεις για έναν φοιτητή που είχε επίδοση 7 στην πρόοδο, έχει 2 απουσίες και είναι αγόρι (gender = 1) είναι:

Εξετάσεις = 8,466 + 0,052 * 7+ 0,389 * 1 - 1,018 *2 ≈ 7,2

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

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

(β) να εφαρμόσει ανάλυση σε κύριες συνιστώσες (Principal Component Analysis - PCA) και να ορίσει με τον τρόπο αυτό έναν ή περισσότερους παράγοντες (factors) από γραμμικούς συνδυασμούς των ισχυρά συσχετισμένων ανεξαρτήτων μεταβλητών οι οποίοι παράγοντες θα είναι εκ κατασκευής ασυσχέτιστοι και θα συμμετάσχουν στο γραμμικό μοντέλο στη θέση των αρχικών μεταβλητών.

5.1 Προϋποθέσεις εγκυρότητας

Το ιστόγραμμα των υπολοίπων παράγεται με την εξής εντολή:

my_histFrequency_with_normal_curve(my.model$fit$residuals, main = 'Ιστόγραμμα υπολοίπων')

Κάποια επιπλέον διαγράμματα με τα οποία αξιολογείται η ποιότητα του γραμμικού μοντέλου και που είναι απαραίτητα στην περίπτωση όπου υπάρχουν περισσότερες από μία επεξηγηματικές μεταβλητές δίνονται με την εντολή plot(fit) η οποία παράγει 4 διαγράμματα. Εναλλακτικά, το σύνολο των 4 διαγραμμάτων είναι διαθέσιμο από το αποτέλεσμα της συνάρτησης my_lm, ως εξής:

my.model$diagnostic.plots

5.1.1 Διάγραμμα Normal Q – Q

Το διάγραμμα Quantile - Quantile φέρνει σε αντιπαράθεση τις τιμές των υπολοίπων με αυτές που θα περιμέναμε να υπάρχουν αν τα υπόλοιπα ακολουθούν κανονική κατανομή. Επιθυμητό είναι τα σημεία του διαγράμματος να βρίσκονται πάνω στην ευθεία y = x, δηλαδή την ευθεία όπου κάθε της σημείο έχει ίδιες συντεταγμένες. Σημαντική απόκλιση των σημείων από την ευθεία υποδεικνύει μη κανονική κατανομή των υπολοίπων.

Μία καλύτερη ιδέα για τον τρόπο που “διαβάζεται” το διάγραμμα Quantile - Quantile δίνεται στο παρακάτω διάγραμμα: Πηγή: Glen_b (https://stats.stackexchange.com/users/805/glen-b), How to interpret a QQ plot, URL (version: 2021-09-04): https://stats.stackexchange.com/q/101290)

Για την καλύτερη κατανόηση του τρόπου που μετασχηματίζεται το Q-Q plot καθώς η κατανομή απομακρύνεται από την κανονική δίνεται εδώ.

Στην περίπτωσή μας, παρατηρούμε πως υπάρχουν τρεις παρατηρήσεις (στις γραμμές των δεδομένων 7, 8 και 10) οι οποίες δημιουργούν ένα μικρό πρόβλημα στην κανονικότητα των υπολοίπων.

5.1.2 Διάγραμμα Residuals vs Fitted

Το διάγραμμα Residuals vs Fitted απεικονίζει την διασπορά των υπολοίπων σε όλο το εύρος τιμών που προβλέπει το γραμμικό μοντέλο. Για να είναι έγκυρη η χρήση του γραμμικού μοντέλου πρέπει τα υπόλοιπα να έχουν την ίδια διασπορά σε όλο το μήκος τιμών του οριζόντιου άξονα. Μία καλύτερη ιδέα για τον τρόπο που “διαβάζεται” το διάγραμμα Residuals vs Fitted δίνεται στο παρακάτω διάγραμμα: Ενδεικτικά διαγράμματα Residuals vs Fitted Πηγή: Linear Models with R (Chapman & Hall/CRC Texts in Statistical Science)

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

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

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

Στην περίπτωσή μας, η προϋπόθεση της ομοσκεδαστικότητας μπορεί να γίνει αποδεκτή για τα δεδομένα μας καθώς όπως παρατηρείται από το αντίστοιχο διάγραμμα διασποράς (Residuals vs Fitted) τα σημεία έχουν ανάλογη κύμανση για όλο το εύρος τιμών της τιμής που προβλέπεται από το γραμμικό μοντέλο.

5.1.3 Διάγραμμα Residuals vs Leverage

Η απόσταση Cook δείχνει το πόσο πολύ θα μεταβαλλόταν τα υπόλοιπα των άλλων παρατηρήσεων αν απομακρυνόταν η ίδια η παρατήρηση. Όσο μεγαλύτερη είναι η τιμή της τόσο πιο ισχυρή είναι η επιρροή του σημείου στον υπολογισμό των συντελεστών του μοντέλου. Από το διάγραμμα Residuals vs Leverage (υπόλοιπα παρατήρησης vs ικανότητα μόχλευσης) παρατηρούμε ότι υπάρχουν δεν υπάρχουν σημεία με σημαντική ικανότητα μόχλευσης, δηλαδή σημεία με απόσταση Cook μεγαλύτερη από το 1.

5.1.4 Διάγραμμα Scale-Location

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

Στην περίπτωσή μας, από το διάγραμμα Scale-Location υποδεικνύεται ένα πρόβλημα ετεροσκεδαστικότητας που οφείλεται στις παρατηρήσεις 8, 9 και 10.

Περισσότερες λεπτομέρειες για το είδος των διαγραμμάτων που προκύπτει από την plot(fit) παρουσιάζονται εδώ ή εδώ

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

Ανοχή (Tolerance) και συγγραμικότητα.

Από τις στήλες t, p, φανερώνεται μία μη επιθυμητή πραγματικότητα. Οι συντελεστές των μεταβλητών Φύλο και Πρόοδος δεν προκύπτουν να είναι στατιστικώς σημαντικοί (t(12) = 0,827, p = 0,428 και t(12) = 0,417, p = 0,686 αντίστοιχα), γεγονός που οφείλεται στην έντονη συγγραμικότητα μεταξύ των μεταβλητών Πρόοδος και Απουσίες που παρατηρήσαμε στην αρχή της διαδικασίας με το συντελεστή συσχέτισης του Pearson.

Η συγγραμικότητα των μεταβλητών μπορεί να μετρηθεί με το παράγοντα μεγέθυνσης της κύμανσης (Variance Inflation Factor - VIF). Στην R, αυτός είναι εύκολο να υπολογιστεί με την εντολή vif της βιβλιοθήκης συναρτήσεων car

library(car)
thefit = my.model$fit
vif(thefit)
##     test   gender absences 
## 3.994111 1.862616 4.799177

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

Ο παράγοντας πληθωρισμού της κύμανσης ορίζεται ως ο αντίστροφος της ανοχής (VIF = 1 / Tolerance) και ερμηνεύεται με αντίστροφο τρόπο. Δεν υπάρχει κάποιο γενικά αποδεκτό όριο για το μέγεθος του VIF που να καταδεικνύει σημαντικό πρόβλημα συγγραμικότητας, με πολλούς ερευνητές να χρησιμοποιούν το 10 ως το μέγιστό όριο του παράγοντα αυτού που αντιστοιχεί σε ελάχιστη τιμή 0,1 για την ανοχή του ίδιου παράγοντα.

5.2 Μέθοδος Stepwise σταδιακής εισαγωγής – εξαγωγής

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

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

fit = lm(exam ~ test + gender + absences, data = exam.df)
step(fit, direction="both") # ή direction = “forward” ή direction = “backward”
## Start:  AIC=-9.29
## exam ~ test + gender + absences
## 
##            Df Sum of Sq     RSS      AIC
## - test      1    0.0708  4.1432 -11.0462
## - gender    1    0.2782  4.3506 -10.3623
## <none>                   4.0724  -9.2875
## - absences  1   13.2123 17.2848   8.9507
## 
## Step:  AIC=-11.05
## exam ~ gender + absences
## 
##            Df Sum of Sq    RSS      AIC
## - gender    1     0.273  4.417 -12.1517
## <none>                   4.143 -11.0462
## + test      1     0.071  4.072  -9.2875
## - absences  1    38.065 42.208  19.4498
## 
## Step:  AIC=-12.15
## exam ~ absences
## 
##            Df Sum of Sq    RSS     AIC
## <none>                   4.417 -12.152
## + gender    1     0.273  4.143 -11.046
## + test      1     0.066  4.351 -10.362
## - absences  1    62.941 67.357  23.993
## 
## Call:
## lm(formula = exam ~ absences, data = exam.df)
## 
## Coefficients:
## (Intercept)     absences  
##       8.966       -1.014

Ερμηνεία του output
Η εντολή step, υπολογίζει το δείκτη AIC (Akaike Information Criterion) για όλα τα πιθανά γραμμικά μοντέλα, επιλέγοντας τελικά το μοντέλο που μειώνει στο μεγαλύτερο βαθμό την τιμή του. Με το κριτήριο αυτό, η διαδικασία τελικά προτείνει μοντέλο πρόβλεψης με μοναδική ανεξάρτητη μεταβλητή τις απουσίες (absences). Η εξίσωση πρόβλεψης είναι η:

Εξετάσεις = 8,966 - 1,014 * Απουσίες,

και η ερμηνεία των συντελεστών είναι άμεση: Ένας φοιτητής δίχως απουσίες στη διάρκεια των μαθημάτων αναμένεται να έχει βαθμό 8,966 ≈ 9 στην τελική εξέταση ενώ για κάθε απουσία που πραγματοποιείται αναμένεται μείωση στο βαθμό κατά 1,014 ≈ 1 μονάδα.

6 Μερικός Συντελεστής Συσχέτισης (Partial Correlation Coefficient)

Ο μερικός συντελεστής συσχέτισης των δύο μεταβλητών Χ, Y ελέγχοντας ως προς την τρίτη μεταβλητή Ζ είναι ένας αριθμός μεταξύ -1 και +1 ο οποίος δείχνει τη συσχέτιση των τιμών των Χ και Y μετά την διόρθωση ως προς την γραμμική επιρροή της μεταβλητής Z πάνω σε αυτές. Πιο συγκεκριμένα αν:

Χ = α * Ζ + β + υ1
και
Υ = γ * Ζ + δ + υ2,

τότε, ο μερικός συντελεστής συσχέτισης των Χ, Υ ελέγχοντας ως προς τη επίδραση της Ζ, ορίζεται ως ο απλός συντελεστής συσχέτισης Pearson των υπολοίπων υ1 και υ2 των παραπάνω γραμμικών εξισώσεων.

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

Παράδειγμα
Ένας ερευνητής κατέγραψε τις επενδύσεις στο χώρο της δημόσιας υγείας 50 πόλεων (μεταβλητή funding), τις αναφορές των ασθενειών ανά 10.000 κατοίκους (μεταβλητή disease), και τις μέσες ημερήσιες επισκέψεις στα νοσοκομεία (μεταβλητή visits). Τα στοιχεία καταχωρούνται στις αντίστοιχες μεταβλητές με τις παρακάτω εντολές.

funding = c(155.33, 177.34, 165.09, 154.28, 185.56, 186.96, 198.63, 172.14, 198.25, 193.9, 157.03, 158.07, 182.09, 173.2, 151.84, 182.9, 197.09, 163.86, 192.3, 172.03, 166.97, 187.86, 184.23, 166.86, 174.64, 159.5, 153.18, 169.46, 169.9, 161.57, 183.99, 187.98, 195.88, 154.76, 177.6, 163.02, 190.34, 179.27, 192.2, 181.23, 159.41, 167.47, 151.38, 195.56, 162.69, 173.17, 176.52, 200.73, 190.53, 180.12)

disease = c(158.34, 157.23, 162.92, 130.59, 202.81, 221.43, 189.22, 166.42, 203.07, 198.57, 161.79, 168.82, 180.41, 178.52, 157.77, 212.44, 198.9, 174.25, 212.07, 149.45, 190.35, 198.38, 172.99, 186.81, 179.42, 148.46, 133.64, 152.11, 183.13, 163.32, 176.52, 178.6, 207.62, 143.88, 159.74, 172.06, 188.55, 202.31, 190.16, 201.54, 153.99, 153.99, 122.39, 184.35, 154.44, 179.7, 167.49, 164.54, 184.12, 178.94)

visits = c(152.13, 167.84, 162.21, 146.69, 186.93, 188.16, 195.39, 175.5, 197.04, 190.16, 156.47, 159.21, 185.65, 176.5, 156.97, 181.62, 205.51, 167.09, 192.71, 165.02, 162.07, 186.16, 189, 166.93, 171.45, 157.3, 150.41, 171.18, 168.88, 151.86, 185.03, 184.52, 197.07, 153.74, 184.18, 155.93, 197.91, 187.9, 194.07, 186.68, 150.64, 164.11, 148.32, 193.66, 162.36, 172.7, 172.51, 196.14, 191.1, 174.83)

my.df = data.frame(funding, disease, visits)
Hmisc::label(my.df$funding) = "Δαπάνη"
Hmisc::label(my.df$disease) = "Αναφερόμενες ασθένειες"
Hmisc::label(my.df$visits) = "Επισκέψεις"

Αρχικά, υπολογίζουμε τους συντελεστές συσχέτισης Pearson μεταξύ των τριών μεταβλητών.

my_cor_table(my.df, c('funding', 'disease', 'visits'))
Pearson Correlation Coefficients (*:p<.05, **:p<.01, ***:p<.001)
Δαπάνη Αναφερόμενες ασθένειες Επισκέψεις
Δαπάνη - .737*** .964***
Αναφερόμενες ασθένειες .737*** - .762***
Επισκέψεις .964*** .762*** -

Παρατηρούμε πως ο συντελεστής Pearson μεταξύ της δαπάνης μίας πόλης για το τοπικό σύστημα υγείας (funding) και του πλήθους αναφορών για ασθένειες (disease) στην πόλη αυτή είναι r(50) = 0,737, p < 0,001, στατιστικό που καταγράφεται ως ισχυρή θετική συσχέτιση του κόστους με το πλήθος των ασθενειών που αναφέρονται στην πόλη αυτή.

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

sjPlot::plot_scatter(my.df, funding, disease)

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

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

library(ppcor)
pcor.test(funding, disease, visits)
##     estimate   p.value  statistic  n gp  Method
## 1 0.01334953 0.9274621 0.09152795 50  1 pearson

Προκύπτει πως ο μερικός συντελεστής συσχέτισης μεταξύ των funding και disease ελέγχοντας ως προς τις τιμές της μεταβλητής visit είναι ίσος με 0,013 και μη στατιστικά σημαντικός p = 0,928, παρόλο που ο συντελεστής συσχέτισης του Pearson μεταξύ των funding και disease είναι ίσος με r(50) = 0,737, p < 0,001.

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

Πράγματι, η μεταβλητή visits έχει ισχυρή θετική γραμμική συσχέτιση τόσο με τη μεταβλητή funding (r(50) = 0,964, p < 0,001) όσο και με την disease (r(50) = 0,762, p < 0,001) πληροφορία που εξηγεί και την παρατηρούμενη θετική σχέση μεταξύ των funding και disease.

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

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

7 Λογιστική Παλινδρόμηση (Binomial Logistic Regression - BLR).

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

Παράδειγμα
Μία τράπεζα επιθυμεί να αναπτύξει ένα μοντέλο πρόβλεψης των πελατών που δηλώνουν αδυναμία πληρωμής μετά τη λήψη ενός δανείου.

Για το σκοπό αυτό θα αξιοποιήσει ιστορικά στοιχεία για 20 πελάτες για τους οποίους κατέγραψε το οικογενειακό εισόδημα σε ευρώ (μεταβλητή income), το μέγεθος του χρέους ως ποσοστό του οικογενειακού εισοδήματος (μεταβλητή deptinc) και το γεγονός της αδυναμίας αποπληρωμής (μεταβλητή default).

income = c(31, 55, 120, 25, 38, 16, 23, 64, 29, 100, 72, 61, 26, 176, 49, 25, 67, 28, 19, 41)
debtinc = c(17, 6, 3, 10, 4, 2, 5, 10, 16, 9, 8, 6, 2, 9, 9, 20, 31, 17, 24, 16)
default = factor(c(1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1), levels = c(0, 1), labels = c('Όχι', "Ναι"))

lr.data.frame = data.frame(income, debtinc, default)
Hmisc::label(lr.data.frame$income) = "Οικογενειακό εισόδημα"
Hmisc::label(lr.data.frame$debtinc) = "Χρέος ως ποσοστό του εισοδήματος"
Hmisc::label(lr.data.frame$default) = "Αδυναμία αποπληρωμής"

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

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

my_cor_table(lr.data.frame, c("income", "debtinc"))
Pearson Correlation Coefficients (*:p<.05, **:p<.01, ***:p<.001)
Οικογενειακό εισόδημα Χρέος ως ποσοστό του εισοδήματος
Οικογενειακό εισόδημα - -.164
Χρέος ως ποσοστό του εισοδήματος -.164 -

Βρίσκουμε πως Pearson r = -0,164, δηλαδή οι δύο μεταβλητές είναι γραμμικά ασυσχέτιστες, γεγονός που επιβεβαιώνεται οπτικά από το διάγραμμα διασποράς.

sjPlot::plot_scatter(lr.data.frame, income, debtinc)

Ο έλεγχος του λογιστικού μοντέλου γίνεται με τη συνάρτηση my.logistic.regression.

my.model = my_logistic_regression(lr.data.frame, default ~ income + debtinc)
my.model$htmlreport
Coding: 0, 1 = Όχι, Ναι
If you need reverse order please run the function by setting reverse.levels = TRUE
Logistic regression results for Αδυναμία αποπληρωμής
Model: Αδυναμία αποπληρωμής ~ Οικογενειακό εισόδημα + Χρέος ως ποσοστό του εισοδήματος
  95% C.I.
B SE z p Exp(B)   Lower Upper
Constant 0.69 1.839 0.375 0.708 1.993   0.054 73.239
Οικογενειακό εισόδημα -0.089 0.053 -1.668 0.095 0.915   0.824 1.016
Χρέος ως ποσοστό του εισοδήματος 0.264 0.124 2.127 0.033 1.302   1.021 1.66
Nagelkerke Coefficient of determination R2 = 0.721
Area Under Curve AUC = 0.939
Prediction quality of Αδυναμία αποπληρωμής at threshold 0.5
Observed  
Όχι Ναι   Sum
Predicted
  Όχι 10 2   12
  Ναι 1 7   8
Sum 11 9   20
Sensitivity (SNS): 0.778
Specificity (SPC): 0.909
You may also consider other thresholds by setting e.g. threshold = c(0.3, 0.7, 0.9)

Ερμηνεία του output
Αν p είναι η πιθανότητα για έναν πελάτη της τράπεζας να δηλώσει αδυναμία πληρωμής του δανείου, τότε, η διαδικασία glm προτείνει το μοντέλο πρόβλεψης:

log(p / 1- p) = 0,69 - 0,089*income + 0,264*debtinc

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

log(p / 1- p) = 0, 689 - 0, 089 * 31 + 0, 264 * 17 = 2,418

από όπου βρίσκουμε

p / 1- p = e2,418 = 11,22

ή

p = 11,22 - 11,22*p

ή

12,22*p = 11,22

ή

p = 0,92 = 92%.

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

print(round(fitted(my.model$fit), 3))
##     1     2     3     4     5     6     7     8     9    10    11    12    13 
## 0.918 0.067 0.000 0.750 0.162 0.448 0.490 0.085 0.911 0.003 0.026 0.041 0.250 
##    14    15    16    17    18    19    20 
## 0.000 0.214 0.977 0.948 0.936 0.995 0.779

Έχουμε τη δυνατότητα να προσθέσουμε τις προβλέψεις ως νέα στήλη στο data.frame με την εντολή:

lr.data.frame$predicted = round(fitted(my.model$fit), 3)

Το data.frame τώρα έχει αποκτήσει τη μορφή

print(lr.data.frame)
##    income debtinc default predicted
## 1      31      17     Ναι     0.918
## 2      55       6     Όχι     0.067
## 3     120       3     Όχι     0.000
## 4      25      10     Όχι     0.750
## 5      38       4     Ναι     0.162
## 6      16       2     Ναι     0.448
## 7      23       5     Όχι     0.490
## 8      64      10     Όχι     0.085
## 9      29      16     Ναι     0.911
## 10    100       9     Όχι     0.003
## 11     72       8     Όχι     0.026
## 12     61       6     Όχι     0.041
## 13     26       2     Όχι     0.250
## 14    176       9     Όχι     0.000
## 15     49       9     Όχι     0.214
## 16     25      20     Ναι     0.977
## 17     67      31     Ναι     0.948
## 18     28      17     Ναι     0.936
## 19     19      24     Ναι     0.995
## 20     41      16     Ναι     0.779

Ψεύδο-Συντελεστής Προσδιορισμού Nagelkerke R2
Ανάλογα με το συντελεστή προσδιορισμού της γραμμικής παλινδρόμησης, είναι δυνατό να υπολογιστεί ο ψεύδο-συντελεστής προσδιορισμού Nagelkerke R2 που ποσοτικοποιεί την ικανότητα πρόβλεψης του μοντέλου στο διάστημα 0 έως 1. Μία τιμή μικρότερη από 0,2 χαρακτηρίζει το μοντέλο ως φτωχό ως προς τη προβλεπτική του ικανότητα, μεταξύ 0,2 και 0,4 είναι αποδεκτό, ενώ μεγαλύτερη από 0,4 είναι πολύ καλή. Στην πετίπτωσή μας, η τιμή του Nagelkerke R2 = 0,721, κάτι που σημαίνει πως το μοντέλο μας προσαρμόζεται ικανοποιητικά στα δεδομένα.

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

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

lr.data.frame = lr.data.frame[order(lr.data.frame$default, lr.data.frame$predicted, decreasing = TRUE), ] 
print(lr.data.frame)
##    income debtinc default predicted
## 19     19      24     Ναι     0.995
## 16     25      20     Ναι     0.977
## 17     67      31     Ναι     0.948
## 18     28      17     Ναι     0.936
## 1      31      17     Ναι     0.918
## 9      29      16     Ναι     0.911
## 20     41      16     Ναι     0.779
## 6      16       2     Ναι     0.448
## 5      38       4     Ναι     0.162
## 4      25      10     Όχι     0.750
## 7      23       5     Όχι     0.490
## 13     26       2     Όχι     0.250
## 15     49       9     Όχι     0.214
## 8      64      10     Όχι     0.085
## 2      55       6     Όχι     0.067
## 12     61       6     Όχι     0.041
## 11     72       8     Όχι     0.026
## 10    100       9     Όχι     0.003
## 3     120       3     Όχι     0.000
## 14    176       9     Όχι     0.000

Ευαισθησία, προσδιοριστικότητα και άλλα μέτρα απόδοσης
Από το ενημερωμένο με την πιθανότητα πρόβλεψης dataframe, ορίζοντας κάποιο κατώφλι πιθανότητας μεταξύ 0 και 1, έχουμε τη δυνατότητα να συμπληρώσουμε έναν πίνακα συμπτώσεων από τον οποίο μπορούν να υπολογιστούν διάφορα μέτρα απόδοσης του λογιστικού μοντέλου:
• το ποσοστό ορθής ταξινόμησης (percentage accuracy in classification - PAC),
• η ευαισθησία (sensitivity) του μοντέλου (αληθώς ταξινομημένα ως θετικά - true positives),
• η προσδιοριστικότητα (specificity) του μοντέλου (αληθώς ταξινομημένα ως αρνητικά - true negatives),
• η θετική διαγνωστική του αξία δηλαδή το ποσοστό των ταξινομημένων ως θετικά προς το σύνολο των θετικών προβλέψεων του μοντέλου (positive predictive value - PPV) και
• η αρνητική διαγνωστική του αξία δηλαδή το ποσοστό των ταξινομημένων ως αρνητικά προς το σύνολο των αρνητικών συμβάντων (negative predictive value - NPV).

Στον επόμενο πίνακα παρουσιάζονται οι ορισμοί των παραπάνω δεικτών απόδοσης ενός μοντέλου πρόβλεψης μίας δίτιμης μεταβλητής.

Παρατήρηση: Όχι Παρατήρηση: Ναι Σύνολο
Πρόβλεψη: Όχι α γ P1 = α + β
Πρόβλεψη: Ναι β δ P2 = γ + δ
Σύνολο O1 = α + β O2 = γ + δ Σ = α + β + γ + δ

Βάσει του παραπάνω πίνακα είναι:

• Ποσοστό ορθής ταξινόμησης (PAC) = (α + δ) / Σ.
• Ευαισθησία (sensitivity) (true positives) = δ / Ο2.
• Προσδιοριστικότητα (specificity) (true negatives) = α / Ο1.
• Θετική διαγνωστική αξία (PPV) (positive predictive value) είναι ίση με δ / P2.
• Αρνητική διαγνωστική αξία (NPV) (negative predictive value) είναι ίση με α / P1.

Εφαρμογή για όριο (threshold) 0.5
Σχετικά με το προτεινόμενο μοντέλο πρόβλεψης, αν θεωρήσουμε ότι μία παρατήρηση θα προβλέπεται ως θετική (αδυναμία αποπληρωμής) όταν η predicted πιθανότητα ξεπερνάει το 0.5, τότε α = 10, β = 1, γ = 2, δ = 7, Ο1 = 11, Ο2 = 11, P1 = 12, P2 = 8, και:

• Το ποσοστό ορθής ταξινόμησης είναι (PAC) = (10 + 7) / 20 = 0,85 = 85%.
• Η ευαισθησία του (sensitivity) είναι 7 / 9 = 0,778 = 77,8%.
• Η προσδιοριστικότητα του (specificity) είναι 10 / 11 = 0,909 = 90,9%.
• Η θετική διαγνωστική αξία (PPV) είναι 7 / 8 = 0,875 = 87,5%.
• Η αρνητική διαγνωστική αξία (NPV) είναι 10 / 12 = 0,833 =83,3%.

Από τα παραπάνω ποσοστά γίνεται αντιληπτό πως το λογιστικό μοντέλο που προέκυψε κάνει αρκετά καλή πρόβλεψη της αδυναμίας πληρωμής ενός πελάτη της τράπεζας.

Στο σημείο αυτό αναδεικνύεται το ερώτημα: Ποιο είναι το ιδανικό cut-off που μας εξασφαλίζει ένα καλό συνδυασμό ευαισθησίας και ειδικότητας;

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

Εμβόλιμη Δραστηριότητα

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

α) Μοντέλο πρόβλεψης της πιθανότητας ένας δανειολήπτης να δηλώσει αδυναμία αποπληρωμής.

β) Για μία σοβαρή ασθένεια αναπτύχθηκε ένα ισχυρό φάρμακο με πολλές παρενέργειες. Η πιθανότητα για έναν άνθρωπο να πάσχει από αυτήν την ασθένεια υπολογίζεται από ένα λογιστικό μοντέλο πρόβλεψης που παίρνει υπόψη του πολλές βιομετρικές παραμέτρους.

7.1 Καμπύλη ROC

Το σύνολο των ζευγών (ευαισθησία, 1 – ειδικότητα) είναι δυνατόν να τοποθετηθούν σε ένα διάγραμμα που θα διευκολύνει την απόφαση του ερευνητή σχετικά με το όριο (cut - off) από το οποίο και πάνω θα δέχεται ως σωστή την πρόβλεψη του λογιστικού μοντέλου. Το διάγραμμα αυτό είναι γνωστό ως Καμπύλη λειτουργικού χαρακτηριστικού δέκτη - receiver operating characteristic curve ή πιο απλά καμπύλη ROC.

Η ονομασία “Καμπύλη ROC” χρονολογείται από τον Δεύτερο Παγκόσμιο Πόλεμο και τους ηλεκτρολόγους μηχανικούς που υποστήριζαν τη λειτουργία των ραντάρ.

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

Η καμπύλη ROC επινοήθηκε ως ένα απλό στη δημιουργία του διάγραμμα που θα βοηθούσε το χειριστή να αποφασίσει ποιο επίπεδο ενίσχυσης θα προσέφερε ικανοποιητική ευαισθησία (True Positive Rate, TPR) και σχετικά μεγάλη ειδικότητα ή ισοδύναμα μικρή πιθανότητα εσφαλμένης ανίχνευσης (False Positive Rate FPR = 1 – ειδικότητα).

Πηγή: Streiner, D. L., & Cairney, J. (2007). What’s under the ROC? An introduction to receiver operating characteristics curves. Canadian journal of psychiatry. Revue canadienne de psychiatrie, 52(2), 121–128. https://doi.org/10.1177/070674370705200210

Μία παρουσίαση της καμπύλης ROC είναι διαθέσιμη εδώ.

Η συνάρτηση my.logistic.regression, δημιουργεί το διάγραμμα ROC και το αποδίδει ως δεύτερο μέρος στη λίστα με τα output. Είναι προσβάσιμο ως:

my.model$ROCplot

Από την παρατήρηση του συνόλου των ζευγών (Sensitivity, Specificity) και τα αντίστοιχα όρια πιθανότητας μπορεί να επιλεγεί το όριο που αντιστοιχεί στο συνδυασμό που ικανοποιεί περισσότερο τις ανάγκες του ερευνητή. Αυτά τα ζεύγη μπορούν να παρουσιαστούν με τον κώδικα:

round(my.model$ROCplot$data, 3)
##    threshold specificity sensitivity 1-specificity
## 21       Inf       1.000       0.000         0.000
## 20     0.986       1.000       0.111         0.000
## 19     0.962       1.000       0.222         0.000
## 18     0.942       1.000       0.333         0.000
## 17     0.927       1.000       0.444         0.000
## 16     0.914       1.000       0.556         0.000
## 15     0.845       1.000       0.667         0.000
## 14     0.764       1.000       0.778         0.000
## 13     0.620       0.909       0.778         0.091
## 12     0.469       0.818       0.778         0.182
## 11     0.349       0.818       0.889         0.182
## 10     0.232       0.727       0.889         0.273
## 9      0.188       0.636       0.889         0.364
## 8      0.124       0.636       1.000         0.364
## 7      0.076       0.545       1.000         0.455
## 6      0.054       0.455       1.000         0.545
## 5      0.033       0.364       1.000         0.636
## 4      0.015       0.273       1.000         0.727
## 3      0.001       0.182       1.000         0.818
## 2      0.000       0.091       1.000         0.909
## 1       -Inf       0.000       1.000         1.000

Μία αυτοματοποιημένη πρόταση βάσει του δείκτη Youden’s J μπορεί να προκύψει με τον παρακάτω κώδικα:

library(OptimalCutpoints)
lr.data.frame$default2 = as.numeric(lr.data.frame$default) - 1
lr.data.frame$predicted = fitted(my.model$fit)
summary(optimal.cutpoints(X = "predicted", status = "default2", tag.healthy = 0, methods = "Youden", data = lr.data.frame, control = control.cutpoints(), ci.fit = TRUE))
## 
## Call:
## optimal.cutpoints.default(X = "predicted", status = "default2", 
##     tag.healthy = 0, methods = "Youden", data = lr.data.frame, 
##     control = control.cutpoints(), ci.fit = TRUE)
## 
## Area under the ROC curve (AUC):  0.455 (0.181, 0.728) 
## 
## CRITERION: Youden
## Number of optimal cutoffs: 1
## 
##                     Estimate 95% CI lower limit 95% CI upper limit
## cutoff            0.06743724                  -                  -
## Se                0.88888889          0.5175035          0.9971909
## Sp                0.36363636          0.1092634          0.6920953
## PPV               0.53333333          0.1970018          0.9806620
## NPV               0.80000000          0.3490757          0.9402425
## DLR.Positive      1.39682540          0.8447628          2.3096675
## DLR.Negative      0.30555556          0.0410873          2.2723370
## FP                7.00000000                  -                  -
## FN                1.00000000                  -                  -
## Optimal criterion 0.25252525                  -                  -

7.2 Διαγνωστικά εργαλεία λογιστικης παλινδρόμησης

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

7.2.1 Αναζήτηση ιδιαζόντων σημείων

Για την ανίχνευση ιδιαζόντων σημείων, παρατηρούμε την απόσταση Cook των διαφορετικών παρατηρήσεων με το εξής διάγραμμα:

plot(my.model$fit, which = 4, id.n = 3)

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

Τα τυποποιημένα υπολοιπα όπως και η απόσταση Cook, μπορούν να ενσωματωθούν στα δεδομένα:

lr.data.frame <- broom::augment(my.model$fit) %>% 
  mutate(index = 1:n())

Στη συνέχεια μπορούμε να παρατηρήσουμε τις 3 περισσότερο ιδιάζουσες παρατηρήσεις με τις εντολές:

lr.data.frame %>% top_n(3, .cooksd)
## # A tibble: 3 × 10
##   dependent_final income  debtinc .fitted .resid  .hat .sigma .cooksd .std.resid
##   <fct>           <label> <label>   <dbl>  <dbl> <dbl>  <dbl>   <dbl>      <dbl>
## 1 Όχι             25      10        1.10   -1.67 0.162  0.738   0.231      -1.82
## 2 Ναι             38       4       -1.64    1.91 0.140  0.699   0.325       2.06
## 3 Ναι             16       2       -0.208   1.27 0.370  0.770   0.383       1.60
## # ℹ 1 more variable: index <int>
lr.data.frame %>% top_n(3, .std.resid)
## # A tibble: 3 × 10
##   dependent_final income  debtinc .fitted .resid  .hat .sigma .cooksd .std.resid
##   <fct>           <label> <label>   <dbl>  <dbl> <dbl>  <dbl>   <dbl>      <dbl>
## 1 Ναι             38       4       -1.64   1.91  0.140  0.699  0.325       2.06 
## 2 Ναι             16       2       -0.208  1.27  0.370  0.770  0.383       1.60 
## 3 Ναι             41      16        1.26   0.708 0.209  0.844  0.0317      0.796
## # ℹ 1 more variable: index <int>

7.2.2 Συγγραμικότητα

Όπως και στην περίπτωση της γραμμικής παλινδρόμησης, η ύπαρξη επεξηγηματικών μεταβλητών με υψηλή συσχέτιση μειώνει την αξιοπιστία του υπολογισμού των αντίστοιχων συντελεστών. Όπως και στη γραμμική περίπτωση, η ανίχνευση του φαινομένου μπορεί να γίνει με τον υπολογισμό του συντελεστή Variance Inflation Factor (VIF):

round(car::vif(my.model$fit), 1)
##  income debtinc 
##     1.2     1.2

Τιμές VIF μεγαλύτερες από 5 θεωρούνται ενοχλητικές και τιμές πάνω από 10 απαγορευτικές για τη συμπερίληψη όλων των αντίστοιχων μεταβλητών στο μοντέλο πρόβλεψης.

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

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

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

parents = create_factors(parents)
  1. Εκτελέστε την εντολή summary(parents) για να εξοικειωθείτε με τα περιεχόμενα του αρχείου
summary(parents)
  1. Μεταβάλετε τις ετικέτες του παράγοντα 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. Υπολογίστε τη συσχέτιση μεταξύ των ψυχομετρικών χαρακτηριστικών των γονέων (δεδομένα parents), με την εντολή:
my_cor_table(parents, c("hdhq.total", "scl_90_gen_deikt_sympt", "oas_total_score", "ess_total", "stai_present", "stai_general"))
  1. Υπολογίστε τη μερική συσχέτιση μεταξύ του άγχους (stai_general) και της ψυχοπαθολογίας (scl_90_gen_deikt_sympt) ελέγχοντας ως προς την επιρροή της επιθετικότητας (hdhq.total), δίνοντας την εντολή:
ppcor::pcor.test(parents$stai_general, parents$scl_90_gen_deikt_sympt, parents$hdhq.total)

Συγκρίνετε το αποτέλεσμα με τη συσχέτιση μεταξύ του άγχους (stai_general) και της ψυχοπαθολογίας (scl_90_gen_deikt_sympt). Πόσο επηρεάζεται η σχέση του άγχους και της ψυχοπαθολογίας από την επιθετικότητα;

  1. Υπολογίστε γραμμικό μοντέλο πρόβλεψης της ψυχοπαθολογίας (scl_90_gen_deikt_sympt) από το άγχος (stai_general), δίνοντας την εντολή:
my.model = my_lm(parents, scl_90_gen_deikt_sympt ~ stai_general)
my.model$htmlreport

Ελέγξτε την ποιότητα του μοντέλου παρατηρώντας τα διαγράμματα που παράγονται από την plot(fit)

layout(matrix(c(1,2,3,4),2,2))
plot(my.model$fit)
  1. Υπολογίστε γραμμικό μοντέλο πρόβλεψης της ψυχοπαθολογίας (scl_90_gen_deikt_sympt) από το άγχος (stai_general) και την ντροπή (ess_total), δίνοντας την εντολή:
my.model = my_lm(parents, scl_90_gen_deikt_sympt ~ stai_general + ess_total)
my.model$htmlreport
  1. Ελέγξτε την ποιότητα προσαρμογής και την ύπαρξη ιδιαζόντων σημείων με την εντολή:
layout(matrix(c(1,2,3,4),2,2))
plot(my.model$fit)
  1. Υπολογίστε γραμμικό μοντέλο πρόβλεψης της ψυχοπαθολογίας (scl_90_gen_deikt_sympt) από το άγχος (stai_general), την ντροπή (ess_total), το φύλο (gender) και την οικογενειακή κατάσταση (marital) δίνοντας την εντολή:
my.model = my_lm(parents, scl_90_gen_deikt_sympt ~ stai_general + ess_total + gender + marital)
my.model$htmlreport
  1. Ελέγξτε την ποιότητα προσαρμογής και την ύπαρξη ιδιαζόντων σημείων με την εντολή:
layout(matrix(c(1,2,3,4),2,2))
plot(my.model$fit)

Θέλοντας να προβλέψουμε την πιθανότητα επίσκεψης σε ψυχολόγο (schico) από την επιθετικότητα (hdhq.total) και την ψυχοπαθολογία (scl_90_gen_deikt_sympt), εφαρμόσαμε λογιστική παλινδρόμηση.

  1. Εκτελέστε τις εντολές:
my.model1 = my_logistic_regression(parents, schico ~ hdhq.total + scl_90_gen_deikt_sympt)
my.model1$htmlreport
parents$predicted1 = fitted(my.model1$fit)
  1. Παρατηρήστε όλους τους συνδυασμούς (Sensitivity, Specificity):
my.model1$ROCplot$data
  1. Προσαρμόστε και ένα λογιστικό μοντέλο πρόβλεψης της επίσκεψης σε ψυχολόγο από τις συνιστώσες της ντροπής:
my.model2 = my_logistic_regression(parents, schico ~ ess_total + oas_total_score)
my.model2$htmlreport
parents$predicted2 = fitted(my.model2$fit)
  1. Δημιουργήστε και την καμπύλη ROC για το δεύτερο μοντέλο:
my.model2$ROCplot
  1. Συγκρίνετε την επίδοση των δύο μοντέλων με τον παρακάτω κώδικα:
roc1 <- roc(parents$schico, parents$predicted1)
roc2 <- roc(parents$schico, parents$predicted2)
roc.test(roc1, roc2)
  1. Επαναλάβατε τη διαδικασία με το μοντέλο πρόβλεψης:
my.model3 = my_logistic_regression(parents, schico ~ hdhq.total + scl_90_gen_deikt_sympt + child_age + gender + child_gender)
my.model3$htmlreport
parents$predicted3 = fitted(my.model3$fit)
roc3 <- roc(parents$schico, parents$predicted3)
  1. Συγκρίνετε την επίδοση των τριών μοντέλων με τον παρακάτω κώδικα:
roc.test(roc1, roc2)
roc.test(roc1, roc3)
roc.test(roc2, roc3)

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

Στα πλαίσια διπλωματικής με αντικείμενο τον τρόπο επεξεργασίας μαθηματικών προβλημάτων από μαθητές Λυκείου, καταγράφηκαν για 200 μαθητές, ορισμένα δημογραφικά χαρακτηριστικά τους όπως και ο τρόπος με τον οποίον αντιμετώπισαν 12 μαθηματικά προβλήματα. Τα αποτελέσματα αποθηκεύτηκαν σε ένα αρχείο SPSS που είναι διαθέσιμα για λήψη εδώ και θα αξιοποιηθούν στη συνέχεια σε δραστηριότητες εξοικείωσης.
1) Να κάνετε λήψη τα δεδομένα.

  1. Να εισάγετε τα δεδομένα στο R Studio ως data.frame με όνομα students

  2. Μετατρέψτε όσες στήλες πρέπει από αλφαριθμητικές (character) σε παράγοντες (factors) με τις εντολές:

cols <- c("taxi", "tmima", "diglosia", "math.diskolies")
students[cols] <- lapply(students[cols], factor)

  1. Εκτελέστε την εντολή summary(students) για να εξοικειωθείτε με τα περιεχόμενα του αρχείου

  2. Προσθέστε ετικέτες στις μεταβλητές taxi, diglosia, math.diskolies, με τις εντολές

library(plyr)
students$taxi = revalue(students$taxi, c("1"="A", "2"="B", "3"="Γ"))
students$diglosia = revalue(students$diglosia, c("1"="Όχι", "2"="Ναι"))
students$math.diskolies = revalue(students$math.diskolies, c("0"="Όχι", "1"="Ναι"))
Hmisc::label(students$taxi) = 'Τάξη'
Hmisc::label(students$tmima) = 'Τμήμα'
Hmisc::label(students$diglosia) = 'Διγλωσία'
Hmisc::label(students$math.diskolies) = 'Μαθησιακές Δυσκολίες'
  1. Υπολογίστε γραμμικό μοντέλο πρόβλεψης του πλήθους των σωστών αποκρίσεων (sostes.apantiseis) από τη διγλωσία (diglosia) και τις μαθησιακές δυσκολίες (math.diskolies), δίνοντας την εντολή:
my.model = my_lm(students, sostes.apantiseis ~ diglosia + math.diskolies)
my.model$htmlreport

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

  1. Εκτελέστε τις εντολές:
my.model = my_logistic_regression(diglosia ~ sostes.apantiseis, data = students)
my.model$htmlreport
  1. Εμφανίστε την καμπύλη ROC με τις εντολές:
my.model$ROCplot
  1. Εκτυπώστε όλα τους συνδυασμούς (Sensitivity, Specificity) με τον κώδικα:
my.model$ROCplot$data

8.3 Εφαρμογή στο data frame patients

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

patients = create_factors(patients)
  1. Εκτελέστε την εντολή summary(patients) για να εξοικειωθείτε με τα περιεχόμενα του αρχείου
summary(patients)
  1. Εισάγετε κατάλληλες ετικέτες σε όσες μεταβλητές χρειάζεται, με τις εντολές
levels(patients$gender) = c('Γυναίκα', 'Ανδρας')
levels(patients$marital) = c('Ελεύθερος/η', 'Παντρεμένος/η', 'Διαζευγμένος/η', 'Χήρος/α)')
levels(patients$livingplace) = c('Χωριό', 'Μικρή Πόλη', 'Μεγάλη Πόλη')
levels(patients$educational) = c('Αναλφάβητος/Πρωτοβάθμια', 'Δευτεροβάθμια/Τριτοβάθμια')
levels(patients$pension) = c('Όχι', 'Ναι')
levels(patients$pensionisgood) = c('Όχι', 'Ναι')
levels(patients$worktodayformoney) = c('Όχι', 'Ναι')
levels(patients$siblings) = c('Όχι', 'Ναι')
levels(patients$sister) = c('Όχι', 'Ναι')
levels(patients$brother) = c('Όχι', 'Ναι')
levels(patients$children) = c('Όχι', 'Ναι')
levels(patients$daughter) = c('Όχι', 'Ναι')
levels(patients$son) = c('Όχι', 'Ναι')
levels(patients$daughterlivenearyou) = c('Όχι', 'Ναι')
levels(patients$sonlivesnearyou) = c('Όχι', 'Ναι')
levels(patients$holidays) = c('Όχι', 'Ναι')
levels(patients$walking) = c('Όχι', 'Ναι')
levels(patients$garden) = c('Όχι', 'Ναι')
levels(patients$difficultgiratia) = c('Όχι', 'Ναι')
levels(patients$friendsinthepast) = c('Όχι', 'Ναι')
levels(patients$friendstoday) = c('Όχι', 'Ναι')
  1. Εισάγετε ετικέτες στις μεταβλητές που θα χρησιμοποιήσουμε στη συνέχεια.
library(Hmisc)
label(patients$age) = 'Ηλικία'
label(patients$gender) = 'Φύλο'
label(patients$marital) = 'Οικογενειακή κατάσταση'
label(patients$livingplace) = 'Τόπος Κατοικίας'
label(patients$children) = 'Παιδιά'
label(patients$siblings) = 'Συγγενείς'
label(patients$sister) = 'Αδελφή'
label(patients$brother) = 'Αδελφός'
label(patients$daughter) = 'Κόρη'
label(patients$son) = 'Γιος'
label(patients$daughterlivenearyou) = 'Η κόρη ζει κοντά'
label(patients$sonlivesnearyou) = 'Ο γιος ζει κοντά'
label(patients$holidays) = 'Διακοπές'
label(patients$walking) = 'Περπάτημα'
label(patients$garden) = 'Κηπουρική'
label(patients$difficultgiratia) = 'Δύσκολα γηρατειά'
label(patients$friendsinthepast) = 'Φίλοι στο παρελθόν'
label(patients$friendstoday) = 'Φίλοι σήμερα'
label(patients$DEP) = 'Depression'
label(patients$ANX) = 'Anxiety'
label(patients$OAS) = 'External Shame'
label(patients$ESS) = 'Internal Shame'
label(patients$SELF_COMPASSION) = 'Self Compassion'
label(patients$CDRISC) = 'Resilience'
label(patients$GOHAI.PF) = 'Physical Functioning'
label(patients$GOHAI.PS) = 'Psyhosocial Functioning'
label(patients$GOHAI.PD) = 'Pain / Discomfort'
  1. Εφαρμόζουμε λογιστική παλινδρόμηση για να βρούμε κατά πόσο επηρεάζει το ψυχοπαθολογικό προφίλ του ηλικιωμένου την πιθανότητα, αυτός να έχει φίλους σήμερα. Εκτελέστε τον παρακάτω κώδικα και σχολιάστε τα αποτελέσματα:
my.model = my_logistic_regression(friendstoday ~ age + gender + DEP + ANX + OAS + ESS + CDRISC, data = patients)
my.model$htmlreport
  1. Επαναλάβατε τη δοκιμασία για τις μεταβλητές holidays, walking και σχολιάστε τα αποτελέσματα.

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