Στο dataset infert.

1. ΔΙΕΥΡΕΥΝΗΣΗ ΤΟΥ DATASET

ΕΙΣΑΓΩΓΗ

Η παρούσα εργασία εξετάζει ένα case study Λογιστικής Παλινδρόμησης με χρήση του έτοιμου dataset infert, το οποίο περιέχει ιατρικά και δημογραφικά δεδομένα από γυναίκες που συμμετείχαν σε μελέτη για τη γονιμότητα.

Ο στόχος είναι να μελετήσουμε ποιες μεταβλητές επηρεάζουν την πιθανότητα στειρότητας, και να δημιουργήσουμε ένα μοντέλο που προβλέπει αν μια γυναίκα έχει ιστορικό στειρότητας (μεταβλητή case).

Το dataset περιέχει 248 εγγραφές και απαρτίζεται από 9 διαφορετικές στήλες/χαρακτηριστικά.

ΑΡΧΙΚΕΣ ΔΙΑΓΝΩΣΤΙΚΕΣ ΕΝΤΟΛΕΣ

str(infert)       # Δομή δεδομένων
## 'data.frame':    248 obs. of  8 variables:
##  $ education     : Factor w/ 3 levels "0-5yrs","6-11yrs",..: 1 1 1 1 2 2 2 2 2 2 ...
##  $ age           : num  26 42 39 34 35 36 23 32 21 28 ...
##  $ parity        : num  6 1 6 4 3 4 1 2 1 2 ...
##  $ induced       : num  1 1 2 2 1 2 0 0 0 0 ...
##  $ case          : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ spontaneous   : num  2 0 0 0 1 1 0 0 1 0 ...
##  $ stratum       : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ pooled.stratum: num  3 1 4 2 32 36 6 22 5 19 ...
summary(infert)   # Περιγραφικά στατιστικά
##    education        age            parity         induced      
##  0-5yrs : 12   Min.   :21.00   Min.   :1.000   Min.   :0.0000  
##  6-11yrs:120   1st Qu.:28.00   1st Qu.:1.000   1st Qu.:0.0000  
##  12+ yrs:116   Median :31.00   Median :2.000   Median :0.0000  
##                Mean   :31.50   Mean   :2.093   Mean   :0.5726  
##                3rd Qu.:35.25   3rd Qu.:3.000   3rd Qu.:1.0000  
##                Max.   :44.00   Max.   :6.000   Max.   :2.0000  
##       case         spontaneous        stratum      pooled.stratum 
##  Min.   :0.0000   Min.   :0.0000   Min.   : 1.00   Min.   : 1.00  
##  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:21.00   1st Qu.:19.00  
##  Median :0.0000   Median :0.0000   Median :42.00   Median :36.00  
##  Mean   :0.3347   Mean   :0.5766   Mean   :41.87   Mean   :33.58  
##  3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:62.25   3rd Qu.:48.25  
##  Max.   :1.0000   Max.   :2.0000   Max.   :83.00   Max.   :63.00
head(infert, 5)  # Πρώτες 5 εγγραφές
##   education age parity induced case spontaneous stratum pooled.stratum
## 1    0-5yrs  26      6       1    1           2       1              3
## 2    0-5yrs  42      1       1    1           0       2              1
## 3    0-5yrs  39      6       2    1           0       3              4
## 4    0-5yrs  34      4       2    1           0       4              2
## 5   6-11yrs  35      3       1    1           1       5             32

ΠΕΡΙΓΡΑΦΗ ΜΕΤΑΒΛΗΤΩΝ

Μεταβλητή Περιγραφή Μονάδα Μέτρησης
case Η Εξαρτημένη μεταβλητή: αν υπάρχει πρόβλημα στειρότητας Δυαδική (0/1)
age Ηλικία της γυναίκας έτη
education Επίπεδο εκπαίδευσης (1 έως 4) Κατηγορική (αριθμός)
parity Πλήθος προηγούμενων γεννήσεων αριθμός
induced Πλήθος τεχνητών αποβολών αριθμός
spontaneous Πλήθος αυτόματων αποβολών αριθμός
stratum Ομάδα κατάταξης (για ανάλυση matched cases) Κατηγορική (αριθμός)
pooled.stratum Συνδυασμένη πληροφορία κατανομής Κατηγορική (αριθμός)

ΕΠΕΞΗΓΗΜΑΤΙΚΑ ΣΧΟΛΙΑ

Σημείωση1: Η μεταβλητή education είναι κατηγορική και παίρνει τιμές 1 έως 4. Αυτές δεν είναι αριθμητικά κλιμακωμένες, αλλά αναπαριστούν διακριτές ομάδες εκπαίδευσης.

