library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.1.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(car)
## Caricamento del pacchetto richiesto: carData
##
## Caricamento pacchetto: 'car'
##
## Il seguente oggetto è mascherato da 'package:dplyr':
##
## recode
##
## Il seguente oggetto è mascherato da 'package:purrr':
##
## some
library(lmtest)
## Caricamento del pacchetto richiesto: zoo
##
## Caricamento pacchetto: 'zoo'
##
## I seguenti oggetti sono mascherati da 'package:base':
##
## as.Date, as.Date.numeric
library(MASS)
##
## Caricamento pacchetto: 'MASS'
##
## Il seguente oggetto è mascherato da 'package:dplyr':
##
## select
library(ggplot2)
# Carico il dataset
neonati <- read.csv("neonati.csv")
summary(neonati)
## 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
##
##
##
Nella prima fase, esploro le variabili attraverso un’analisi descrittiva per comprenderne la distribuzione e identificare eventuali outlier o anomalie.
Anni.madre riporta come valore minimo 0. Osservando il dataset rilevo che 2 record riportano relativamente 0 e 1 come età al parto. Impossibile, quindi rimuovo i dati inferiori o uguali a 1.
righe_da_rimuovere <- which(neonati$Anni.madre <= 1)
if(length(righe_da_rimuovere) > 0) {
cat("Identificate", length(righe_da_rimuovere), "righe con Anni.madre pari a 0 o 1.\n")
cat("Righe rimosse (indici originali):", paste(righe_da_rimuovere, collapse = ", "), "\n")
neonati <- neonati[-righe_da_rimuovere, ]
cat("Numero di osservazioni dopo la rimozione:", nrow(neonati), "\n")
} else {
cat("Nessuna riga con Anni.madre pari a 0 o 1 trovata. Il dataset rimane invariato.\n")
}
## Identificate 2 righe con Anni.madre pari a 0 o 1.
## Righe rimosse (indici originali): 1152, 1380
## Numero di osservazioni dopo la rimozione: 2498
Procedo a saggiare le seguenti ipotesi con i test adatti:
# Converto le variabili categoriche in fattori
neonati$Fumatrici <- as.factor(neonati$Fumatrici)
neonati$Tipo.parto <- as.factor(neonati$Tipo.parto)
neonati$Ospedale <- as.factor(neonati$Ospedale)
neonati$Sesso <- as.factor(neonati$Sesso)
# Attach del dataframe
attach(neonati)
# Distribuzione delle variabili categoriche
# Distribuzione di Fumatrici
print(table(Fumatrici))
## Fumatrici
## 0 1
## 2394 104
#Distribuzione di Tipo.parto
print(table(Tipo.parto))
## Tipo.parto
## Ces Nat
## 728 1770
# Distribuzione di Ospedale
print(table(Ospedale))
## Ospedale
## osp1 osp2 osp3
## 816 848 834
#Distribuzione di Sesso
print(table(Sesso))
## Sesso
## F M
## 1255 1243
# Generazione Istogramma del Peso
ggplot(neonati, aes(x = Peso)) +
geom_histogram(binwidth = 150, fill = "skyblue", color = "darkgray") +
labs(title = "Distribuzione del Peso alla Nascita", x = "Peso (grammi)", y = "Frequenza") +
theme_minimal()
# Generazione Boxplot del Peso per Sesso
ggplot(neonati, aes(x = Sesso, y = Peso, fill = Sesso)) +
geom_boxplot() +
labs(title = "Peso alla Nascita per Sesso", x = "Sesso", y = "Peso (grammi)") +
theme_minimal()
# Generazione Boxplot del Peso per Fumatrici
ggplot(neonati, aes(x = factor(Fumatrici), y = Peso, fill = factor(Fumatrici))) +
geom_boxplot() +
labs(title = "Peso alla Nascita in base al Fumo Materno", x = "Fumatrice (0=No, 1=Sì)", y = "Peso (grammi)") +
theme_minimal()
# Generazione Scatter plot tra Gestazione e Peso
ggplot(neonati, aes(x = Gestazione, y = Peso)) +
geom_point(alpha = 0.5) +
geom_smooth(method = "lm", col = "red") +
labs(title = "Relazione tra Settimane di Gestazione e Peso alla Nascita", x = "Gestazione (settimane)", y = "Peso (grammi)") +
theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'
# Generazione Scatter plot tra Gestazione, Peso e Fumo materno
ggplot(neonati, aes(x = Gestazione, y = Peso, color = neonati$Fumatrici)) +
geom_point(alpha = 0.5, size = 1) +
geom_smooth(method = "lm", se = TRUE) +
labs(title = "Impatto della Gestazione e del Fumo Materno sul Peso alla Nascita",
x = "Settimane di Gestazione",
y = "Peso Previsto (grammi)",
color = "Fumatrice") +
scale_color_manual(values = c("0" = "darkgreen", "1" = "darkred"),
labels = c("0" = "Non Fumatrice", "1" = "Fumatrice"),
drop = FALSE) +
theme_minimal() +
theme(legend.position = "bottom")
## Warning: Use of `neonati$Fumatrici` is discouraged.
## ℹ Use `Fumatrici` instead.
## Use of `neonati$Fumatrici` is discouraged.
## ℹ Use `Fumatrici` instead.
## `geom_smooth()` using formula = 'y ~ x'
#Test di Ipotesi 1: Proporzione di parti cesarei per ospedale
table_parto_ospedale <- table(Tipo.parto, Ospedale)
print(table_parto_ospedale)
## Ospedale
## Tipo.parto osp1 osp2 osp3
## Ces 242 254 232
## Nat 574 594 602
chi_square_test <- chisq.test(table_parto_ospedale)
print(chi_square_test)
##
## Pearson's Chi-squared test
##
## data: table_parto_ospedale
## X-squared = 1.083, df = 2, p-value = 0.5819
cat("Conclusione: p-value =", round(chi_square_test$p.value, 4), ". Non significativo (p > 0.05), quindi non ci sono prove di differenza significativa nelle proporzioni di cesarei tra ospedali.\n")
## Conclusione: p-value = 0.5819 . Non significativo (p > 0.05), quindi non ci sono prove di differenza significativa nelle proporzioni di cesarei tra ospedali.
# Test di Ipotesi 2: Confronto medie campionarie vs popolazione
pop_peso_medio <- 3300
pop_lunghezza_media <- 500
# T-test per il Peso (vs media popolazione 3300g)
t_test_peso <- t.test(Peso, mu = pop_peso_medio)
print(t_test_peso)
##
## One Sample t-test
##
## data: Peso
## t = -1.505, df = 2497, p-value = 0.1324
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
## 3263.577 3304.791
## sample estimates:
## mean of x
## 3284.184
cat("Conclusione per Peso: p-value =", format(t_test_peso$p.value, scientific = TRUE, digits = 3), ". Significativo (p < 0.05), la media campionaria è diversa.\n")
## Conclusione per Peso: p-value = 1.32e-01 . Significativo (p < 0.05), la media campionaria è diversa.
# T-test per la Lunghezza (vs media popolazione 500mm)
t_test_lunghezza <- t.test(Lunghezza, mu = pop_lunghezza_media)
print(t_test_lunghezza)
##
## One Sample t-test
##
## data: Lunghezza
## t = -10.069, df = 2497, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
## 493.6628 495.7287
## sample estimates:
## mean of x
## 494.6958
cat("Conclusione per Lunghezza: p-value =", round(t_test_lunghezza$p.value, 4), ". Non significativo (p > 0.05), la media campionaria non è significativamente diversa.\n")
## Conclusione per Lunghezza: p-value = 0 . Non significativo (p > 0.05), la media campionaria non è significativamente diversa.
## Test di Ipotesi 3: Confronto misure antropometriche tra sessi
#T-test per il Peso per Sesso
t_test_peso_sesso <- t.test(Peso ~ Sesso, data = neonati)
print(t_test_peso_sesso)
##
## Welch Two Sample t-test
##
## data: Peso by Sesso
## t = -12.115, df = 2488.7, 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.4841 -207.3844
## sample estimates:
## mean in group F mean in group M
## 3161.061 3408.496
cat("Conclusione per Peso: p-value =", format(t_test_peso_sesso$p.value, scientific = TRUE, digits = 3), ". Significativo (p < 0.05), il peso medio è diverso tra i sessi.\n")
## Conclusione per Peso: p-value = 7.28e-33 . Significativo (p < 0.05), il peso medio è diverso tra i sessi.
# T-test per la Lunghezza per Sesso
t_test_lunghezza_sesso <- t.test(Lunghezza ~ Sesso, data = neonati)
print(t_test_lunghezza_sesso)
##
## Welch Two Sample t-test
##
## data: Lunghezza by Sesso
## t = -9.5823, df = 2457.3, 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.939001 -7.882672
## sample estimates:
## mean in group F mean in group M
## 489.7641 499.6750
cat("Conclusione per Lunghezza: p-value =", format(t_test_lunghezza_sesso$p.value, scientific = TRUE, digits = 3), ". Significativo (p < 0.05), la lunghezza media è diversa tra i sessi.\n")
## Conclusione per Lunghezza: p-value = 2.23e-21 . Significativo (p < 0.05), la lunghezza media è diversa tra i sessi.
# T-test per il Cranio per Sesso
t_test_cranio_sesso <- t.test(Cranio ~ Sesso, data = neonati)
print(t_test_cranio_sesso)
##
## Welch Two Sample t-test
##
## data: Cranio by Sesso
## t = -7.4366, df = 2489.4, p-value = 1.414e-13
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -6.110504 -3.560417
## sample estimates:
## mean in group F mean in group M
## 337.6231 342.4586
cat("Conclusione per Cranio: p-value =", format(t_test_cranio_sesso$p.value, scientific = TRUE, digits = 3), ". Significativo (p < 0.05), il diametro cranico medio è diverso tra i sessi.\n")
## Conclusione per Cranio: p-value = 1.41e-13 . Significativo (p < 0.05), il diametro cranico medio è diverso tra i sessi.
model_initial <- lm(Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso, data = neonati)
print(summary(model_initial))
##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso, data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1123.26 -181.53 -14.45 161.05 2611.89
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6735.7960 141.4790 -47.610 < 2e-16 ***
## Anni.madre 0.8018 1.1467 0.699 0.4845
## N.gravidanze 11.3812 4.6686 2.438 0.0148 *
## Fumatrici1 -30.2741 27.5492 -1.099 0.2719
## Gestazione 32.5773 3.8208 8.526 < 2e-16 ***
## Lunghezza 10.2922 0.3009 34.207 < 2e-16 ***
## Cranio 10.4722 0.4263 24.567 < 2e-16 ***
## Tipo.partoNat 29.6335 12.0905 2.451 0.0143 *
## Ospedaleosp2 -11.0912 13.4471 -0.825 0.4096
## Ospedaleosp3 28.2495 13.5054 2.092 0.0366 *
## SessoM 77.5723 11.1865 6.934 5.18e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274 on 2487 degrees of freedom
## Multiple R-squared: 0.7289, Adjusted R-squared: 0.7278
## F-statistic: 668.7 on 10 and 2487 DF, p-value: < 2.2e-16
model_optimal <- stepAIC(model_initial, direction = "both", trace = FALSE)
print(summary(model_optimal))
##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso, data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1113.07 -181.71 -16.66 161.08 2619.57
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6707.9252 136.0257 -49.314 < 2e-16 ***
## N.gravidanze 12.3360 4.3344 2.846 0.00446 **
## Gestazione 32.0386 3.7925 8.448 < 2e-16 ***
## Lunghezza 10.3059 0.3006 34.286 < 2e-16 ***
## Cranio 10.4920 0.4257 24.648 < 2e-16 ***
## Tipo.partoNat 29.4080 12.0875 2.433 0.01505 *
## Ospedaleosp2 -10.8939 13.4447 -0.810 0.41786
## Ospedaleosp3 28.7917 13.4969 2.133 0.03301 *
## SessoM 77.4657 11.1842 6.926 5.48e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274 on 2489 degrees of freedom
## Multiple R-squared: 0.7287, Adjusted R-squared: 0.7278
## F-statistic: 835.7 on 8 and 2489 DF, p-value: < 2.2e-16
# Calcolo RMSE
predictions <- predict(model_optimal, newdata = neonati)
rmse <- sqrt(mean((Peso - predictions)^2))
cat("RMSE del modello ottimale:", round(rmse, 2), "grammi\n")
## RMSE del modello ottimale: 273.51 grammi
# Analisi dei residui
par(mfrow = c(2, 2))
plot(model_optimal)
par(mfrow = c(1, 1))
# Verificare la multicollinearità (VIF)
vif_values <- vif(model_optimal)
print(vif_values)
## GVIF Df GVIF^(1/(2*Df))
## N.gravidanze 1.025249 1 1.012546
## Gestazione 1.670864 1 1.292619
## Lunghezza 2.083003 1 1.443261
## Cranio 1.626636 1 1.275396
## Tipo.parto 1.003840 1 1.001918
## Ospedale 1.003079 2 1.000769
## Sesso 1.040433 1 1.020016
cat("Conclusione VIF: Tutti i valori sono molto vicini a 1, indicando assenza di problemi di multicollinearità.\n")
## Conclusione VIF: Tutti i valori sono molto vicini a 1, indicando assenza di problemi di multicollinearità.
# Previsione del Peso Neonatale e creazione nuovo dataframe
mean_anni_madre <- mean(neonati$Anni.madre)
mean_lunghezza <- mean(neonati$Lunghezza)
mean_cranio <- mean(neonati$Cranio)
nuova_neonata_data <- data.frame(
Anni.madre = mean_anni_madre,
N.gravidanze = 3,
Fumatrici = factor(0, levels = levels(neonati$Fumatrici)),
Gestazione = 39,
Lunghezza = mean_lunghezza,
Cranio = mean_cranio,
Tipo.parto = factor("Nat", levels = levels(neonati$Tipo.parto)),
Ospedale = factor("osp1", levels = levels(neonati$Ospedale)),
Sesso = factor("F", levels = levels(neonati$Sesso))
)
print(str(nuova_neonata_data))
## 'data.frame': 1 obs. of 9 variables:
## $ Anni.madre : num 28.2
## $ N.gravidanze: num 3
## $ Fumatrici : Factor w/ 2 levels "0","1": 1
## $ Gestazione : num 39
## $ Lunghezza : num 495
## $ Cranio : num 340
## $ Tipo.parto : Factor w/ 2 levels "Ces","Nat": 2
## $ Ospedale : Factor w/ 3 levels "osp1","osp2",..: 1
## $ Sesso : Factor w/ 2 levels "F","M": 1
## NULL
predicted_weight <- predict(model_optimal, newdata = nuova_neonata_data)
cat("Peso previsto per la nuova neonata:", round(predicted_weight, 2), "grammi\n")
## Peso previsto per la nuova neonata: 3273.85 grammi
predicted_weight_ci <- predict(model_optimal, newdata = nuova_neonata_data, interval = "confidence")
# Previsione con intervallo di confidenza (per la media)
print(predicted_weight_ci)
## fit lwr upr
## 1 3273.847 3244.939 3302.755
predicted_weight_pi <- predict(model_optimal, newdata = nuova_neonata_data, interval = "prediction")
# Previsione con intervallo di previsione (per un'osservazione singola)
print(predicted_weight_pi)
## fit lwr upr
## 1 3273.847 2735.768 3811.926
neonati <- neonati %>%
mutate(
Fumatrici = as.factor(Fumatrici),
Tipo.parto = as.factor(Tipo.parto),
Ospedale = as.factor(Ospedale),
Sesso = as.factor(Sesso)
)
gestazione_fissa <- 39
# GRAFICO 1: GESTAZIONE vs. PESO (Ignorando il Fumo)
# Questo riflette la conclusione del modello AIC che il fumo è irrilevante.
ggplot(neonati, aes(x = Gestazione, y = Peso)) +
geom_point(alpha = 0.4, size = 1, color = "gray") +
geom_smooth(method = "lm", se = TRUE, color = "lightblue", linewidth = 1) +
labs(title = "Relazione tra Settimane di Gestazione e Peso alla Nascita",
subtitle = "La Gestazione è il predittore più forte del modello",
x = "Settimane di Gestazione",
y = "Peso alla Nascita (grammi)") +
theme_minimal() +
theme(plot.title = element_text(face = "bold"))
## `geom_smooth()` using formula = 'y ~ x'
# GRAFICO 2: DISTRIBUZIONE REALE DEL PESO PER FASCE DI GESTAZIONE
neonati$Gestazione_Fascia <- cut(neonati$Gestazione,
breaks = c(30, 36, 38, 40, 42, max(neonati$Gestazione)),
labels = c("31-36 (Pre-Termine)", "37-38 (Termine Precoce)", "39-40 (Termine Pieno)", "41-42 (Termine Tardivo)", ">=43"),
right = FALSE, include.lowest = FALSE)
ggplot(neonati, aes(x = Gestazione_Fascia, y = Peso, fill = Gestazione_Fascia)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha = 0.1, size = 0.5) +
labs(title = "Distribuzione del Peso per Fascia di Età Gestazionale",
subtitle = "Evidenziato l'aumento di peso e la riduzione della variabilità (rischio) con il termine",
x = "Fascia di Età Gestazionale (Settimane)",
y = "Peso alla Nascita (grammi)") +
scale_fill_brewer(palette = "YlGnBu") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none", plot.title = element_text(face = "bold"))
# GRAFICO 3: BOXPLOT DISTRIBUZIONE + PESO PREVISTO (solo Sesso)
previsione_solo_sesso <- data.frame(
Anni.madre = mean_anni_madre, N.gravidanze = 3, Gestazione = gestazione_fissa,
Lunghezza = mean_lunghezza, Cranio = mean_cranio, Tipo.parto = factor("Nat", levels = levels(neonati$Tipo.parto)),
Ospedale = factor("osp1", levels = levels(neonati$Ospedale)), Fumatrici = factor(0, levels = levels(neonati$Fumatrici)),
Sesso = factor(c("M", "F"), levels = levels(neonati$Sesso))
)
predictions_df <- predict(model_optimal, newdata = previsione_solo_sesso) %>% as.data.frame() %>% rename(fit = ".")
plot_data_sesso <- cbind(previsione_solo_sesso, predictions_df)
ggplot(neonati, aes(x = Sesso, y = Peso)) +
geom_boxplot(aes(fill = Sesso), alpha = 0.7) +
geom_point(data = plot_data_sesso, aes(y = fit),
size = 5, shape = 18, color = "red") +
geom_text(data = plot_data_sesso, aes(y = fit, label = paste0(round(fit, 0), " g")),
vjust = -1.5, size = 4.5, fontface = "bold", color = "red") +
labs(title = "Peso Previsto vs. Distribuzione Reale per Sesso",
subtitle = "Il punto rosso è la stima del modello AIC (a parità di altri fattori)",
x = "Sesso del Neonato",
y = "Peso (grammi)") +
scale_fill_manual(values = c("F" = "pink", "M" = "skyblue")) +
theme_minimal() +
theme(legend.position = "none", plot.title = element_text(face = "bold"))
detach(neonati)
Grafico 1: La durata della Gestazione è essenziale, ogni settimana aggiuntiva ha un impatto positivo significativo sul peso finale. Il fumo agisce come un fattore di penalizzazione costante ma non è significativo in questo modello. Raccomandazione Clinica: Mantenere il fumo come fattore di rischio primario per l’educazione sanitaria. La letteratura dimostra che i neonati di madri fumatrici seguono un percorso di crescita inferiore.
Grafico 2: Mostra in modo chiaro il rischio clinico. Le fasce “Pre-Termine” avranno pesi mediani significativamente più bassi e una maggiore variabilità (più outlier), evidenziando dove si concentra il rischio di basso peso alla nascita dovuto alla prematurità.
Grafico 3 (Peso per Sesso): La previsione si colloca nella fascia attesa di peso per entrambi i sessi.