ΕΙΣΑΓΩΓΗ

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

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

ΕΞΕΡΕΥΝΗΣΗ ΤΟΥ DATASET

data(infert)      # Φόρτωση του dataset (περιλαμβάνεται στο R)
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, 10)  # Πρώτες 10 εγγραφές
##    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
## 6    6-11yrs  36      4       2    1           1       6             36
## 7    6-11yrs  23      1       0    1           0       7              6
## 8    6-11yrs  32      2       0    1           0       8             22
## 9    6-11yrs  21      1       0    1           1       9              5
## 10   6-11yrs  28      2       0    1           0      10             19

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

Το dataset περιλαμβάνει 248 παρατηρήσεις και τις παρακάτω βασικές μεταβλητές:

  • case: Η εξαρτημένη μεταβλητή
    • 0 = χωρίς πρόβλημα στειρότητας
    • 1 = με πρόβλημα στειρότητας
  • age: Ηλικία της γυναίκας
  • education: Επίπεδο εκπαίδευσης (π.χ. λυκείου, πανεπιστημίου)
  • parity: Πλήθος γεννήσεων
  • induced: Τεχνητές αποβολές (π.χ. προγραμματισμένες)
  • spontaneous: Αυτόματες αποβολές (μη προγραμματισμένες)
  • stratum: Στρωματοποίηση (για στατιστικό έλεγχο σε ομάδες)
  • pooled.stratum: Συνδυασμένη στρωματοποίηση

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

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

# Επιλογή αριθμητικών μεταβλητών
numeric_vars <- infert[, sapply(infert, is.numeric)]

# Υπολογισμός πίνακα συσχέτισης
cor_matrix <- cor(numeric_vars)

# Μετατροπή σε long format για ggplot
cor_melt <- melt(cor_matrix)

# Heatmap με αριθμούς
ggplot(cor_melt, aes(x = Var1, y = Var2, fill = value)) +
  geom_tile(color = "white") +
  geom_text(aes(label = round(value, 2)), size = 3.5) +  # Τιμές συσχέτισης με 2 δεκαδικά
  scale_fill_gradient2(low = "blue", high = "red", mid = "white",
                       midpoint = 0, limit = c(-1,1), space = "Lab",
                       name = "Συσχέτιση") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  labs(title = "Heatmap Συσχέτισης Αριθμητικών Μεταβλητών",
       x = "", y = "")

ΔΙΑΓΡΑΜΜΑΤΑ

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

ggplot(infert, aes(x = age)) +
  geom_histogram(fill = "steelblue", color = "white", bins = 15) +
  labs(
    title = "Κατανομή Ηλικιών στο Δείγμα",
    x = "Ηλικία",
    y = "Συχνότητα"
  ) +
  theme_minimal()

ggplot(infert, aes(x = age, y = parity)) +
  geom_jitter(width = 0.3, height = 0.1, alpha = 0.6, color = "darkgreen") +
  labs(
    title = "Σχέση Ηλικίας και Πλήθους Γεννήσεων",
    x = "Ηλικία",
    y = "Πλήθος γεννήσεων (parity)"
  ) +
  theme_minimal()

ggplot(infert, aes(x = factor(case))) +
  geom_bar(fill = "steelblue") +
  labs(x = "Περίπτωση στειρότητας", y = "Συχνότητα",
       title = "Κατανομή εξαρτημένης μεταβλητής (case)") +
  scale_x_discrete(labels = c("0" = "Χωρίς πρόβλημα", "1" = "Με πρόβλημα"))

ggplot(infert, aes(x = education, y = age)) +
  geom_boxplot(fill = "skyblue", color = "black") +
  labs(
    title = "Ηλικίες ανά επίπεδο εκπαίδευσης",
    x = "Επίπεδο εκπαίδευσης",
    y = "Ηλικία"
  ) +
  theme_minimal()

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

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

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

ΠΛΗΘΟΣ ΔΕΙΓΜΑΤΟΣ TRAIN VS TEST

set.seed(906)  # Αντικατέστησέ το με το δικό σου όταν έχεις το seed
split <- sample.split(infert$case, SplitRatio = 0.65)

train <- subset(infert, split == TRUE)
test <- subset(infert, split == FALSE)

cat("Πλήθος στο training set:", nrow(train), "\n")
## Πλήθος στο training set: 161
cat("Πλήθος στο test set:", nrow(test), "\n")
## Πλήθος στο 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 ✅ Σημαντική τεχνική μείωση