Τιμή Περιγραφή εκπαίδευσης
1 0–5 χρόνια
2 6–11 χρόνια
3 12+ χρόνια
4 Δε γνωρίζουμε / missing value

Σημείωση2: Αυτές είναι ειδικές μεταβλητές που δεν περιέχουν πρωτογενή πληροφορία, αλλά είναι τεχνητά δημιουργημένες μεταβλητές για το design του πειράματος: stratum και pooled.stratum (θα τις αγνοήσω παρόλο που έχουν υψηλές συσχετίσεις)

  • stratum: Δηλώνει το matching group: κάθε ομάδα περιλαμβάνει μία περίπτωση (case) και δύο ελέγχους (controls) που “ταιριάζουν” ως προς συγκεκριμένα χαρακτηριστικά.

  • pooled.stratum: Είναι αριθμός που ταξινομεί τα matched sets από όλες τις περιπτώσεις.

Σημείωση3: Και η πιο σημαντική! Η age στο dataset αναφέρεται στην ηλικία της γυναίκας τη στιγμή της μελέτης, δηλαδή όταν συμμετείχε στην έρευνα (και όχι απαραίτητα όταν γέννησε). Δεν είναι η ηλικία κατά τον τοκετό. Είναι η ηλικία τη στιγμή που καταγράφηκαν τα δεδομένα για προηγούμενες κυήσεις (γεννήσεις, αποβολές, κλπ). Δηλαδή μια γυναίκα 33 ετών μπορεί να έχει κάνει 3 γέννες στα 24, 26 και 29, αλλά η ηλικία της καταγράφεται ως 33, και το parity είναι 3.

ΠΙΝΑΚΑΣ ΣΥΣΧΕΤΙΣΕΩΝ

Ένα correlation matrix σε βοηθάει ώστε με μια γρήγορη ματιά να ξέρεις ποιες είναι οι δυνατές συσχετίσεις που υπάρχουν στο dataset.

ΔΙΑΓΡΑΜΜΑΤΑ

Παραθέτω κάποια διαγράμματα μεταξύ κάποιων μεταβλητών-ζευγαριών για καλύτερη κατανόηση του dataset

graph1: geom_histogram

graph2: geom_jitter

graph3: geom_bar

graph4: geom_jitter

graph5: geom_col

ΚΑΤΑΝΟΜΗ ΠΕΡΙΠΤΩΣΕΩΝ

table(infert$case)
## 
##   0   1 
## 165  83

Η παραπάνω εντολή μας δείχνει την κατανομή των περιπτώσεων στειρότητας στο δείγμα. Συγκεκριμένα, παρατηρούμε ότι 165 γυναίκες δεν παρουσίασαν πρόβλημα στειρότητας (case = 0), ενώ 83 γυναίκες παρουσίασαν (case = 1). Αυτή η αναλογία είναι σημαντική, καθώς μας δίνει μια πρώτη εικόνα για την ισορροπία των κατηγοριών στην εξαρτημένη μεταβλητή, κάτι που επηρεάζει τόσο την εκπαίδευση του μοντέλου όσο και την αξιολόγηση των αποτελεσμάτων.

2. ΔΗΜΙΟΥΡΓΙΑ ΜΟΝΤΕΛΟΥ ΠΑΛΙΝΔΡΟΜΗΣΗΣ

SPLIT ΤΟ DATASET

set.seed(906)  
split <- sample.split(infert$case, SplitRatio = 0.65)

train <- subset(infert, split == TRUE)
test <- subset(infert, split == FALSE)
## Πλήθος στο training set: 161
## Πλήθος στο test set: 87

ΠΛΗΡΕΣ ΜΟΝΤΕΛΟ

model <- glm(case ~ ., data = train, family = binomial)
summary(model)
## 
## Call:
## glm(formula = case ~ ., family = binomial, data = train)
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -5.05008    2.72136  -1.856 0.063494 .  
## education6-11yrs  3.11528    2.09867   1.484 0.137701    
## education12+ yrs  7.07402    3.82722   1.848 0.064552 .  
## age               0.07911    0.04753   1.664 0.096081 .  
## parity           -0.18540    0.33060  -0.561 0.574931    
## induced           1.48317    0.40778   3.637 0.000276 ***
## spontaneous       2.13860    0.41799   5.116 3.11e-07 ***
## stratum          -0.01749    0.01767  -0.990 0.322133    
## pooled.stratum   -0.12260    0.05596  -2.191 0.028465 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 205.42  on 160  degrees of freedom
## Residual deviance: 168.10  on 152  degrees of freedom
## AIC: 186.1
## 
## Number of Fisher Scoring iterations: 5

