Previsione peso neonatale
Cristian Costache
2025-01-03
data <- read.csv("neonati.csv", stringsAsFactors = FALSE)
head(data)## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto
## 1 26 0 0 42 3380 490 325 Nat
## 2 21 2 0 39 3150 490 345 Nat
## 3 34 3 0 38 3640 500 375 Nat
## 4 28 1 0 41 3690 515 365 Nat
## 5 20 0 0 38 3700 480 335 Nat
## 6 32 0 0 40 3200 495 340 Nat
## Ospedale Sesso
## 1 osp3 M
## 2 osp1 F
## 3 osp2 M
## 4 osp2 M
## 5 osp3 F
## 6 osp2 Fstr(data)## 'data.frame': 2500 obs. of 10 variables:
## $ Anni.madre : int 26 21 34 28 20 32 26 25 22 23 ...
## $ N.gravidanze: int 0 2 3 1 0 0 1 0 1 0 ...
## $ Fumatrici : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Gestazione : int 42 39 38 41 38 40 39 40 40 41 ...
## $ Peso : int 3380 3150 3640 3690 3700 3200 3100 3580 3670 3700 ...
## $ Lunghezza : int 490 490 500 515 480 495 480 510 500 510 ...
## $ Cranio : int 325 345 375 365 335 340 345 349 335 362 ...
## $ Tipo.parto : chr "Nat" "Nat" "Nat" "Nat" ...
## $ Ospedale : chr "osp3" "osp1" "osp2" "osp2" ...
## $ Sesso : chr "M" "F" "M" "M" ...summary(data)## Anni.madre N.gravidanze Fumatrici Gestazione
## Min. : 0.00 Min. : 0.0000 Min. :0.0000 Min. :25.00
## 1st Qu.:25.00 1st Qu.: 0.0000 1st Qu.:0.0000 1st Qu.:38.00
## Median :28.00 Median : 1.0000 Median :0.0000 Median :39.00
## Mean :28.16 Mean : 0.9812 Mean :0.0416 Mean :38.98
## 3rd Qu.:32.00 3rd Qu.: 1.0000 3rd Qu.:0.0000 3rd Qu.:40.00
## Max. :46.00 Max. :12.0000 Max. :1.0000 Max. :43.00
## Peso Lunghezza Cranio Tipo.parto
## Min. : 830 Min. :310.0 Min. :235 Length:2500
## 1st Qu.:2990 1st Qu.:480.0 1st Qu.:330 Class :character
## Median :3300 Median :500.0 Median :340 Mode :character
## Mean :3284 Mean :494.7 Mean :340
## 3rd Qu.:3620 3rd Qu.:510.0 3rd Qu.:350
## Max. :4930 Max. :565.0 Max. :390
## Ospedale Sesso
## Length:2500 Length:2500
## Class :character Class :character
## Mode :character Mode :character
##
##
## colSums(is.na(data))## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza
## 0 0 0 0 0 0
## Cranio Tipo.parto Ospedale Sesso
## 0 0 0 0sum(duplicated(data))## [1] 0library(corrplot)## Warning: il pacchetto 'corrplot' è stato creato con R versione 4.4.2## corrplot 0.95 loadedboxplot(data$Anni.madre, main = "Boxplot Età della Madre")
hist(data$Peso, main = "Istogramma Peso Neonati", xlab = "Peso (g)", breaks = 30)
corrplot(cor(data[, sapply(data, is.numeric)]), method = "color")
cor(data[, sapply(data, is.numeric)])## Anni.madre N.gravidanze Fumatrici Gestazione Peso
## Anni.madre 1.000000000 0.38063184 0.006057852 -0.13566714 -0.02247017
## N.gravidanze 0.380631844 1.00000000 0.046869466 -0.10149194 0.00240730
## Fumatrici 0.006057852 0.04686947 1.000000000 0.03220907 -0.01894534
## Gestazione -0.135667143 -0.10149194 0.032209066 1.00000000 0.59176872
## Peso -0.022470167 0.00240730 -0.018945344 0.59176872 1.00000000
## Lunghezza -0.063491571 -0.06040371 -0.020781411 0.61892045 0.79603676
## Cranio 0.016076703 0.03900707 -0.008665549 0.46082894 0.70480151
## Lunghezza Cranio
## Anni.madre -0.06349157 0.016076703
## N.gravidanze -0.06040371 0.039007074
## Fumatrici -0.02078141 -0.008665549
## Gestazione 0.61892045 0.460828941
## Peso 0.79603676 0.704801507
## Lunghezza 1.00000000 0.603340974
## Cranio 0.60334097 1.000000000table(data$Fumatrici)##
## 0 1
## 2396 104table(data$Tipo.parto)##
## Ces Nat
## 728 1772table(data$Ospedale)##
## osp1 osp2 osp3
## 816 849 835table(data$Sesso)##
## F M
## 1256 1244data <- data[data$Anni.madre > 0, ]
data$Fumatrici <- factor(data$Fumatrici)
data$Tipo.parto <- factor(data$Tipo.parto)
data$Ospedale <- factor(data$Ospedale)
data$Sesso <- factor(data$Sesso)
data$Prematuro <- ifelse(data$Gestazione < 37, 1, 0)Dopo aver importato il set ed aver effettuato l’EDA, ho applicato alcune modifiche per rendere il dataset più attendibile. Ho rimosso i dati delle madri con età = 0 e ho convertito le variabili categoriche in fattori, ho anche creato una nuova colonna che ci indichi i parti prematuri (quelli avvenuti prima della 37esima settimana)
boxplot(data$Peso ~ data$Fumatrici,
main = "Peso Neonati vs Fumo Materno",
xlab = "Fumatrici", ylab = "Peso (g)")
plot(data$Gestazione, data$Peso,
main = "Peso Neonati vs Settimane di Gestazione",
xlab = "Settimane di Gestazione", ylab = "Peso (g)", pch = 19)
Dai grafici possiamo dedurre che: Il peso mediano dei neonati è leggermente inferiore nel gruppo delle madri fumatrici; Nel gruppo delle madri fumatrici ci sono più outlier sul lato inferiore, ad indicare la presenza di più neonati con un peso più basso; La variabilità del peso è simile tra i due gruppi, ma la distribuzione per le madri fumatrici sembra più compressa verso i valori inferiori.
Per le settimane di gestazione si nota una relazione positiva e il peso del neonato. Più lunga è la durata della gestazione e maggiore è il peso mediano del neonato; I neonati prematuri tendono ad avere un peso significativamente inferirore rispetto ai neonati a termine; C’è una variabilità marcata nei pesi dei neonati anche a durate di gestazione simili, ciò potrebbe indicare che altre variabili possono influenzare il peso dei neonati.
cont_table_ospedale_parto <- table(data$Ospedale, data$Tipo.parto)
test_chi_ospedale_parto <- chisq.test(cont_table_ospedale_parto)
test_chi_ospedale_parto##
## Pearson's Chi-squared test
##
## data: cont_table_ospedale_parto
## X-squared = 1.0604, df = 2, p-value = 0.5885Poiché il p-value è maggiore di 0.05, non rifiutiamo l’ipotesi nulla. Questo significa che non ci sono prove sufficienti per affermare che ci siano differenze significative nel numero di parti cesarei tra i diversi ospedali.
In base ai dati a disposizione, possiamo concludere che la distribuzione dei parti cesarei sembra indipendente dall’ospedale. Non ci sono evidenze che in alcuni ospedali si facciano significativamente più parti cesarei rispetto ad altri.
media_peso_campione <- mean(data$Peso)
media_lunghezza_campione <- mean(data$Lunghezza)
media_peso_popolazione <- 3300
media_lunghezza_popolazione <- 50
test_peso <- t.test(data$Peso, mu = media_peso_popolazione)
test_lunghezza <- t.test(data$Lunghezza, mu = media_lunghezza_popolazione)
test_peso##
## One Sample t-test
##
## data: data$Peso
## t = -1.5069, df = 2498, p-value = 0.132
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
## 3263.572 3304.769
## sample estimates:
## mean of x
## 3284.17test_lunghezza##
## One Sample t-test
##
## data: data$Lunghezza
## t = 844.49, df = 2498, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 50
## 95 percent confidence interval:
## 493.6613 495.7265
## sample estimates:
## mean of x
## 494.6939Per il peso: Non c’è una differenza significativa rispetto alla popolazione. Per la lunghezza: La differenza rispetto alla popolazione è altamente significativa.
t.test(data$Peso ~ data$Sesso)##
## Welch Two Sample t-test
##
## data: data$Peso by data$Sesso
## t = -12.116, df = 2490, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -287.3966 -207.3302
## sample estimates:
## mean in group F mean in group M
## 3161.132 3408.496t.test(data$Lunghezza ~ data$Sesso)##
## Welch Two Sample t-test
##
## data: data$Lunghezza by data$Sesso
## t = -9.5864, df = 2459, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -11.937899 -7.883398
## sample estimates:
## mean in group F mean in group M
## 489.7643 499.6750t.test(data$Cranio ~ data$Sesso)##
## Welch Two Sample t-test
##
## data: data$Cranio by data$Sesso
## t = -7.4237, df = 2490.6, p-value = 1.555e-13
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -6.100260 -3.550953
## sample estimates:
## mean in group F mean in group M
## 337.6330 342.4586In tutti e tre i test (Peso, Lunghezza, e Cranio), la differenza tra i sessi è statisticamente significativa, con le femmine che tendono ad avere valori medi inferiori per ciascuna delle misure rispetto ai maschi. Pertanto, possiamo rifiutare l’ipotesi nulla di assenza di differenze tra i gruppi per tutte e tre le variabili.
library(car)## Warning: il pacchetto 'car' è stato creato con R versione 4.4.2## Caricamento del pacchetto richiesto: carData## Warning: il pacchetto 'carData' è stato creato con R versione 4.4.2modello_base <- lm(Peso ~ ., data = data)
vif(modello_base)## GVIF Df GVIF^(1/(2*Df))
## Anni.madre 1.188698 1 1.090274
## N.gravidanze 1.187004 1 1.089497
## Fumatrici 1.008263 1 1.004123
## Gestazione 2.475330 1 1.573318
## Lunghezza 2.097481 1 1.448268
## Cranio 1.637117 1 1.279499
## Tipo.parto 1.005089 1 1.002541
## Ospedale 1.004573 2 1.001141
## Sesso 1.043444 1 1.021491
## Prematuro 1.978516 1 1.406597library(caret)## Warning: il pacchetto 'caret' è stato creato con R versione 4.4.2## Caricamento del pacchetto richiesto: ggplot2## Caricamento del pacchetto richiesto: latticeset.seed(123)
index <- createDataPartition(data$Peso, p = 0.8, list = FALSE)
trainData <- data[index, ]
testData <- data[-index, ]
dim(trainData)## [1] 2002 11dim(testData)## [1] 497 11modello_base <- lm(Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso + Prematuro,
data = trainData)
summary(modello_base)##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso + Prematuro,
## data = trainData)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1054.8 -183.0 -15.2 163.1 2538.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6485.7414 208.7668 -31.067 < 2e-16 ***
## Anni.madre 1.4693 1.2964 1.133 0.257
## N.gravidanze 8.7847 5.3436 1.644 0.100
## Fumatrici1 -47.1189 30.2905 -1.556 0.120
## Gestazione 29.6231 5.1926 5.705 1.34e-08 ***
## Lunghezza 9.9393 0.3424 29.025 < 2e-16 ***
## Cranio 10.5968 0.4822 21.977 < 2e-16 ***
## Tipo.partoNat 21.2799 13.6579 1.558 0.119
## Ospedaleosp2 -17.7828 15.2651 -1.165 0.244
## Ospedaleosp3 16.5139 15.3140 1.078 0.281
## SessoM 81.4633 12.6628 6.433 1.56e-10 ***
## Prematuro -42.4878 35.1832 -1.208 0.227
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 278.1 on 1990 degrees of freedom
## Multiple R-squared: 0.723, Adjusted R-squared: 0.7215
## F-statistic: 472.3 on 11 and 1990 DF, p-value: < 2.2e-16library(MASS)
modello_migliorato <- stepAIC(lm(Peso ~ Anni.madre + N.gravidanze + Fumatrici +
Gestazione + Lunghezza + Cranio +
Tipo.parto + Ospedale + Sesso + Prematuro,
data = trainData),
direction = "both")## Start: AIC=22545.89
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso + Prematuro
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 99323 153965952 22545
## - Prematuro 1 112759 153979387 22545
## <none> 153866629 22546
## - Fumatrici 1 187099 154053727 22546
## - Tipo.parto 1 187700 154054329 22546
## - N.gravidanze 1 208961 154075590 22547
## - Ospedale 2 396537 154263166 22547
## - Gestazione 1 2516398 156383026 22576
## - Sesso 1 3200050 157066679 22585
## - Cranio 1 37343233 191209862 22979
## - Lunghezza 1 65138056 219004685 23251
##
## Step: AIC=22545.18
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso + Prematuro
##
## Df Sum of Sq RSS AIC
## - Prematuro 1 115766 154081718 22545
## <none> 153965952 22545
## - Tipo.parto 1 187202 154153155 22546
## - Fumatrici 1 189259 154155211 22546
## + Anni.madre 1 99323 153866629 22546
## - Ospedale 2 401307 154367259 22546
## - N.gravidanze 1 396645 154362597 22548
## - Gestazione 1 2446080 156412033 22575
## - Sesso 1 3206093 157172045 22584
## - Cranio 1 37738255 191704207 22982
## - Lunghezza 1 65111070 219077022 23249
##
## Step: AIC=22544.69
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## <none> 154081718 22545
## - Tipo.parto 1 178746 154260463 22545
## - Fumatrici 1 180941 154262659 22545
## + Prematuro 1 115766 153965952 22545
## + Anni.madre 1 102330 153979387 22545
## - Ospedale 2 400916 154482634 22546
## - N.gravidanze 1 402861 154484579 22548
## - Sesso 1 3164929 157246647 22583
## - Gestazione 1 4524189 158605907 22601
## - Cranio 1 38010568 192092286 22984
## - Lunghezza 1 66004807 220086525 23257summary(modello_migliorato)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso, data = trainData)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1048.85 -184.66 -14.21 162.65 2557.36
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6603.9460 151.7161 -43.528 < 2e-16 ***
## N.gravidanze 11.2314 4.9214 2.282 0.0226 *
## Fumatrici1 -46.3162 30.2828 -1.529 0.1263
## Gestazione 32.6533 4.2696 7.648 3.16e-14 ***
## Lunghezza 9.9714 0.3413 29.212 < 2e-16 ***
## Cranio 10.6584 0.4808 22.168 < 2e-16 ***
## Tipo.partoNat 20.7570 13.6545 1.520 0.1286
## Ospedaleosp2 -17.9422 15.2578 -1.176 0.2398
## Ospedaleosp3 16.5392 15.3008 1.081 0.2799
## SessoM 80.9567 12.6561 6.397 1.97e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 278.1 on 1992 degrees of freedom
## Multiple R-squared: 0.7226, Adjusted R-squared: 0.7214
## F-statistic: 576.7 on 9 and 1992 DF, p-value: < 2.2e-16predictions <- predict(modello_migliorato, newdata = testData)
rmse <- sqrt(mean((testData$Peso - predictions)^2))
mae <- mean(abs(testData$Peso - predictions))
sst <- sum((testData$Peso - mean(testData$Peso))^2)
sse <- sum((testData$Peso - predictions)^2)
r2 <- 1 - (sse / sst)
cat("RMSE:", rmse, "\n")## RMSE: 257.7244cat("MAE:", mae, "\n")## MAE: 200.7319cat("R-squared:", r2, "\n")## R-squared: 0.752335residui <- testData$Peso - predictions
par(mfrow = c(2, 2))
plot(predictions, residui, main = "Residui vs Predetti", xlab = "Predetti", ylab = "Residui")
abline(h = 0, col = "red")
hist(residui, main = "Istogramma dei Residui", xlab = "Residui", breaks = 30)
qqnorm(residui, main = "QQ-Plot dei Residui")
qqline(residui, col = "red")
plot(testData$Peso, predictions, main = "Valori Osservati vs Predetti",
xlab = "Osservati", ylab = "Predetti")
abline(0, 1, col = "blue")
par(mfrow = c(1, 1))In questa fase del progetto ho creato il modello base e poi quello migliorato con l’AIC. ###
new_data <- data.frame(
Anni.madre = c(30, 28, 35, 32, 40),
N.gravidanze = c(3, 1, 2, 4, 2),
Fumatrici = factor(c("0", "1", "0", "0", "1"), levels = levels(trainData$Fumatrici)),
Gestazione = c(39, 40, 38, 37, 41),
Lunghezza = c(480, 472, 518, 490, 505),
Cranio = c(335, 364, 403, 384, 346),
Tipo.parto = factor(c("Nat", "Ces", "Nat", "Nat", "Ces"), levels = levels(trainData$Tipo.parto)),
Ospedale = factor(c("osp2", "osp1", "osp3", "osp2", "osp3"), levels = levels(trainData$Ospedale)),
Sesso = factor(c("F", "M", "M", "F", "M"), levels = levels(trainData$Sesso))
)
previsioni <- predict(modello_migliorato, newdata = new_data)
new_data$Predizione_Peso <- previsioni
print(previsioni)## 1 2 3 4 5
## 3062.870 3334.209 4238.108 3630.771 3531.837print(new_data)## Anni.madre N.gravidanze Fumatrici Gestazione Lunghezza Cranio Tipo.parto
## 1 30 3 0 39 480 335 Nat
## 2 28 1 1 40 472 364 Ces
## 3 35 2 0 38 518 403 Nat
## 4 32 4 0 37 490 384 Nat
## 5 40 2 1 41 505 346 Ces
## Ospedale Sesso Predizione_Peso
## 1 osp2 F 3062.870
## 2 osp1 M 3334.209
## 3 osp3 M 4238.108
## 4 osp2 F 3630.771
## 5 osp3 M 3531.837Con il modello migliorato sono state effettuate alcune previsioni su dati casuali ottendo risultati attendibili.
library(ggplot2)
ggplot(trainData, aes(x = Gestazione, y = Peso)) +
geom_point(aes(color = "Real Data"), alpha = 0.6) +
geom_smooth(method = "lm", aes(color = "Fitted Line"), se = FALSE, linetype = "dashed") +
labs(title = "Relazione tra Gestazione e Peso", x = "Settimane di Gestazione", y = "Peso (grammi)") +
theme_minimal() +
scale_color_manual(values = c("Real Data" = "blue", "Fitted Line" = "red"))## `geom_smooth()` using formula = 'y ~ x'
ggplot(trainData, aes(x = factor(Fumatrici), y = Peso, fill = factor(Fumatrici))) +
geom_boxplot() +
labs(title = "Distribuzione del Peso in base al Fumo", x = "Fumatrici (0 = No, 1 = Sì)", y = "Peso (grammi)") +
scale_fill_manual(values = c("0" = "lightblue", "1" = "lightcoral")) +
theme_minimal()
trainData$predictions <- predict(modello_migliorato, newdata = trainData)
ggplot(trainData, aes(x = Peso, y = predictions)) +
geom_point(alpha = 0.6, color = "blue") +
geom_smooth(method = "lm", color = "red", se = FALSE) +
labs(title = "Confronto tra Valori Reali e Previsioni", x = "Peso Reale (grammi)", y = "Peso Predetto (grammi)") +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
importance <- varImp(modello_migliorato, scale = FALSE)
ggplot(importance, aes(x = reorder(rownames(importance), Overall), y = Overall)) +
geom_bar(stat = "identity", fill = "steelblue") +
coord_flip() +
labs(title = "Importanza delle Variabili nel Modello", x = "Variabile", y = "Importanza") +
theme_minimal()