Λογιστική Παλινδρόμιση στο framingham dataset
Δουλγερίδης Αναστάσης iis22022
2025-05-05
Περιγραφή του Dataset
Το Framingham Heart Study είναι μια μακροχρόνια επιδημιολογική μελέτη που επικεντρώνεται στην πρόληψη καρδιαγγειακών παθήσεων. Το dataset περιλαμβάνει δεδομένα από συμμετέχοντες, με κάθε γραμμή να αντιπροσωπεύει ένα άτομο και τις αντίστοιχες μεταβλητές που σχετίζονται με τον καρδιαγγειακό κίνδυνο. Ο στόχος της ανάλυσης είναι η πρόβλεψη της πιθανότητας εμφάνισης καρδιοπάθειας εντός 10 ετών (TenYearCHD).
Το dataset περιέχει 15 μεταβλητές, οι κυριότερες από τις οποίες είναι:
Ηλικία (age): Η ηλικία του συμμετέχοντα σε έτη. Φύλο (male): Δυαδική μεταβλητή (0 = γυναίκα, 1 = άνδρας). Συνολική χοληστερόλη (totChol): Επίπεδα χοληστερόλης στο αίμα (mg/dL). Συστολική αρτηριακή πίεση (sysBP): Μετρήσεις σε mmHg. Κάπνισμα (cigsPerDay): Αριθμός τσιγάρων που καπνίζονται ημερησίως. Δείκτης Μάζας Σώματος (BMI): Υπολογίζεται ως βάρος (kg) / ύψος^2 (m). Επίπεδα γλυκόζης (glucose): Επίπεδα γλυκόζης στο αίμα (mg/dL). Υπέρταση (prevalentHyp): Δυαδική μεταβλητή (0 = όχι, 1 = ναι). Καρδιοπάθεια εντός 10 ετών (TenYearCHD): Η μεταβλητή απόκρισης (0 = όχι, 1 = ναι).
Επεξεργασία και Διαχωρισμός Δεδομένων
Εισαγωγή Βιβλιοθηκών και Δεδομένων Πρώτα, εισάγουμε τις απαραίτητες βιβλιοθήκες για την ανάλυση και φορτώνουμε το dataset.
library(caTools)## Warning: package 'caTools' was built under R version 4.4.3library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.4.3library(dplyr)## Warning: package 'dplyr' was built under R version 4.4.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(knitr)## Warning: package 'knitr' was built under R version 4.4.3library(ROCR)## Warning: package 'ROCR' was built under R version 4.4.3library(kableExtra)## Warning: package 'kableExtra' was built under R version 4.4.3##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rowsΦόρτωση δεδομένων
fram <- read.csv("C:/Users/Tasos/Downloads/framingham.csv")
str(fram)## 'data.frame': 4240 obs. of 16 variables:
## $ male : int 1 0 1 0 0 0 0 0 1 1 ...
## $ age : int 39 46 48 61 46 43 63 45 52 43 ...
## $ education : int 4 2 1 3 3 2 1 2 1 1 ...
## $ currentSmoker : int 0 0 1 1 1 0 0 1 0 1 ...
## $ cigsPerDay : int 0 0 20 30 23 0 0 20 0 30 ...
## $ BPMeds : int 0 0 0 0 0 0 0 0 0 0 ...
## $ prevalentStroke: int 0 0 0 0 0 0 0 0 0 0 ...
## $ prevalentHyp : int 0 0 0 1 0 1 0 0 1 1 ...
## $ diabetes : int 0 0 0 0 0 0 0 0 0 0 ...
## $ totChol : int 195 250 245 225 285 228 205 313 260 225 ...
## $ sysBP : num 106 121 128 150 130 ...
## $ diaBP : num 70 81 80 95 84 110 71 71 89 107 ...
## $ BMI : num 27 28.7 25.3 28.6 23.1 ...
## $ heartRate : int 80 95 75 65 85 77 60 79 76 93 ...
## $ glucose : int 77 76 70 103 85 99 85 78 79 88 ...
## $ TenYearCHD : int 0 0 0 1 0 0 1 0 0 0 ...summary(fram)## male age education currentSmoker
## Min. :0.0000 Min. :32.00 Min. :1.000 Min. :0.0000
## 1st Qu.:0.0000 1st Qu.:42.00 1st Qu.:1.000 1st Qu.:0.0000
## Median :0.0000 Median :49.00 Median :2.000 Median :0.0000
## Mean :0.4292 Mean :49.58 Mean :1.979 Mean :0.4941
## 3rd Qu.:1.0000 3rd Qu.:56.00 3rd Qu.:3.000 3rd Qu.:1.0000
## Max. :1.0000 Max. :70.00 Max. :4.000 Max. :1.0000
## NA's :105
## cigsPerDay BPMeds prevalentStroke prevalentHyp
## Min. : 0.000 Min. :0.00000 Min. :0.000000 Min. :0.0000
## 1st Qu.: 0.000 1st Qu.:0.00000 1st Qu.:0.000000 1st Qu.:0.0000
## Median : 0.000 Median :0.00000 Median :0.000000 Median :0.0000
## Mean : 9.006 Mean :0.02962 Mean :0.005896 Mean :0.3106
## 3rd Qu.:20.000 3rd Qu.:0.00000 3rd Qu.:0.000000 3rd Qu.:1.0000
## Max. :70.000 Max. :1.00000 Max. :1.000000 Max. :1.0000
## NA's :29 NA's :53
## diabetes totChol sysBP diaBP
## Min. :0.00000 Min. :107.0 Min. : 83.5 Min. : 48.0
## 1st Qu.:0.00000 1st Qu.:206.0 1st Qu.:117.0 1st Qu.: 75.0
## Median :0.00000 Median :234.0 Median :128.0 Median : 82.0
## Mean :0.02571 Mean :236.7 Mean :132.4 Mean : 82.9
## 3rd Qu.:0.00000 3rd Qu.:263.0 3rd Qu.:144.0 3rd Qu.: 90.0
## Max. :1.00000 Max. :696.0 Max. :295.0 Max. :142.5
## NA's :50
## BMI heartRate glucose TenYearCHD
## Min. :15.54 Min. : 44.00 Min. : 40.00 Min. :0.0000
## 1st Qu.:23.07 1st Qu.: 68.00 1st Qu.: 71.00 1st Qu.:0.0000
## Median :25.40 Median : 75.00 Median : 78.00 Median :0.0000
## Mean :25.80 Mean : 75.88 Mean : 81.96 Mean :0.1519
## 3rd Qu.:28.04 3rd Qu.: 83.00 3rd Qu.: 87.00 3rd Qu.:0.0000
## Max. :56.80 Max. :143.00 Max. :394.00 Max. :1.0000
## NA's :19 NA's :1 NA's :388Υπολογισμός ελλιπών τιμών
Καθαρισμός Δεδομένων Ελέγχουμε για ελλιπή δεδομένα και τα αφαιρούμε για να εξασφαλίσουμε την ποιότητα του dataset.
cat("Συνολικός αριθμός ελλιπών τιμών: ", sum(is.na(fram)), "\n")## Συνολικός αριθμός ελλιπών τιμών: 645data <- na.omit(fram)
cat("Διαστάσεις dataset μετά τον καθαρισμό: ", dim(data), "\n")## Διαστάσεις dataset μετά τον καθαρισμό: 3658 16Διαχωρισμός σε Train και Test Set Χωρίζουμε το dataset σε εκπαιδευτικό σύνολο (train, 65%) και δοκιμαστικό σύνολο (test, 35%) για την εκπαίδευση και αξιολόγηση του μοντέλου.
set.seed(922)
split <- sample.split(data$TenYearCHD, SplitRatio = 0.65)
train <- subset(data, split == TRUE)
test <- subset(data, split == FALSE)
cat("Μέγεθος εκπαιδευτικού συνόλου: ", nrow(train), "\n")## Μέγεθος εκπαιδευτικού συνόλου: 2378 cat("Μέγεθος δοκιμαστικού συνόλου: ", nrow(test), "\n")## Μέγεθος δοκιμαστικού συνόλου: 1280Δημιουργία Μοντέλων Λογιστικής Παλινδρόμισης
Αρχικό Μοντέλο Δημιουργούμε ένα αρχικό μοντέλο λογιστικής παλινδρόμισης που περιλαμβάνει όλες τις ανεξάρτητες μεταβλητές.
model_full <- glm(TenYearCHD ~ ., data = train, family = "binomial")
summary(model_full)##
## Call:
## glm(formula = TenYearCHD ~ ., family = "binomial", data = train)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.431688 0.897972 -9.390 < 2e-16 ***
## male 0.658538 0.137250 4.798 1.60e-06 ***
## age 0.063789 0.008443 7.555 4.17e-14 ***
## education -0.089975 0.062549 -1.438 0.15030
## currentSmoker 0.132674 0.193670 0.685 0.49331
## cigsPerDay 0.014864 0.007836 1.897 0.05783 .
## BPMeds 0.442929 0.306589 1.445 0.14854
## prevalentStroke 0.070240 0.585284 0.120 0.90448
## prevalentHyp 0.261106 0.169814 1.538 0.12415
## diabetes -0.374707 0.404630 -0.926 0.35442
## totChol 0.003439 0.001379 2.493 0.01266 *
## sysBP 0.013446 0.004773 2.817 0.00485 **
## diaBP -0.004478 0.008143 -0.550 0.58236
## BMI 0.008531 0.015936 0.535 0.59244
## heartRate -0.003907 0.005295 -0.738 0.46055
## glucose 0.008699 0.002869 3.032 0.00243 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2028.7 on 2377 degrees of freedom
## Residual deviance: 1782.6 on 2362 degrees of freedom
## AIC: 1814.6
##
## Number of Fisher Scoring iterations: 5Από την περίληψη του μοντέλου, παρατηρούμε ότι οι μεταβλητές male, age, sysBP, cigsPerDay, prevalentHyp, totChol και glucose είναι στατιστικά σημαντικές (p-value < 0.05).
Μοντέλο 1: Μείωση Μεταβλητών Δημιουργούμε ένα δεύτερο μοντέλο, διατηρώντας μόνο τις σημαντικές μεταβλητές από το αρχικό μοντέλο.
model1 <- glm(TenYearCHD ~ male + age + sysBP + cigsPerDay + prevalentHyp + totChol + glucose,
data = train, family = "binomial")
summary(model1)##
## Call:
## glm(formula = TenYearCHD ~ male + age + sysBP + cigsPerDay +
## prevalentHyp + totChol + glucose, family = "binomial", data = train)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.856858 0.651565 -13.593 < 2e-16 ***
## male 0.644170 0.133998 4.807 1.53e-06 ***
## age 0.065967 0.008067 8.178 2.89e-16 ***
## sysBP 0.012982 0.003482 3.728 0.000193 ***
## cigsPerDay 0.018643 0.005293 3.522 0.000428 ***
## prevalentHyp 0.252438 0.165988 1.521 0.128304
## totChol 0.003197 0.001370 2.333 0.019642 *
## glucose 0.007002 0.002122 3.300 0.000967 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2028.7 on 2377 degrees of freedom
## Residual deviance: 1789.4 on 2370 degrees of freedom
## AIC: 1805.4
##
## Number of Fisher Scoring iterations: 5Μοντέλο 2: Περαιτέρω Μείωση Αφαιρούμε επιπλέον μεταβλητές με χαμηλότερη στατιστική σημασία (p-value > 0.01), διατηρώντας μόνο τις male, age και sysBP.
model2 <- glm(TenYearCHD ~ male + age + sysBP, data = train, family = "binomial")
summary(model2)##
## Call:
## glm(formula = TenYearCHD ~ male + age + sysBP, family = "binomial",
## data = train)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.810395 0.479090 -16.303 < 2e-16 ***
## male 0.749820 0.122492 6.121 9.28e-10 ***
## age 0.062970 0.007637 8.245 < 2e-16 ***
## sysBP 0.018053 0.002616 6.902 5.12e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2028.7 on 2377 degrees of freedom
## Residual deviance: 1819.3 on 2374 degrees of freedom
## AIC: 1827.3
##
## Number of Fisher Scoring iterations: 5Μοντέλο 3: Ελάχιστο Μοντέλο Δημιουργούμε ένα τελικό μοντέλο με τις δύο πιο σημαντικές μεταβλητές (male, age), καθώς το μοντέλο αυτό εμφανίζει υψηλό AIC, υποδεικνύοντας καλύτερη ισορροπία μεταξύ προσαρμογής και απλότητας.
model3 <- glm(TenYearCHD ~ male + age, data = train, family = "binomial")
summary(model3)##
## Call:
## glm(formula = TenYearCHD ~ male + age, family = "binomial", data = train)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.133167 0.395942 -15.490 < 2e-16 ***
## male 0.644000 0.119022 5.411 6.28e-08 ***
## age 0.079543 0.007164 11.104 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2028.7 on 2377 degrees of freedom
## Residual deviance: 1867.2 on 2375 degrees of freedom
## AIC: 1873.2
##
## Number of Fisher Scoring iterations: 5Επιλογή Μοντέλου: Επιλέγουμε το μοντέλο 2, καθώς περιλαμβάνει τις μεταβλητές male, age και sysBP, οι οποίες είναι στατιστικά σημαντικές, και προσφέρει καλύτερη προγνωστική ικανότητα σε σχέση με το μοντέλο 3, ενώ παραμένει πιο απλό από το μοντέλο 1.
Προβλέψεις και Confusion Matrix
Προβλέψεις στο Test και Train Set Χρησιμοποιούμε το επιλεγμένο μοντέλο (model2) για να κάνουμε προβλέψεις στα εκπαιδευτικά και δοκιμαστικά σύνολα.
Προβλέψεις στο test set
predictTest <- predict(model2, newdata = test, type = "response")
cat("Μέση πιθανότητα (Test): ", mean(predictTest), "\n")## Μέση πιθανότητα (Test): 0.1483496 cat("Εύρος πιθανοτήτων (Test): ", range(predictTest), "\n")## Εύρος πιθανοτήτων (Test): 0.02091697 0.7353176Προβλέψεις στο train set
predictTrain <- predict(model2, newdata = train, type = "response")Κατηγοριοποίηση προβλέψεων
pred_class <- ifelse(predictTest >= 0.5, 1, 0)
actual_class <- test$TenYearCHD
conf_mat <- table(Predicted = pred_class, Actual = actual_class)Δημιουργία του επεξηγηματικού πίνακα
Παρακάτω παρουσιάζεται ο πίνακας σύγχυσης με τις πραγματικές τιμές που προκύπτουν από τις προβλέψεις του μοντέλου:
Actual = 0, Predicted = 0: Περιπτώσεις που το μοντέλο προέβλεψε ότι
δεν έχουν την πάθηση και όντως δεν την έχουν (True
Negatives)./n Actual = 0, Predicted = 1: Περιπτώσεις που το
μοντέλο προέβλεψε ότι έχουν την πάθηση, αλλά δεν την έχουν (False
Positives)./n Actual = 1, Predicted = 0: Περιπτώσεις που το
μοντέλο προέβλεψε ότι δεν έχουν την πάθηση, αλλά την έχουν (False
Negatives)./n Actual = 1, Predicted = 1: Περιπτώσεις που το
μοντέλο προέβλεψε ότι έχουν την πάθηση και όντως την έχουν (True
Positives)./n
conf_matrix_values <- matrix(c(conf_mat[1,1], conf_mat[1,2],
conf_mat[2,1], conf_mat[2,2]),
nrow = 2, byrow = TRUE)
colnames(conf_matrix_values) <- c("Predicted = 0", "Predicted = 1")
rownames(conf_matrix_values) <- c("Actual = 0", "Actual = 1")Μορφοποίηση με kableExtra
kable(conf_matrix_values, caption = "Πίνακας Σύγχυσης") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) %>%
add_header_above(c(" " = 1, "Προβλεπόμενη Κατηγορία" = 2)) %>%
column_spec(1, bold = TRUE, background = "#f2dede") %>%
row_spec(0, background = "#a94442", color = "white", bold = TRUE) %>%
column_spec(2:3, background = "#f5f5f5")| Predicted = 0 | Predicted = 1 | |
|---|---|---|
| Actual = 0 | 1081 | 185 |
| Actual = 1 | 4 | 10 |
ROC Καμπύλη και AUC
ROC Καμπύλη Η καμπύλη ROC (Receiver Operating Characteristic) απεικονίζει την ικανότητα του μοντέλου να διακρίνει μεταξύ των κλάσεων.
ROCRpred <- prediction(predictTest, test$TenYearCHD)
ROCRperf <- performance(ROCRpred, "tpr", "fpr")
plot(ROCRperf, colorize = TRUE, print.cutoffs.at = seq(0, 1, 0.1),
text.adj = c(-0.2, 1.7), main = "ROC Καμπύλη")
AUC (Area Under the Curve)
Η τιμή AUC μετρά τη συνολική απόδοση του μοντέλου. Τιμές κοντά στο 1 υποδηλώνουν εξαιρετική διακριτική ικανότητα.
auc_value <- as.numeric(performance(ROCRpred, "auc")@y.values)
cat("AUC: ", auc_value, "\n")## AUC: 0.7169703Συμπεράσματα
Η ανάλυση έδειξε ότι το μοντέλο λογιστικής παλινδρόμισης με τις μεταβλητές male, age και sysBP (model2) παρέχει μια ισορροπημένη προσέγγιση για την πρόβλεψη του κινδύνου καρδιοπάθειας εντός 10 ετών. Ο πίνακας σύγχυσης παρουσιάζει τις πραγματικές τιμές των προβλέψεων, βοηθώντας στην κατανόηση της απόδοσης του μοντέλου. Η ROC καμπύλη και η τιμή AUC υποδεικνύουν ικανοποιητική διακριτική ικανότητα. Το μοντέλο μπορεί να βελτιωθεί περαιτέρω με την προσθήκη αλληλεπιδράσεων μεταβλητών ή τη χρήση πιο πολύπλοκων αλγορίθμων μηχανικής μάθησης.