Σημαντικο σχόλιο

Το . είναι shorthand που σημαίνει: “Βάλε όλες τις ανεξάρτητες μεταβλητές του dataset, εκτός από την εξαρτημένη (case)”

Ανάλυση Συντελεστών του Μοντέλου

Από την έξοδο του summary(model) παρατηρούμε ότι οι μεταβλητές induced και spontaneous έχουν στατιστικά σημαντική επίδραση (p < 0.001) στην πιθανότητα εμφάνισης προβλήματος στειρότητας (case = 1). Και οι δύο έχουν θετικούς συντελεστές, γεγονός που σημαίνει ότι όσο αυξάνονται οι τεχνητές ή αυτόματες αποβολές, τόσο αυξάνεται και η πιθανότητα στειρότητας.

ΕΡΜΗΝΕΙΑ ΣΥΝΤΕΛΕΣΤΩΝ

Ο παρακάτω πίνακας συνοψίζει τους συντελεστές του μοντέλου λογιστικής παλινδρόμησης, την στατιστική τους σημαντικότητα και μια σύντομη ερμηνεία της επίδρασής τους στην πιθανότητα εμφάνισης στειρότητας (case = 1).

Μεταβλητή Estimate p-value Ερμηνεία
(Intercept) –5.05 0.063 .
age +0.08 0.096 . αυξάνεται ελάχιστα η πιθανότητα όσο αυξάνεται και η ηλικία
education6–11 yrs +3.13 0.138 Όχι σημαντική αύξηση
education12+ yrs +7.07 0.065 . Οριακά σημαντική αύξηση
parity –0.23 0.575 Μη σημαντική μείωση
induced +1.48 0.000 *** Σημαντική αύξηση
spontaneous +2.14 <0.001 *** Πολύ σημαντική αύξηση
stratum –0.01 0.322 Μη σημαντική τεχνική
pooled.stratum –0.12 0.029 * Σημαντική τεχνική μείωση

3. ΠΡΟΒΛΕΨΗ

ΕΦΑΡΜΟΓΗ ΠΡΟΒΛΕΨΗΣ

predictTest <- predict(model, type = "response", newdata = test)
head(predictTest)
##         1         2         6        15        18        20 
## 0.7806339 0.3572198 0.6809057 0.7173836 0.7113758 0.7436504
summary(predictTest)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.01123 0.12788 0.34212 0.34950 0.55436 0.81663

Η συνάρτηση predict() με το όρισμα type = "response" επιστρέφει, για κάθε παρατήρηση στο test set, την εκτιμώμενη πιθανότητα να ανήκει στην κατηγορία case = 1, δηλαδή να έχει πρόβλημα στειρότητας. Οι τιμές αυτές κυμαίνονται μεταξύ 0 και 1 και αντιπροσωπεύουν το πόσο “σίγουρο” είναι το μοντέλο για την πρόβλεψή του.

Για παράδειγμα, μια τιμή predictTest = 0.78 σημαίνει ότι το μοντέλο εκτιμά με 78% πιθανότητα ότι η συγκεκριμένη γυναίκα ανήκει στην κατηγορία με ιστορικό στειρότητας (case = 1). Αντίστοιχα, μια τιμή predictTest = 0.35 υποδηλώνει ότι η πιθανότητα είναι μόνο 35%, άρα μάλλον δεν έχει πρόβλημα.

THRESHOLD VALUE

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

table(test$case, predictTest > 0.5)
##    
##     FALSE TRUE
##   0    46   12
##   1    13   16

ΥΠΟΛΟΓΙΣΜΟΣ ΜΕΤΡΙΚΩΝ

Φτιάχνω τον πίνακα σύγχησης:

predictedClass <- ifelse(predictTest > 0.5, 1, 0)

conf_matrix <- table(Actual = test$case, Predicted = predictedClass)
conf_matrix
##       Predicted
## Actual  0  1
##      0 46 12
##      1 13 16
TN <- conf_matrix[1, 1]
FP <- conf_matrix[1, 2]
FN <- conf_matrix[2, 1]
TP <- conf_matrix[2, 2]

Υπολογίζω τις μετρικές:

# N = συνολικά παραδείγματα
N <- sum(conf_matrix)

# Υπολογισμοί:
accuracy  <- (TP + TN) / N
sensitivity <- TP / (TP + FN)
specificity <- TN / (TN + FP)
error_rate <- (FP + FN) / N
false_negative_rate <- FN / (TP + FN)
false_positive_rate <- FP / (TN + FP)