ΤΙ ΔΕΙΧΝΕΙ Η PREDICT

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

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

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

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

ΚΑΜΠΥΛΗ ROC ΜΕ ΠΛΗΡΕΣ ΜΟΝΤΕΛΟ

# Προβλέψεις με το πλήρες μοντέλο
predictTest_full <- predict(model, type = "response", newdata = test)

# Δημιουργία ROCR αντικειμένου
ROCRpred_full <- prediction(predictTest_full, 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 Curve για το Πλήρες Μοντέλο")

ΚΑΜΠΥΛΗ ROC ΜΕ ΠΛΗΡΕΣ & ΕΝΔΙΑΜΕΣΟ ΜΟΝΤΕΛΟ

model_sig <- glm(case ~ induced + spontaneous + pooled.stratum, data = train, family = binomial)
summary(model_sig)
## 
## Call:
## glm(formula = case ~ induced + spontaneous + pooled.stratum, 
##     family = binomial, data = train)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -1.21361    0.42783  -2.837  0.00456 ** 
## induced         0.57740    0.25787   2.239  0.02515 *  
## spontaneous     1.24221    0.27298   4.551 5.35e-06 ***
## pooled.stratum -0.01889    0.01166  -1.620  0.10534    
## ---
## 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: 181.35  on 157  degrees of freedom
## AIC: 189.35
## 
## Number of Fisher Scoring iterations: 4
predictTest <- predict(model_sig, type = "response", newdata = test)
head(predictTest)
##         1         2         6        15        18        20 
## 0.8571171 0.3418419 0.6232785 0.5521554 0.4741213 0.4273988
table(test$case, predictTest > 0.5)
##    
##     FALSE TRUE
##   0    54    4
##   1    19   10
# Προβλέψεις για τα δύο μοντέλα
predictTest_full <- predict(model, type = "response", newdata = test)
predictTest_sig  <- predict(model_sig, type = "response", newdata = test)

# ROCR objects
ROCRpred_full <- prediction(predictTest_full, test$case)
ROCRpred_sig  <- prediction(predictTest_sig, test$case)

# ROC performance objects
ROCRperf_full <- performance(ROCRpred_full, "tpr", "fpr")
ROCRperf_sig  <- performance(ROCRpred_sig,  "tpr", "fpr")

# Plot both 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 ΜΕ ΠΛΗΡΕΣ, ΕΝΔΙΑΜΕΣΟ & ΑΔΥΝΑΜΟ ΜΟΝΤΕΛΟ

model_weak <- glm(case ~ education + parity + stratum, data = train, family = binomial)
summary(model_weak)
## 
## Call:
## glm(formula = case ~ education + parity + stratum, family = binomial, 
##     data = train)
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)
## (Intercept)      -0.90915    1.02118  -0.890    0.373
## education6-11yrs  0.36403    0.95937   0.379    0.704
## education12+ yrs  0.82618    1.26303   0.654    0.513
## parity            0.06753    0.14007   0.482    0.630
## stratum          -0.01153    0.01490  -0.773    0.439
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 205.42  on 160  degrees of freedom
## Residual deviance: 204.36  on 156  degrees of freedom
## AIC: 214.36
## 
## Number of Fisher Scoring iterations: 4
# Προβλέψεις
predictTest_weak <- predict(model_weak, type = "response", newdata = test)

# ROCR objects
ROCRpred_weak <- prediction(predictTest_weak, test$case)
ROCRperf_weak <- performance(ROCRpred_weak, "tpr", "fpr")

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)

ΥΠΟΛΟΓΙΣΜΟΣ AUC

# AUC τιμές για όλα τα μοντέλα
auc_full <- performance(ROCRpred_full, measure = "auc")@y.values[[1]]
auc_sig  <- performance(ROCRpred_sig,  measure = "auc")@y.values[[1]]
auc_weak <- performance(ROCRpred_weak, measure = "auc")@y.values[[1]]

# Δημιουργία data frame με τις τιμές
auc_table <- data.frame(
  Μοντέλο = c("Πλήρες", "Απλοποιημένο", "Αδύναμο"),
  AUC = round(c(auc_full, auc_sig, auc_weak), 3)
)

print(auc_table)
##        Μοντέλο   AUC
## 1       Πλήρες 0.798
## 2 Απλοποιημένο 0.740
## 3      Αδύναμο 0.375

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

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

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