# Προβολή
results <- data.frame(
  Metric = c("Accuracy", "Sensitivity", "Specificity", 
             "Error Rate", "False Negative Rate", "False Positive Rate"),
  Value = round(c(accuracy, sensitivity, specificity, 
                  error_rate, false_negative_rate, false_positive_rate), 4)
)
knitr::kable(results, caption = "Αποτίμηση απόδοσης μοντέλου")
Αποτίμηση απόδοσης μοντέλου
Metric Value
Accuracy 0.7126
Sensitivity 0.5517
Specificity 0.7931
Error Rate 0.2874
False Negative Rate 0.4483
False Positive Rate 0.2069

AREA UNDER THE CURCE (AUC)

Καμπύλη roc 1 με το πλήρες μοντέλο

# Δημιουργία ROCR αντικειμένου
ROCRpred_full <- prediction(predictTest, test$case)

# ROC καμπύλη
ROCRperf_full <- performance(ROCRpred_full, "tpr", "fpr")

plot(ROCRperf_full, colorize = TRUE, print.cutoffs.at = seq(0.1, 0.9, by = 0.1),
     text.adj = c(-0.2, 1.2),
     main = "Σύγκριση ROC Καμπυλών Μοντέλων")

Καμπύλη roc 2 με πλήρες μοντέλο vs το ενδιάμεσο

# φτιάχνουμε το μοντέλο απλά με τις πιο δυνατές μεταβλητές και όχι με όλες
model_sig <- glm(case ~ induced + spontaneous + pooled.stratum, data = train, family = binomial)

# κάνουμε την πρόβλεψη με το ενδιάμεσο μοντέλο
predictTest_sig <- predict(model_sig, type = "response", newdata = test)

# βλέπουμε τα αποτελέσματα
table(test$case, predictTest_sig > 0.5)
##    
##     FALSE TRUE
##   0    54    4
##   1    19   10
# φτιάχνουμε την καμπύλη roc για το ενδιάμεσο
ROCRpred_sig  <- prediction(predictTest_sig, test$case)

# ROC performance object
ROCRperf_sig  <- performance(ROCRpred_sig,  "tpr", "fpr")
# κάνουμε Plot και τις 2 ROC curves που έχουμε μέχρι τώρα
plot(ROCRperf_full, col = "blue", lwd = 2, main = "Σύγκριση ROC Καμπυλών Μοντέλων")
plot(ROCRperf_sig,  col = "red", lwd = 2, add = TRUE)
abline(a = 0, b = 1, lty = 2, col = "gray")  # γραμμή τύχης

legend("bottomright",
       legend = c("Πλήρες Μοντέλο", "Ενδιάμεσο Μοντέλο"),
       col = c("blue", "red"), lwd = 2)

Καμπύλη roc 3 με πλήρες και ενδιάμεσο μοντέλο vs το αδύναμο

# φτιάχνουμε το μοντέλο μόνο με τις πιο αδύναμες μεταβλητές
model_weak <- glm(case ~ education + parity + stratum, data = train, family = binomial)

# κάνουμε την πρόβλεψη με το αδύναμο μοντέλο
predictTest_weak <- predict(model_weak, type = "response", newdata = test)

# βλέπουμε τα αποτελέσματα
table(test$case, predictTest_weak > 0.5)
##    
##     FALSE
##   0    58
##   1    29
# φτιάχνουμε την καμπύλη roc για το ενδιάμεσο
ROCRpred_weak <- prediction(predictTest_weak, test$case)

# ROC performance object
ROCRperf_weak <- performance(ROCRpred_weak, "tpr", "fpr")

# κάνουμε Plot και τις 2 ROC curves που έχουμε μέχρι τώρα
plot(ROCRperf_full, col = "blue", lwd = 2, main = "Σύγκριση ROC Καμπυλών Μοντέλων")
plot(ROCRperf_sig, col = "red", lwd = 2, add = TRUE)
plot(ROCRperf_weak, col = "darkgreen", lwd = 2, add = TRUE)
abline(a = 0, b = 1, lty = 2, col = "gray")

legend("bottomright",
       legend = c("Πλήρες Μοντέλο", "Ενδιάμεσο Μοντέλο", "Αδύναμο Μοντέλο"),
       col = c("blue", "red", "darkgreen"),
       lwd = 2)

4. ΑΠΟΤΕΛΕΣΜΑΤΑ - ΠΑΡΑΔΕΙΓΜΑΤΑ

ΣΥΓΚΡΙΣΗ ΜΕΤΡΙΚΩΝ

Παρακάτω κάνουμε σύγκριση των 3 μοντέλων λογιστικής παλινδρόμησης που φτιάξαμε. Το AIC που βλέπετε είναι ο δείκτης Akaike Information Criterion τον οποίο βρίσκουμε άνα κάνουμε summary τα μοντέλα μας κάτω κάτω. Είναι ένα μέτρο που αξιολογεί πόσο καλά ταιριάζει το μοντέλο στα δεδομένα, αλλά επιβάλει ποινή για την πολυπλοκότητα. Από το semi προς το full βελτιώνεται λίγο αλλά από το semi στο weak χάνει πολύ.
Model AIC AUC
full 186.100 0.798
semi 189.351 0.740
weak 214.358 0.375

ΣΥΜΠΕΡΑΣΜΑΤΑ

  • Στην παρούσα εργασία αναλύσαμε το dataset infert (το οποίο είναι από τα datasets που υπάρχουν προεγκατεστημένα στην R οπως το mtcars κλπ κλπ) με στόχο τη δημιουργία μοντέλου λογιστικής παλινδρόμησης για την πρόβλεψη πιθανότητας στειρότητας.

  • Μέσα από την ανάλυση, εντοπίστηκαν ως στατιστικά σημαντικοί παράγοντες οι μεταβλητές induced και spontaneous, που αφορούν τεχνητές και αυτόματες αποβολές αντίστοιχα.

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

  • Είδαμε πόσο πιο σημαντικό είναι να έχουμε τις πιο επιδραστικές μεταβλητές για την δημιουργία του μοντέλου παρά να τις έχουμε όλες. Δείτε πόσο λίγο βελτιώνεται από το semi προς το full αλλά πόσο πολύ χάνει από το semi στο weak.

  • Το default κατώφλι 0.5 μπορεί να χάνει πολλά πραγματικά θετικά (TP) σε datasets με ανισορροπία (class imbalance).

ΠΡΟΒΛΗΜΑΤΙΣΜΟΙ

  • Νόμιζα πώς το Sensitivity και Specificity είναι συμπληρωματικά (δηλ. 1 - sensitivity = specificity) αλλά αντίθετα, μετρούν διαφορετικά πράγματα. Το sensitivity κοιτάει τη γραμμή των πραγματικών θετικών ενώ το specificity κοιτάει τη γραμμή των πραγματικών αρνητικών.

  • Όταν αρχικά είδα το AIC στο summary() φαντάστηκα ότι πρόκειται για το AUC αλλά όχι. Το AIC επιβραβεύει καλή προσαρμογή και απλότητα, ενώ το AUC μετράει την ικανότητα διάκρισης θετικών και αρνητικών. Είδαμε ότι μοντέλο με χαμηλό AIC δεν σημαίνει απαραίτητα ότι έχει και υψηλό AUC, και η αξιολόγηση πρέπει να λαμβάνει υπόψη πολλούς δείκτες ταυτόχρονα.

  • Δεν μπορούσα να καταλάβω πως γίνεται να υπάρχει κανονική κατανομή στων αριθμό γεννήσεων ανάλογα την ηλικία όταν θεωρητικά όσο πιο μεγάλη ηλικιακά η γυναίκα, τόσες περισσότερες γέννες θα είχε προλάβει να κάνει. Αλλά η κατανομή έχει κορυφή γύρω στα 30 επειδή οι περισσότερες γυναίκες στο δείγμα είναι ηλικίας 30–34, και άρα η πλειοψηφία των καταγεγραμμένων γεννήσεων συνδέεται με αυτές τις ηλικίες όχι γιατί γέννησαν τότε, αλλά γιατί αυτές έχουν ήδη τεκνοποιήσει μέχρι την ηλικία τους. Παρόλα αυτά φαίνεται μια γενική τάση: Μέχρι ~25 οι περισσότερες έχουν λίγα παιδιά, Στα ~30 έχουν ήδη αποκτήσει παιδιά, Στα 40+ δεν αυξάνεται πολύ (ή δεν υπάρχουν πολλές γυναίκες στο δείγμα). Αυτό δεν λέει πότε γέννησαν, αλλά πόσα παιδιά είχαν ήδη όταν έγιναν π.χ. 30.

  • Επίσης ακόμα και τώρα δεν έχω καταλάβει γιατί το μοντέλο επέλεξε να σπάσει την μεταβλητή education στα 2 από τα 3 σκαλοπάτια τιμών του. Άμα δείτε σε κάποιο summary() βγάζει education6-11yrs και education12+ yrs και δεν ξέρω γιατί λείπει το πρώτο. Επίσης δεν έκανα πολλά summary() για να κερδίσουμε σε ευκολία ανάγνωσης γιατί πιάνει πολύ χώρο.